From e645888e51768f46c0f70d7b4825c612a5b2c94c Mon Sep 17 00:00:00 2001 From: Michael Reed <18372368+chakravala@users.noreply.github.com> Date: Fri, 28 Jun 2019 01:12:47 -0400 Subject: [PATCH] created imag,real,even,odd,angular,radial methods --- src/algebra.jl | 22 ++++++----- src/multivectors.jl | 5 ++- src/parity.jl | 93 +++++++++++++++++++++++++++++++++++++++++++++ src/utilities.jl | 4 +- 4 files changed, 112 insertions(+), 12 deletions(-) diff --git a/src/algebra.jl b/src/algebra.jl index fcef0b0..efdd08b 100644 --- a/src/algebra.jl +++ b/src/algebra.jl @@ -175,7 +175,8 @@ for Blade ∈ MSB end end -export ⊛ +export ⊛, ⊖ +const ⊖ = * ## exterior product @@ -300,7 +301,9 @@ end ## sandwich product ->>>(x::TensorAlgebra{V},y::TensorAlgebra{V}) where V = x * y * ~x +>>>(x::TensorAlgebra{V},y::TensorAlgebra{V}) where V = (~x) * y * x +const ⊘ = >>> +export ⊘ ### Product Algebra Constructor @@ -980,9 +983,10 @@ for (nv,d) ∈ ((:inv,:/),(:inv_rat,://)) rm = ~m d = rm*m fd = norm(d) - for k ∈ 0:ndims(V) - dk = d(k) - norm(dk) ≈ fd && (return $d(rm,dk)) + sd = scalar(d) + value(sd) ≈ fd && (return $d(rm,sd)) + for k ∈ 1:ndims(V) + norm(d[k]) ≈ fd && (return $d(rm,d(k))) end throw(error("inv($m) is undefined")) end @@ -1006,11 +1010,11 @@ for op ∈ (:div,:rem,:mod,:mod1,:fld,:fld1,:cld,:ldexp) @eval Base.$op(b::$Value{V,G,B,T},m) where {V,G,B,T} = $Value{V,G,B}($op(value(b),m)) end for Blade ∈ MSB - @eval Base.$op(a::$Blade{T,V,G},m) where {T,V,G} = $Blade{T,V,G}($op.(value(a),m)) + @eval Base.$op(a::$Blade{T,V,G},m::S) where {T,V,G,S} = $Blade{promote_type(T,S),V,G}($op.(value(a),m)) end @eval begin Base.$op(a::Basis{V,G},m) where {V,G} = Basis{V,G}($op(value(a),m)) - Base.$op(a::MultiVector{T,V},m) where {T,V} = MultiVector{T,V}($op.(value(a),m)) + Base.$op(a::MultiVector{T,V},m::S) where {T,V,S} = MultiVector{promote_type(T,S),V}($op.(value(a),m)) end end for op ∈ (:mod2pi,:rem2pi,:rad2deg,:deg2rad) @@ -1018,11 +1022,11 @@ for op ∈ (:mod2pi,:rem2pi,:rad2deg,:deg2rad) @eval Base.$op(b::$Value{V,G,B,T}) where {V,G,B,T} = $Value{V,G,B}($op(value(b))) end for Blade ∈ MSB - @eval Base.$op(a::$Blade{T,V,G}) where {T,V,G} = $Blade{T,V,G}($op.(value(a))) + @eval Base.$op(a::$Blade{T,V,G}) where {T,V,G} = $Blade{promote_type(T,Float64),V,G}($op.(value(a))) end @eval begin Base.$op(a::Basis{V,G}) where {V,G} = Basis{V,G}($op(value(a))) - Base.$op(a::MultiVector{T,V}) where {T,V} = MultiVector{T,V}($op.(value(a))) + Base.$op(a::MultiVector{T,V}) where {T,V} = MultiVector{promote_type(T,Float64),V}($op.(value(a))) end end for Value ∈ MSV diff --git a/src/multivectors.jl b/src/multivectors.jl index 069251c..5886fb9 100644 --- a/src/multivectors.jl +++ b/src/multivectors.jl @@ -213,7 +213,8 @@ setindex!(m::MultiVector{T},k::T,i::Int,j::Int) where T = (m[i][j] = k) Base.firstindex(m::MultiVector) = 0 Base.lastindex(m::MultiVector{T,V} where T) where V = ndims(V) -function (m::MultiVector{T,V})(g::Int,::Type{B}=SBlade) where {T,V,B} +(m::MultiVector{T,V})(g::Int,b::Type{B}=SBlade) where {T,V,B} = m(Dim{g}(),b) +function (m::MultiVector{T,V})(::Dim{g},::Type{B}=SBlade) where {T,V,g,B} B ≠ SBlade ? MBlade{T,V,g}(m[g]) : SBlade{T,V,g}(m[g]) end function (m::MultiVector{T,V})(g::Int,i::Int,::Type{B}=SValue) where {T,V,B} @@ -366,7 +367,7 @@ end @inline vector(t::T) where T<:TensorTerm{V,1} where V = t @inline vector(t::T) where T<:TensorTerm{V} where V = zero(V) -@inline vector(t::MultiVector{T,V}) where {T,V} = SBlade{T,V,0}(t[1]) +@inline vector(t::MultiVector{T,V}) where {T,V} = SBlade{T,V,1}(t[1]) for Blade ∈ MSB @eval begin @inline vector(t::$Blade{T,V,1} where {T,V}) = t diff --git a/src/parity.jl b/src/parity.jl index 4a105d3..5d46e26 100644 --- a/src/parity.jl +++ b/src/parity.jl @@ -197,3 +197,96 @@ for par ∈ (:conformal,:regressive,:interior,:crossprod) @eval @pure $par(a::Basis{V,G,B},b::Basis{V,L,C}) where {V,G,B,L,C} = $par(bits(a),bits(b),V) end +import Base: signbit, imag, real +export odd, even, angular, radial + +@pure signbit(V::T) where T<:VectorSpace{N} where N = (ib=indexbasis(N); parity.(ib,ib,Ref(V))) +@pure signbit(V::T,G) where T<:VectorSpace{N} where N = (ib=indexbasis(N,G); parity.(ib,ib,Ref(V))) +@pure angular(V::T) where T<:VectorSpace = SVector(findall(signbit(V))...) +@pure radial(V::T) where T<:VectorSpace = SVector(findall(.!signbit(V))...) +@pure angular(V::T,G) where T<:VectorSpace = findall(signbit(V,G)) +@pure radial(V::T,G) where T<:VectorSpace = findall(.!signbit(V,G)) + +for (op,other) ∈ ((:angular,:radial),(:radial,:angular)) + @eval $op(t::T) where T<:TensorTerm{V,G} where {V,G} = basisindex(ndims(V),bits(basis(t))) ∈ $op(V,G) ? t : zero(V) + for Blade ∈ MSB + @eval function $op(t::$Blade{T,V,G}) where {T,V,G} + out = copy(value(t,mvec(ndims(V),G,T))) + for k ∈ $other(V,G) + @inbounds out[k]≠0 && (out[k] = zero(T)) + end + MBlade{T,V,G}(out) + end + end + @eval function $op(t::MultiVector{T,V}) where {T,V} + out = copy(value(t,mvec(ndims(V),T))) + for k ∈ $other(V) + @inbounds out[k]≠0 && (out[k] = zero(T)) + end + MultiVector{T,V}(out) + end +end + +odd(t::T) where T<:TensorTerm{V,G} where {V,G} = parityinvolute(G) ? t : zero(V) +even(t::T) where T<:TensorTerm{V,G} where {V,G} = parityinvolute(G) ? zero(V) : t +for Blade ∈ MSB + @eval begin + odd(t::$Blade{V,G}) where {V,G} = parityinvolute(G) ? t : zero(V) + even(t::$Blade{V,G}) where {V,G} = parityinvolute(G) ? zero(V) : t + end +end +function odd(t::MultiVector{T,V}) where {T,V} + N = ndims(V) + out = copy(value(t,mvec(N,T))) + bs = binomsum_set(N) + out[1]≠0 && (out[1] = zero(T)) + for g ∈ 3:2:N+1 + for k ∈ bs[g]+1:bs[g+1] + @inbounds out[k]≠0 && (out[k] = zero(T)) + end + end + MultiVector{T,V}(out) +end +function even(t::MultiVector{T,V}) where {T,V} + N = ndims(V) + out = copy(value(t,mvec(N,T))) + bs = binomsum_set(N) + for g ∈ 2:2:N+1 + for k ∈ bs[g]+1:bs[g+1] + @inbounds out[k]≠0 && (out[k] = zero(T)) + end + end + MultiVector{T,V}(out) +end + +imag(t::T) where T<:TensorTerm{V,G} where {V,G} = parityreverse(G) ? t : zero(V) +real(t::T) where T<:TensorTerm{V,G} where {V,G} = parityreverse(G) ? zero(V) : t +for Blade ∈ MSB + @eval begin + imag(t::$Blade{V,G}) where {V,G} = parityreverse(G) ? t : zero(V) + real(t::$Blade{V,G}) where {V,G} = parityreverse(G) ? zero(V) : t + end +end +function imag(t::MultiVector{T,V}) where {T,V} + N = ndims(V) + out = copy(value(t,mvec(N,T))) + bs = binomsum_set(N) + out[1]≠0 && (out[1] = zero(T)) + for g ∈ 3:N+1 + !parityreverse(g-1) && for k ∈ bs[g]+1:bs[g+1] + out[k]≠0 && (out[k] = zero(T)) + end + end + MultiVector{T,V}(out) +end +function real(t::MultiVector{T,V}) where {T,V} + N = ndims(V) + out = copy(value(t,mvec(N,T))) + bs = binomsum_set(N) + for g ∈ 3:N+1 + parityreverse(g-1) && for k ∈ bs[g]+1:bs[g+1] + out[k]≠0 && (out[k] = zero(T)) + end + end + MultiVector{T,V}(out) +end diff --git a/src/utilities.jl b/src/utilities.jl index e2eb39e..fd910eb 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -3,7 +3,7 @@ # Grassmann Copyright (C) 2019 Michael Reed import Base: @pure, print, show, getindex, setindex!, promote_rule, ==, convert, ndims -import DirectSum: Bits, bit2int, doc2m, indexbits, indices, diffmode +import DirectSum: Bits, bit2int, doc2m, indexbits, indices, diffmode, Dim, Grade bcast(op,arg) = op ∈ (:(Reduce.Algebra.:+),:(Reduce.Algebra.:-)) ? Expr(:.,op,arg) : Expr(:call,op,arg.args...) @@ -201,6 +201,8 @@ indexbasis(Int((sparse_limit+cache_limit)/2),1) @pure indexbasis_set(N) = SVector(((N≠0 && N