Merge 0646f8c into 1bd9920

README.md

@@ -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
 
@@ -131,8 +131,11 @@ v₂₃
 ```
 Effort of computation depends on the sparsity. Then it is possible to assign the **quaternion** generators `i,j,k` with
 ```Julia
-julia> i,j,k = complementright.((-Λ(3).v1,-Λ(3).v2,-Λ(3).v3))
-(-1v₂₃, 1v₁₃, -1v₁₂)
+julia> i,j,k = hyperplanes(ℝ^3)
+3-element Array{SValue{⟨+++⟩,2,B,Int64} where B,1}:
+ -1v₂₃
+ 1v₁₃
+ -1v₁₂
 
 julia> @btime i^2, j^2, k^2, i*j*k
   158.925 ns (5 allocations: 112 bytes)

src/Grassmann.jl

@@ -17,6 +17,12 @@ include("forms.jl")
 include("symbolic.jl")
 include("generators.jl")
 
+## fundamentals
+
+export hyperplanes
+
+hyperplanes(V::VectorSpace{N}) where N = map(n->I*getbasis(V,1<<n),0:N-1)
+
 abstract type SubAlgebra{V} <: TensorAlgebra{V} end
 
 adjoint(G::A) where A<:SubAlgebra{V} where V = Λ(dual(V))

src/algebra.jl

@@ -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
@@ -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}())
@@ -137,6 +151,19 @@ end
 #*(a::Basis{V},b::MultiGrade{V}) where V = MultiGrade{V}(b.v,a*basis(b))
 #*(a::MultiGrade{V},b::MultiGrade{V}) where V = MultiGrade{V}(a.v*b.v,basis(a)*basis(b))
 
+for Value ∈ MSV
+    @eval begin
+        *(a::UniformScaling,b::$Value{V}) where V = V(a)*b
+        *(a::$Value{V},b::UniformScaling) where V = a*V(b)
+    end
+end
+for Blade ∈ MSB
+    @eval begin
+        *(a::UniformScaling,b::$Blade{T,V} where T) where V = V(a)*b
+        *(a::$Blade{T,V} where T,b::UniformScaling) where V = a*V(b)
+    end
+end
+
 ## exterior product
 
 export ∧, ∨
@@ -159,7 +186,10 @@ function ∧(a::X,b::Y) where {X<:TensorTerm{V},Y<:TensorTerm{V}} where V
     return SValue{V}(parity(x,y) ? -v : v,Basis{V}(A⊻B))
 end
 
-∧(a::X,b::Y) where {X<:TensorAlgebra,Y<:TensorAlgebra} = interop(∧,a,b)
+@inline ∧(a::X,b::Y) where {X<:TensorAlgebra,Y<:TensorAlgebra} = interop(∧,a,b)
+@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
+
 
 @inline exterior_product!(V::VectorSpace,out,α,β,γ) = (count_ones(α&β)==0) && joinaddmulti!(V,out,α,β,γ)
 
@@ -171,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
@@ -196,25 +228,28 @@ 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
         @eval $r(b::$Value) = value(b) ≠ 0 ? value(b) * $r(basis(b)) : zero(vectorspace(b))
     end
 end
 
+reverse(a::UniformScaling{Bool}) = UniformScaling(!a.λ)
+reverse(a::UniformScaling{T}) where T<:Field = UniformScaling(-a.λ)
+
 ## inner product: a ∨ ⋆(b)
 
 import LinearAlgebra: dot, ⋅
@@ -237,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)))
@@ -256,20 +291,52 @@ function ∨(a::X,b::Y) where {X<:TensorTerm{V},Y<:TensorTerm{V}} where V
     return SValue{V}(p ? -v : v,Basis{V}(C))
 end
 
-∨(a::X,b::Y) where {X<:TensorAlgebra,Y<:TensorAlgebra} = interop(∨,a,b)
+@inline ∨(a::X,b::Y) where {X<:TensorAlgebra,Y<:TensorAlgebra} = interop(∨,a,b)
+@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
     end
 end
 
+## cross product
+
+import LinearAlgebra: cross
+export ×
+
+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
 
 @eval begin
@@ -306,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)
@@ -552,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
@@ -583,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
@@ -617,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)...)

src/multivectors.jl

@@ -11,6 +11,12 @@ abstract type TensorMixed{T,V} <: TensorAlgebra{V} end
 
 import DirectSum: shift_indices, printindex, printindices, VTI
 
+## pseudoscalar
+
+using LinearAlgebra
+import LinearAlgebra: I
+export UniformScaling, I
+
 ## MultiBasis{N}
 
 struct Basis{V,G,B} <: TensorTerm{V,G}
@@ -382,6 +388,8 @@ end
 
 ## conversions
 
+@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)
     !(V⊆W) && throw(error("cannot convert from $(V) to $(W)"))

src/symbolic.jl

@@ -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}())