Skip to content

Commit

Permalink
separated definition of hodge and complement, crossprod methods
Browse files Browse the repository at this point in the history
  • Loading branch information
@chakravala
chakravala committed Mar 2, 2019
1 parent 9536ca1 commit 0646f8c
Show file tree
Hide file tree
Showing 4 changed files with 84 additions and 29 deletions.
4 changes: 2 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,9 @@ This package is a work in progress providing the necessary tools to work with ar

It is currently possible to do both high-performance numerical computations with `Grassmann` and it is also currently possible to use symbolic scalar coefficients when the `Reduce` package is also loaded (compatibility instructions at bottom).

Products available for high-performance and sparse computation include `∧,∨,⋅,*` (exterior, regressive, interior, geometric).
Products available for sparse & high-performance ops include `∧,∨,⋅,*,×` (exterior, regressive, interior, geometric, cross).

Some unary operations include `complementleft`,`complementright`,`reverse,`,`involve`,`conj`, and `adjoint`.
Some unary operations include `complementleft`,`complementright`,`reverse`,`involute`,`conj`, and `adjoint`.

#### Design, code generation

Expand Down
76 changes: 58 additions & 18 deletions src/algebra.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ for (op,set) ∈ ((:add,:(+=)),(:set,:(=)))
return m
end
@inline function $(Symbol(:join,sm))(m::MArray{Tuple{M},T,1,M},v::T,A::Basis{V},B::Basis{V}) where {V,T<:Field,M}
if v ≠ 0 !(hasdual(V) && hasdual(A) && hasdual(B))
if v ≠ 0 && !(hasdual(V) && hasdual(A) && hasdual(B))
$sm(m,parity(A,B) ? -(v) : v,bits(A) ⊻ bits(B),Dimension{ndims(V)}())
end
return m
Expand Down Expand Up @@ -71,6 +71,20 @@ for (op,set) ∈ ((:add,:(+=)),(:set,:(=)))
end
return m
end
@inline function $(Symbol(:cross,sm))(V::VectorSpace{N,D},m::MArray{Tuple{M},T,1,M},A::Bits,B::Bits,v::T) where {N,D,T<:Field,M}
if v ≠ 0
p,C,t = crossprod(A,B,V)
t && $sm(m,p ? -(v) : v,C,Dimension{N}())
end
return m
end
@inline function $(Symbol(:cross,sm))(m::MArray{Tuple{M},T,1,M},A::Basis{V},B::Basis{V},v::T) where {V,T<:Field,M}
if v ≠ 0
p,C,t = crossprod(bits(A),bits(B),V)
t && $sm(m,p ? -(v) : v,C,Dimension{N}())
end
return m
end
@inline function $(Symbol(:join,sb))(V::VectorSpace{N,D},m::MArray{Tuple{M},T,1,M},A::Bits,B::Bits,v::T) where {N,D,T<:Field,M}
if v ≠ 0 && !(hasdual(V) && isodd(A) && isodd(B))
$sb(m,parity(A,B,V) ? -(v) : v,A ⊻ B,Dimension{N}())
Expand Down Expand Up @@ -187,15 +201,17 @@ end

## complement

export complementleft, complementright
export complementleft, complementright, complementhodge

@pure complement(N::Int,B::UInt) = (~B)&(one(Bits)<<N-1)

@pure parityright(V::Bits,B::Bits) = parityright(count_ones(V&B),sum(indices(B)),count_ones(B))
@pure parityright(V::Int,B,G,N=nothing) = isodd(V+B+Int((G+1)*G/2))
@pure parityhodge(V::Bits,B::Bits) = parityhodge(count_ones(V&B),sum(indices(B)),count_ones(B))

@pure parityright(V::Int,B,G,N=nothing) = isodd(B+Int((G+1)*G/2))
@pure parityleft(V::Int,B,G,N) = (isodd(G) && iseven(N)) ⊻ parityright(V,B,G,N)
@pure parityhodge(V::Int,B,G,N=nothing) = isodd(V)⊻parityright(V,B,G)

for side ∈ (:left,:right)
for side ∈ (:left,:right,:hodge)
c = Symbol(:complement,side)
p = Symbol(:parity,side)
@eval begin
Expand All @@ -212,18 +228,18 @@ for side ∈ (:left,:right)
end

export ⋆
const ⋆ = complementright
const ⋆ = complementhodge

## reverse

import Base: reverse, conj
export involve
export involute

@pure parityreverse(G) = isodd(Int((G-1)*G/2))
@pure parityinvolve(G) = isodd(G)
@pure parityconj(G) = parityreverse(G)⊻parityinvolve(G)
@pure parityinvolute(G) = isodd(G)
@pure parityconj(G) = parityreverse(G)⊻parityinvolute(G)

