text:
int_rules_4_5_1_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.5.1.3*(d*sin)^n*(a+b*sec)^m =#
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => Antiderivative(((~g * cos(~e + ~f * ~x)) ^ ~p * (~b + ~a * sin(~e + ~f * ~x)) ^ ~m) / sin(~e + ~f * ~x) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && IntegerQ(~m)
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-1 / (~f * (~b) ^ (~p - 1))) * Subst(Antiderivative(((-(~a) + ~b * ~x) ^ ((~p - 1) / 2) * (~a + ~b * ~x) ^ (~m + (~p - 1) / 2)) / (~x) ^ (~p + 1), ~x), ~x, csc(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (IntegerQ((~p - 1) / 2) && EqQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative(cos(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (-1 / ~f) * Subst(Antiderivative(((-1 + ~x) ^ ((~p - 1) / 2) * (1 + ~x) ^ ((~p - 1) / 2) * (~a + ~b * ~x) ^ ~m) / (~x) ^ (~p + 1), ~x), ~x, csc(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (IntegerQ((~p - 1) / 2) && NeQ((~a) ^ 2 - (~b) ^ 2, 0))
    @apply_utils Antiderivative((~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m / cos(~(e') + ~(f') * ~x) ^ 2, ~x) => (tan(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m) / ~f + ~b * ~m * Antiderivative(csc(~e + ~f * ~x) * (~a + ~b * csc(~e + ~f * ~x)) ^ (~m - 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x)
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~m, ~x) => ((sin(~e + ~f * ~x) ^ FracPart(~m) * (~a + ~b * csc(~e + ~f * ~x)) ^ FracPart(~m)) / (~b + ~a * sin(~e + ~f * ~x)) ^ FracPart(~m)) * Antiderivative(((~g * cos(~e + ~f * ~x)) ^ ~p * (~b + ~a * sin(~e + ~f * ~x)) ^ ~m) / sin(~e + ~f * ~x) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) || IntegersQ(2 * ~m, ~p))
    @apply_utils Antiderivative((~(g') * cos(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => Unintegrable((~g * cos(~e + ~f * ~x)) ^ ~p * (~a + ~b * csc(~e + ~f * ~x)) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x)
    #= Antiderivative(((~g')*sec((~e')+(~f')*(~x)))^(~p)*((~a)+(~b')*csc((~e')+(~f')*(~x)))^(~m'),~x)  := Antiderivative((g*sec(e+f*x))^p*(b+a*sin(e+f*x))^m/sin(e+f*x)^m,x) <--  FreeQ([a,b,e,f,g,p],x) && Not(IntegerQ(p)) && IntegerQ(m) =#
    @apply_utils Antiderivative((~(g') * sec(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * csc(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (~g) ^ IntPart(~p) * (~g * sec(~e + ~f * ~x)) ^ FracPart(~p) * cos(~e + ~f * ~x) ^ FracPart(~p) * Antiderivative((~a + ~b * csc(~e + ~f * ~x)) ^ ~m / cos(~e + ~f * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && Not(IntegerQ(~p))
end

max_stars_repo_path: src/rules/4 Trig functions/4.5 Secant/4.5.1.3 (d sin)^n (a+b sec)^m.jl


text:
module Compat

include("jacobi.jl")

end

ext: jl


text:
"""
    _qrE!(A, E, Q, B; withQ = true)

Reduce the regular matrix pencil `A - λE` to an equivalent form `A1 - λE1 = Q1'*(A - λE)` using an
orthogonal or unitary transformation matrix `Q1` such that the transformed matrix `E1` is upper triangular.
The reduction is performed using the QR-decomposition of E.

The performed left orthogonal or unitary transformations Q1 are accumulated in the matrix `Q <- Q*Q1`
if `withQ = true`. Otherwise, `Q` is unchanged.

`Q1'*B` is returned in `B` unless `B = missing`.
"""
function _qrE!(A::AbstractMatrix{T}, E::AbstractMatrix{T}, Q::Union{AbstractMatrix{T},Nothing},
               B::Union{AbstractVecOrMat{T},Missing} = missing; withQ::Bool = true) where T <: BlasFloat

   # fast return for dimensions 0 or 1
   size(A,1) <= 1 && return
   T <: Complex ? tran = 'C' : tran = 'T'
   # compute in-place the QR-decomposition E = Q1*E1
   _, τ = LinearAlgebra.LAPACK.geqrf!(E)
   # A <- Q1'*A
   LinearAlgebra.LAPACK.ormqr!('L',tran,E,τ,A)
   # Q <- Q*Q1
   withQ && LinearAlgebra.LAPACK.ormqr!('R','N',E,τ,Q)
   # B <- Q1'*B
   ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,E,τ,B)
   triu!(E)
end
"""
    _svdlikeAE!(A, E, Q, Z, B, C; svdA = true, fast = true, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (rE, rA22)

Reduce the regular matrix pencil `A - λE` to an equivalent form `A1 - λE1 = Q1'*(A - λE)*Z1` using
orthogonal or unitary transformation matrices `Q1` and `Z1` such that the transformed matrices `A1` and `E1`
are, for `svdA = true`, in the following SVD-like coordinate form

                   | A11-λE11 |  A12  |  A13  |
                   |----------|-------|-------|
        A1 - λE1 = |    A21   |  A22  |   0   | ,
                   |----------|-------|-------|
                   |    A31   |   0   |   0   |

where the `rE x rE` matrix `E11` and `rA22 x rA22` matrix `A22` are nosingular, and `E11` and `A22` are upper triangular,
if `fast = true`, and diagonal, if `fast = false`.

If `svdA = false`, only `E` is reduced to SVD-like form and `A1 - λE1` has the form

                   | A11-λE11 |  A12  |
        A1 - λE1 = |----------|-------| ,
                   |    A21   |  A22  |

where the `rE x rE` matrix `E11` is nonsingular upper triangular, if `fast = true`,
and diagonal, if `fast = false`, and `A22` is unreduced and has rank `rA22`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations Q1 are accumulated in the matrix `Q <- Q*Q1`
if `withQ = true`. Otherwise, `Q` is unchanged.
The performed right orthogonal or unitary transformations Z1 are accumulated in the matrix `Z <- Z*Z1` if `withZ = true`.
Otherwise, `Z` is unchanged.

`Q1'*B` is returned in `B` unless `B = missing` and `C*Z1` is returned in `C` unless `C = missing` .

"""
function _svdlikeAE!(A::AbstractMatrix{T}, E::AbstractMatrix{T},
                     Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing},
                     B::Union{AbstractVecOrMat{T},Missing} = missing, C::Union{AbstractMatrix{T},Missing} = missing;
                     svdA::Bool = true, fast::Bool = true, atol1::Real = zero(real(T)), atol2::Real = zero(real(T)),
                     rtol::Real = (size(A,1)*eps(real(float(one(T)))))*iszero(min(atol1,atol2)),
                     withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat

   # fast returns for null dimensions
   n = size(A,1)
   if n == 0
      return 0, 0
   end
   T <: Complex ? tran = 'C' : tran = 'T'

   if fast
      # compute in-place the QR-decomposition E*P1 = Q1*[E1;0] with column pivoting
      _, τ, jpvt = LinearAlgebra.LAPACK.geqp3!(E)
      tol = max(atol2, rtol*abs(E[1,1]))
      rE = count(x -> x > tol, abs.(diag(E)))
      n2 = n-rE
      i1 = 1:rE
      # A <- Q1'*A
      LinearAlgebra.LAPACK.ormqr!('L',tran,E,τ,A)
      # A <- A*P1
      A[:,:] = A[:,jpvt]
      # Q <- Q*Q1
      withQ && LinearAlgebra.LAPACK.ormqr!('R','N',E,τ,Q)
      # Z <- Z*P1
      withZ && (Z[:,:] = Z[:,jpvt])
      # B <- Q1'*B
      ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,E,τ,B)
      # E = [E1;0]
      E[:,:] = [ triu(E[i1,:]); zeros(T,n2,n) ]
      # C <- C*P1
      ismissing(C) || (C[:,:] = C[:,jpvt])

      # compute in-place the complete orthogonal decomposition E*Z1 = [E11 0; 0 0] with E11 nonsingular and UT
      E1 = view(E,i1,:)
      _, tau = LinearAlgebra.LAPACK.tzrzf!(E1)
      # A <- A*Z1
      LinearAlgebra.LAPACK.ormrz!('R',tran,E1,tau,A)
      withZ && LinearAlgebra.LAPACK.ormrz!('R',tran,E1,tau,Z)
      # C <- C*Z1
      ismissing(C) || LinearAlgebra.LAPACK.ormrz!('R',tran,E1,tau,C);
      E1[:,:] = [ triu(E[i1,i1]) zeros(T,rE,n2)  ]
      n2 == 0 && (return rE, 0)
      i22 = rE+1:n
      tolA = max(atol1, rtol*opnorm(A,1))
      svdA || (return rE, rank(view(A,i22,i22), atol = tolA))
      # assume
      #    A = [A11 A12]
      #        [A21 A22]
      # compute in-place the QR-decomposition A22*P2 = Q2*[R2;0] with column pivoting
      A22 = view(A,i22,i22)
      _, τ, jpvt = LinearAlgebra.LAPACK.geqp3!(A22)
      rA22 = count(x -> x > tolA, abs.(diag(A22)))
      n3 = n2-rA22
      i2 = rE+1:rE+rA22
      i3 = rE+rA22+1:n
      # A21 <- Q2'*A21
      LinearAlgebra.LAPACK.ormqr!('L',tran,A22,τ,view(A,i22,i1))
      # A12 <- A12*P1
      A[i1,i22] = A[i1,i22[jpvt]]
      # Q <- Q*Q2
      withQ && LinearAlgebra.LAPACK.ormqr!('R','N',A22,τ,view(Q,:,i22))
      # Z <- Z*P2
      withZ && (Z[:,i22] = Z[:,i22[jpvt]])
      # B <- Q2'*B
      ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,A22,τ,view(B,i22,:))
      # A22 = [R2;0]
      A22[:,:] = [ triu(A[i2,i22]); zeros(T,n3,n2) ]
      # R <- R*P1
      ismissing(C) || (C[:,i22] = C[:,i22[jpvt]])

      # compute in-place the complete orthogonal decomposition R2*Z2 = [A2 0; 0 0] with A22 nonsingular and UT
      A2 = view(A,i2,i22)
      _, tau = LinearAlgebra.LAPACK.tzrzf!(A2)
      # A12 <- A12*Z2
      LinearAlgebra.LAPACK.ormrz!('R',tran,A2,tau,view(A,i1,i22))
      withZ && LinearAlgebra.LAPACK.ormrz!('R',tran,A2,tau,view(Z,:,i22))
      # C <- C*Z2
      ismissing(C) || LinearAlgebra.LAPACK.ormrz!('R',tran,A2,tau,view(C,:,i22));
      A2[:,:] = [ triu(A[i2,i2]) zeros(T,rA22,n3)  ]
   else
      # compute the complete orthogonal decomposition of E using the SVD-decomposition
      U, S, Vt = LinearAlgebra.LAPACK.gesdd!('A',E)
      tolE = max(atol2, rtol*S[1])
      rE = count(x -> x > tolE, S)
      n2 = n-rE
      # A <- A*V
      A[:,:] = A[:,:]*Vt'
      # A <- U'*A
      A[:,:] = U'*A[:,:]
      # Q <- Q*U
      withQ && (Q[:,:] = Q[:,:]*U)
      # Z <- Q*V
      withZ && (Z[:,:] = Z[:,:]*Vt')
      # B <- U'*B
      ismissing(B) || (B[:,:] = U'*B[:,:])
      # C <- C*V
      ismissing(C) || (C[:,:] = C[:,:]*Vt')
      E[:,:] = [ Diagonal(S[1:rE]) zeros(T,rE,n2) ; zeros(T,n2,n) ]
      n2 == 0 && (return rE, 0)
      i22 = rE+1:n
      tolA = max(atol1, rtol*opnorm(A,1))
      svdA || (return rE, rank(view(A,i22,i22), atol = tolA))
      # assume
      #    A = [A11 A12]
      #        [A21 A22]
      # compute the complete orthogonal decomposition of A22 using the SVD-decomposition
      A22 = view(A,i22,i22)
      U, S, Vt = LinearAlgebra.LAPACK.gesdd!('A',A22)
      rA22 = count(x -> x > tolA, S)
      n3 = n2-rA22
      i1 = 1:rE
      i2 = rE+1:rE+rA22
      i3 = rE+rA22+1:n
      # A12 <- A12*V
      A[i1,i22] = A[i1,i22]*Vt'
      # A21 <- U'*A21
      A[i22,i1] = U'*A[i22,i1]
      # Q <- Q*U
      withQ && (Q[:,i22] = Q[:,i22]*U)
      # Z <- Q*V
      withZ && (Z[:,i22] = Z[:,i22]*Vt')
      # B <- U'*B
      ismissing(B) || (B[i22,:] = U'*B[i22,:])
      # C <- C*V
      ismissing(C) || (C[:,i22] = C[:,i22]*Vt')
      A22[:,:] = [ Diagonal(S[1:rA22]) zeros(T,rA22,n3); zeros(T,n3,n2) ]
   end
   return rE, rA22
end
"""
    isregular(A, E, γ; atol::Real = 0,  rtol::Real = atol > 0 ? 0 : n*ϵ) -> Bool

Test whether the matrix pencil `A-λE` is regular at `λ = γ` (i.e., `A-λE` is square and ``{\\small\\det(A-λE) \\neq 0}``).
The underlying computational procedure checks the maximal rank of `A-γE` if `γ` is finite and of `E` if
`γ` is infinite .

The keyword arguements `atol` and `rtol` specify the absolute and relative tolerances for the nonzero
elements of `A-γE`, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of  `A`, and `ϵ` is the
machine epsilon of the element type of `A`.
"""
function isregular(A::AbstractMatrix, E::Union{AbstractMatrix,Nothing}, γ::Number; atol::Real = zero(real(eltype(A))),
                   rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(atol))

   m, n = size(A)
   m == n || (return false)
   E === nothing && (return rank(A, atol = atol, rtol = rtol) == n )
   (m,n) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
   if isinf(γ)
      return rank(E,atol = atol,rtol=rtol) == n
   else
      return rank(A-γ*E,atol = atol,rtol=rtol) == n
   end
end
"""
    isregular(A, E; atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> Bool

Test whether the matrix pencil `A-λE` is regular (i.e., `A-λE` is square and ``{\\small\\det(A-λE) \\not\\equiv 0}``).
The underlying computational procedure reduces the pencil `A-λE` to an appropriate Kronecker-like form,
which provides information on the rank of `A-λE`.

The keyword arguements `atol1`, `atol2` and `rtol` specify the absolute tolerance for the nonzero
elements of `A`, the absolute tolerance for the nonzero elements of `E`, and the relative tolerance
for the nonzero elements of `A` and `E`, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of  `A`, and `ϵ` is the
machine epsilon of the element type of `A`.
"""
function isregular(A::AbstractMatrix, E::Union{AbstractMatrix,Nothing}; atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))),
                   rtol::Real = (min(size(A)...)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)))
   mA, nA = size(A)
   mA == nA || (return false)
   E === nothing && (return rank(A, atol = atol1, rtol = rtol) == mA )
   (mA,nA) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
   T = promote_type(eltype(A), eltype(E))
   T <: BlasFloat || (T = promote_type(Float64,T))
   A1 = copy_oftype(A,T)
   E1 = copy_oftype(E,T)

   Q = nothing
   Z = nothing

   # Step 0: Reduce to the standard form
   n, m, p = _preduceBF!(A1, E1, Q, Z; atol = atol2, rtol = rtol, fast = false, withQ = false, withZ = false)

   mrinf = 0
   tol1 = max(atol1, rtol*opnorm(A1,1))
   while m > 0
      # Steps 1 & 2: Standard algorithm PREDUCE
      i1 = mrinf+1:mA
      τ, ρ = _preduce1!(n, m, m, view(A1,i1,i1), view(E1,i1,i1), Q, Z, tol1; fast = false, withQ = false, withZ = false)
      ρ+τ == m || (return false)
      mrinf += m
      n -= ρ
      m = ρ
   end
   return true
end
"""
    isunimodular(A, E; atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> Bool

Test whether the matrix pencil `A-λE` is unimodular (i.e., `A-λE` is square, regular and `det(A-λE) == constant`).
The underlying computational procedure reduces the pencil `A-λE` to an appropriate Kronecker-like form,
which provides information to check the full rank of `A-λE` and the lack of finite eigenvalues.

The keyword arguements `atol1`, `atol2` and `rtol` specify the absolute tolerance for the nonzero
elements of `A`, the absolute tolerance for the nonzero elements of `E`, and the relative tolerance
for the nonzero elements of `A` and `E`, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of  `A`, and `ϵ` is the
machine epsilon of the element type of `A`.
"""
function isunimodular(A::AbstractMatrix, E::Union{AbstractMatrix,Nothing}; atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(A))),
                      rtol::Real = (size(A,1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)))

   mA, nA = size(A)
   mA == nA || (return false)
   E === nothing && (return rank(A, atol = atol1, rtol = rtol) == mA )
   (mA,nA) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
   mA == 0 && (return true)
   T = promote_type(eltype(A), eltype(E))
   T <: BlasFloat || (T = promote_type(Float64,T))
   A1 = copy_oftype(A,T)
   E1 = copy_oftype(E,T)

   Q = nothing
   Z = nothing

   # Step 0: Reduce to the standard form
   n, m, p = _preduceBF!(A1, E1, Q, Z; atol = atol2, rtol = rtol, fast = false, withQ = false, withZ = false)
   n == 0 && (return true)
   mrinf = 0
   tol1 = max(atol1, rtol*opnorm(A1,1))
   while m > 0
      # Steps 1 & 2: Standard algorithm PREDUCE
      i1 = mrinf+1:mA
      τ, ρ = _preduce1!(n, m, m, view(A1,i1,i1), view(E1,i1,i1), Q, Z, tol1; fast = false, withQ = false, withZ = false)
      ρ+τ == m || (return false)
      mrinf += m
      n -= ρ
      m = ρ
   end
   return n == 0
end
"""
    fisplit(A, E, B, C; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (At, Et, Bt, Ct, Q, Z, ν, blkdims)

Reduce the regular matrix pencil `A - λE` to an equivalent form `At - λEt = Q'*(A - λE)*Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `At` and `Et` are in one of the
following block upper-triangular forms:

(1) if `finite_infinite = false`, then

                   | Ai-λEi |   *     |
        At - λEt = |--------|---------|,
                   |    O   | Af-λEf  |

where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary
divisors and the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, contains the finite eigenvalues of the pencil `A-λE`.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni, nf)`.

(2) if `finite_infinite = true`, then

                   | Af-λEf |   *    |
        At - λEt = |--------|--------|,
                   |   O    | Ai-λEi |

where the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular,
contains the finite eigenvalues of the pencil `A-λE` and the `ni x ni` subpencil `Ai-λEi`
contains the infinite elementary divisors.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular  and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `j` (with `ν[0] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (nf, ni)`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.

`Bt = Q'*B`, unless `B = missing`, in which case `Bt = missing` is returned, and `Ct = C*Z`,
unless `C = missing`, in which case `Ct = missing` is returned .

"""
function fisplit(A::AbstractMatrix, E::AbstractMatrix, B::Union{AbstractVecOrMat,Missing}, C::Union{AbstractMatrix,Missing};
   fast::Bool = true, finite_infinite::Bool = false,
   atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
   rtol::Real = (size(A,1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)),
   withQ::Bool = true, withZ::Bool = true)

   n = LinearAlgebra.checksquare(A)
   n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
   (!ismissing(B) && n != size(B,1)) && throw(DimensionMismatch("A and B must have the same number of rows"))
   (!ismissing(C) && n != size(C,2)) && throw(DimensionMismatch("A and C must have the same number of columns"))
   T = promote_type(eltype(A), eltype(E))
   ismissing(B) || (T = promote_type(T,eltype(B)))
   ismissing(C) || (T = promote_type(T,eltype(C)))
   T <: BlasFloat || (T = promote_type(Float64,T))

   A1 = copy_oftype(A,T)
   E1 = copy_oftype(E,T)
   ismissing(B) ? B1 = missing : B1 = copy_oftype(B,T)
   ismissing(C) ? C1 = missing : C1 = copy_oftype(C,T)

   withQ ? (Q = Matrix{T}(I,n,n)) : (Q = nothing)
   withZ ? (Z = Matrix{T}(I,n,n)) : (Z = nothing)

   ν, blkdims = fisplit!(A1, E1, Q, Z, B1, C1;
                        fast  = fast, finite_infinite = finite_infinite,
                        atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = withQ, withZ = withZ)

   return A1, E1, B1, C1, Q, Z, ν, blkdims
end
"""
    fisplit!(A, E, Q, Z, B, C; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (ν, blkdims)

Reduce the regular matrix pencil `A - λE` to an equivalent form `At - λEt = Q1'*(A - λE)*Z1` using
orthogonal or unitary transformation matrices `Q1` and `Z1` such that the transformed matrices `At` and `Et` are in one of the
following block upper-triangular forms:

(1) if `finite_infinite = false`, then

                   | Ai-λEi |   *     |
        At - λEt = |--------|---------|,
                   |    O   | Af-λEf  |

where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary
divisors and the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, contains the finite eigenvalues of the pencil `A-λE`.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni, nf)`.

(2) if `finite_infinite = true`, then

                   | Af-λEf |   *    |
        At - λEt = |--------|--------|,
                   |   O    | Ai-λEi |

where the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular,
contains the finite eigenvalues of the pencil `A-λE` and the `ni x ni` subpencil `Ai-λEi`
contains the infinite elementary divisors.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `j` (with `ν[0] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (nf, ni)`.

The reduced matrices `At` and `Et` are returned in `A` and `E`, respectively,
while `Q1'*B` is returned in `B`, unless `B = missing`, and `C*Z1`, is returned in `C`,
unless `C = missing`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations Q1 are accumulated in the matrix `Q <- Q*Q1`
if `withQ = true`. Otherwise, `Q` is unchanged.
The performed right orthogonal or unitary transformations Z1 are accumulated in the matrix `Z <- Z*Z1`
if `withZ = true`. Otherwise, `Z` is unchanged.

"""
function fisplit!(A::AbstractMatrix{T}, E::AbstractMatrix{T},
                  Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing},
                  B::Union{AbstractVecOrMat{T},Missing} = missing, C::Union{AbstractMatrix{T},Missing} = missing;
                  fast::Bool = true, finite_infinite::Bool = false,
                  atol1::Real = zero(real(T)), atol2::Real = zero(real(T)),
                  rtol::Real = (size(A,1)*eps(real(float(one(T)))))*iszero(min(atol1,atol2)),
                  withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat

   # fast returns for null dimensions
   n = size(A,1)
   ν = Vector{Int}(undef,n)
   if n == 0
      return ν, (0, 0)
   end
   T <: Complex ? tran = 'C' : tran = 'T'

   # Step 0: Reduce to the standard form
   nf, m1, p1 = _preduceBF!(A, E, Q, Z, B, C; atol = atol2, rtol = rtol, fast = fast)

   tolA = max(atol1, rtol*opnorm(A,1))

   if finite_infinite

      # Reduce A-λE to the Kronecker-like form by splitting the finite-infinite structures
      #
      #                  | Af - λ Ef |     *      |
      #      At - λ Et = |-----------|------------|,
      #                  |    0      | Ai -  λ Ei |
      #
      # where Ai - λ Ei is in a staircase form.

      i = 0
      ni = 0
      while p1 > 0
         # Step 1 & 2: Dual algorithm PREDUCE
         τ, ρ  = _preduce2!(nf, m1, p1, A, E, Q, Z, tolA, B, C; fast = fast,
                            rtrail = ni, ctrail = ni, withQ = withQ, withZ = withZ)
         ρ+τ == p1 || error("A-λE is not regular")
         ni += p1
         nf -= ρ
         i += 1
         ν[i] = p1
         p1 = ρ
         m1 -= τ
      end
      reverse!(view(ν,1:i))
      klf_left_refineinf!(view(ν,1:i), A, E, Q, B; roff = nf, coff = nf, withQ = withQ)
      return ν[1:i], (nf, ni)
   else

      # Reduce A-λE to the Kronecker-like form by splitting the infinite-finite structures
      #
      #                  | Ai - λ Ei |    *      |
      #      At - λ Et = |-----------|-----------|
      #                  |    0      | Af - λ Ef |
      #
      # where Ai - λ Ei is in a staircase form.

      i = 0
      ni = 0
      while m1 > 0
         # Steps 1 & 2: Standard algorithm PREDUCE
         τ, ρ = _preduce1!(nf, m1, p1, A, E, Q, Z, tolA, B, C; fast = fast,
                           roff = ni, coff = ni, withQ = withQ, withZ = withZ)
         ρ+τ == m1 || error("A-λE is not regular")
         ni += m1
         nf -= ρ
         i += 1
         ν[i] = m1
         m1 = ρ
         p1 -= τ
      end

      klf_right_refineinf!(view(ν,1:i), A, E, Z, C; withZ = withZ)
      # the following code can be used to make the supradiagonal blocks of E upper triangular
      # but it is not necessary in this context
      # k2 = ni
      # for k = i:-1:1
      #     nk = ν[k]
      #     k1 = k2-nk+1
      #     kk = k1:k2
      #     if nk > 1
      #        Ak = view(A,kk,kk)
      #        tau = similar(A,nk)
      #        LinearAlgebra.LAPACK.gerqf!(Ak,tau)
      #        LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(A,1:k1-1,kk))
      #        withZ && LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(Z,:,kk))
      #        LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(E,1:k1-1,kk))
      #        ismissing(C) || LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(C,:,kk))
      #        triu!(Ak)
      #     end
      #     if k > 1
      #        nk1 = ν[k-1]
      #        k1e = k1 - nk1
      #        k2e = k1 - 1
      #        kke = k1e:k2e
      #        if nk1 > 1
      #           Ek = view(E,kke,kk)
      #           tau = similar(A,nk)
      #           LinearAlgebra.LAPACK.geqrf!(Ek,tau)
      #           LinearAlgebra.LAPACK.ormqr!('L',tran,Ek,tau,view(A,kke,k1e:n))
      #           withQ && LinearAlgebra.LAPACK.ormqr!('R','N',Ek,tau,view(Q,:,kke))
      #           LinearAlgebra.LAPACK.ormqr!('L',tran,Ek,tau,view(E,kke,k2+1:n))
      #           ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,Ek,tau,view(B,kke,:))
      #           triu!(Ek)
      #        end
      #     end
      #     k2 = k1-1
      # end

      return ν[1:i], (ni, nf)
   end
end
"""
    sfisplit(A, E, B, C; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (At, Et, Bt, Ct, Q, Z, ν, blkdims)

Reduce the regular matrix pencil `A - λE` to an equivalent form `At - λEt = Q'*(A - λE)*Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `At` and `Et` are in one of the
following block upper-triangular forms:

(1) if `finite_infinite = true`, then

                   | Ai1    *        *     |
        At - λEt = | O    Af-λEf     *     |
                   | O      0     Ai2-λEi2 |

where the `ni1 x ni1` upper triangular nonsingular matrix `Ai1` and the `ni2 x ni2` subpencil `Ai2-λEi2` contain the infinite elementary
divisors and the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, contains the finite eigenvalues of the pencil `A-λE`.

The subpencil `Ai2-λEi2` is in a staircase form, with `Ai2` nonsingular and upper triangular and `Ei2` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai2-λEi2` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+2]-ν[nb-j+1]` for `j = 1, 2, ..., nb+1` is the number of infinite elementary
divisors of degree `j` (with `ν[0] := 0` and `ν[nb+1] := ni1`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni1, nf, ni2)`.

(2) if `finite_infinite = false`, then

                   | Ai1-λEi1    *      *  |
        At - λEt = | O         Af-λEf   *  |
                   | O           0     Ai2 |

where the `ni1 x ni1` subpencil `Ai1-λEi1` and the upper triangular nonsingular matrix `Ai2` contain the infinite elementary
divisors and the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, contains the finite eigenvalues of the pencil `A-λE`.

The subpencil `Ai1-λEi1` is in a staircase form, with `Ai1` nonsingular and upper triangular and `Ei1` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai1-λEi1` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 0, 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0` and `ν[0] = ni2`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni1, nf, ni2)`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.

`Bt = Q'*B`, unless `B = missing`, in which case `Bt = missing` is returned, and `Ct = C*Z`,
unless `C = missing`, in which case `Ct = missing` is returned .
"""
function sfisplit(A::AbstractMatrix, E::AbstractMatrix, B::Union{AbstractVecOrMat,Missing}, C::Union{AbstractMatrix,Missing};
   fast::Bool = true, finite_infinite::Bool = false,
   atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
   rtol::Real = (size(A,1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)),
   withQ::Bool = true, withZ::Bool = true)

   n = LinearAlgebra.checksquare(A)
   n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
   (!ismissing(B) && n != size(B,1)) && throw(DimensionMismatch("A and B must have the same number of rows"))
   (!ismissing(C) && n != size(C,2)) && throw(DimensionMismatch("A and C must have the same number of columns"))
   T = promote_type(eltype(A), eltype(E))
   ismissing(B) || (T = promote_type(T,eltype(B)))
   ismissing(C) || (T = promote_type(T,eltype(C)))
   T <: BlasFloat || (T = promote_type(Float64,T))

   A1 = copy_oftype(A,T)
   E1 = copy_oftype(E,T)
   ismissing(B) ? B1 = missing : B1 = copy_oftype(B,T)
   ismissing(C) ? C1 = missing : C1 = copy_oftype(C,T)

   withQ ? (Q = Matrix{T}(I,n,n)) : (Q = nothing)
   withZ ? (Z = Matrix{T}(I,n,n)) : (Z = nothing)

   ν, blkdims = sfisplit!(A1, E1, Q, Z, B1, C1;
                        fast  = fast, finite_infinite = finite_infinite,
                        atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = withQ, withZ = withZ)

   return A1, E1, B1, C1, Q, Z, ν, blkdims
end
"""
    sfisplit!(A, E, Q, Z, B, C; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (ν, blkdims)

Reduce the regular matrix pencil `A - λE` to an equivalent form `At - λEt = Q1'*(A - λE)*Z1` using
orthogonal or unitary transformation matrices `Q1` and `Z1` such that the transformed matrices `At` and `Et` are in one of the
following block upper-triangular forms:

(1) if `finite_infinite = true`, then

                   | Ai1    *        *     |
        At - λEt = | O    Af-λEf     *     |
                   | O      0     Ai2-λEi2 |

where the `ni1 x ni1` upper triangular nonsingular matrix `Ai1` and the `ni2 x ni2` subpencil `Ai2-λEi2` contain the infinite elementary
divisors and the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, contains the finite eigenvalues of the pencil `A-λE`.

The subpencil `Ai2-λEi2` is in a staircase form, with `Ai2` nonsingular and upper triangular and `Ei2` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai2-λEi2` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+2]-ν[nb-j+1]` for `j = 1, 2, ..., nb+1` is the number of infinite elementary
divisors of degree `j` (with `ν[0] := 0` and `ν[nb+1] := ni1`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni1, nf, ni2)`.

(2) if `finite_infinite = false`, then

                   | Ai1-λEi1    *      *  |
        At - λEt = | O         Af-λEf   *  |
                   | O           0     Ai2 |

where the `ni1 x ni1` subpencil `Ai1-λEi1` and the upper triangular nonsingular matrix `Ai2` contain the infinite elementary
divisors and the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, contains the finite eigenvalues of the pencil `A-λE`.

The subpencil `Ai1-λEi1` is in a staircase form, with `Ai1` nonsingular and upper triangular and `Ei1` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai1-λEi1` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 0, 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0` and `ν[0] = ni2`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni1, nf, ni2)`.

The reduced matrices `At` and `Et` are returned in `A` and `E`, respectively,
while `Q1'*B` is returned in `B`, unless `B = missing`, and `C*Z1`, is returned in `C`,
unless `C = missing`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations Q1 are accumulated in the matrix `Q <- Q*Q1`
if `withQ = true`. Otherwise, `Q` is unchanged.
The performed right orthogonal or unitary transformations Z1 are accumulated in the matrix `Z <- Z*Z1`
if `withZ = true`. Otherwise, `Z` is unchanged.
"""
function sfisplit!(A::AbstractMatrix{T}, E::AbstractMatrix{T},
                  Q::Union{AbstractMatrix{T},Nothing}, Z::Union{AbstractMatrix{T},Nothing},
                  B::Union{AbstractVecOrMat{T},Missing} = missing, C::Union{AbstractMatrix{T},Missing} = missing;
                  fast::Bool = true, finite_infinite::Bool = false,
                  atol1::Real = zero(real(T)), atol2::Real = zero(real(T)),
                  rtol::Real = (size(A,1)*eps(real(float(one(T)))))*iszero(min(atol1,atol2)),
                  withQ::Bool = true, withZ::Bool = true) where T <: BlasFloat

   # fast returns for null dimensions
   n = size(A,1)
   ν = Vector{Int}(undef,n)
   if n == 0
      return ν, (0, 0, 0)
   end
   T <: Complex ? tran = 'C' : tran = 'T'

   # Step 0: Reduce to the standard form
   nf, m1, p1 = _preduceBF!(A, E, Q, Z, B, C; atol = atol2, rtol = rtol, fast = fast, withQ = withQ, withZ = withZ)

   tolA = max(atol1, rtol*opnorm(A,1))
   Q === nothing && (withQ = false)
   Z === nothing && (withZ = false)

   if finite_infinite
      # Reduce A-λE to the form
      #
      #                  | Ai1  |    *      |
      #      At - λ Et = |------|-----------|
      #                  |  0   | A2 - λ E2 |
      #
      # where Ai1 is upper triangular and nonsingular and E2 upper triangular.
      τ, ρ = _preduce1!(nf, m1, p1, A, E, Q, Z, tolA, B, C; fast = fast,
                        roff = 0, coff = 0, withQ = withQ, withZ = withZ)
      ρ+τ == m1 || error("A-λE is not regular")
      nf -= ρ
      ni1 = m1
      m1 = ρ
      p1 -= τ
      if ni1 > 1
         kk = 1:ni1
         Ak = view(A,kk,kk)
         tau = similar(A,ni1)
         LinearAlgebra.LAPACK.gerqf!(Ak,tau)
         withZ && LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(Z,:,kk))
         ismissing(C) || LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(C,:,kk))
         triu!(Ak)
      end

      # Reduce A2-λE2 to the Kronecker-like form by splitting the finite-infinite structures
      #
      #                  | Af - λ Ef |     *       |
      #      A2 - λ E2 = |-----------|-------------|,
      #                  |    0      | Ai2 - λ Ei2 |
      #
      # where Ai2 - λ Ei2 is in a staircase form with Ai2 nonsingular and upper triangular and
      # Ei2 upper triangular and nilpotent.

      i = 0
      ni = 0
      while p1 > 0
         # Step 1 & 2: Dual algorithm PREDUCE
         τ, ρ  = _preduce2!(nf, m1, p1, A, E, Q, Z, tolA, B, C; fast = fast,
                            roff = ni1, coff = ni1, rtrail = ni, ctrail = ni, withQ = withQ, withZ = withZ)
         ρ+τ == p1 || error("A-λE is not regular")
         ni += p1
         nf -= ρ
         i += 1
         ν[i] = p1
         p1 = ρ
         m1 -= τ
      end
      k1 = ni1+nf+1
      reverse!(view(ν,1:i))
      for k = 1:i
          nk = ν[k]
          k2 = k1+nk-1
          kk = k1:k2
          if nk > 1
             Ak = view(A,kk,kk)
             tau = similar(A,nk)
             LinearAlgebra.LAPACK.geqrf!(Ak,tau)
             LinearAlgebra.LAPACK.ormqr!('L',tran,Ak,tau,view(A,kk,k2+1:n))
             withQ && LinearAlgebra.LAPACK.ormqr!('R','N',Ak,tau,view(Q,:,kk))
             LinearAlgebra.LAPACK.ormqr!('L',tran,Ak,tau,view(E,kk,k2+1:n))
             ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,Ak,tau,view(B,kk,:))
             triu!(Ak)
          end
          k1 = k2+1
      end
      return ν[1:i], (ni1, nf, ni)
   else
      # Reduce A-λE to the form
      #
      #                  | A1 - λ E1 |  *   |
      #      At - λ Et = |-----------|------|,
      #                  |    0      | Ai2  |
      #
      # where Ai2 is upper triangular and nonsingular and E1 is upper triangular.
      τ, ρ  = _preduce2!(nf, m1, p1, A, E, Q, Z, tolA, B, C; fast = fast,
                         rtrail = 0, ctrail = 0, withQ = withQ, withZ = withZ)
      ρ+τ == p1 || error("A-λE is not regular")
      ni2 = p1
      nf -= ρ
      p1 = ρ
      m1 -= τ
      if ni2 > 1
         kk = n-ni2+1:n
         Ak = view(A,kk,kk)
         tau = similar(A,ni2)
         LinearAlgebra.LAPACK.geqrf!(Ak,tau)
         #LinearAlgebra.LAPACK.ormqr!('L',tran,Ak,tau,view(A,kk,k2+1:n))
         withQ && LinearAlgebra.LAPACK.ormqr!('R','N',Ak,tau,view(Q,:,kk))
         #LinearAlgebra.LAPACK.ormqr!('L',tran,Ak,tau,view(E,kk,k2+1:n))
         ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,Ak,tau,view(B,kk,:))
         triu!(Ak)
      end

      # Reduce A1-λE1 to the Kronecker-like form by splitting the infinite-finite structures
      #
      #                  | Ai1 - λ Ei1 |    *      |
      #      A1 - λ E1 = |-------------|-----------|
      #                  |    0        | Af - λ Ef |
      #
      # where Ai1 - λ Ei1 is in a staircase form with Ai1 nonsingular and upper triangular and
      # Ei1 upper triangular and nilpotent.

      i = 0
      ni = 0
      while m1 > 0
         # Steps 1 & 2: Standard algorithm PREDUCE
         τ, ρ = _preduce1!(nf, m1, p1, A, E, Q, Z, tolA, B, C; fast = fast,
                           roff = ni, coff = ni, rtrail = ni2, ctrail = ni2, withQ = withQ, withZ = withZ)
         ρ+τ == m1 || error("A-λE is not regular")
         ni += m1
         nf -= ρ
         i += 1
         ν[i] = m1
         m1 = ρ
         p1 -= τ
      end

      k2 = ni
      for k = i:-1:1
          nk = ν[k]
          k1 = k2-nk+1
          kk = k1:k2
          if nk > 1
             Ak = view(A,kk,kk)
             tau = similar(A,nk)
             LinearAlgebra.LAPACK.gerqf!(Ak,tau)
             LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(A,1:k1-1,kk))
             withZ && LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(Z,:,kk))
             LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(E,1:k1-1,kk))
             ismissing(C) || LinearAlgebra.LAPACK.ormrq!('R',tran,Ak,tau,view(C,:,kk))
             triu!(Ak)
          end
         #  the following code can be used to make the supradiagonal blocks of E upper triangular
         #  but it is not necessary in this context
         #  if k > 1
         #     nk1 = ν[k-1]
         #     k1e = k1 - nk1
         #     k2e = k1 - 1
         #     kke = k1e:k2e
         #     if nk1 > 1
         #        Ek = view(E,kke,kk)
         #        tau = similar(A,nk)
         #        LinearAlgebra.LAPACK.geqrf!(Ek,tau)
         #        LinearAlgebra.LAPACK.ormqr!('L',tran,Ek,tau,view(A,kke,k1e:n))
         #        withQ && LinearAlgebra.LAPACK.ormqr!('R','N',Ek,tau,view(Q,:,kke))
         #        LinearAlgebra.LAPACK.ormqr!('L',tran,Ek,tau,view(E,kke,k2+1:n))
         #        ismissing(B) || LinearAlgebra.LAPACK.ormqr!('L',tran,Ek,tau,view(B,kke,:))
         #        triu!(Ek)
         #     end
         k2 = k1-1
      end

      return ν[1:i], (ni, nf, ni2)
   end
end


text:
module CommunalHelperConnectedSwapBlock

using ..Ahorn, Maple
using Ahorn.CommunalHelper

function swapFinalizer(entity)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    entity.data["nodes"] = [(x + width, y)]
end

@mapdef Entity "CommunalHelper/ConnectedSwapBlock" ConnectedSwapBlock(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockWidth,
    theme::String="Normal",
    customGreenBlockTexture::String="",
    customRedBlockTexture::String="",
)

const placements = Ahorn.PlacementDict(
    "Connected Swap Block ($theme) (Communal Helper)" => Ahorn.EntityPlacement(
        ConnectedSwapBlock,
        "rectangle",
        Dict{String,Any}(
            "theme" => theme,
        ),
        swapFinalizer,
    ) for theme in Maple.swap_block_themes
)

placements["Connected Swap Block (Reskinnable) (Communal Helper)"] = Ahorn.EntityPlacement(
    ConnectedSwapBlock,
    "rectangle",
    Dict{String,Any}(
        "customGreenBlockTexture" => "CommunalHelper/customConnectedBlock/customConnectedBlock",
        "customRedBlockTexture" => "CommunalHelper/customConnectedBlock/customConnectedBlock",
    ),
    swapFinalizer,
)

Ahorn.editingOptions(entity::ConnectedSwapBlock) = Dict{String,Any}(
    "theme" => Maple.swap_block_themes,
)

Ahorn.nodeLimits(entity::ConnectedSwapBlock) = 1, 1

Ahorn.minimumSize(entity::ConnectedSwapBlock) = 16, 16
Ahorn.resizable(entity::ConnectedSwapBlock) = true, true

function Ahorn.selection(entity::ConnectedSwapBlock)
    x, y = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return [
        Ahorn.Rectangle(x, y, width, height),
        Ahorn.Rectangle(stopX, stopY, width, height),
    ]
end

function getTextures(entity::ConnectedSwapBlock)
    theme = lowercase(get(entity, "theme", "normal"))
    themePath = (theme == "normal") ? "" : string(theme, "/")

    return (
        "objects/swapblock/$(themePath)blockRed",
        "objects/swapblock/$(themePath)target",
        "objects/swapblock/$(themePath)midBlockRed00",
        "objects/CommunalHelper/connectedSwapBlock/$(themePath)innerCornersRed",
    )
end

function renderTrail(ctx, x::Number, y::Number, width::Number, height::Number, trail::String)
    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + 2, 6, 0, 8, 6)
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + height - 8, 6, 14, 8, 6)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, trail, x + 2, y + (i - 1) * 8, 0, 6, 6, 8)
        Ahorn.drawImage(ctx, trail, x + width - 8, y + (i - 1) * 8, 14, 6, 6, 8)
    end

    for i in 2:tilesWidth - 1, j in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, trail, x + (i - 1) * 8, y + (j - 1) * 8, 6, 6, 8, 8)
    end

    Ahorn.drawImage(ctx, trail, x + width - 8, y + 2, 14, 0, 6, 6)
    Ahorn.drawImage(ctx, trail, x + width - 8, y + height - 8, 14, 14, 6, 6)
    Ahorn.drawImage(ctx, trail, x + 2, y + 2, 0, 0, 6, 6)
    Ahorn.drawImage(ctx, trail, x + 2, y + height - 8, 0, 14, 6, 6)
