refined utility methods

src/composite.jl

@@ -197,7 +197,7 @@ end
 @inline Base.expm1(A::Chain{V,G,<:Chain{V,G}},_=nothing) where {V,G} = exp(A)-I
 @inline Base.exp(A::Chain{V,G,<:Chain{V,G},1},_=nothing) where {V,G} = Chain{V,G}(Values(Chain{V,G}(exp(A[1][1]))))
 
-@inline function Base.exp(A::Chain{V,G,Chain{V,G,<:Real,2},2},_=nothing) where {V,G}
+@inline function Base.exp(A::Chain{V,G,<:Chain{V,G,<:Real,2},2},_=nothing) where {V,G}
     T = typeof(exp(zero(valuetype(A))))
     @inbounds a = A[1][1]
     @inbounds c = A[1][2]
@@ -224,7 +224,7 @@ end
     Chain{V,G}(Chain{V,G}(m11, m21), Chain{V,G}(m12, m22))
 end
 
-@inline function Base.exp(A::Chain{V,G,Chain{V,G,<:Complex,2},2},_=nothing) where {V,G}
+@inline function Base.exp(A::Chain{V,G,<:Chain{V,G,<:Complex,2},2},_=nothing) where {V,G}
     T = typeof(exp(zero(valuetype(A))))
     @inbounds a = A[1][1]
     @inbounds c = A[1][2]
@@ -244,7 +244,7 @@ end
 
 # Adapted from implementation in Base; algorithm from
 # Higham, "Functions of Matrices: Theory and Computation", SIAM, 2008
-function Base.exp(_A::Chain{W,G,Chain{W,G,T,N},N},_=nothing) where {W,G,T,N}
+function Base.exp(_A::Chain{W,G,<:Chain{W,G,T,N},N},_=nothing) where {W,G,T,N}
     S = typeof((zero(T)*zero(T) + zero(T)*zero(T))/one(T))
     A = Chain{W,G}(map.(S,value(_A)))
     # omitted: matrix balancing, i.e., LAPACK.gebal!
@@ -289,7 +289,7 @@ function Base.exp(_A::Chain{W,G,Chain{W,G,T,N},N},_=nothing) where {W,G,T,N}
             S(64764752532480000)*I
         expA = (V - U) \ (V + U)
         if s > 0            # squaring to reverse dividing by power of 2
-            for t=1:si
+            for t=Base.OneTo(si)
                 expA = expA*expA
             end
         end
@@ -821,8 +821,9 @@ function inv_approx(t::Chain{M,1,<:Chain{V,1}}) where {M,V}
     mdims(M) < mdims(V) ? (inv(tt⋅t))⋅tt : tt⋅inv(t⋅tt)
 end
 
-Base.:\(t::LinearAlgebra.UniformScaling,v::Chain{V,G,<:Chain}) where {V,G} = inv(v)#value(Chain{V,G}(t))\v
-Base.:\(t::Chain{M,1,<:Chain{W,1}},v::Chain{V,1,<:Chain}) where {M,W,V} = t*inv(v)#Chain{V,1}(t.\value(v))
+Base.:\(v::Chain{V,G,<:Chain},t::LinearAlgebra.UniformScaling) where {V,G} = inv(v)#value(Chain{V,G}(t))\v
+Base.:\(t::LinearAlgebra.UniformScaling,v::Chain{V,G,<:Chain}) where {V,G} = v
+Base.:\(t::Chain{M,1,<:Chain{W,1}},v::Chain{V,1,<:Chain}) where {M,W,V} = inv(t)*v#Chain{V,1}(t.\value(v))
 Base.:\(t::Chain{M,1,<:Chain{W,1}},v::Chain{V,1}) where {M,W,V} = value(t)\v
 Base.in(v::Chain{V,1},t::Chain{W,1,<:Chain{V,1}}) where {V,W} = v ∈ value(t)
 #Base.inv(t::Chain{V,1,<:Chain{V,G}}) where {V,G} = value(Chain{V,G}(I))\t