for r ∈ (:reverse,:involve,:conj)
for r ∈ (:reverse,:involute,:conj)
p = Symbol(:parity,r)
@eval @pure $r(b::Basis{V,G,B}) where {V,G,B} =$p(G) ? SValue{V}(-value(b),b) : b
for Value ∈ MSV
Expand Down Expand Up @@ -256,7 +272,7 @@ end
@pure function interior_calc(N,M,S,A,B)
γ = complement(N,B)
p,C,t = regressive_calc(N,M,S,A,γ)
return t ? p⊻parityright(S,B) : p, C, t
return t ? p⊻parityhodge(S,B) : p, C, t
end

## regressive product: (L = grade(a) + grade(b); (-1)^(L*(L-ndims(V)))*⋆(⋆(a)∧⋆(b)))
Expand All @@ -279,13 +295,12 @@ end
@inline ∨(a::TensorAlgebra{V},b::UniformScaling{T}) where {V,T<:Field} = a∨V(b)
@inline ∨(a::UniformScaling{T},b::TensorAlgebra{V}) where {V,T<:Field} = V(a)∨b


@pure function regressive_calc(N,M,S,A,B)
α,β = complement(N,A),complement(N,B)
if (count_ones(α&β)==0) && !(M ∈ (1,3,5,7,9,11) && isodd(α) && isodd(β))
C = α ⊻ β
L = count_ones(A)+count_ones(B)
pa,pb,pc = parityright(S,A),parityright(S,B),parityright(S,C)
pa,pb,pc = parityhodge(S,A),parityhodge(S,B),parityhodge(S,C)
return isodd(L*(L-N))⊻pa⊻pb⊻parity(N,S,α,β)⊻pc, complement(N,C), true
else
return false, zero(Bits), false
Expand All @@ -297,7 +312,30 @@ end
import LinearAlgebra: cross
export ×

cross(a::TensorAlgebra{V},b::TensorAlgebra{V}) where V = I*(a∧b)
cross(a::TensorAlgebra{V},b::TensorAlgebra{V}) where V = ⋆(a∧b)

@pure function cross(a::Basis{V},b::Basis{V}) where V
p,C,t = crossprod(a,b)
!t && (return zero(V))
d = Basis{V}(C)
return p ? SValue{V}(-1,d) : d
end

function cross(a::X,b::Y) where {X<:TensorTerm{V},Y<:TensorTerm{V}} where V
p,C,t = crossprod(bits(basis(a)),bits(basis(b)),V)
!t && (return zero(V))
v = value(a)*value(b)
return SValue{V}(p ? -v : v,Basis{V}(C))
end

@pure function crossprod_calc(N,M,S,A,B)
if (count_ones(A&B)==0) && !(M ∈ (1,3,5,7,9,11) && isodd(A) && isodd(B))
C = A ⊻ B
return parity(N,S,A,B)⊻parityhodge(S,C), complement(N,C), true
else
return false, zero(Bits), false
end
end

### parity cache

Expand Down Expand Up @@ -335,7 +373,8 @@ end

### parity cache 2