end

function renderSwapBlock(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, width::Number, height::Number, midResource::String, block::String, innerCorners::String, txOffset::Integer, entity::ConnectedSwapBlock, room::Maple.Room)
    midSprite = Ahorn.getSprite(midResource, "Gameplay")

    tileWidth = div(width, 8)
    tileHeight = div(height, 8)

    rects = getExtensionRectangles(room)
    rx, ry = Ahorn.position(entity)
    rect = Ahorn.Rectangle(rx, ry, width, height)

    if !(rect in rects)
        push!(rects, rect)
    end

    for i in 1:tileWidth, j in 1:tileHeight
        drawX, drawY = (i - 1) * 8, (j - 1) * 8

        closedLeft = !notAdjacent(entity, drawX - 8, drawY, rects)
        closedRight = !notAdjacent(entity, drawX + 8, drawY, rects)
        closedUp = !notAdjacent(entity, drawX, drawY - 8, rects)
        closedDown = !notAdjacent(entity, drawX, drawY + 8, rects)
        completelyClosed = closedLeft && closedRight && closedUp && closedDown

        if completelyClosed
            if notAdjacent(entity, drawX + 8, drawY - 8, rects)
                # up right
                Ahorn.drawImage(ctx, innerCorners, x + drawX, y + drawY, 8 + txOffset, 0, 8, 8)

            elseif notAdjacent(entity, drawX - 8, drawY - 8, rects)
                # up left
                Ahorn.drawImage(ctx, innerCorners, x + drawX, y + drawY, 0 + txOffset, 0, 8, 8)

            elseif notAdjacent(entity, drawX + 8, drawY + 8, rects)
                # down right
                Ahorn.drawImage(ctx, innerCorners, x + drawX, y + drawY, 8 + txOffset, 8, 8, 8)

            elseif notAdjacent(entity, drawX - 8, drawY + 8, rects)
                # down left
                Ahorn.drawImage(ctx, innerCorners, x + drawX, y + drawY, 0 + txOffset, 8, 8, 8)

            else
                # entirely surrounded, fill tile
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 8, 8, 8, 8)
            end
        else
            if closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 8, 0, 8, 8)

            elseif closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 8, 16, 8, 8)

            elseif closedLeft && !closedRight && closedUp && closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 16, 8, 8, 8)

            elseif !closedLeft && closedRight && closedUp && closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 0, 8, 8, 8)

            elseif closedLeft && !closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 16, 0, 8, 8)

            elseif !closedLeft && closedRight && !closedUp && closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 0, 0, 8, 8)

            elseif !closedLeft && closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 0, 16, 8, 8)

            elseif closedLeft && !closedRight && closedUp && !closedDown
                Ahorn.drawImage(ctx, block, x + drawX, y + drawY, 16, 16, 8, 8)
            end
        end
    end

    Ahorn.drawImage(ctx, midSprite, x + div(width - midSprite.width, 2), y + div(height - midSprite.height, 2))
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::ConnectedSwapBlock, room::Maple.Room)
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    frame, _, mid, innerCorners = getTextures(entity)
    customBlockTexture = String(get(entity.data, "customRedBlockTexture", ""))
    hasCustomTexture = customBlockTexture != ""
    txOffset = 0

    if hasCustomTexture
        frame = innerCorners = "objects/" * customBlockTexture
        txOffset = 24
    end

    renderSwapBlock(ctx, stopX, stopY, width, height, mid, frame, innerCorners, txOffset, entity, room)
    Ahorn.drawArrow(ctx, startX + width / 2, startY + height / 2, stopX + width / 2, stopY + height / 2, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ConnectedSwapBlock, room::Maple.Room)
    startX, startY = Int(entity.data["x"]), Int(entity.data["y"])
    stopX, stopY = Int.(entity.data["nodes"][1])

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    frame, trail, mid, innerCorners = getTextures(entity)
    customBlockTexture = String(get(entity.data, "customRedBlockTexture", ""))
    hasCustomTexture = customBlockTexture != ""
    txOffset = 0

    if hasCustomTexture
        frame = innerCorners = "objects/" * customBlockTexture
        txOffset = 24
    end

    renderTrail(ctx, min(startX, stopX), min(startY, stopY), abs(startX - stopX) + width, abs(startY - stopY) + height, trail)
    renderSwapBlock(ctx, startX, startY, width, height, mid, frame, innerCorners, txOffset, entity, room)
