Skip to content

Commit

Permalink
refined some source code
Browse files Browse the repository at this point in the history
  • Loading branch information
@chakravala
chakravala committed May 8, 2024
1 parent 0e2bf7b commit 5a4c30a
Show file tree
Hide file tree
Showing 6 changed files with 563 additions and 414 deletions.
7 changes: 0 additions & 7 deletions src/Grassmann.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,13 +59,6 @@ include("forms.jl")

export cayley, hyperplanes, points, TensorAlgebra

cayley(V) = cayley(Vector(Λ(V).b))
cayley(V,op) = cayley(Vector(Λ(V).b),op)
cayley(b::AbstractVector) = cayley(b,b)
cayley(b::AbstractVector,op) = cayley(b,b,op)
cayley(a::AbstractVector,b::AbstractVector) = a*transpose(b)
cayley(a::AbstractVector,b::AbstractVector,op) = TensorAlgebra{Manifold(Manifold(eltype(a)))}[op(x,y) for x ∈ a, y ∈ b]

@pure hyperplanes(V::Manifold) = map(n->UniformScaling{Bool}(false)*getbasis(V,1<<n),0:rank(V)-1-diffvars(V))

for M ∈ (:Signature,:DiagonalForm)
Expand Down
54 changes: 24 additions & 30 deletions src/algebra.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ end end
Geometric algebraic product: ω⊖η = (-1)ᵖdet(ω∩η)⊗(Λ(ω⊖η)∪L(ω⊕η))
"""
@pure times(a::Submanifold{V},b::Submanifold{V}) where V = mul(a,b)
@pure ⟑(a::Submanifold{V},b::Submanifold{V}) where V = mul(a,b)
*(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
Expand All @@ -46,12 +46,12 @@ Geometric algebraic product: ω⊖η = (-1)ᵖdet(ω∩η)⊗(Λ(ω⊖η)∪L(ω
return cc ? (v=value(out);out+Single{V}(hasinforigin(V,A,B) ? -(v) : v,getbasis(V,conformalmask(V)⊻UInt(d)))) : out
end

function times(a::Single{V},b::Submanifold{V}) where V
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)
end
function times(a::Submanifold{V},b::Single{V}) where V
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)
Expand Down Expand Up @@ -251,13 +251,6 @@ export ⊘, sandwich, pseudosandwich, antisandwich
@generated ⊘(a::TensorGraded{V,G},b::Couple{V}) where {V,G} = product_sandwich(a,b)
@generated ⊘(a::TensorGraded{V,G},b::PseudoCouple{V}) where {V,G} = product_sandwich(a,b)
@generated ⊘(a::TensorGraded{V,G},b::TensorGraded{V,L}) where {V,G,L} = product_sandwich(a,b)
#=for t ∈ (:Spinor,:AntiSpinor)
@eval quote
@generated ⊘(a::Spinor{V,G},b::$$t{V}) where {V,G} = product_sandwich(a,b)
@generated ⊘(a::AntiSpinor{V,G},b::$$t{V}) where {V,G} = product_sandwich(a,b)
@generated ⊘(a::Multivector{V,G},b::$$t{V}) where {V,G} = product_sandwich(a,b)
end
end=#

@doc """
⊘(ω::TensorAlgebra,η::TensorAlgebra)
Expand Down Expand Up @@ -427,36 +420,37 @@ for (nv,d) ∈ ((:inv,:/),(:inv_rat,://))
end
end

/(a::TensorTerm{V,0},b::Couple{V,B,S}) where {V,B,S} = (T = promote_type(valuetype(a),S); Couple{V,B}(value(a)*inv(T(b.v))))
/(a::Couple{V,B},b::TensorTerm{V,0}) where {V,B} = Couple{V,B}(Complex(a.v.re/value(b),a.v.im/value(b)))
/(a::TensorTerm{V,0},b::Couple{V,B,S}) where {V,B,S} = a*inv(b)
/(a::Couple{V,B},b::TensorTerm{V,0}) where {V,B} = Couple{V,B}(realvalue(a)/value(b),imagvalue(a)/value(b))

function /(a::Couple{V,B,T}, b::Couple{V,B,T}) where {V,B,T<:Real}
are,aim = reim(a); bre,bim = reim(b)
B2 = value(abs2_inv(B))
Couple{V,B}(if abs(bre) <= abs(bim)
(rout,iout) = if abs(bre) <= abs(bim)
if isinf(bre) && isinf(bim)
r = sign(bre)/sign(bim)
else
r = bre / bim
end
den = bim*B2 + r*bre
Complex((are*r + aim*B2)/den, (aim*r - are)/den)
((are*r + aim*B2)/den, (aim*r - are)/den)
else
if isinf(bre) && isinf(bim)
r = sign(bim)/sign(bre)
else
r = bim / bre
end
den = bre + (r*bim)*B2
Complex((are + (aim*r)*B2)/den, (aim - are*r)/den)
end)
((are + (aim*r)*B2)/den, (aim - are*r)/den)
end
Chain{V,B}(rout,iout)
end

inv(z::Couple{V,B,<:Union{Float16,Float32}}) where {V,B} =
(w = inv(widen(z)); Couple{V,B}(oftype(z.v,w.v)))
inv(z::Couple{V,B,T}) where {V,B,T<:Union{Float16,Float32}} =
(w = inv(widen(z)); Couple{V,B}(convert(T,realvalue(w)),convert(T,imagvalue(w))))

/(z::Couple{V,B,T}, w::Couple{V,B,T}) where {V,B,T<:Union{Float16,Float32}} =
(w = widen(z)*inv(widen(w)); Couple{V,B}(oftype(z.v, w.v)))
(w = widen(z)*inv(widen(w)); Couple{V,B}(convert(T,realvalue(w)),convert(T,imagvalue(w))))

# robust complex division for double precision
# variables are scaled & unscaled to avoid over/underflow, if necessary
Expand All @@ -477,7 +471,7 @@ function /(z::Couple{V,B,Float64}, w::Couple{V,B,Float64}) where {V,B}
else
p,q = cdiv(a,b,c,d,e)
end
return Couple{V,B}(ComplexF64(p,q))
return Couple{V,B,Float64}(p,q)
end

# sub-functionality for /(z::ComplexF64, w::ComplexF64)
Expand Down Expand Up @@ -622,16 +616,16 @@ adder(a,b,op=:+) = adder(typeof(a),typeof(b),op)
:(Single{V,L}($bop(value(a),value(b)),basis(a)))
elseif !istangent(V) && !hasconformal(V) && L == 0 &&
valuetype(a)<:Real && valuetype(b)<:Real
:(Couple{V,basis(b)}(Complex(value(a),$bop(value(b)))))
:(Couple{V,basis(b)}(value(a),$bop(value(b))))
elseif !istangent(V) && !hasconformal(V) && G == 0 &&
valuetype(a)<:Real && valuetype(b)<:Real
:(Couple{V,basis(a)}(Complex($bop(value(b)),value(a))))
:(Couple{V,basis(a)}($bop(value(b)),value(a)))
elseif !istangent(V) && !hasconformal(V) && L == grade(V) &&
valuetype(a)<:Real && valuetype(b)<:Real
:(PseudoCouple{V,basis(b)}(Complex($bop(value(b)),value(a))))
:(PseudoCouple{V,basis(b)}($bop(value(b)),value(a)))
elseif !istangent(V) && !hasconformal(V) && G == grade(V) &&
valuetype(a)<:Real && valuetype(b)<:Real
:(PseudoCouple{V,basis(a)}(Complex(value(a),$bop(value(b)))))
:(PseudoCouple{V,basis(a)}(value(a),$bop(value(b))))
elseif L == G
if binomial(mdims(V),G)<(1<<cache_limit)
$(insert_expr((:N,:ib,:t),:mvec)...)
Expand Down Expand Up @@ -734,23 +728,23 @@ adder(a,b,op=:+) = adder(typeof(a),typeof(b),op)
if !istangent(V) && !hasconformal(V) && L == 0 && G == mdims(V) &&
valuetype(a)<:Real && valuetype(b)<:Real
if swap
:(Couple{V,basis(V)}(Complex($right(value(a)),b.v[1])))
:(Couple{V,basis(V)}($right(value(a)),b.v[1]))
else
:(Couple{V,basis(V)}(Complex(value(a),$right(b.v[1]))))
:(Couple{V,basis(V)}(value(a),$right(b.v[1])))
end
elseif !istangent(V) && !hasconformal(V) && G == 0 &&
valuetype(a)<:Real && valuetype(b)<:Real
if swap
:(Couple{V,basis(a)}(Complex(b.v[1],$right(value(a)))))
:(Couple{V,basis(a)}(b.v[1],$right(value(a))))
else
:(Couple{V,basis(a)}(Complex($right(b.v[1]),value(a))))
:(Couple{V,basis(a)}($right(b.v[1]),value(a)))
end
elseif !istangent(V) && !hasconformal(V) && G == grade(V) &&
valuetype(a)<:Real && valuetype(b)<:Real
if swap
:(PseudoCouple{V,basis(a)}(Complex($right(value(a)),b.v[1])))
:(PseudoCouple{V,basis(a)}($right(value(a)),b.v[1]))
else
:(PseudoCouple{V,basis(a)}(Complex(value(a),$right(b.v[1]))))
:(PseudoCouple{V,basis(a)}(value(a),$right(b.v[1])))
end
elseif iseven(L) && iseven(G)
if mdims(V)-1<cache_limit
Expand Down
57 changes: 16 additions & 41 deletions src/composite.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ end
Base.expm1(t::Phasor{V}) where V = exp(t)-One(V)
function Base.exp(t::Phasor{V,B}) where {V,B}
z = exp(angle(t))
Phasor{V,B}(exp(radius(t)+z.v.re),z.v.im)
Phasor{V,B}(exp(radius(t)+realvalue(z)),imagvalue(z))
end

function Base.exp(t::Couple{V,B}) where {V,B}
Expand All @@ -166,7 +166,8 @@ function Base.expm1(t::PseudoCouple{V,B}) where {V,B}
end
function Base.exp(t::PseudoCouple{V,B}) where {V,B}
if isscalar(B)
PseudoCouple{V,B}(value(exp(Couple{V,Submanifold(V)}(value(t)))))
out = exp(Couple{V,Submanifold(V)}(realvalue(t),imagvalue(t)))
PseudoCouple{V,B}(realvalue(out),imagvalue(out))
else
exp(multispin(t))
end
Expand Down Expand Up @@ -404,26 +405,29 @@ end
end
end

Base.log(A::Chain{V,G,<:Chain{V,G}}) where {V,G} = Chain{V,G,Chain{V,G}}(log(Matrix(A)))
Base.log(t::TensorTerm) = log(Couple(t))
Base.log(t::Phasor) = (r=radius(t); log(r)+angle(t))
Base.log1p(t::Phasor{V}) where V = log(One(V)+t)
Base.log(t::Couple{V,B}) where {V,B} = value(B*B)==-1 ? Couple{V,B}(log(t.v)) : log(radius(t))+angle(t)
Base.log1p(t::Couple{V,B}) where {V,B} = value(B*B)==-1 ? Couple{V,B}(log1p(t.v)) : log(One(V)+t)
Base.log(t::Couple{V,B}) where {V,B} = value(B*B)==-1 ? Couple{V,B}(log(Complex(t))) : log(radius(t))+angle(t)
Base.log1p(t::Couple{V,B}) where {V,B} = value(B*B)==-1 ? Couple{V,B}(log1p(Complex(t))) : log(One(V)+t)
Base.log(t::Quaternion{V}) where V = iszero(metric(V)) ? log(radius(t))+angle(t,r) : qlog((t-One(V))/(t+One(V)))
Base.log1p(t::Quaternion{V}) where V = iszero(metric(V)) ? log(One(V)+t) : qlog(t/(t+2))
@inline Base.log(t::T) where T<:TensorAlgebra{V} where V = qlog((t-One(V))/(t+One(V)))
@inline Base.log1p(t::T) where T<:TensorAlgebra = qlog(t/(t+2))

function Base.log(t::PseudoCouple{V,B}) where {V,B}
if isscalar(B)
PseudoCouple{V,B}(value(log(Couple{V,Submanifold(V)}(value(t)))))
out = log(Couple{V,Submanifold(V)}(realvalue(t),imagvalue(t)))
PseudoCouple{V,B}(realvalue(out),imagvalue(out))
else
log(multispin(t))
end
end
function Base.log1p(t::PseudoCouple{V,B}) where {V,B}
if isscalar(B)
PseudoCouple{V,B}(value(log1p(Couple{V,Submanifold(V)}(value(t)))))
out = log1p(Couple{V,Submanifold(V)}(realvalue(t),imagvalue(t)))
PseudoCouple{V,B}(realvalue(out),imagvalue(out))
else
log1p(multispin(t))
end
Expand Down Expand Up @@ -480,7 +484,7 @@ for (qrt,n) ∈ ((:sqrt,2),(:cbrt,3))
iszero(metric(V)) ? $qrt(radius(t))*exp(angle(t)/$n) : exp(log(t)/$n)
end
@inline function Base.$qrt(t::Couple{V,B}) where {V,B}
value(B*B)==-1 ? Couple{V,B}($qrt(t.v)) :
value(B*B)==-1 ? Couple{V,B}($qrt(Complex(t))) :
$qrt(radius(t))*exp(angle(t)/$n)
end
@inline Base.$qrt(t::Phasor) = Phasor($qrt(radius(t)),angle(t)/$n)
Expand Down Expand Up @@ -705,11 +709,10 @@ end
end=#

function Base.angle(z::Couple{V,B}) where {V,B}
c = value(z)
if value(B^2) == -1
atan(imag(c),real(c))*B
atan(imagvalue(z),realvalue(z))*B
elseif value(B^2) == 1
atanh(imag(c),real(c))*B
atanh(imagvalue(z),realvalue(z))*B
else
error("Unsupported trigonometric angle")
end
Expand Down Expand Up @@ -907,9 +910,11 @@ 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::Chain{M,1,<:Chain{W,1}},v::Chain{V,1,<:Chain}) where {M,W,V} = Chain{V,1}(t.\value(v))
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.:\(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
Base.inv(t::Chain{V,1,<:Chain{W,1}}) where {W,V} = inv(value(t))
grad(t::Chain{V,1,<:Chain{W,1}}) where {V,W} = grad(value(t))

Expand Down Expand Up @@ -1186,33 +1191,3 @@ 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})