for (par,T) ∈ ((:interior,Tuple{Bool,Bits,Bool}),(:regressive,Tuple{Bool,Bits,Bool}))
for par ∈ (:interior,:regressive,:crossprod)
T = Tuple{Bool,Bits,Bool}
extra = Symbol(par,:_extra)
cache = Symbol(par,:_cache)
calc = Symbol(par,:_calc)
Expand Down Expand Up @@ -581,7 +620,7 @@ function generate_product_algebra(Field=Field,MUL=:*,ADD=:+,SUB=:-,VEC=:mvec,CON
end
end
end
for side ∈ (:left,:right)
for side ∈ (:left,:right,:hodge)
c = Symbol(:complement,side)
p = Symbol(:parity,side)
for Blade ∈ MSB
Expand Down Expand Up @@ -612,7 +651,7 @@ function generate_product_algebra(Field=Field,MUL=:*,ADD=:+,SUB=:-,VEC=:mvec,CON
end
end
end
for reverse ∈ (:reverse,:involve,:conj)
for reverse ∈ (:reverse,:involute,:conj)
p = Symbol(:parity,reverse)
for Blade ∈ MSB
@eval begin
Expand Down Expand Up @@ -646,7 +685,8 @@ function generate_product_algebra(Field=Field,MUL=:*,ADD=:+,SUB=:-,VEC=:mvec,CON
end

for (op,product!) ∈ ((:∧,:exterior_product!),(:*,:joinaddmulti!),
(:∨,:meetaddmulti!),(:dot,:skewaddmulti!))
(:∨,:meetaddmulti!),(:dot,:skewaddmulti!),
(:cross,:crossaddmulti!))
@eval begin
function $op(a::MultiVector{T,V},b::Basis{V,G}) where {T<:$Field,V,G}
$(insert_expr((:N,:t,:out,:bs,:bn),VEC)...)
Expand Down
2 changes: 1 addition & 1 deletion src/multivectors.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -388,7 +388,7 @@ end

## conversions

@inline (V::VectorSpace)(s::UniformScaling{T}) where T = SValue{V}(T<:Bool ? (s.λ ? -(one(T)) : one(T)) : s.λ,getbasis(V,(one(Bits)<<ndims(V))-1))
@inline (V::VectorSpace)(s::UniformScaling{T}) where T = SValue{V}(T<:Bool ? (s.λ ? one(Int) : -(one(Int))) : s.λ,getbasis(V,(one(T)<<ndims(V))-1))

@pure function (W::VectorSpace)(b::Basis{V}) where V
V==W && (return b)
Expand Down
31 changes: 23 additions & 8 deletions src/symbolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,33 +32,48 @@ for (op,set) ∈ ((:add,:(+=)),(:set,:(=)))
return m
end
@inline function $(Symbol(:meet,sm))(V::VectorSpace{N,D},m::SizedArray{Tuple{M},T,1,1},A::Bits,B::Bits,v::T) where {N,D,T<:SymField,M}
if v ≠ 0
p,C,t = regressive(N,value(V),A,B)
if v ≠ 0 && !(hasdual(V) && hasdual(A) && hasdual(B))
p,C,t = regressive(A,B,V)
t && $sm(m,p ? $Sym.:-(v) : v,C,Dimension{N}())
end
return m
end
@inline function $(Symbol(:meet,sm))(m::SizedArray{Tuple{M},T,1,1},A::Basis{V},B::Basis{V},v::T) where {V,T<:SymField,M}
if v ≠ 0
p,C,t = regressive(N,value(V),bits(A),bits(B))
if v ≠ 0 && !(hasdual(V) && hasdual(A) && hasdual(B))
p,C,t = regressive(bits(A),bits(B),V)
t && $sm(m,p ? $Sym.:-(v) : v,C,Dimension{N}())
end
return m
end
@inline function $(Symbol(:skew,sm))(V::VectorSpace{N,D},m::SizedArray{Tuple{M},T,1,1},A::Bits,B::Bits,v::T) where {N,D,T<:SymField,M}
if v ≠ 0
p,C,t = interior(N,value(V),A,B)
if v ≠ 0 && !(hasdual(V) && hasdual(A) && hasdual(B))
p,C,t = interior(A,B,V)
t && $sm(m,p ? $Sym.:-(v) : v,C,Dimension{N}())
end
return m
end
@inline function $(Symbol(:skew,sm))(m::SizedArray{Tuple{M},T,1,1},A::Basis{V},B::Basis{V},v::T) where {V,T<:SymField,M}
if v ≠ 0
p,C,t = interior(N,value(V),bits(A),bits(B))
if v ≠ 0 && !(hasdual(V) && hasdual(A) && hasdual(B))
p,C,t = interior(bits(A),bits(B),V)
t && $sm(m,p ? $Sym.:-(v) : v,C,Dimension{N}())
end
return m
end
@inline function $(Symbol(:cross,sm))(V::VectorSpace{N,D},m::SizedArray{Tuple{M},T,1,1},A::Bits,B::Bits,v::T) where {N,D,T<:SymField,M}
if v ≠ 0 && !(hasdual(V) && hasdual(A) && hasdual(B))
p,C,t = crossprod(A,B,V)
t && $sm(m,p ? $Sym.:-(v) : v,C,Dimension{N}())
end
return m
end
@inline function $(Symbol(:cross,sm))(m::SizedArray{Tuple{M},T,1,1},A::Basis{V},B::Basis{V},v::T) where {V,T<:SymField,M}
if v ≠ 0 && !(hasdual(V) && hasdual(A) && hasdual(B))
p,C,t = crossprod(bits(A),bits(B),V)
t && $sm(m,p ? $Sym.:-(v) : v,C,Dimension{N}())
end
return m
end

@inline function $(Symbol(:join,sb))(V::VectorSpace{N,D},m::SizedArray{Tuple{M},T,1,1},A::Bits,B::Bits,v::T) where {N,D,T<:SymField,M}
if v ≠ 0 && !(hasdual(V) && isodd(A) && isodd(B))
$sb(m,parity(A,B,V) ? $Sym.:-(v) : v,A ⊻ B,Dimension{N}())
Expand Down

0 comments on commit 0646f8c

Please sign in to comment.