end

end


text:
"""
    size(A) -> Tuple
    size(A, dim) -> Union{Int,StaticInt}

Returns the size of each dimension of `A` or along dimension `dim` of `A`. If the size of
any axes are known at compile time, these should be returned as `Static` numbers. Otherwise,
`ArrayInterface.size(A)` is identical to `Base.size(A)`

```julia
julia> using StaticArrays, ArrayInterface

julia> A = @SMatrix rand(3,4);

julia> ArrayInterface.size(A)
(static(3), static(4))
```
"""
@inline function size(a::A) where {A}
    if is_forwarding_wrapper(A)
        return size(parent(a))
    else
        return _maybe_size(Base.IteratorSize(A), a)
    end
end
size(a::Base.Broadcast.Broadcasted) = map(length, axes(a))

_maybe_size(::Base.HasShape{N}, a::A) where {N,A} = map(length, axes(a))
_maybe_size(::Base.HasLength, a::A) where {A} = (length(a),)
@inline size(x::SubArray) = flatten_tuples(map(size, x.indices))
@inline size(B::VecAdjTrans) = (One(), length(parent(B)))
@inline size(B::MatAdjTrans) = permute(size(parent(B)), to_parent_dims(B))
@inline function size(B::PermutedDimsArray{T,N,I1}) where {T,N,I1}
    permute(size(parent(B)), static(I1))
