From 2523f2548287733dbb60df960add153a11b1908f Mon Sep 17 00:00:00 2001 From: Michael Reed <18372368+chakravala@users.noreply.github.com> Date: Tue, 21 May 2024 16:00:10 -0400 Subject: [PATCH] improved robustness of algebra --- Project.toml | 2 +- src/Grassmann.jl | 4 +- src/algebra.jl | 73 ++++++++----- src/composite.jl | 4 +- src/forms.jl | 128 +++++++++++++++++++++-- src/multivectors.jl | 36 ++++--- src/parity.jl | 220 +++++++++++++++++++++++++++++++++++++-- src/products.jl | 237 ++++++++++++++++++++++++++++++------------- test/generictests.jl | 2 +- 9 files changed, 579 insertions(+), 127 deletions(-) diff --git a/Project.toml b/Project.toml index 7ef2514..1150ba8 100644 --- a/Project.toml +++ b/Project.toml @@ -16,7 +16,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] julia = "1" Leibniz = "0.2" -DirectSum = "0.8.1" +DirectSum = "0.8.12" AbstractTensors = "0.8" ComputedFieldTypes = "1" Requires = "1" diff --git a/src/Grassmann.jl b/src/Grassmann.jl index 59151cf..c3ef353 100644 --- a/src/Grassmann.jl +++ b/src/Grassmann.jl @@ -26,7 +26,7 @@ import AbstractTensors: Values, Variables, FixedVector, clifford, hodge, wedge, export ⊕, ℝ, @V_str, @S_str, @D_str, Manifold, Submanifold, Signature, DiagonalForm, value export @basis, @basis_str, @dualbasis, @dualbasis_str, @mixedbasis, @mixedbasis_str, Λ export ℝ0, ℝ1, ℝ2, ℝ3, ℝ4, ℝ5, ℝ6, ℝ7, ℝ8, ℝ9, mdims, tangent, metric, antimetric -export hodge, wedge, vee, complement, dot, antidot, istangent +export hodge, wedge, vee, complement, dot, antidot, istangent, Values import Base: @pure, ==, isapprox import Base: print, show, getindex, setindex!, promote_rule, convert, adjoint @@ -34,7 +34,7 @@ import DirectSum: V0, ⊕, generate, basis, getalgebra, getbasis, dual, Zero, On import Leibniz: hasinf, hasorigin, dyadmode, value, pre, vsn, metric, mdims, gdims import Leibniz: bit2int, indexbits, indices, diffvars, diffmask import Leibniz: symmetricmask, indexstring, indexsymbol, combo, digits_fast -import DirectSum: antimetric +import DirectSum: antimetric, signbool import Leibniz: hasconformal, hasinf2origin, hasorigin2inf import AbstractTensors: valuetype, scalar, isscalar, trivector, istrivector, ⊗, complement diff --git a/src/algebra.jl b/src/algebra.jl index fae49c5..4603046 100644 --- a/src/algebra.jl +++ b/src/algebra.jl @@ -36,25 +36,33 @@ Geometric algebraic product: ω⊖η = (-1)ᵖdet(ω∩η)⊗(Λ(ω⊖η)∪L(ω *(a::X,b::Y,c::Z...) where {X<:TensorAlgebra,Y<:TensorAlgebra,Z<:TensorAlgebra} = *(a*b,c...) @pure function mul(a::Submanifold{V},b::Submanifold{V},der=derive_mul(V,UInt(a),UInt(b),1,true)) where V - ba,bb = UInt(a),UInt(b) - (diffcheck(V,ba,bb) || iszero(der)) && (return Zero(V)) - A,B,Q,Z = symmetricmask(V,ba,bb) - pcc,bas,cc = (hasinf(V) && hasorigin(V)) ? conformal(V,A,B) : (false,A⊻B,false) - d = getbasis(V,bas|Q) - out = (typeof(V)<:Signature || count_ones(A&B)==0) ? (parity(a,b)⊻pcc ? Single{V}(-1,d) : d) : Single{V}((pcc ? -1 : 1)*parityinner(V,A,B),d) - diffvars(V)≠0 && !iszero(Z) && (out = Single{V}(getbasis(loworder(V),Z),out)) - return cc ? (v=value(out);out+Single{V}(hasinforigin(V,A,B) ? -(v) : v,getbasis(V,conformalmask(V)⊻UInt(d)))) : out + if isdiag(V) + ba,bb = UInt(a),UInt(b) + (diffcheck(V,ba,bb) || iszero(der)) && (return Zero(V)) + A,B,Q,Z = symmetricmask(V,ba,bb) + d = getbasis(V,(A⊻B)|Q) + out = if typeof(V)<:Signature || count_ones(A&B)==0 + (parity(a,b) ? Single{V}(-1,d) : d) + else + Single{V}(signbool(parityinner(V,A,B)),d) + end + diffvars(V)≠0 && !iszero(Z) && (out = Single{V}(getbasis(loworder(V),Z),out)) + return out + else + out = paritygeometric(V,UInt(a),UInt(b)) + isempty(out) ? Zero(V) : +(Single{V}.(out)...) + end end function ⟑(a::Single{V},b::Submanifold{V}) where V v = derive_mul(V,UInt(basis(a)),UInt(b),a.v,true) bas = mul(basis(a),b,v) - order(a.v)+order(bas)>diffmode(V) ? Zero(V) : Single{V}(v,bas) + order(a.v)+order(bas)>diffmode(V) ? Zero(V) : v*bas end function ⟑(a::Submanifold{V},b::Single{V}) where V v = derive_mul(V,UInt(a),UInt(basis(b)),b.v,false) bas = mul(a,basis(b),v) - order(b.v)+order(bas)>diffmode(V) ? Zero(V) : Single{V}(v,bas) + order(b.v)+order(bas)>diffmode(V) ? Zero(V) : v*bas end export ∗, ⊛, ⊖ @@ -173,24 +181,43 @@ export ⋅ Interior (right) contraction product: ω⋅η = ω∨⋆η """ @pure function contraction(a::Submanifold{V},b::Submanifold{V}) where V - g,C,t,Z = interior(a,b) - (!t || iszero(derive_mul(V,UInt(a),UInt(b),1,true))) && (return Zero(V)) - d = getbasis(V,C) - istangent(V) && !iszero(Z) && (d = Single{V}(getbasis(loworder(V),Z),d)) - return isone(g) ? d : Single{V}(g,d) + iszero(derive_mul(V,UInt(a),UInt(b),1,true)) && (return Zero(V)) + if isdiag(V) || hasconformal(V) + g,C,t,Z = interior(a,b) + (!t) && (return Zero(V)) + d = getbasis(V,C) + istangent(V) && !iszero(Z) && (d = Single{V}(getbasis(loworder(V),Z),d)) + return isone(g) ? d : Single{V}(g,d) + else + Cg,Z = parityinterior(V,UInt(a),UInt(b),Val(true)) + #d = getbasis(V,C) + #istangent(V) && !iszero(Z) && (d = Single{V}(getbasis(loworder(V),Z),d)) + return isempty(Cg) ? Zero(V) : +(Single{V}.(Cg)...) + end end function contraction(a::X,b::Y) where {X<:TensorTerm{V},Y<:TensorTerm{V}} where V ba,bb = UInt(basis(a)),UInt(basis(b)) - g,C,t,Z = interior(V,ba,bb) - !t && (return Zero(V)) - v = derive_mul(V,ba,bb,value(a),value(b),AbstractTensors.dot) - if istangent(V) && !iszero(Z) - _,_,Q,_ = symmetricmask(V,ba,bb) - v = !(typeof(v)<:TensorTerm) ? Single{V}(v,getbasis(V,Z)) : Single{V}(v,getbasis(loworder(V),Z)) - count_ones(Q)+order(v)>diffmode(V) && (return Zero(V)) + if isdiag(V) || hasconformal(V) + g,C,t,Z = interior(V,ba,bb) + !t && (return Zero(V)) + v = derive_mul(V,ba,bb,value(a),value(b),AbstractTensors.dot) + if istangent(V) && !iszero(Z) + _,_,Q,_ = symmetricmask(V,ba,bb) + v = !(typeof(v)<:TensorTerm) ? Single{V}(v,getbasis(V,Z)) : Single{V}(v,getbasis(loworder(V),Z)) + count_ones(Q)+order(v)>diffmode(V) && (return Zero(V)) + end + return Single{V}(g*v,getbasis(V,C)) + else + Cg,Z = parityinterior(V,ba,bb,Val(true)) + v = derive_mul(V,ba,bb,value(a),value(b),AbstractTensors.dot) + if istangent(V) && !iszero(Z) + _,_,Q,_ = symmetricmask(V,ba,bb) + v = !(typeof(v)<:TensorTerm) ? Single{V}(v,getbasis(V,Z)) : Single{V}(v,getbasis(loworder(V),Z)) + count_ones(Q)+order(v)>diffmode(V) && (return Zero(V)) + end + return isempty(Cg) ? Zero(V) : (+(Single{V}.(Cg)...))*v end - return Single{V}(g*v,getbasis(V,C)) end export ⨼, ⨽ diff --git a/src/composite.jl b/src/composite.jl index c0267af..735812c 100644 --- a/src/composite.jl +++ b/src/composite.jl @@ -183,7 +183,8 @@ end function Base.exp(t::T) where T<:TensorGraded S,B,V = T<:Submanifold,T<:TensorTerm,Manifold(t) if B && isnull(t) - return One(V) + vt = valuetype(t) + return Couple{V,basis(t)}(one(vt),zero(vt)) elseif isR301(V) && grade(t)==2 # && abs(t[0])<1e-9 && !options.over u = sqrt(abs(abs2(t)[1])) u<1e-5 && (return One(V)+t) @@ -804,6 +805,7 @@ end DirectSum.Λ(x::Chain{V,1,<:Chain{V,1}},G) where V = compound(x,G) compound(x,G::T) where T<:Integer = compound(x,Val(G)) +compound(x::Chain{V,1,<:Chain{V,1}},::Val{0}) where V = Chain{V,0}(Values(Chain{V,0}(1))) @generated function compound(x::Chain{V,1,<:Chain{V,1}},::Val{G}) where {V,G} Expr(:call,:(Chain{V,G}),Expr(:call,:Values,[Expr(:call,:∧,[:(@inbounds x[$i]) for i ∈ indices(j)]...) for j ∈ indexbasis(mdims(V),G)]...)) end diff --git a/src/forms.jl b/src/forms.jl index c492bc6..5301104 100644 --- a/src/forms.jl +++ b/src/forms.jl @@ -339,12 +339,12 @@ for pinor ∈ (:Spinor,:AntiSpinor)#,:Multivector) end end -display_matrix(m::Chain{V,G,<:TensorGraded{W,G}}) where {V,G,W} = vcat(transpose([V,chainbasis(V,G)...]),hcat(chainbasis(W,G),matrix(m))) -display_matrix(m::TensorGraded{V,G,<:Spinor{W}}) where {V,G,W} = vcat(transpose([V,evenbasis(V)...]),hcat(evenbasis(W),matrix(m))) -display_matrix(m::TensorGraded{V,G,<:AntiSpinor{W}}) where {V,G,W} = vcat(transpose([V,oddbasis(V)...]),hcat(oddbasis(W),matrix(m))) -display_matrix(m::TensorAlgebra{V,<:Multivector{W}}) where {V,W} = vcat(transpose([V,fullbasis(V)...]),hcat(fullbasis(W),matrix(m))) +display_matrix(m::Chain{V,G,<:TensorGraded{W,G}}) where {V,G,W} = vcat(transpose([Submanifold(V),chainbasis(V,G)...]),hcat(chainbasis(W,G),matrix(m))) +display_matrix(m::TensorGraded{V,G,<:Spinor{W}}) where {V,G,W} = vcat(transpose([Submanifold(V),evenbasis(V)...]),hcat(evenbasis(W),matrix(m))) +display_matrix(m::TensorGraded{V,G,<:AntiSpinor{W}}) where {V,G,W} = vcat(transpose([Submanifold(V),oddbasis(V)...]),hcat(oddbasis(W),matrix(m))) +display_matrix(m::TensorAlgebra{V,<:Multivector{W}}) where {V,W} = vcat(transpose([Submanifold(V),fullbasis(V)...]),hcat(fullbasis(W),matrix(m))) for (pinor,bas) ∈ ((:Spinor,:evenbasis),(:AntiSpinor,:oddbasis),(:Multivector,:fullbasis)) - @eval display_matrix(m::$pinor{V,<:TensorAlgebra{W}}) where {V,W} = vcat(transpose([V,$bas(V)...]),hcat($bas(W),matrix(m))) + @eval display_matrix(m::$pinor{V,<:TensorAlgebra{W}}) where {V,W} = vcat(transpose([Submanifold(V),$bas(V)...]),hcat($bas(W),matrix(m))) end export Projector, Dyadic, Proj @@ -401,6 +401,15 @@ Chain(P::Dyadic{V}) where V = Chain{V}(P) export TensorOperator, Endomorphism +#=struct TensorOperator{V,W,S<:TensorAlgebra{W},T<:TensorAlgebra{V,S}} <: TensorNested{V,S} + v::T + TensorOperator{V,W}(t::T) where {V,W,S<:TensorAlgebra{W},T<:TensorAlgebra{V,S}} = new{V,W,S,T}(t) + TensorOperator{V}(t::T) where {V,W,S<:TensorAlgebra{W},T<:TensorAlgebra{V,S}} = new{V,W,S,T}(t) + TensorOperator(t::T) where {V,W,S<:TensorAlgebra{W},T<:TensorAlgebra{V,S}} = new{V,W,S,T}(t) +end + +Endomorphism{V,S<:TensorAlgebra{V},T<:TensorAlgebra{V,S}} = TensorOperator{V,V,S,T}=# + struct TensorOperator{V,W,T<:TensorAlgebra{V,<:TensorAlgebra{W}}} <: TensorNested{V,T} v::T TensorOperator{V,W}(t::T) where {V,W,T<:TensorAlgebra{V,<:TensorAlgebra{W}}} = new{V,W,T}(t) @@ -412,6 +421,8 @@ Endomorphism{V,T<:TensorAlgebra{V,<:TensorAlgebra{V}}} = TensorOperator{V,V,T} value(t::TensorOperator) = t.v matrix(m::TensorOperator) = matrix(value(m)) +compound(m::TensorOperator,g) = TensorOperator(compound(value(m),g)) +compound(m::TensorOperator,g::Integer) = TensorOperator(compound(value(m),g)) getindex(t::TensorOperator,i::Int,j::Int) = value(value(t.v)[j])[i] for op ∈ (:(Base.inv),) @@ -432,6 +443,9 @@ function show(io::IO, ::MIME"text/plain", t::TensorOperator) if isempty(X) && get(io, :compact, false)::Bool return show(io, X) end + show_matrix(io, t, X) +end +function show_matrix(io::IO, t, X) # 0) show summary before setting :compact summary(io, t) isempty(X) && return @@ -544,6 +558,15 @@ end @inline @generated function matmul(x::Values{N,<:Multivector{V}},y::Values{N}) where {N,V} Expr(:call,:Values,[Expr(:call,:+,[:(@inbounds y[$i]*value(x[$i])[$j]) for i ∈ list(1,N)]...) for j ∈ list(1,1<value(contraction(x,y))) + end +end + +metricfull(V) = TensorOperator(Multivector{V}(vcat(value.(map.(Multivector,value.(metrictensor.(V,list(0,mdims(V))))))...))) +metriceven(V) = TensorOperator(Spinor{V}(vcat(value.(map.(Spinor,value.(metrictensor.(V,evens(0,mdims(V))))))...))) +metricodd(V) = TensorOperator(AntiSpinor{V}(vcat(value.(map.(AntiSpinor,value.(metrictensor.(V,evens(1,mdims(V))))))...))) + +struct MetricTensor{n,ℙ,g,Vars,Diff,Name} <: TensorBundle{n,ℙ,g,Vars,Diff,Name} + @pure MetricTensor{N,M,S,F,D,L}() where {N,M,S,F,D,L} = new{N,M,S,F,D,L}() +end + +export MetricTensor, metrictensor +@pure MetricTensor{N,M,S,F,D}() where {N,M,S,F,D} = MetricTensor{N,M,S,F,D,1}() +@pure MetricTensor{N,M,S}() where {N,M,S} = MetricTensor{N,M,S,0,0}() +@pure MetricTensor{N,M}(b::Values{N,<:Tuple}) where {N,M} = MetricTensor{N,M,metricsig(M,Values.(b))}() +@pure MetricTensor{N,M}(b::Values{N,<:Values}) where {N,M} = MetricTensor{N,M,metricsig(M,b)}() +@pure MetricTensor{N,M}(b::Values{N,<:Chain}) where {N,M} = MetricTensor{N,M,metricsig(M,value.(b))}() +@pure MetricTensor{N,M}(b::Values{N,<:AbstractVector}) where {N,M} = MetricTensor{N,M,metricsig(M,Values{N}.(b))}() +MetricTensor{N,M}(b::Vector) where {N,M} = MetricTensor{N,M}(Values(b...)) +@pure MetricTensor(b::Tuple) = MetricTensor(Values(b)) +@pure MetricTensor(b::Values{N}) where N = MetricTensor{N,0}(b) +MetricTensor(b::Values{N,<:Real}) where N = DiagonalForm(b) +MetricTensor(b::AbstractVector{<:Real}) = DiagonalForm(b) +MetricTensor(b::AbstractVector) = MetricTensor{length(b),0}(b) +MetricTensor(b::AbstractMatrix) = MetricTensor(getindex.(Ref(b),:,1:(size(b)[1]))) +MetricTensor(b::Chain{V,G,<:Chain} where {V,G}) = MetricTensor(value(b)) +MetricTensor(b::Chain) = DiagonalForm(value(b)) +MetricTensor(b::TensorOperator) = MetricTensor(value(b)) +MetricTensor(b...) = MetricTensor(b) + +@pure Manifold(::Type{T}) where T<:MetricTensor = T() + +@pure metrictensor(b::Submanifold{V}) where V = isbasis(b) ? metrictensor(V) : TensorOperator(map(b,b(value(metrictensor(V))))) +@pure metrictensor(V,G) = compound(metrictensor(V),G) +@pure function metrictensor(V::MetricTensor{N,M,S} where N) where {M,S} + out = Chain{Submanifold(V)}.(metrictensor_cache[S]) + TensorOperator(Chain{Submanifold(V)}(isdual(V) ? SUB(out) : out)) +end +const metrictensor_cache = Values[] +@pure function metricsig(M,b::Values) + a = dyadmode(M)>0 ? SUB(b) : b + if a ∈ metrictensor_cache + findfirst(x->x==a,metrictensor_cache) + else + push!(metrictensor_cache,a) + length(metrictensor_cache) + end +end + +@pure DirectSum.getalgebra(V::MetricTensor) = DirectSum.getalgebra(Submanifold(V)) + +for t ∈ (Any,Integer) + @eval @inline getindex(s::MetricTensor{N,M,S} where {N,M},i::T) where {S,T<:$t} = value(value(metrictensor(s))[i]) +end +@inline getindex(vs::MetricTensor,i::Vector) = [getindex(vs,j) for j ∈ i] +@inline getindex(vs::MetricTensor,i::UnitRange{Int}) = [getindex(vs,j) for j ∈ i] +@inline getindex(vs::MetricTensor{N,M,S} where M,i::Colon) where {N,S} = Vector(value(metrictensor(vs))) + +@pure Signature(V::MetricTensor{N,M}) where {N,M} = Signature{N,M,UInt(0)}() + +# anything array-like gets summarized e.g. 10-element Array{Int64,1} +Base.summary(io::IO, a::MetricTensor) = Base.array_summary(io, a, _axes(metrictensor(a))) + +show(io::IO,M::MetricTensor) = Base.show(io,Submanifold(M)) +function show(io::IO, ::MIME"text/plain", M::MetricTensor) + X = display_matrix(value(metrictensor(M))) + if isempty(X) && get(io, :compact, false)::Bool + return show(io, X) + end + show_matrix(io,M,X) +end + +(M::MetricTensor)(b::Int...) = Submanifold{M}(b) +(M::MetricTensor)(b::T) where T<:AbstractVector{Int} = Submanifold{M}(b) +(M::MetricTensor)(b::T) where T<:AbstractRange{Int} = Submanifold{M}(b) + +import LinearAlgebra: isdiag; export isdiag +isdiag(::MetricTensor) = false diff --git a/src/multivectors.jl b/src/multivectors.jl index 9c1cb44..bb622c7 100644 --- a/src/multivectors.jl +++ b/src/multivectors.jl @@ -189,6 +189,11 @@ Base.ones(::Type{Chain{V,G,T,X}}) where {V,G,T<:Chain,X} = Chain{V,G,T}(ones.(nt ⊗(a::Type{<:Chain{V,1}},b::Type{<:Chain{W,1}}) where {V,W} = Chain{V,1,Chain{W,1,Float64,mdims(W)},mdims(V)} ⊗(a::Type{<:Chain{V,1}},b::Type{<:Chain{W,1,T}}) where {V,W,T} = Chain{V,1,Chain{W,1,T,mdims(W)},mdims(V)} +function single(t::Chain) + i = findfirst(x->!isnull(x),value(t)) + norm(t[i]) == norm(t) ? t(i) : t +end + """ ChainBundle{V,G,T,P} <: Manifold{V,T} <: TensorAlgebra{V,T} @@ -297,6 +302,7 @@ end end) end +Multivector(::Zero{V}) where V = Multivector{V}(zeros(tvec(mdims(V),Int))) Multivector(val::TensorAlgebra{V}) where V = Multivector{V}(val) Multivector{V}(val::NTuple{N,T}) where {V,N,T} = Multivector{V}(Values{N,T}(val)) Multivector{V}(val::NTuple{N,Any}) where {V,N} = Multivector{V}(Values{N}(val)) @@ -442,18 +448,16 @@ abstract type AbstractSpinor{V,T} <: TensorMixed{V,T} end @pure log2sub2(N) = log2sub(2N) for pinor ∈ (:Spinor,:AntiSpinor) - dpinor = Symbol(:Dyadic,pinor) @eval begin @computed struct $pinor{V,T} <: AbstractSpinor{V,T} v::Values{1<<(mdims(V)-1),T} - $pinor{V,T}(v::Values{N,T}) where {N,V,T} = new{DirectSum.submanifold(V),T}(v) + $pinor{V,T}(v) where {V,T} = new{DirectSum.submanifold(V),T}(v) end - $pinor{V,T}(v::AbstractVector{T}) where {V,T} = $pinor{V,T}(Values{1<<(mdims(V)-1),T}(v)) + #$pinor{V,T}(v::AbstractVector{T}) where {V,T} = $pinor{V,T}(Values{1<<(mdims(V)-1),T}(v) $pinor{V}(v::AbstractVector{T}) where {V,T} = $pinor{V,T}(v) $pinor{V,𝕂}(val::Single{V,G,B,𝕂}) where {V,G,B,𝕂} = $pinor{V,𝕂}(val.v,B) $pinor{V}(val::Single{V,G,B,𝕂}) where {V,G,B,𝕂} = $pinor{V,𝕂}(val) $pinor(val::TensorAlgebra{V}) where V = $pinor{V}(val) - $pinor(t::Chain{V}) where V = $pinor{V}(t) $pinor{V}(val::NTuple{N,T}) where {V,N,T} = $pinor{V}(Values{N,T}(val)) $pinor{V}(val::NTuple{N,Any}) where {V,N} = $pinor{V}(Values{N}(val)) $pinor(val::NTuple{N,T}) where {N,T} = $pinor{log2sub2(N)}(Values{N,T}(val)) @@ -476,7 +480,6 @@ for pinor ∈ (:Spinor,:AntiSpinor) equal(a::Multivector{V,T},b::$pinor{V,S}) where {V,T,S} = equal(a,Multivector(b)) equal(a::Chain{V,G,T},b::$pinor{V,S}) where {V,S,G,T} = b == a equal(a::T,b::$pinor{V,S} where S) where T<:TensorTerm{V} where V = b==a - $dpinor{V,T,N} = $pinor{V,$pinor{V,T,N},N} end end @@ -516,15 +519,17 @@ AntiSpinor{V,T}(v::Submanifold{V,G}) where {V,G,T} = AntiSpinor{V,T}(one(T),v) end end -@generated function Spinor{V}(t::Chain{V,G,T}) where {V,G,T} - isodd(G) && error("$t is not expressible as a Spinor") - N = mdims(V) - chain_src(N,G,T,evens(0,N),:Spinor) -end -@generated function AntiSpinor{V}(t::Chain{V,G,T}) where {V,G,T} - iseven(G) && error("$t is not expressible as a AntiSpinor") - N = mdims(V) - chain_src(N,G,T,evens(1,N),:AntiSpinor) +for var ∈ ((:V,:T),(:T,),()) + @eval @generated function Spinor{$(var...)}(t::Chain{V,G,T}) where {V,G,T} + isodd(G) && error("$t is not expressible as a Spinor") + N = mdims(V) + chain_src(N,G,T,evens(0,N),:Spinor) + end + @eval @generated function AntiSpinor{$(var...)}(t::Chain{V,G,T}) where {V,G,T} + iseven(G) && error("$t is not expressible as a AntiSpinor") + N = mdims(V) + chain_src(N,G,T,evens(1,N),:AntiSpinor) + end end @generated function (t::Spinor{V,T})(G::Int) where {V,T} @@ -713,6 +718,9 @@ Couple(m::TensorTerm{V}) where V = Couple{V,basis(m)}(Complex(m)) Couple(m::TensorTerm{V,0}) where V = Couple{V,Submanifold(V)}(Complex(m)) PseudoCouple(m::TensorTerm{V,G}) where {V,G} = G==grade(V) ? PseudoCouple{V,One(V)}(value(m)*im) : PseudoCouple{V,basis(m)}(Complex(value(m))) +Couple{V,B}(m::TensorTerm{V}) where {V,B} = Couple{V,B}(Complex(m)) +Couple{V,B}(m::TensorTerm{V,0}) where {V,B} = Couple{V,B}(Complex(m)) + @generated Multivector{V}(a::Single{V,L},b::Single{V,G}) where {V,L,G} = addermulti(a,b,:+) Multivector{V,T}(z::Couple{V,B,T}) where {V,B,T} = Multivector{V}(scalar(z), imaginary(z)) Multivector{V,T}(z::PseudoCouple{V,B,T}) where {V,B,T} = Multivector{V}(imaginary(z), volume(z)) diff --git a/src/parity.jl b/src/parity.jl index bcb76a6..f9edc83 100644 --- a/src/parity.jl +++ b/src/parity.jl @@ -64,11 +64,11 @@ end A,B,Q,Z = symmetricmask(V,a,b) α,β = complement(N,A,D),complement(N,B,D) cc = skew && (hasinforigin(V,A,β) || hasorigininf(V,A,β)) - if ((count_ones(α&β)==0) && !diffcheck(V,α,β)) || cc + if ((count_ones(α&β)==0) && !diffcheck(V,α,β)) #|| cc C,L = α ⊻ β, count_ones(A)+count_ones(B) bas = complement(N,C,D) pcc,bas = if skew - A3,β3,i2o,o2i,xor = conformalcheck(V,A,β) + A3,β3,i2o,o2i,xor = UInt(0),UInt(0),false,false,false#conformalcheck(V,A,β) cx = cc || xor cx && parity(V,A3,β3)⊻(i2o || o2i)⊻(xor&!i2o), cx ? (A3|β3)⊻bas : bas else @@ -82,7 +82,7 @@ end end @pure parityregressive(V::Signature,a,b,skew=Val(false)) = _parityregressive(V,a,b,skew) -@pure function parityregressive(V::M,A,B) where M<:Manifold +@pure function parityregressive(::M,A,B) where M<:Manifold{V} where V p,C,t,Z = parityregressive(Signature(V),A,B) return p ? -1 : 1, C, t, Z end @@ -101,12 +101,47 @@ end return t ? p⊻parityrighthodge(S,B,N) : p, C|Q, t, Z end=# -@pure function parityinterior(V::M,a,b) where M<:Manifold +@pure function diffcheck2(V,A::UInt,B::UInt) + d,db = diffvars(V),diffmask(V) + v = isdyadic(V) ? db[1]|db[2] : db + d≠0 && count_ones(A&v)+count_ones(B&v)>diffmode(V) +end + +@pure function parityinterior(V::M,a,b,::Val{lim}=Val(false)) where {W,M<:Manifold{W},lim} A,B,Q,Z = symmetricmask(V,a,b); N = rank(V) - (diffcheck(V,A,B) || cga(V,A,B)) && (return false,Zero(UInt),false,Z) - p,C,t = parityregressive(Signature(V),A,complement(N,B,diffvars(V)),Val{true}()) - ind = indices(B,N); g = prod(V[ind]) - return t ? (p⊻parityright(0,sum(ind),count_ones(B)) ? -(g) : g) : g, C|Q, t, Z + diffcheck(V,A,B) && (return lim ? (Values{0,Tuple{UInt,Int}}(),Z) : (1,Zero(UInt),false,Z)) + gout,tout,G,diag = 0,false,count_ones(B),isdiag(V) + bas = diag ? Values((B,)) : indexbasis(grade(V),G) + g = if diag + gg = V[indices(B,N)] + Values((prod(isempty(gg) ? 1 : signbool.(gg)),)) + else + value(metrictensor(V,G))[bladeindex(grade(V),B)] + end + bs,bg = (),() + for i ∈ list(1,diag ? 1 : gdims(grade(V),G)) + if (!isnull(g[i])) && !(diffcheck2(V,A,bas[i]) || cga(V,A,bas[i])) + p,C,t = parityregressive(Signature(W),A,complement(N,bas[i],diffvars(V)),Val{true}()) + CQ,tout = C|Q,tout|t + if t + ggg = (p⊻parityright(0,sum(indices(bas[i],N)),G)) ? -(g[i]) : g[i] + if CQ ∉ bs + bs = (bs...,CQ) + bg = (bg...,(CQ,ggg)) + else + j = findfirst(x->x==CQ,bs) + bg = @inbounds (bg[1:j-1]...,(CQ,(bg[j][2]+ggg)),bg[j+1:length(bs)]...) + end + gout += ggg + end + end + end + if lim + return (isempty(bg) ? Values{0,Tuple{UInt,Int}}() : Values(bg)), Z + else + length(bs) > 1 && throw("this is a limited variant of interior product") + return @inbounds gout, isempty(bs) ? UInt(0) : bs[1], tout, Z + end end @pure function parityinner(V::Int,a::UInt,b::UInt) @@ -114,10 +149,173 @@ end parity(V,A,B) ? -1 : 1 end -@pure function parityinner(V::M,a::UInt,b::UInt) where M<:Manifold +@pure function parityinner(V::M,a::UInt,b::UInt) where {W,M<:Manifold{W}} A,B = symmetricmask(V,a,b) - g = abs(prod(V[indices(A&B,mdims(V))])) - parity(Signature(V),A,B) ? -(g) : g + C = A&B; G = count_ones(C) + g = if isdiag(V) || hasconformal(V) + gg = V[indices(C,mdims(V))] + abs(prod(isempty(gg) ? 1 : signbool.(gg))) + else + value(metrictensor(V,G))[bladeindex(mdims(V),C)] + end + parity(Signature(W),A,B) ? -(g) : g +end + +function parityseq(V,B) + l,out = length(B),false + (isdiag(V) || isone(l)) && (return 1) + for i ∈ list(2,length(B)) + out = out⊻parity(grade(V),B[i-1],B[i]) + end + return out ? -1 : 1 +end + +maxempty(x) = isempty(x) ? 0 : maximum(x) + +function paritygeometric(V,A::UInt,B::UInt) + #if iszero(B) + # Values(((A,1),)) + a,b = splitbasis(V,A),splitbasis(V,B) + ga,gb = maxempty(count_ones.(a)),maxempty(count_ones.(b)) + if (ga≤1 && gb≤1) ? count_ones(A) ≥ count_ones(B) : ga ≥ gb + paritygeometric(V,A,b) + else + paritygeometric(V,a,B) + end +end +function paritygeometric(V,A::UInt,b::Tuple) + i = length(b) + vals = if iszero(i) + Values((((UInt(0),1),(A,1)),)) + else + p = parityseq(V,b) + out = paritygeometricright(V,((UInt(0),p),(A,1)),b[1]) + isone(i) ? out : paritygeometricright(V,out,b,2) + end + combinebasis(V,vals) +end +function paritygeometric(V,a::Tuple,B::UInt) + i = length(a) + vals = if iszero(i) + Values((((UInt(0),1),(B,1)),)) + else + p = parityseq(V,a) + out = paritygeometricleft(V,a[i],((UInt(0),p),(B,1))) + isone(i) ? out : paritygeometricleft(V,a,out,i-1) + end + combinebasis(V,vals) +end +function paritygeometricright(V,a,b::Tuple,i) + out = vcat(paritygeometricright.(Ref(V),a,Ref(b[i]))...) + i==length(b) ? out : paritygeometricright(V,out,b,i+1) +end +function paritygeometricright(V,aeai,B::UInt) + ae,ai = @inbounds (aeai[1],aeai[2]) + Ae,Aeg = @inbounds (ae[1],ae[2]) + Ai,Aig = @inbounds (ai[1],ai[2]) + G = count_ones(B) + CCg = if isdiag(V) || hasconformal(V) + g,C,t,Z = interior(V,Ai,B) + Cg = parityclifford(G)⊻isodd(G*count_ones(Ai)) ? -Aig*g : Aig*g + t ? ((ae,(C,Cg)),) : () + else + Cg,Z = parityinterior(V,Ai,B,Val(true)) + g = parityclifford(G)⊻isodd(G*count_ones(Ai)) ? -Aig*Aeg : Aig*Aeg + ([((Ae,g),cg) for cg ∈ Cg]...,) + end + out = if iszero(Ai&B) + p = parity(grade(V),Ai⊻Ae,B) ? -Aeg : Aeg + (combinegeometric(V,(Ae⊻B,p),ai),CCg...) + else + CCg + end + isempty(out) ? Values{0,Tuple{Tuple{UInt,Int},Tuple{UInt,Int}}}() : Values(out) +end +function paritygeometricleft(V,a::Tuple,b,i) + out = vcat(paritygeometricleft.(Ref(V),Ref(a[i]),b)...) + isone(i) ? out : paritygeometricleft(V,a,out,i-1) +end +function paritygeometricleft(V,A::UInt,bebi) + be,bi = @inbounds (bebi[1],bebi[2]) + Be,Beg = @inbounds (be[1],be[2]) + Bi,Big = @inbounds (bi[1],bi[2]) + G = count_ones(A) + CCg = if isdiag(V) || hasconformal(V) + g,C,t,Z = interior(V,Bi,A) + Cg = parityreverse(G) ? -Big*g : Big*g + t ? ((be,(C,Cg)),) : () + else + Cg,Z = parityinterior(V,Bi,A,Val(true)) + g = parityreverse(G) ? -Big*Beg : Big*Beg + ([((Be,g),cg) for cg ∈ Cg]...,) + end + out = if iszero(A&Bi) + p = parity(grade(V),A,Be⊻Bi) ? -Beg : Beg + (combinegeometric(V,(A⊻Be,p),bi),CCg...) + else + CCg + end + isempty(out) ? Values{0,Tuple{Tuple{UInt,Int},Tuple{UInt,Int}}}() : Values(out) +end + +function splitbasis(V,ind::AbstractVector) + g = metrictensor(V) + f = [findall(.!iszero.(value(value(g)[i]))) for i ∈ 1:mdims(V)] + for i ∈ list(1,length(f)) + i ∉ f[i] && push!(f[i],i) + end + j = 2 + while j ≤ length(f) + t = false + for k ∈ list(1,j-1) + for q ∈ f[j] + if q ∈ f[k] + t = true + for p ∈ f[j] + (p ∉ f[k]) && push!(f[k],p) + end + popat!(f,j) + break + end + end + t && (break) + end + !t && (j += 1) + end + return getindex.(Ref(ind),f) +end + +function splitbasis(V,B) + iszero(B) && (return ()) + ind = indices(B,mdims(V)) + bs = ((UInt(1).<<(ind.-1))...,) + if isdiag(V) + return bs + else + f = splitbasis(V(ind...),ind) + ([|((UInt(1).<<(j.-1))...) for j ∈ f]...,) + end +end + +function combinebasis(V,vals) + result = () + for val ∈ vals + cg = combinegeometric(V,val) + !isnothing(cg) && (result = (result...,cg)) + end + isempty(result) ? Values{0,Tuple{UInt,Int}}() : Values(result) +end +function combinegeometric(V,ei) + e,i = @inbounds ei[1],ei[2] + E,Eg = @inbounds e[1],e[2] + I,Ig = @inbounds i[1],i[2] + iszero(E&I) ? (E⊻I,Eg*Ig) : nothing +end +function combinegeometric(V,e,i) + E,Eg = @inbounds e[1],e[2] + I,Ig = @inbounds i[1],i[2] + #p = parity(grade(V),E,I) ? -Eg*Ig : Eg*Ig + return ((E,Eg*Ig),(I,1)) end ### parity cache diff --git a/src/products.jl b/src/products.jl index a25edf4..4d0eef9 100644 --- a/src/products.jl +++ b/src/products.jl @@ -187,7 +187,7 @@ function generate_mutators(M,F,set_val,SUB,MUL) @inline function $(Symbol(:join,s))(V,m::$M,a::UInt,b::UInt,v::S) where {T<:$F,S<:$F,M} if v ≠ 0 && !diffcheck(V,a,b) A,B,Q,Z = symmetricmask(V,a,b) - val = $MUL(parityinner(V,A,B),v) + val = $MUL(parityinner(grade(V),A,B),v) if diffvars(V)≠0 !iszero(Z) && (T≠Any ? (return true) : (val *= getbasis(loworder(V),Z))) count_ones(Q)+order(val)>diffmode(V) && (return false) @@ -199,7 +199,7 @@ function generate_mutators(M,F,set_val,SUB,MUL) @inline function $(Symbol(:join,spre))(V,m::$M,a::UInt,b::UInt,v::S) where {T<:$F,S<:$F,M} if v ≠ 0 && !diffcheck(V,a,b) A,B,Q,Z = symmetricmask(V,a,b) - val = :($$MUL($(parityinner(V,A,B)),$v)) + val = :($$MUL($(parityinner(grade(V),A,B)),$v)) if diffvars(V)≠0 !iszero(Z) && (val = Expr(:call,:*,val,getbasis(loworder(V),Z))) val = :(h=$val;iszero(h)||$(count_ones(Q))+order(h)>$(diffmode(V)) ? 0 : h) @@ -211,28 +211,52 @@ function generate_mutators(M,F,set_val,SUB,MUL) @inline function $(Symbol(:geom,s))(V,m::$M,a::UInt,b::UInt,v::S) where {T<:$F,S<:$F,M} if v ≠ 0 && !diffcheck(V,a,b) A,B,Q,Z = symmetricmask(V,a,b) - pcc,bas,cc = (hasinf(V) && hasorigin(V)) ? conformal(V,A,B) : (false,A⊻B,false) - val = $MUL(parityinner(V,A,B),pcc ? $SUB(v) : v) - if istangent(V) - !iszero(Z) && (T≠Any ? (return true) : (val *= getbasis(loworder(V),Z))) - count_ones(Q)+order(val)>diffmode(V) && (return false) + if isdiag(V) + pcc,bas,cc = (hasinf(V) && hasorigin(V)) ? conformal(V,A,B) : (false,A⊻B,false) + g = parityinner(V,A,B) + val = $MUL(g,pcc ? $SUB(v) : v) + if istangent(V) + !iszero(Z) && (T≠Any ? (return true) : (val *= getbasis(loworder(V),Z))) + count_ones(Q)+order(val)>diffmode(V) && (return false) + end + $s(m,val,bas|Q,Val(mdims(V))) + cc && $s(m,hasinforigin(V,A,B) ? $SUB(val) : val,(conformalmask(V)⊻bas)|Q,Val(mdims(V))) + else + for (bas,g) ∈ paritygeometric(V,A,B) + val = $MUL(g,v) + if istangent(V) + !iszero(Z) && (T≠Any ? (return true) : (val *= getbasis(loworder(V),Z))) + count_ones(Q)+order(val)>diffmode(V) && (return false) + end + $s(m,val,bas|Q,Val(mdims(V))) + end end - $s(m,val,bas|Q,Val(mdims(V))) - cc && $s(m,hasinforigin(V,A,B) ? $SUB(val) : val,(conformalmask(V)⊻bas)|Q,Val(mdims(V))) end return false end @inline function $(Symbol(:geom,spre))(V,m::$M,a::UInt,b::UInt,v::S) where {T<:$F,S<:$F,M} if v ≠ 0 && !diffcheck(V,a,b) A,B,Q,Z = symmetricmask(V,a,b) - pcc,bas,cc = (hasinf(V) && hasorigin(V)) ? conformal(V,A,B) : (false,A⊻B,false) - val = :($$MUL($(parityinner(V,A,B)),$(pcc ? :($$SUB($v)) : v))) - if istangent(V) - !iszero(Z) && (val = Expr(:call,:*,val,getbasis(loworder(V),Z))) - val = :(h=$val;iszero(h)||$(count_ones(Q))+order(h)>$(diffmode(V)) ? 0 : h) + if isdiag(V) + pcc,bas,cc = (hasinf(V) && hasorigin(V)) ? conformal(V,A,B) : (false,A⊻B,false) + g = parityinner(V,A,B) + val = :($$MUL($g,$(pcc ? :($$SUB($v)) : v))) + if istangent(V) + !iszero(Z) && (val = Expr(:call,:*,val,getbasis(loworder(V),Z))) + val = :(h=$val;iszero(h)||$(count_ones(Q))+order(h)>$(diffmode(V)) ? 0 : h) + end + $spre(m,val,bas|Q,Val(mdims(V))) + cc && $spre(m,hasinforigin(V,A,B) ? :($$SUB($val)) : val,(conformalmask(V)⊻bas)|Q,Val(mdims(V))) + else + for (bas,g) ∈ paritygeometric(V,A,B) + val = :($$MUL($g,$v)) + if istangent(V) + !iszero(Z) && (val = Expr(:call,:*,val,getbasis(loworder(V),Z))) + val = :(h=$val;iszero(h)||$(count_ones(Q))+order(h)>$(diffmode(V)) ? 0 : h) + end + $spre(m,val,bas|Q,Val(mdims(V))) + end end - $spre(m,val,bas|Q,Val(mdims(V))) - cc && $spre(m,hasinforigin(V,A,B) ? :($$SUB($val)) : val,(conformalmask(V)⊻bas)|Q,Val(mdims(V))) end return false end @@ -246,28 +270,55 @@ function generate_mutators(M,F,set_val,SUB,MUL) @eval begin @inline function $(Symbol(prod,s))(V,m::$M,A::UInt,B::UInt,val::T) where {T,M} if val ≠ 0 - g,C,t,Z = $uct(V,A,B) - v = val - if istangent(V) && !iszero(Z) - T≠Any && (return true) - _,_,Q,_ = symmetricmask(V,A,B) - v *= getbasis(loworder(V),Z) - count_ones(Q)+order(v)>diffmode(V) && (return false) + if isdiag(V) + g,C,t,Z = $uct(V,A,B) + v = val + if istangent(V) && !iszero(Z) + T≠Any && (return true) + _,_,Q,_ = symmetricmask(V,A,B) + v *= getbasis(loworder(V),Z) + count_ones(Q)+order(v)>diffmode(V) && (return false) + end + t && $s(m,$MUL(g,v),C,Val(mdims(V))) + else + Cg,Z = $(Symbol(:parity,uct))(V,A,B,Val(true)) + v = val + if istangent(V) && !iszero(Z) + T≠Any && (return true) + _,_,Q,_ = symmetricmask(V,A,B) + v *= getbasis(loworder(V),Z) + count_ones(Q)+order(v)>diffmode(V) && (return false) + end + for (C,g) ∈ Cg + $s(m,$MUL(g,v),C,Val(mdims(V))) + end end - t && $s(m,$MUL(g,v),C,Val(mdims(V))) end return false end @inline function $(Symbol(prod,spre))(V,m::$M,A::UInt,B::UInt,val::T) where {T,M} if val ≠ 0 - g,C,t,Z = $uct(V,A,B) - v = val - if istangent(V) && !iszero(Z) - _,_,Q,_ = symmetricmask(V,A,B) - v = Expr(:call,:*,v,getbasis(loworder(V),Z)) - v = :(h=$v;iszero(h)||$(count_ones(Q))+order(h)>$(diffmode(V)) ? 0 : h) + if isdiag(V) + g,C,t,Z = $uct(V,A,B) + v = val + if istangent(V) && !iszero(Z) + _,_,Q,_ = symmetricmask(V,A,B) + v = Expr(:call,:*,v,getbasis(loworder(V),Z)) + v = :(h=$v;iszero(h)||$(count_ones(Q))+order(h)>$(diffmode(V)) ? 0 : h) + end + t && $spre(m,Expr(:call,$(QuoteNode(MUL)),g,v),C,Val(mdims(V))) + else + Cg,Z = $(Symbol(:parity,uct))(V,A,B,Val(true)) + v = val + if istangent(V) && !iszero(Z) + _,_,Q,_ = symmetricmask(V,A,B) + v = Expr(:call,:*,v,getbasis(loworder(V),Z)) + v = :(h=$v;iszero(h)||$(count_ones(Q))+order(h)>$(diffmode(V)) ? 0 : h) + end + for (C,g) ∈ Cg + $spre(m,Expr(:call,$(QuoteNode(MUL)),g,v),C,Val(mdims(V))) + end end - t && $spre(m,Expr(:call,$(QuoteNode(MUL)),g,v),C,Val(mdims(V))) end return false end @@ -277,7 +328,7 @@ function generate_mutators(M,F,set_val,SUB,MUL) end end -@inline exterbits(V,α,β) = diffvars(V)≠0 ? ((a,b)=symmetricmask(V,α,β);count_ones(a&b)==0) : count_ones(α&β)==0 +@inline exterbits(V,α,β) = diffvars(V)≠0 ? ((a,b)=symmetricmask(V,α,β);iszero(a&b)) : iszero(α&β) @inline exteraddmulti!(V,out,α,β,γ) = exterbits(V,α,β) && joinaddmulti!(V,out,α,β,γ) @inline exteraddblade!(V,out,α,β,γ) = exterbits(V,α,β) && joinaddblade!(V,out,α,β,γ) @inline exteraddspin!(V,out,α,β,γ) = exterbits(V,α,β) && joinaddspin!(V,out,α,β,γ) @@ -921,7 +972,7 @@ function generate_products(Field=Field,VEC=:mvec,MUL=:*,ADD=:+,SUB=:-,CONJ=:conj function ⟑(a::Single{V,G,A,T} where {G,A},b::Single{V,L,B,S} where {L,B}) where {V,T<:$Field,S<:$Field} ba,bb = basis(a),basis(b) v = derive_mul(V,UInt(ba),UInt(bb),a.v,b.v,$MUL) - Single(v,mul(ba,bb,v)) + v*mul(ba,bb,v) end ∧(a::$Field,b::$Field) = $MUL(a,b) ∧(a::F,b::B) where B<:TensorTerm{V,G} where {F<:$EF,V,G} = Single{V,G}(a,b) @@ -1255,10 +1306,10 @@ end end for side ∈ (:left,:right) - c,p = Symbol(:complement,side),Symbol(:parity,side) - h,pg,pn = Symbol(c,:hodge),Symbol(p,:hodge),Symbol(p,:null) + cc,p = Symbol(:complement,side),Symbol(:parity,side) + h,pg,pn = Symbol(cc,:hodge),Symbol(p,:hodge),Symbol(p,:null) pnp = :(Leibniz.$(Symbol(pn,:pre))) - for (c,p) ∈ ((c,p),(h,pg)) + for (c,p) ∈ ((cc,p),(h,pg)) @eval begin $c(z::Phasor) = $c(Couple(z)) function $c(z::Couple{V}) where V @@ -1269,30 +1320,31 @@ for side ∈ (:left,:right) @generated function $c(b::Chain{V,G,T}) where {V,G,T} isdyadic(V) && throw(error("Complement for dyadic tensors is undefined")) istangent(V) && (return :($$c(Multivector(b)))) + $(c≠h ? nothing : :((!isdiag(V)) && (return :($$cc(metric(b)))) )) SUB,VEC,MUL = subvec(b) if binomial(mdims(V),G)<(1<