@@ -1098,3 +1099,146 @@ Base.rand(::AbstractRNG,::SamplerType{PseudoCouple{V}}) where V = rand(PseudoCou
 Base.rand(::AbstractRNG,::SamplerType{PseudoCouple{V,B}}) where {V,B} = PseudoCouple{V,B}(rand(Complex{Float64}))
 Base.rand(::AbstractRNG,::SamplerType{PseudoCouple{V,B,T}}) where {V,B,T} = PseudoCouple{V,B}(rand(Complex{T}))
 Base.rand(::AbstractRNG,::SamplerType{PseudoCouple{V,B,T} where B}) where {V,T} = rand(PseudoCouple{V,Submanifold{V}(UInt(rand(0:(1<<mdims(V)-1)-1))),T})
+
+zero!(x::T) where T = x≈0 ? zero(T) : x
+zero!(x::Complex) = Complex(zero!(real(x)),zero!(imag(x)))
+subzero(a,b) = (ab = a/b; ab ≈ 1 ? zero(ab) : a-b)
+subsqrt(a,b) = (ab = a/b; ab ≈ 1 ? zero(ab) : sqrt(a-b))
+subsqrtcomplex(a,b) = (ab = a/b; ab ≈ 1 ? zero(ab) : (amb = a-b; sqrt(amb < 0 ? Complex(amb) : amb)))
+
+export roots, rootsreal, rootscomplex, monicroots, monicrootsreal, monicrootscomplex
+
+roots(a) = roots(value(a)...)
+roots(a0::Real) = zero(a0)
+roots(a0::Complex) = zero(a0)
+roots(a0,a1) = monicroots(a0/a1)
+roots(a0,a1,a2) = monicroots(a0/a2,a1/a2)
+roots(a::Vararg{T,N}) where {N,T<:Real} = monicroots((@inbounds a[list(1,N-1)]./a[N])...)
+
+rootsreal(a) = rootsreal(value(a)...)
+rootsreal(a0::Real) = zero(a0)
+rootsreal(a::Vararg{T,N}) where {N,T<:Real} = monicrootsreal((@inbounds a[list(1,N-1)]./a[N])...)
+
+rootscomplex(a) = rootscomplex(value(a)...)
+rootscomplex(a0::Real) = zero(Complex{typeof(a0)})
+rootscomplex(a0::Complex) = zero(Complex{typeof(a0)})
+rootscomplex(a0,a1) = monicrootscomplex(a0/a1)
+rootscomplex(a0,a1,a2) = monicrootscomplex(a0/a2,a1/a2)
+rootscomplex(a::Vararg{T,N}) where {N,T<:Real} = monicrootscomplex((@inbounds a[list(1,N-1)]./a[N])...)
+
+monicroots(a) = monicroots(value(a)...)
+monicroots(a...) = value(eigvals(companion(Values(a...))))
+monicroots(a0::Real) = -a0 # z+a0
+monicroots(a0::Complex) = -a0 # z+a0
+function monicroots(a0,a1) # z^2+a1*z+a0
+    sq = a1*a1-4a0
+    quadratic(a0,a1,sqrt(sq < 0 ? Complex(sq) : sq))
+end
+function quadratic(a0,a1,rt)
+    if a1 < 0
+        Values(2a0/(-a1+rt),(-a1+rt)/2)
+    else
+        Values((-a1-rt)/2,2a0/(-a1-rt))
+    end
+end
+function monicroots(a0::Real,a1::Real,a2::Real,::Val{C}=Val(false)) where C
+    a22,a23 = a2*a2,a2/3 # z^3+a2*z^2+a1*z+a0
+    q,r = subzero((a1/3),a22/9),(a1*a2/6)-(a0/2)-(a22*a2/27)
+    r2,q3 = r*r,q*q*q
+    if r2+q3 > 0
+        A = cbrt(abs(r)+sqrt(r2+q3))
+        qA = q/A
+        t = r < 0 ? qA-A : A-qA
+        x,y,z = zero!(-((t/2)+a23)),zero!((sqrt(3)/2)*(A+qA)),t-a23
+        if z < x
+            Values(z,Complex(x,-y),Complex(x,y))
+        else
+            Values(Complex(x,-y),Complex(x,y),z)
+        end
+    else # Real (Viete)
+        sq = sqrt(-q)
+        ϕ1,sq2 = if q<0
+            c = r/(q>-1 ? sqrt(-q3) : sq*sq*sq)
+            acos(abs(c)≈1 ? float(sign(c)) : c)/3,2sq
+        else
+            sq/3,2sq
+        end
+        ϕ2,ϕ3 = ϕ1-(2π/3),ϕ1+(2π/3)
+        out = Values(sq2*cos(ϕ3)-a23,sq2*cos(ϕ2)-a23,sq2*cos(ϕ1)-a23)
+        C ? Complex.(out) : out
+    end
+end
+cubicmax(a) = cubicmax(value(a)...)
+function cubicmax(a0::Real,a1::Real,a2::Real) # z^3+a2*z^2+a1*z+a0
+    a22,a23 = a2*a2,a2/3
+    q,r = subzero((a1/3),a22/9),(a1*a2/6)-(a0/2)-(a22*a2/27)
+    r2,q3 = r*r,q*q*q
+    if r2+q3 > 0
+        A = cbrt(abs(r)+sqrt(r2+q3))
+        (r < 0 ? (q/A)-A : A-(q/A))-a23
+    else # Real (Viete)
+        sq = sqrt(-q)
+        θ = if q<0
+            c = r/(q>-1 ? sqrt(-q3) : sq*sq*sq)
+            acos(abs(c)≈1 ? float(sign(c)) : c)
+        else
+            sq
+        end
+        2sq*cos(θ/3)-a23
+    end
+end
+function quartic(a0::Real,a1::Real,a2::Real,a3::Real) # z^4+a3*z^3+a2*z^2+a1*z+a0
+    a04 = 4a0
+    u = cubicmax(a04*a2-a1*a1-a0*a3*a3,a1*a3-a04,-a2)
+    a32,u2 = a3/2,u/2
+    z1 = zero!(a32*a32+u-a2)
+    psq,qsq = z1≤0 ? zero(typeof(z1)) : sqrt(z1),subsqrt(u2*u2,a0)
+    p1,p2,qsqpm = a32-psq,a32+psq,(a1-(a3*u)/2 > 0 ? qsq : -qsq)
+    q1,q2,p12,p22 = u2+qsqpm,u2-qsqpm,(p1/-2),(p2/-2)
+end
+function monicroots(a0::Real,a1::Real,a2::Real,a3::Real) # z^4+a3*z^3+a2*z^2+a1*z+a0
+    q1,q2,p12,p22 = quartic(a0,a1,a2,a3)
+    sq1,sq2 = p12*p12-q1,p22*p22-q2
+    rt1,rt2 = sqrt(sq1 < 0 ? Complex(sq1) : sq1),sqrt(sq2 < 0 ? Complex(sq2) : sq2)
+    if isreal(rt1) && isreal(rt2)
+        Values{4,typeof(q1)}(p22-rt2,p22+rt2,p12-rt1,p12+rt1)
+    else
+        Values{4,Complex{typeof(q1)}}(p22-rt2,p22+rt2,p12-rt1,p12+rt1)
+    end
+end
+
+monicrootsreal(a) = monicrootsreal(value(a)...)
+monicrootsreal(a...) = value(eigvalsreal(companion(Values(a...))))
+monicrootsreal(a0::Real) = -a0 # z+z0
+monicrootsreal(a0::Real,a1::Real) = quadratic(a0,a1,sqrt(a1*a1-4a0)) # z^2+a1*z+a0
+function monicrootsreal(a0::Real,a1::Real,a2::Real) # z^3+a2*z^2+a1*z+a0
+    a22,a23 = a2*a2,a2/3
+    q,r = subzero((a1/3),a22/9),(a1*a2/6)-(a0/2)-(a22*a2/27)
+    sq = sqrt(-q)
+    ϕ1,sq2 = if q<0
+        c = r/(q>-1 ? sqrt(-q*q*q) : sq*sq*sq)
+        acos(abs(c)≈1 ? float(sign(c)) : c)/3,2sq
+    else
+        sq/3,2sq
+    end
+    ϕ2,ϕ3 = ϕ1-(2π/3),ϕ1+(2π/3)
+    Values(sq2*cos(ϕ3)-a23,sq2*cos(ϕ2)-a23,sq2*cos(ϕ1)-a23)
+end
+function monicrootsreal(a0::Real,a1::Real,a2::Real,a3::Real) # z^4+a3*z^3+a2*z^2+a1*z+a0
+    q1,q2,p12,p22 = quartic(a0,a1,a2,a3)
+    rt1,rt2 = sqrt(zero!(p12*p12-q1)),sqrt(zero!(p22*p22-q2))
+    Values(p22-rt2,p22+rt2,p12-rt1,p12+rt1)
+end
+
+monicrootscomplex(a) = monicrootscomplex(value(a)...)
+monicrootscomplex(a...) = value(eigvalscomplex(companion(Values(a...))))
+monicrootscomplex(a0::Real) = Complex(-a0) # z+a0
+monicrootscomplex(a0::Complex) = -a0 # z+a0
+monicrootscomplex(a0,a1) = quadratic(a0,a1,sqrt(Complex(a1*a1-4a0))) # z^2+a1*z+a0
+monicrootscomplex(a0::Real,a1::Real,a2::Real) = monicroots(a0,a1,a2,Val(true))
+function monicrootscomplex(a0::Real,a1::Real,a2::Real,a3::Real)# z^4+a3*z^3+a2*z^2+a1*z+a0
+    q1,q2,p12,p22 = quartic(a0,a1,a2,a3)
+    rt1,rt2 = sqrt(Complex(subzero(p12*p12,q1))),sqrt(Complex(subzero(p22*p22,q2)))
+    Values(p22-rt2,p22+rt2,p12-rt1,p12+rt1)
+end
+

src/forms.jl

@@ -145,7 +145,7 @@ function (W::Submanifold{Q,M,S})(m::Multivector{V,T}) where {Q,M,V,S,T}
     end
 end
 
-# need separate storage for m and F for caching
+#= need separate storage for m and F for caching
 
 const dualform_cache = Vector{Tuple{Int,Bool}}[]
 const dualformC_cache = Vector{Tuple{Int,Bool}}[]
@@ -292,7 +292,10 @@ end
         end
         return Single{V}(out::t,Submanifold{V}())
     end
-end
+end=#
+
+(t::TensorGraded)(y::TensorGraded...) = contraction(t,∧(y...))
+(t::TensorMixed)(y::TensorGraded...) = contraction(t,∧(y...))
 
 # Dyadic
 
@@ -1344,7 +1347,7 @@ bracket(X,Y) = X(Y) - Y(X)
 bracket(X,Y,Z) = X(bracket(Y,Z)) + Y(bracket(Z,X)) + Z(bracket(X,Y))
 bracket(W,X,Y,Z) = W(bracket(X,Y,Z)) + X(bracket(W,Z,Y)) + Y(bracket(W,X,Z)) + Z(bracket(W,Y,X))
 bracket(V,W,X,Y,Z) = V(bracket(W,X,Y,Z)) + W(bracket(V,X,Z,Y)) + X(bracket(V,W,Y,Z)) + Y(bracket(V,W,Z,X)) + Z(bracket(V,W,X,Y))
-@generated function bracket(X::Vararg{T,N} where T) where N
+@generated function bracket(X::Vararg{T,N}) where {N,T}
     Expr(:call,:+,[:($(isodd(i) ? :+ : :-)(X[$i](bracket($(vcat([j≠i ? [:(X[$j])] : [] for j ∈ list(1,N)]...)...))))) for i ∈ list(1,N)]...)
 end
 

src/parity.jl

@@ -440,7 +440,7 @@ construct_cache(:DiagonalForm)
 
 @pure signbit(V::T) where T<:Manifold = (ib=indexbasis(rank(V)); parity.(Ref(V),ib,ib))
 @pure signbit(V::T,G) where T<:Manifold = (ib=indexbasis(rank(V),G); parity.(Ref(V),ib,ib))
-@pure angular(V::T) where T<:Manifold = Values(findall(signbit(V))...)
+#=@pure angular(V::T) where T<:Manifold = Values(findall(signbit(V))...)
 @pure radial(V::T) where T<:Manifold = Values(findall(.!signbit(V))...)
 @pure angular(V::T,G) where T<:Manifold = findall(signbit(V,G))
 @pure radial(V::T,G) where T<:Manifold = findall(.!signbit(V,G))
@@ -456,7 +456,7 @@ for (op,other) ∈ ((:angular,:radial),(:radial,:angular))
             Chain{V,G}(out)
         end
     end
-end
+end=#
 
 Base.iseven(t::Zero) = true
 Base.isodd(t::Zero) = true
@@ -481,8 +481,35 @@ even(t::Couple{V,B}) where {V,B} = iseven(grade(B)) ? t : scalar(t)
 odd(t::Couple{V,B}) where {V,B} = isodd(grade(B)) ? imaginary(t) : Zero{V}()
 even(t::PseudoCouple{V,B}) where {V,B} = even(imaginary(t)) + even(volume(t))
 odd(t::PseudoCouple{V,B}) where {V,B} = odd(imaginary(t)) + odd(volume(t))
-
-function imag(t::Multivector{V,T}) where {V,T}
+even(t::Multivector{V,T,4} where {V,T}) = scalar(t) + volume(t)
+odd(t::Multivector{V,T,4} where {V,T}) = vector(t)
+
+real(t::TensorTerm{V,G}) where {V,G} = parityreverse(G) ? Zero{V}() : t
+imag(t::TensorTerm{V,G}) where {V,G} = !parityreverse(G) ? Zero{V}() : t
+real(t::Couple) = scalar(t) + real(imaginary(t))
+imag(t::Couple) = imag(imaginary(t))
+real(t::PseudoCouple) = real(volume(t)) + real(imaginary(t))
+imag(t::PseudoCouple) = imag(volume(t)) + imag(imaginary(t))
+real(t::Multivector{V,T,2} where {V,T}) = t
+imag(t::Multivector{V,T,2}) where {V,T} = Zero{V}()
+imag(t::Multivector{V,T,4} where {V,T}) = volume(t)
+real(t::Spinor{V,T,2} where {V,T}) = scalar(t)
+real(t::Quaternion) = scalar(t)
+real(t::Spinor{V,T,8} where {V,T}) = scalar(t) + volume(t)
+imag(t::Spinor{V,T,2} where {V,T}) = volume(t)
+imag(t::Quaternion) = imaginary(t)
+imag(t::Spinor{V,T,8} where {V,T}) = bivector(t)
+imag(t::Spinor{V,T,16} where {V,T}) = bivector(t)
+real(t::CoSpinor{V,T,2} where {V,T}) = t
+real(t::CoSpinor{V,T,4} where {V,T}) = vector(t)
+real(t::CoSpinor{V,T,8} where {V,T}) = vector(t)
+imag(t::CoSpinor{V,T,2}) where {V,T} = Zero{V}()
+imag(t::CoSpinor{V,T,4} where {V,T}) = volume(t)
+imag(t::CoSpinor{V,T,8} where {V,T}) = trivector(t)
+imag(t::CoSpinor{V,T,16} where {V,T}) = trivector(t)
+imag(t::CoSpinor{V,T,32} where {V,T}) = trivector(t)
+
+#=function imag(t::Multivector{V,T}) where {V,T}
     N = mdims(V)
     out = copy(value(t,mvec(N,T)))
     bs = binomsum_set(N)
@@ -504,27 +531,49 @@ function real(t::Multivector{V,T}) where {V,T}
         end
     end
     Multivector{V}(out)
-end
-function imag(t::Spinor{V,T}) where {V,T}
+end=#
+#=function imag(t::Spinor{V,T}) where {V,T}
     N = mdims(V)
     out = copy(value(t,mvecs(N,T)))
-    bs = spinsum_set(N)
+    bs = spincumsum(N)
     @inbounds out[1]≠0 && (out[1] = zero(T))
-    for g ∈ evens(2,N+1)
+    for g ∈ evens(3,N+1)
         @inbounds !parityreverse(g-1) && for k ∈ bs[g]+1:bs[g+1]
             @inbounds out[k]≠0 && (out[k] = zero(T))
         end
     end
     Spinor{V}(out)
-end
-function real(t::Spinor{V,T}) where {V,T}
+end=#
+#=function real(t::Spinor{V,T}) where {V,T}
     N = mdims(V)
     out = copy(value(t,mvecs(N,T)))
-    bs = spinsum_set(N)
+    bs = spincumsum(N)
     for g ∈ evens(3,N+1)
         @inbounds parityreverse(g-1) && for k ∈ bs[g]+1:bs[g+1]
             @inbounds out[k]≠0 && (out[k] = zero(T))
         end
     end
     Spinor{V}(out)
-end
+end=#
+#=function imag(t::CoSpinor{V,T}) where {V,T}
+    N = mdims(V)
+    out = copy(value(t,mvecs(N,T)))
+    bs = anticumsum(N)
+    for g ∈ evens(2,N+1)
+        @inbounds !parityreverse(g-1) && for k ∈ bs[g]+1:bs[g+1]
+            @inbounds out[k]≠0 && (out[k] = zero(T))
+        end
+    end
+    CoSpinor{V}(out)
+end=#
+#=function real(t::CoSpinor{V,T}) where {V,T}
+    N = mdims(V)
+    out = copy(value(t,mvecs(N,T)))
+    bs = anticumsum(N)
+    for g ∈ evens(2,N+1)
+        @inbounds parityreverse(g-1) && for k ∈ bs[g]+1:bs[g+1]
+            @inbounds out[k]≠0 && (out[k] = zero(T))
+        end
+    end
+    CoSpinor{V}(out)
+end=#