end
function size(a::ReinterpretArray{T,N,S,A,IsReshaped}) where {T,N,S,A,IsReshaped}
    psize = size(parent(a))
    if IsReshaped
        if sizeof(S) === sizeof(T)
            return psize
        elseif sizeof(S) > sizeof(T)
            return (static(div(sizeof(S), sizeof(T))), psize...)
        else
            return tail(psize)
        end
    else
        return (div(first(psize) * static(sizeof(S)), static(sizeof(T))), tail(psize)...,)
    end
end
size(A::ReshapedArray) = Base.size(A)
size(A::AbstractRange) = (length(A),)
size(x::Base.Generator) = size(getfield(x, :iter))
size(x::Iterators.Reverse) = size(getfield(x, :itr))
size(x::Iterators.Enumerate) = size(getfield(x, :itr))
size(x::Iterators.Accumulate) = size(getfield(x, :itr))
size(x::Iterators.Pairs) = size(getfield(x, :itr))
@inline function size(x::Iterators.ProductIterator)
    eachop(_sub_size, ntuple(static, StaticInt(ndims(x))), getfield(x, :iterators))
end
_sub_size(x::Tuple, ::StaticInt{dim}) where {dim} = length(getfield(x, dim))

size(a, dim) = size(a, to_dims(a, dim))
size(a::Array, dim::CanonicalInt) = Base.arraysize(a, convert(Int, dim))
function size(a::A, dim::CanonicalInt) where {A}
    if is_forwarding_wrapper(A)
        return size(parent(a), dim)
    else
        len = known_size(A, dim)
        if len === nothing
            return Int(length(axes(a, dim)))
        else
            return StaticInt(len)
        end
    end