# Dyadic

export operator, gradedoperator, evenoperator, oddoperator

@generated function operator(t::TensorAlgebra{V},::Val{G}=Val(1)) where {V,G}
N = mdims(V)
r,b = binomsum(N,G),binomial(N,G)
bas = Λ(V).b[list(r+1,r+b)]
:(Chain{V,G}($bas .⊘ Ref(t)))
end
operator(t::TensorAlgebra,G::Int) = operator(t,Val(G))
gradedoperator(t::TensorAlgebra{V}) where V = Multivector{V}(Λ(V).b .⊘ Ref(t))

@generated function operator(fun,V,::Val{G}=Val(1)) where G
N = mdims(V)
r,b = binomsum(N,G),binomial(N,G)
bas = Λ(V()).b[list(r+1,r+b)]
:(Chain{V,G}(fun.($bas)))
end
operator(fun,V,G::Int) = operator(fun,V,Val(G))
gradedoperator(fun,V) = Multivector{V}(fun.(Λ(V).b))

@pure function evenbasis(V,even=true)
N = mdims(V)
r,b = binomsum_set(N),binomial_set(N)
vcat([Λ(V).b[list(r[g]+1,r[g]+b[g])] for g ∈ evens(even ? 1 : 2,N+1)]...)
end
evenoperator(t::TensorAlgebra{V}) where V = Spinor{V}(evenbasis(V) .⊘ Ref(t))
oddoperator(t::TensorAlgebra{V}) where V = AntiSpinor{V}(evenbasis(V,false) .⊘ Ref(t))
Loading

0 comments on commit 5a4c30a

Please sign in to comment.