end
size(x::Iterators.Zip) = Static.reduce_tup(promote_shape, map(size, getfield(x, :is)))

"""
    known_size(::Type{T}) -> Tuple
    known_size(::Type{T}, dim) -> Union{Int,Nothing}

Returns the size of each dimension of `A` or along dimension `dim` of `A` that is known at
compile time. If a dimension does not have a known size along a dimension then `nothing` is
returned in its position.
"""
known_size(x) = known_size(typeof(x))
@inline known_size(@nospecialize T::Type{<:Number}) = ()
@inline known_size(@nospecialize T::Type{<:VecAdjTrans}) = (1, known_length(parent_type(T)))
@inline function known_size(@nospecialize T::Type{<:MatAdjTrans})
    s1, s2 = known_size(parent_type(T))
    (s2, s1)
end
function known_size(@nospecialize T::Type{<:Diagonal})
    s = known_length(parent_type(T))
    (s, s)
end
known_size(@nospecialize T::Type{<:Union{Symmetric,Hermitian}}) = known_size(parent_type(T))
@inline function known_size(::Type{<:Base.ReinterpretArray{T,N,S,A,IsReshaped}}) where {T,N,S,A,IsReshaped}
    psize = known_size(A)
    if IsReshaped
        if sizeof(S) > sizeof(T)
            return (div(sizeof(S), sizeof(T)), psize...)
        elseif sizeof(S) < sizeof(T)
            return Base.tail(psize)
        else
            return psize
        end
    else
        if Base.issingletontype(T) || first(psize) === nothing
            return psize
        else
            return (div(first(psize) * sizeof(S), sizeof(T)), Base.tail(psize)...)
        end
    end
end

@inline function known_size(::Type{<:PermutedDimsArray{<:Any,N,I1,<:Any,P}}) where {N,I1,P}
    sz = known_size(P)
    ntuple(i -> getfield(sz, getfield(I1, i)), Val{N}())
end

@inline function known_size(::Type{T}) where {T}
    if is_forwarding_wrapper(T)
        return known_size(parent_type(T))
    else
        return _maybe_known_size(Base.IteratorSize(T), T)
    end
end
function _maybe_known_size(::Base.HasShape{N}, ::Type{T}) where {N,T}
    eachop(_known_size, ntuple(static, StaticInt(N)), axes_types(T))
end
_maybe_known_size(::Base.IteratorSize, ::Type{T}) where {T} = (known_length(T),)
known_size(::Type{<:Base.IdentityUnitRange{I}}) where {I} = known_size(I)
known_size(::Type{<:Base.Generator{I}}) where {I} = known_size(I)
known_size(::Type{<:Iterators.Reverse{I}}) where {I} = known_size(I)
known_size(::Type{<:Iterators.Enumerate{I}}) where {I} = known_size(I)
known_size(::Type{<:Iterators.Accumulate{<:Any,I}}) where {I} = known_size(I)
known_size(::Type{<:Iterators.Pairs{<:Any,<:Any,I}}) where {I} = known_size(I)
@inline function known_size(::Type{<:Iterators.ProductIterator{T}}) where {T}
    ntuple(i -> known_length(T.parameters[i]), Val(known_length(T)))
end
@inline function known_size(@nospecialize T::Type{<:AbstractRange})
    if is_forwarding_wrapper(T)
        return known_size(parent_type(T))
    else
        return (_range_length(known_first(T), known_step(T), known_last(T)),)
    end
end
@inline function known_size(@nospecialize T::Type{<:Union{LinearIndices,CartesianIndices}})
    I = fieldtype(T, :indices)
    ntuple(i -> known_length(I.parameters[i]), Val(ndims(T)))
end
@inline function known_size(@nospecialize T::Type{<:SubArray})
    _known_sub_sizes(fieldtype(T, :indices))
end
@inline function _known_sub_sizes(T::Type{<:Tuple})
    flatten_tuples(ntuple(i -> known_size(T.parameters[i]), Val(known_length(T))))
end

# 1. `Zip` doesn't check that its collections are compatible (same size) at construction,
#   but we assume as much b/c otherwise it will error while iterating. So we promote to the
#   known size if matching a `Nothing` and `Int` size.
# 2. `promote_shape(::Tuple{Vararg{CanonicalInt}}, ::Tuple{Vararg{CanonicalInt}})` promotes
#   trailing dimensions (which must be of size 1), to `static(1)`. We want to stick to
#   `Nothing` and `Int` types, so we do one last pass to ensure everything is dynamic
@inline function known_size(::Type{<:Iterators.Zip{T}}) where {T}
    dynamic(reduce_tup(Static._promote_shape, eachop(_unzip_size, ntuple(static, StaticInt(known_length(T))), T)))
end
_unzip_size(::Type{T}, n::StaticInt{N}) where {T,N} = known_size(field_type(T, n))
_known_size(::Type{T}, dim::StaticInt) where {T} = known_length(field_type(T, dim))
@inline known_size(x, dim) = known_size(typeof(x), dim)
@inline known_size(::Type{T}, dim) where {T} = known_size(T, to_dims(T, dim))
@inline function known_size(::Type{T}, dim::CanonicalInt) where {T}
    if ndims(T) < dim
        return 1
    else
        return known_size(T)[dim]
    end
end

"""
    length(A) -> Union{Int,StaticInt}

Returns the length of `A`.  If the length is known at compile time, it is
returned as `Static` number.  Otherwise, `ArrayInterface.length(A)` is identical
to `Base.length(A)`.

```julia
julia> using StaticArrays, ArrayInterface

julia> A = @SMatrix rand(3,4);

julia> ArrayInterface.length(A)
static(12)
```
"""
@inline length(a::UnitRange{T}) where {T} = last(a) - first(a) + oneunit(T)
@inline length(x) = Static.maybe_static(known_length, Base.length, x)

# Alias to to-be-depreciated internal function
const static_length = length

"""
    known_length(::Type{T}) -> Union{Int,Nothing}

If `length` of an instance of type `T` is known at compile time, return it.
Otherwise, return `nothing`.
"""
known_length(x) = known_length(typeof(x))
known_length(::Type{<:NamedTuple{L}}) where {L} = length(L)
known_length(::Type{T}) where {T<:Slice} = known_length(parent_type(T))
known_length(::Type{<:Tuple{Vararg{Any,N}}}) where {N} = N
known_length(::Type{<:Number}) = 1
known_length(::Type{<:AbstractCartesianIndex{N}}) where {N} = N
known_length(::Type{T}) where {T} = _maybe_known_length(Base.IteratorSize(T), T)
function known_length(::Type{<:Iterators.Flatten{I}}) where {I}
  _prod_or_nothing((known_length(I),known_length(eltype(I))))
end

_prod_or_nothing(x::Tuple{Vararg{Int}}) = prod(x)
_prod_or_nothing(_) = nothing

_maybe_known_length(::Base.HasShape, ::Type{T}) where {T} = _prod_or_nothing(known_size(T))
_maybe_known_length(::Base.IteratorSize, ::Type) = nothing


text:
export set_mode, get_mode, set_pull_up_down

exceptions = true
"""
A class to store socket and lock.
"""
mutable struct SockLock
    s::TCPSocket
    l::ReentrantLock
end

"""
A class to store pulse information.

gpio_on: the GPIO to switch on at the start of the pulse.
gpio_off: the GPIO to switch off at the start of the pulse.
delay: the delay in microseconds before the next pulse.
"""
mutable struct Pulse
    gpio_on::Int
    gpio_off::Int
    delay::Int
end

"""
Returns a text description of a pigpio error.

errnum:= <0, the error number

...
print(pigpio.error_text(-5))
level not 0-1
...
"""
function error_text(errnum)
    for e in _errors
        if e[1] == errnum
           return e[2]
       end
    end
    return "unknown error ($errnum)"
end

"""
Returns the microsecond difference between two ticks.

t1:= the earlier tick
t2:= the later tick

...
print(pigpio.tickDiff(4294967272, 12))
36
...
"""
function tickDiff(t1, t2)
   tDiff = t2 - t1
   if tDiff < 0
      tDiff += (1 << 32)
   end
   return tDiff
end

"""
Converts a 32 bit unsigned number to signed.  If the number
is negative it indicates an error.  On error a pigpio
exception will be raised if exceptions is true.
"""
function _u2i(x::UInt32)
   v = convert(Int32, x)
   if v < 0
      if exceptions
          error(error_text(v))
     end
   end
   return v
end

struct InMsg
    cmd::Cuint # a bits type
    p1::Cuint # an array of bits types
    p2::Cuint # a string with a fixed number of bytes
    d::Cuint
end

struct OutMsg
    dummy::Array{UInt8,1} # a bits type
    res::Cuint # an array of bits types
end

"""
Runs a pigpio socket command.

sl:= command socket and lock.
cmd:= the command to be executed.
p1:= command parameter 1 (if applicable).
 p2:=  command parameter 2 (if applicable).
"""
function _pigpio_command(sl::SockLock, cmd::Integer, p1::Integer, p2::Integer, rl=true)
    lock(sl.l)
    Base.write(sl.s, UInt32.([cmd, p1, p2, 0]))
    out = IOBuffer(Base.read(sl.s, 16))
    msg = reinterpret(Cuint, take!(out))[4]
    if rl
        unlock(sl.l)
    end
   return msg
end

"""
Runs an extended pigpio socket command.

    sl:= command socket and lock.
   cmd:= the command to be executed.
    p1:= command parameter 1 (if applicable).
    p2:= command parameter 2 (if applicable).
    p3:= total size in bytes of following extents
extents:= additional data blocks
"""
function _pigpio_command_ext(sl, cmd, p1, p2, p3, extents, rl=true)
    ext = IOBuffer()
    Base.write(ext, Array(reinterpret(UInt8, [cmd, p1, p2, p3])))
    for x in extents
       Base.write(ext, string(x))
    end
    lock(sl.l)
    Base.write(sl.s, ext)
    msg = reinterpret(Cuint, sl.s)[4]
    if rl
         unlock(sl.l)
    end
    return msg
end

"""An ADT class to hold callback information

   gpio:= Broadcom GPIO number.
   edge:= EITHER_EDGE, RISING_EDGE, or FALLING_EDGE.
   func:= a user function taking three arguments (GPIO, level, tick).
"""
mutable struct Callback_ADT
    gpio::Int
    edge::Int
    func::Function
    bit::Int

    Callback_ADT(gpio, edge, func) = new(gpio, edge, func, 1<<gpio)
end

"""A class to encapsulate pigpio notification callbacks."""
mutable struct CallbackThread #(threading.Thread)
    control::SockLock
    sl::SockLock
    go::Bool
    daemon::Bool
    monitor::Int
    handle::Cuint
    callbacks::Array{Any, 1}
end

"""Initialises notifications."""
function CallbackThread(control, host, port)
    socket = connect(host, port)
    sl = SockLock(socket, ReentrantLock())
    self = CallbackThread(control, sl, false, true, 0, 0, Any[])
    self.handle = _pigpio_command(sl, _PI_CMD_NOIB, 0, 0)
    self.go = true
    return self
    #self.start()  #TODO
end

"""Stops notifications."""
function stop(self::CallbackThread)
    if self.go
        self.go = false
        Base.write(self.sl.s, _PI_CMD_NC, self.handle, 0, 0)
    end
end

"""Adds a callback to the notification thread."""
function append(self::CallbackThread, callb)
    push!(self.callbacks, callb)
    self.monitor = self.monitor | callb.bit
    _pigpio_command(self.control, _PI_CMD_NB, self.handle, self.monitor)
end

"""Removes a callback from the notification thread."""
function remove(self::CallbackThread, callb)
    if callb in self.callbacks
        self.callbacks.remove(callb)
        newMonitor = 0
        for c in self.callbacks
            newMonitor |= c.bit
        end

        if newMonitor != self.monitor
            self.monitor = newMonitor
            _pigpio_command(
                self.control, _PI_CMD_NB, self.handle, self.monitor)
        end
    end
end

struct CallbMsg
    seq::Cushort
    flags::Cushort
    tick::Cuint
    level::Cuint
end

"""Runs the notification thread."""
function run(self::CallbackThread)
    lastLevel = _pigpio_command(self.control,  _PI_CMD_BR1, 0, 0)
    MSG_SIZ = 12
    while self.go
        buf = readbytes!(self.sl.s, MSG_SIZ, all=true)
        if self.go
            msg = reinterpret(CallbMsg, buf)
            seq = msg.seq
            seq, flags, tick, level = (msg.seq, msg.flags, msg.tick, msg.level)
            if flags == 0
                changed = level ^ lastLevel
                lastLevel = level
                for  cb in self.callbacks
                    if cb.bit && changed
                        newLevel = 0
                    elseif cb.bit & level
                        newLevel = 1
                    end
                    if (cb.edge ^ newLevel)
                        cb.func(cb.gpio, newLevel, tick)
                    end
                end
            else
                if flags & NTFY_FLAGS_WDOG
                    gpio = flags & NTFY_FLAGS_GPIO
                    for cb in self.callbacks
                        if cb.gpio == gpio
                            cb.func(cb.gpio, TIMEOUT, tick)
                        end
                    end
                end
            end
        end
    end
    close(self.sl.s)
end

"""A class to provide GPIO level change callbacks."""
mutable struct Callback
    notify::Array{Callback_ADT, 1}
    count::Int
    reset::Bool
    callb::Callback_ADT
end

function Callback(notify::CallbackThread, user_gpio::Int, edge::Int=RISING_EDGE, func=nothing)
    # Pi.notify is passed here which is a CallbackThread, and so the array CallbackThread.callbacks must be passed to the Callback constructor since it is type Array{Callback_ADT, 1}
    #  self = Callback(notify.callbacks, 0, false, nothing)
    if func === nothing
        func = _tally
    end
    self = Callback(notify.callbacks, 0, false, Callback_ADT(user_gpio, edge, func))
    #  self.callb = Callback_ADT(user_gpio, edge, func)
    push!(self.notify, self.callb)
    return self
end

"""Cancels a callback by removing it from the notification thread."""
function cancel(self)
    filter(x->x!=self.callb, self.notify)
end

"""Increment the callback called count."""
function _tally(self::Callback, user_gpio, level, tick)
    if self.reset
        self.reset = false
        self.count = 0
    end
    self.count += 1
end

"""
Provides a count of how many times the default tally
callback has triggered.

The count will be zero if the user has supplied their own
callback function.
"""
function tally(self::Callback)
    return self.count
end

"""
Resets the tally count to zero.
"""
function reset_tally(self::Callback)
    self.reset = true
    self.count = 0
end

"""Encapsulates waiting for GPIO edges."""
mutable struct WaitForEdge
    notify
    callb::Function
    trigger
    start
end

"""Initialises a wait_for_edge."""
function WaitForEdge( notify, gpio::Int, edge, timeout)
    callb = Callback_ADT(gpio, edge, func)
    self = WaitForEdge(notify, callb, false, time())
    push!(self.notify, self.callb)
    while (self.trigger == false) && ((time()-self.start) < timeout)
      time.sleep(0.05)
    end
    self.notify.remove(self.callb)
end

"""Sets wait_for_edge triggered."""
function func(self::WaitForEdge, gpio, level, tick)
    self.trigger = true
end

mutable struct Pi
    host::String
    port::Int
    connected::Bool
    sl::SockLock
    notify::CallbackThread
end

"""Returns count bytes from the command socket."""
function rxbuf(self::Pi, count)
    ext = readbytes!(self.sl.s, count, all)
    return ext
end

"""
Sets the GPIO mode.

gpio:= 0-53.
mode:= INPUT, OUTPUT, ALT0, ALT1, ALT2, ALT3, ALT4, ALT5.

...
set_mode(pi,  4, pigpio.INPUT)  # GPIO  4 as input
set_mode(pi, 17, pigpio.OUTPUT) # GPIO 17 as output
set_mode(pi, 24, pigpio.ALT2)   # GPIO 24 as ALT2
...
"""
function set_mode(self::Pi, gpio, mode)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_MODES, gpio, mode))
end

"""
Returns the GPIO mode.

gpio:= 0-53.

Returns a value as follows

. .
0 = INPUT
1 = OUTPUT
2 = ALT5
3 = ALT4
4 = ALT0
5 = ALT1
6 = ALT2
7 = ALT3
. .

...
print(get_mode(pi, 0))
4
...
"""
function get_mode(self::Pi, gpio)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_MODEG, gpio, 0))
end

"""
Sets or clears the internal GPIO pull-up/down resistor.

gpio:= 0-53.
pud:= PUD_UP, PUD_DOWN, PUD_OFF.

...
set_pull_up_down(pi, 17, pigpio.PUD_OFF)
set_pull_up_down(pi, 23, pigpio.PUD_UP)
set_pull_up_down(pi, 24, pigpio.PUD_DOWN)
...
"""
function set_pull_up_down(self::Pi, gpio, pud)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_PUD, gpio, pud))
end

"""
Returns the GPIO level.

gpio:= 0-53.

...
set_mode(pi, 23, pigpio.INPUT)

set_pull_up_down(pi, 23, pigpio.PUD_DOWN)
print(read(pi, 23))
0

set_pull_up_down(pi, 23, pigpio.PUD_UP)
print(read(pi, 23))
1
...
"""
function read(self::Pi, gpio)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_READ, gpio, 0))
end

"""
Sets the GPIO level.

GPIO:= 0-53.
level:= 0, 1.

If PWM or servo pulses are active on the GPIO they are
switched off.

...
set_mode(pi, 17, pigpio.OUTPUT)

write(pi, 17,0)
print(read(pi, 17))
0

write(pi, 17,1)
print(read(pi, 17))
1
...
"""
function write(self::Pi, gpio, level)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_WRITE, gpio, level))
end

"""
Starts (non-zero dutycycle) or stops (0) PWM pulses on the GPIO.

user_gpio:= 0-31.
dutycycle:= 0-range (range defaults to 255).

The [*set_PWM_range*] function can change the default range of 255.

...
set_PWM_dutycycle(pi, 4,   0) # PWM off
set_PWM_dutycycle(pi, 4,  64) # PWM 1/4 on
set_PWM_dutycycle(pi, 4, 128) # PWM 1/2 on
set_PWM_dutycycle(pi, 4, 192) # PWM 3/4 on
set_PWM_dutycycle(pi, 4, 255) # PWM full on
...
"""
function set_PWM_dutycycle(self::Pi, user_gpio, dutycycle)
    return _u2i(_pigpio_command(
        self.sl, _PI_CMD_PWM, user_gpio, Int(dutycycle)))
end

"""
Returns the PWM dutycycle being used on the GPIO.

user_gpio:= 0-31.

Returns the PWM dutycycle.

For normal PWM the dutycycle will be out of the defined range
for the GPIO (see [*get_PWM_range*]).

If a hardware clock is active on the GPIO the reported
dutycycle will be 500000 (500k) out of 1000000 (1M).

If hardware PWM is active on the GPIO the reported dutycycle
will be out of a 1000000 (1M).

...
set_PWM_dutycycle(pi, 4, 25)
print(get_PWM_dutycycle(pi, 4))
25

set_PWM_dutycycle(pi, 4, 203)
print(get_PWM_dutycycle(pi, 4))
203
...
"""
function get_PWM_dutycycle(self::Pi, user_gpio)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_GDC, user_gpio, 0))
end

"""
Sets the range of PWM values to be used on the GPIO.

user_gpio:= 0-31.
 range_:= 25-40000.

...
set_PWM_range(pi, 9, 100)  # now  25 1/4,   50 1/2,   75 3/4 on
set_PWM_range(pi, 9, 500)  # now 125 1/4,  250 1/2,  375 3/4 on
set_PWM_range(pi, 9, 3000) # now 750 1/4, 1500 1/2, 2250 3/4 on
...
"""
function set_PWM_range(self::Pi, user_gpio, range_)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_PRS, user_gpio, range_))
end

"""
Returns the range of PWM values being used on the GPIO.

user_gpio:= 0-31.

If a hardware clock or hardware PWM is active on the GPIO
the reported range will be 1000000 (1M).

...
set_PWM_range(pi, 9, 500)
print(get_PWM_range(pi, 9))
500
...
"""
function get_PWM_range(self::Pi, user_gpio)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_PRG, user_gpio, 0))
end

"""
Returns the real (underlying) range of PWM values being
used on the GPIO.

user_gpio:= 0-31.

If a hardware clock is active on the GPIO the reported
real range will be 1000000 (1M).

If hardware PWM is active on the GPIO the reported real range
will be approximately 250M divided by the set PWM frequency.

...
set_PWM_frequency(pi, 4, 800)
print(get_PWM_real_range(pi, 4))
250
...
"""
function get_PWM_real_range(self::Pi, user_gpio)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_PRRG, user_gpio, 0))
end

"""
Sets the frequency (in Hz) of the PWM to be used on the GPIO.

user_gpio:= 0-31.
frequency:= >=0 Hz

Returns the numerically closest frequency if OK, otherwise
PI_BAD_USER_GPIO or PI_NOT_PERMITTED.

If PWM is currently active on the GPIO it will be switched
off and then back on at the new frequency.

Each GPIO can be independently set to one of 18 different
PWM frequencies.

The selectable frequencies depend upon the sample rate which
may be 1, 2, 4, 5, 8, or 10 microseconds (default 5).  The
sample rate is set when the pigpio daemon is started.

The frequencies for each sample rate are

. .
                     Hertz

     1: 40000 20000 10000 8000 5000 4000 2500 2000 1600
         1250  1000   800  500  400  250  200  100   50

     2: 20000 10000  5000 4000 2500 2000 1250 1000  800
          625   500   400  250  200  125  100   50   25

     4: 10000  5000  2500 2000 1250 1000  625  500  400
          313   250   200  125  100   63   50   25   13
sample
rate
(us)  5:  8000  4000  2000 1600 1000  800  500  400  320
          250   200   160  100   80   50   40   20   10

     8:  5000  2500  1250 1000  625  500  313  250  200
          156   125   100   63   50   31   25   13    6

    10:  4000  2000  1000  800  500  400  250  200  160
          125   100    80   50   40   25   20   10    5
. .

...
set_PWM_frequency(pi, 4,0)
print(get_PWM_frequency(pi, 4))
10

set_PWM_frequency(pi, 4,100000)
print(get_PWM_frequency(pi, 4))
8000
...
"""
function set_PWM_frequency(self::Pi, user_gpio, frequency)
    return _u2i(
        _pigpio_command(self.sl, _PI_CMD_PFS, user_gpio, frequency))
end

"""
Returns the frequency of PWM being used on the GPIO.

user_gpio:= 0-31.

Returns the frequency (in Hz) used for the GPIO.

For normal PWM the frequency will be that defined for the GPIO
by [*set_PWM_frequency*].

If a hardware clock is active on the GPIO the reported frequency
will be that set by [*hardware_clock*].

If hardware PWM is active on the GPIO the reported frequency
will be that set by [*hardware_PWM*].

...
set_PWM_frequency(pi, 4,0)
print(get_PWM_frequency(pi, 4))
10

set_PWM_frequency(pi, 4, 800)
print(get_PWM_frequency(pi, 4))
800
...
"""
function get_PWM_frequency(self::Pi, user_gpio)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_PFG, user_gpio, 0))
end

"""
Starts (500-2500) or stops (0) servo pulses on the GPIO.

user_gpio:= 0-31.
pulsewidth:= 0 (off),
           500 (most anti-clockwise) - 2500 (most clockwise).

The selected pulsewidth will continue to be transmitted until
changed by a subsequent call to set_servo_pulsewidth.

The pulsewidths supported by servos varies and should probably
be determined by experiment. A value of 1500 should always be
safe and represents the mid-point of rotation.

You can DAMAGE a servo if you command it to move beyond its
limits.

...
set_servo_pulsewidth(pi, 17, 0)    # off
set_servo_pulsewidth(pi, 17, 1000) # safe anti-clockwise
set_servo_pulsewidth(pi, 17, 1500) # centre
set_servo_pulsewidth(pi, 17, 2000) # safe clockwise
...
"""
function set_servo_pulsewidth(self::Pi, user_gpio, pulsewidth)
    return _u2i(_pigpio_command(
        self.sl, _PI_CMD_SERVO, user_gpio, Int(pulsewidth)))
end

"""
Returns the servo pulsewidth being used on the GPIO.

user_gpio:= 0-31.

Returns the servo pulsewidth.

...
set_servo_pulsewidth(pi, 4, 525)
print(get_servo_pulsewidth(pi, 4))
525

set_servo_pulsewidth(pi, 4, 2130)
print(get_servo_pulsewidth(pi, 4))
2130
...
"""
function get_servo_pulsewidth(self::Pi, user_gpio)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_GPW, user_gpio, 0))
end

"""
Returns a notification handle (>=0).

A notification is a method for being notified of GPIO state
changes via a pipe.

Pipes are only accessible from the local machine so this
function serves no purpose if you are using Python from a
remote machine.  The in-built (socket) notifications
provided by [*callback*] should be used instead.

Notifications for handle x will be available at the pipe
named /dev/pigpiox (where x is the handle number).

E.g. if the function returns 15 then the notifications must be
read from /dev/pigpio15.

Notifications have the following structure.

. .
I seqno
I flags
I tick
I level
. .

seqno: starts at 0 each time the handle is opened and then
increments by one for each report.

flags: two flags are defined, PI_NTFY_FLAGS_WDOG and
PI_NTFY_FLAGS_ALIVE.  If bit 5 is set (PI_NTFY_FLAGS_WDOG)
then bits 0-4 of the flags indicate a GPIO which has had a
watchdog timeout; if bit 6 is set (PI_NTFY_FLAGS_ALIVE) this
indicates a keep alive signal on the pipe/socket and is sent
once a minute in the absence of other notification activity.

tick: the number of microseconds since system boot.  It wraps
around after 1h12m.

level: indicates the level of each GPIO.  If bit 1<<x is set
then GPIO x is high.

...
h = notify_open(pi, )
if h >= 0
    notify_begin(pi, h, 1234)
...
"""
function notify_open(self::Pi)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_NO, 0, 0))
end

"""
Starts notifications on a handle.

handle:= >=0 (as returned by a prior call to [*notify_open*])
bits:= a 32 bit mask indicating the GPIO to be notified.

The notification sends state changes for each GPIO whose
corresponding bit in bits is set.

The following code starts notifications for GPIO 1, 4,
6, 7, and 10 (1234 = 0x04D2 = 0b0000010011010010).

...
h = notify_open(pi, )
if h >= 0
    notify_begin(pi, h, 1234)
...
"""
function notify_begin(self::Pi, handle, bits)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_NB, handle, bits))
end

"""
Pauses notifications on a handle.

handle:= >=0 (as returned by a prior call to [*notify_open*])

Notifications for the handle are suspended until
[*notify_begin*] is called again.

...
h = notify_open(pi, )
if h >= 0
    notify_begin(pi, h, 1234)
    ...
    notify_pause(pi, h)
    ...
    notify_begin(pi, h, 1234)
    ...
...
"""
function notify_pause(self::Pi, handle)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_NB, handle, 0))
end

"""
Stops notifications on a handle and releases the handle for reuse.

handle:= >=0 (as returned by a prior call to [*notify_open*])

...
h = notify_open(pi, )
if h >= 0
    notify_begin(pi, h, 1234)
    ...
    notify_close(pi, h)
    ...
...
"""
function notify_close(self::Pi, handle)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_NC, handle, 0))
end

"""
Sets a watchdog timeout for a GPIO.

 user_gpio:= 0-31.
wdog_timeout:= 0-60000.

The watchdog is nominally in milliseconds.

Only one watchdog may be registered per GPIO.

The watchdog may be cancelled by setting timeout to 0.

If no level change has been detected for the GPIO for timeout
milliseconds any notification for the GPIO has a report written
to the fifo with the flags set to indicate a watchdog timeout.

The callback class interprets the flags and will
call registered callbacks for the GPIO with level TIMEOUT.

...
set_watchdog(pi, 23, 1000) # 1000 ms watchdog on GPIO 23
set_watchdog(pi, 23, 0)    # cancel watchdog on GPIO 23
...
"""
function set_watchdog(self::Pi, user_gpio, wdog_timeout)
    return _u2i(_pigpio_command(
        self.sl, _PI_CMD_WDOG, user_gpio, Int(wdog_timeout)))
end

"""
Returns the levels of the bank 1 GPIO (GPIO 0-31).

The returned 32 bit integer has a bit set if the corresponding
GPIO is high.  GPIO n has bit value (1<<n).

...
print(bin(read_bank_1(pi, )))
0b10010100000011100100001001111
...
"""
function read_bank_1(self::Pi)
	return _pigpio_command(self.sl, _PI_CMD_BR1, 0, 0)
end

"""
Returns the levels of the bank 2 GPIO (GPIO 32-53).

The returned 32 bit integer has a bit set if the corresponding
GPIO is high.  GPIO n has bit value (1<<(n-32)).

...
print(bin(read_bank_2(pi, )))
0b1111110000000000000000
...
"""
function read_bank_2(self::Pi)
    return _pigpio_command(self.sl, _PI_CMD_BR2, 0, 0)
end

"""
Clears GPIO 0-31 if the corresponding bit in bits is set.

bits:= a 32 bit mask with 1 set if the corresponding GPIO is
 to be cleared.

A returned status of PI_SOME_PERMITTED indicates that the user
is not allowed to write to one or more of the GPIO.

...
clear_bank_1(int(pi, "111110010000",2))
...
"""
function clear_bank_1(self::Pi, bits)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_BC1, bits, 0))
end

"""
Clears GPIO 32-53 if the corresponding bit (0-21) in bits is set.

bits:= a 32 bit mask with 1 set if the corresponding GPIO is
to be cleared.

A returned status of PI_SOME_PERMITTED indicates that the user
is not allowed to write to one or more of the GPIO.

...
clear_bank_2(pi, 0x1010)
...
"""
function clear_bank_2(self::Pi, bits)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_BC2, bits, 0))
end

"""
Sets GPIO 0-31 if the corresponding bit in bits is set.

bits:= a 32 bit mask with 1 set if the corresponding GPIO is
 to be set.

A returned status of PI_SOME_PERMITTED indicates that the user
is not allowed to write to one or more of the GPIO.

...
set_bank_1(int(pi, "111110010000",2))
...
"""
function set_bank_1(self::Pi, bits)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_BS1, bits, 0))
end

"""
Sets GPIO 32-53 if the corresponding bit (0-21) in bits is set.

bits:= a 32 bit mask with 1 set if the corresponding GPIO is
 to be set.

A returned status of PI_SOME_PERMITTED indicates that the user
is not allowed to write to one or more of the GPIO.

...
set_bank_2(pi, 0x303)
...
"""
function set_bank_2(self::Pi, bits)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_BS2, bits, 0))
end

"""
Starts a hardware clock on a GPIO at the specified frequency.
Frequencies above 30MHz are unlikely to work.

gpio:= see description
clkfreq:= 0 (off) or 4689-250000000 (250M)

Returns 0 if OK, otherwise PI_NOT_PERMITTED, PI_BAD_GPIO,
PI_NOT_HCLK_GPIO, PI_BAD_HCLK_FREQ,or PI_BAD_HCLK_PASS.

The same clock is available on multiple GPIO.  The latest
frequency setting will be used by all GPIO which share a clock.

The GPIO must be one of the following.

. .
4   clock 0  All models
5   clock 1  All models but A and B (reserved for system use)
6   clock 2  All models but A and B
20  clock 0  All models but A and B
21  clock 1  All models but A and Rev.2 B (reserved for system use)

32  clock 0  Compute module only
34  clock 0  Compute module only
42  clock 1  Compute module only (reserved for system use)
43  clock 2  Compute module only
44  clock 1  Compute module only (reserved for system use)
. .

Access to clock 1 is protected by a password as its use will
likely crash the Pi.  The password is given by or'ing 0x5A000000
with the GPIO number.

...
hardware_clock(pi, 4, 5000) # 5 KHz clock on GPIO 4

hardware_clock(pi, 4, 40000000) # 40 MHz clock on GPIO 4
...
"""
function hardware_clock(self::Pi, gpio, clkfreq)
    return _u2i(_pigpio_command(self.sl, _PI_CMD_HC, gpio, clkfreq))
end

"""
Starts hardware PWM on a GPIO at the specified frequency
and dutycycle. Frequencies above 30MHz are unlikely to work.

NOTE: Any waveform started by [*wave_send_once*],
[*wave_send_repeat*], or [*wave_chain*] will be cancelled.

This function is only valid if the pigpio main clock is PCM.
The main clock defaults to PCM but may be overridden when the
pigpio daemon is started (option -t).

gpio:= see descripton
PWMfreq:= 0 (off) or 1-125000000 (125M).
PWMduty:= 0 (off) to 1000000 (1M)(fully on).

Returns 0 if OK, otherwise PI_NOT_PERMITTED, PI_BAD_GPIO,
PI_NOT_HPWM_GPIO, PI_BAD_HPWM_DUTY, PI_BAD_HPWM_FREQ.

The same PWM channel is available on multiple GPIO.
The latest frequency and dutycycle setting will be used
by all GPIO which share a PWM channel.

The GPIO must be one of the following.

. .
12  PWM channel 0  All models but A and B
13  PWM channel 1  All models but A and B
18  PWM channel 0  All models
19  PWM channel 1  All models but A and B

40  PWM channel 0  Compute module only
41  PWM channel 1  Compute module only
45  PWM channel 1  Compute module only
52  PWM channel 0  Compute module only
53  PWM channel 1  Compute module only
. .

The actual number of steps beween off and fully on is the
integral part of 250 million divided by PWMfreq.

The actual frequency set is 250 million / steps.

There will only be a million steps for a PWMfreq of 250.
Lower frequencies will have more steps and higher
frequencies will have fewer steps.  PWMduty is
automatically scaled to take this into account.

...
hardware_PWM(pi, 18, 800, 250000) # 800Hz 25% dutycycle

hardware_PWM(pi, 18, 2000, 750000) # 2000Hz 75% dutycycle
...
"""

function hardware_PWM(self::Pi, gpio, PWMfreq, PWMduty)
# pigpio message format

# I p1 gpio
# I p2 PWMfreq
# I p3 4
## extension ##
# I PWMdutycycle
    extents = IOBuffer()
    extents = Base.write(extents, 10)
    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_HP, gpio, PWMfreq, 4, extents))
end

"""
Returns the current system tick.

Tick is the number of microseconds since system boot.  As an
unsigned 32 bit quantity tick wraps around approximately
every 71.6 minutes.

...
t1 = get_current_tick(pi, )
time.sleep(1)
t2 = get_current_tick(pi, )
...
"""
function get_current_tick(self::Pi)
    return _pigpio_command(self.sl, _PI_CMD_TICK, 0, 0)
end

"""
Returns the Pi's hardware revision number.

The hardware revision is the last few characters on the
Revision line of /proc/cpuinfo.

The revision number can be used to determine the assignment
of GPIO to pins (see [*gpio*]).

There are at least three types of board.

Type 1 boards have hardware revision numbers of 2 and 3.

Type 2 boards have hardware revision numbers of 4, 5, 6, and 15.

Type 3 boards have hardware revision numbers of 16 or greater.

If the hardware revision can not be found or is not a valid
hexadecimal number the function returns 0.

...
print(get_hardware_revision(pi, ))
2
...
"""
function get_hardware_revision(self::Pi)
    return _pigpio_command(self.sl, _PI_CMD_HWVER, 0, 0)
end

"""
Returns the pigpio software version.

...
v = get_pigpio_version(pi, )
...
"""
function get_pigpio_version(self::Pi)
    return _pigpio_command(self.sl, _PI_CMD_PIGPV, 0, 0)
end

"""
Calls a pigpio function customised by the user.

arg1:= >=0, default 0.
arg2:= >=0, default 0.
argx:= extra arguments (each 0-255), default empty.

The returned value is an integer which by convention
should be >=0 for OK and <0 for error.

...
value = custom_1(pi, )

value = custom_1(pi, 23)

value = custom_1(pi, 0, 55)

value = custom_1(pi, 23, 56, [1, 5, 7])

value = custom_1(pi, 23, 56, b"hello")

value = custom_1(pi, 23, 56, "hello")
...
"""
function custom_1(self, arg1=0, arg2=0, argx=[])
    # I p1 arg1
    # I p2 arg2
    # I p3 len
    ## extension ##
    # s len argx bytes

    return _u2i(_pigpio_command_ext(
        self.sl, _PI_CMD_CF1, arg1, arg2, length(argx), [argx]))
end

"""
Calls a pigpio function customised by the user.

arg1:= >=0, default 0.
argx:= extra arguments (each 0-255), default empty.
retMax:= >=0, maximum number of bytes to return, default 8192.

The returned value is a tuple of the number of bytes
returned and a bytearray containing the bytes.  If
there was an error the number of bytes read will be
less than zero (and will contain the error code).

...
(count, data) = custom_2(pi, )

(count, data) = custom_2(pi, 23)

(count, data) = custom_2(pi, 23, [1, 5, 7])

(count, data) = custom_2(pi, 23, b"hello")

(count, data) = custom_2(pi, 23, "hello", 128)
...
"""
function custom_2(self, arg1=0, argx=[], retMax=8192)
    # I p1 arg1
    # I p2 retMax
    # I p3 len
    ## extension ##
    # s len argx bytes

    # Don't raise exception.  Must release lock.
    bytes = _u2i(_pigpio_command_ext(
    self.sl, _PI_CMD_CF2, arg1, retMax, length(argx), [argx], false))
    if bytes > 0
        data = rxbuf(self, bytes)
    else
        data = ""
    end
    unlock(self.sl.l)
    return bytes, data
end

"""
Calls a user supplied function (a callback) whenever the
specified GPIO edge is detected.

user_gpio:= 0-31.
edge:= EITHER_EDGE, RISING_EDGE (default), or FALLING_EDGE.
func:= user supplied callback function.

The user supplied callback receives three parameters, the GPIO,
the level, and the tick.

If a user callback is not specified a default tally callback is
provided which simply counts edges.  The count may be retrieved
by calling the tally function.  The count may be reset to zero
by calling the reset_tally function.

The callback may be cancelled by calling the cancel function.

A GPIO may have multiple callbacks (although I can't think of
a reason to do so).

...
end

function cbf(gpio, level, tick)
print(gpio, level, tick)

cb1 = callback(pi, 22, pigpio.EITHER_EDGE, cbf)

cb2 = callback(pi, 4, pigpio.EITHER_EDGE)

cb3 = callback(pi, 17)

print(cb3.tally())

cb3.reset_tally()

cb1.cancel() # To cancel callback cb1.
...
"""
function callback(self::Pi, user_gpio, edge=RISING_EDGE, func=nothing)
    return Callback(self.notify, user_gpio, edge, func)
end

"""
Wait for an edge event on a GPIO.

user_gpio:= 0-31.
  edge:= EITHER_EDGE, RISING_EDGE (default), or
         FALLING_EDGE.
wait_timeout:= >=0.0 (default 60.0).

The function returns when the edge is detected or after
the number of seconds specified by timeout has expired.

Do not use this function for precise timing purposes,
the edge is only checked 20 times a second. Whenever
you need to know the accurate time of GPIO events use
a [*callback*] function.

The function returns true if the edge is detected,
otherwise false.

...
if wait_for_edge(pi, 23)
print("Rising edge detected")
else
print("wait for edge timed out")

if wait_for_edge(pi, 23, pigpio.FALLING_EDGE, 5.0)
print("Falling edge detected")
else
print("wait for falling edge timed out")
...
"""
function wait_for_edge(self::Pi, user_gpio, edge=RISING_EDGE, wait_timeout=60.0)
    a = WaitForEdge(self.notify, user_gpio, edge, wait_timeout)
    return a.trigger
end

"""
Grants access to a Pi's GPIO.

host:= the host name of the Pi on which the pigpio daemon is
 running.  The default is localhost unless overridden by
 the PIGPIO_ADDR environment variable.

port:= the port number on which the pigpio daemon is listening.
 The default is 8888 unless overridden by the PIGPIO_PORT
 environment variable.  The pigpio daemon must have been
 started with the same port number.

This connects to the pigpio daemon and reserves resources
to be used for sending commands and receiving notifications.

An instance attribute [*connected*] may be used to check the
success of the connection.  If the connection is established
successfully [*connected*] will be true, otherwise false.

...
pi = pigpio.pi()              # use defaults
pi = pigpio.pi('mypi')       # specify host, default port
pi = pigpio.pi('mypi', 7777) # specify host and port

pi = pigpio.pi()             # exit script if no connection
if not pi.connected
exit()
...
"""
function Pi(; host = get(ENV, "PIGPIO_ADDR", ""), port = get(ENV, "PIGPIO_PORT", 8888))
    port = Int(port)
    if host == "" || host === nothing
        host = "localhost"
    end

    try
        sock = connect(host, port)
        ccall(:uv_tcp_nodelay, Cint, (Ptr{Cvoid}, Cuint), sock, 1) # Disable Nagle's Algorithm
        sl = SockLock(sock, ReentrantLock())
        notify = CallbackThread(sl, host, port)
        self = Pi(host, port, true, sl, notify)
        @info "Successfully connected!"
        return self
        #atexit.register(self.stop) #TODO
    catch error
        println("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
        println("Can't connect to pigpio at $host:$port")
        println("Did you start the pigpio daemon? E.g. sudo pigpiod\n")
        println("Did you specify the correct Pi host/port in the environment")
        println("variables PIGPIO_ADDR/PIGPIO_PORT?")
        println("E.g. export PIGPIO_ADDR=soft, export PIGPIO_PORT=8888\n")
        println("Did you specify the correct Pi host/port in the")
        println("Pi() function? E.g. Pi('soft', 8888))")
        println("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
        throw(error)
    end
end

"""Release pigpio resources.
...
stop(pi)
...
"""
function stop(self::Pi)
    self.connected = false
    stop(self.notify)
    close(self.sl.s)
end


text:
using HorizonSideRobots

function alldraw(r)
    y = 0
    z = 0
    y, z = goinSWcorner(r)
    draw(r)
    vozvrat(r,y,z)
end

function obhod(r,side)
    x = 0
    sign = true
    while isborder(r,HorizonSide(side)) && !isborder(r,HorizonSide((side + 3) % 4))
        move!(r,HorizonSide((side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 3) % 4))
        sign = false
    end
    if  !isborder(r,HorizonSide(side))
        move!(r,HorizonSide(side))
    end
    while !isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 1) % 4))
        move!(r,HorizonSide(side))
    end
    for i = 1:x
        move!(r,HorizonSide((side+1) % 4))
    end
    return sign
end

function obhod2(r,side,y)
    x = 0
    sign = true
    while isborder(r,HorizonSide(side)) && !isborder(r,HorizonSide((side + 3) % 4))
        move!(r,HorizonSide((side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 3) % 4))
        sign = false
    end
    if  !isborder(r,HorizonSide(side))
        move!(r,HorizonSide(side))
        y+=1
    end
    while !isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 1) % 4))
        move!(r,HorizonSide(side))
        y+= 1
    end
    for i = 1:x
        move!(r,HorizonSide((side+1) % 4))
    end
    return sign,y
end

function obhod3(r,side,y)
    x = 0
    sign = true
    while isborder(r,HorizonSide(side)) && !isborder(r,HorizonSide((side + 3) % 4))
        move!(r,HorizonSide((side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 3) % 4))
        sign = false
    end
    if  !isborder(r,HorizonSide(side))
        move!(r,HorizonSide(side))
        y-=1
    end
    while !isborder(r,HorizonSide(side)) && isborder(r,HorizonSide((side + 1) % 4))
        move!(r,HorizonSide(side))
        y-= 1
    end
    for i = 1:x
        move!(r,HorizonSide((side+1) % 4))
    end
    return sign,y
end

function goinSWcorner(r)
    sign = true
    y = 0
    z = 0
    while !isborder(r,Sud) || sign
        if !isborder(r,Sud)
            move!(r,Sud)
            y += 1
        end
        if isborder(r,Sud)
            sign,y =  obhod2(r,2,y)
        end
    end
    while !isborder(r,West)
        move!(r,West)
        z += 1
    end
    return y,z
end

function draw(r)
    sign2 = true
    while !isborder(r,Nord) || sign2
        sign1 = true
        while !isborder(r,Ost) || sign1
            putmarker!(r)
            if !isborder(r,Ost)
                move!(r,Ost)
            else
                sign1 = obhod(r,3)
            end
        end
        sign1 = true
        while !isborder(r,West) || sign1
            if !isborder(r,West)
                move!(r,West)
            else
                sign1 = obhod(r,1)
            end
        end
        if !isborder(r,Nord)
            move!(r,Nord)
        else
            sign2 = obhod(r,0)
        end
    end
end

function vozvrat(r,y,z)
    sign = true
    while !isborder(r,Sud)
        move!(r,Sud)
    end
    while !isborder(r,West)
        move!(r,West)
    end
    for i = 1:z
        move!(r,Ost)
    end
    while y > 0
        if !isborder(r,Nord)
            move!(r,Nord)
            y -= 1
        else
            sign, y = obhod3(r,0,y)
        end
    end
end

alldraw(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs