added recursive Shirokov inverse algorithm for general multivectors. …

…renames sortbasis to sort_basis.

README.md

@@ -185,7 +185,7 @@ Transform a 1-vector with the sandwich product.
 
     julia> v = reverse(q)*(1.0e₁+1.0e₂+1.0e₃)*q
 
-    julia> grade(v, 1) |> prune∘sortbasis
+    julia> grade(v, 1) |> prune∘sort_basis
     2-element KVector{Float64,1,2}:
      1.414213562373095e₁
                    1.0e₃

src/KVectors.jl

@@ -29,7 +29,7 @@ KVectors are containers for Blades of the same grade.  New Blade elements are ad
 export 
 KVector,
 NullKVector,
-sortbasis,
+sort_basis,
 isnull,
 coords,
 prune,
@@ -173,13 +173,13 @@ Hodge star operator mapping k to it's Hodge dual.
 ⋆(k::K) where {K<:KVector} = mapreduce(⋆,+,k)
 
 """
-    sortbasis(b)
+    sort_basis(b)
 
 sort blades by bases indices in ascending order within the k-vector
 """
-sortbasis(B::BT) where {BT<:KVector} = BT(sort(Vector(B.k); by=subspace∘typeof))
-sortbasis(B::BT) where {BT<:Blade} = B
-Base.:(==)(B::BT, B2::BT) where {BT<:KVector} = sortbasis(B).k == sortbasis(B2).k
+sort_basis(B::BT) where {BT<:KVector} = BT(sort(Vector(B.k); by=subspace∘typeof))
+sort_basis(B::BT) where {BT<:Blade} = B
+Base.:(==)(B::BT, B2::BT) where {BT<:KVector} = sort_basis(B).k == sort_basis(B2).k
 Base.:(==)(B::BT, B2::BT2) where {BT<:KVector, BT2<:KVector} = false
 #Base.:(==)(B::BT, B2::BT2) where {T,K,K2,BT<:KVector{T,K}, BT2<:KVector{T,K2}} = false
 
@@ -240,8 +240,8 @@ for length(a) = 2, gram( a, u ) = [ a₁⋅u₁ a₁⋅u₂ ;
                                     a₂⋅u₁ a₂⋅u₂ ]
 """
 function gram(a::K,u::L) where {K<:KVector, L<:KVector}
-  a = sortbasis(prune(a))
-  u = sortbasis(prune(u))
+  a = sort_basis(prune(a))
+  u = sort_basis(prune(u))
   n = length(a)
   if n != length(u)
     zero(eltype(eltype(a)))
@@ -300,7 +300,7 @@ Base.in(be::BK, bs::BK2) where {BK<:Union{Blade,KVector}, BK2<:Union{Blade,KVect
 fieldtype(k::K) where {F<:Number, K<:KVector{F}} = F
 
 function Base.isapprox(k::K, l::K2; kwargs...) where {K<:KVector, K2<:KVector} 
-  mapreduce( (b,c)->isapprox(b,c; kwargs...), (acc,e)->acc && e, (sortbasis∘prune)(k), (sortbasis∘prune)(l))
+  mapreduce( (b,c)->isapprox(b,c; kwargs...), (acc,e)->acc && e, (sort_basis∘prune)(k), (sort_basis∘prune)(l))
 end
 
 """

src/Multivectors.jl

@@ -37,7 +37,8 @@ involute,
 cayley_table,
 cayley_matrix_description,
 matrix_representation,
-newton_inv
+newton_inv,
+shirokov_inv
 
 
 using Base.Iterators
@@ -159,6 +160,15 @@ import Base: map, mapreduce
 map(f, m::Multivector) = (f(k) for k in m)
 mapreduce(f, r, m::Multivector) = reduce(r, (f(k) for k in m))
 
+function pseudoscalar(U::M) where {T, M<:Multivector{T}} 
+  for i in 1:length(U.B) 
+    if length(U.B[i]) > 0
+      return pseudoscalar(first(U.B[i]))
+    end
+  end
+  pseudoscalar(one(T)*e₁)
+end
+
 Base.conj(m::Multivector) = mapreduce(conj, +, m)
 
 Base.:+(b::KVector{T,K}, c::Multivector{T,L}) where {T,K,L} = Multivector(b) + c
@@ -301,11 +311,11 @@ Base.:*(s::T, M::MT) where {T<:Real, MT<:Multivector} = M*s
 Base.:/(M::MT, s::T) where {T<:Real, MT<:Multivector} = M*(one(T)/s)
 Base.:/(A::M,B::N) where {M<:CliffordNumber,N<:CliffordNumber} = A*inv(B) 
 
-sortbasis(s::T) where T<:Real = s
+sort_basis(s::T) where T<:Real = s
 
-Base.:(==)(a::M,b::M) where {M<:Multivector} = (a.s==b.s) && sortbasis.(kvectors(prune(a)))==sortbasis.(kvectors(prune(b)))
+Base.:(==)(a::M,b::M) where {M<:Multivector} = (a.s==b.s) && sort_basis.(kvectors(prune(a)))==sort_basis.(kvectors(prune(b)))
 Base.:(==)(a::M,b::N) where {M<:Multivector, N<:Multivector} = (a.s==b.s) && 
-  sortbasis.(kvectors(prune(a)))==sortbasis.(kvectors(prune(b)))
+  sort_basis.(kvectors(prune(a)))==sort_basis.(kvectors(prune(b)))
 
 for M in [Real, Multivector, Blade, KVector]
   for N in [Real, Multivector, Blade, KVector]
@@ -318,7 +328,9 @@ end
 
 Base.reverse(a::M) where M<:Multivector = a.s+mapreduce(reverse, +, a.B; init=M())
 Base.reverse(s::T) where T<:Real = s
+Base.:~(k::K) where {K<:CliffordNumber} = reverse(k)
 
+involute(v::CliffordNumber) = mapreduce(vᵢ->(-1)^grade(vᵢ)*vᵢ, +, v)
 involute(A::M) where {T, M<:Multivector{T}} = M(A.s, map(((k,b),)->(-one(T))^k*b, enumerate(kvectors(A))))
 
 scalarprod(A::M,B::N) where {M<:Multivector,N<:Multivector} = grade(A*B,0)
@@ -628,7 +640,6 @@ end
 
 fieldtype(M::MT) where {F<:Number, MT<:Multivector{F}} = F
 
-Base.:~(M::CliffordNumber) = reverse(M)
 
 """
     newton_inv(m)
@@ -653,6 +664,43 @@ function newton_inv(m,
   xᵢ*n
 end
 
+
+"""
+conjugation from Shirokov inverse paper
+"""
+Δⱼ(v::CliffordNumber, j::Int) = mapreduce(vᵢ->vᵢ*(-1)^binomial(grade(vᵢ), 2^(j-1)), +, v)
+
+# helpers for skirokov_inv
+C(U::CliffordNumber, k, N) = grade(U, 0)*N/k
+
+"""
+  shirokov_inv(U)
+
+  algebraic inverse of multivector of arbitrary grade in a non-degenerate algebra
+
+  Ref: On determinant, other characteristic polynomial coefficients, and inverses in Clifford algebras of arbitrary dimension
+D. S. Shirokov
+
+"""
+function shirokov_inv(U::CliffordNumber)
+  n = grade(pseudoscalar(U))
+  N = 2^((n+1)/2)
+
+  Ukee = Multivector[U]
+  for k in 2:N
+    Uprev = last(Ukee)
+    Unext = U*(Uprev - C(Uprev, k-1, N))
+    push!(Ukee, Unext)
+  end
+
+  U_N_1 = Ukee[end-1]
+  AdjU = C(U_N_1, N-1, N) - U_N_1
+  DetU = -last(Ukee)
+
+  AdjU/DetU
+
+end
+
 function Base.isapprox(M::MT, N::MT2; kwargs...) where {MT<:CliffordNumber, MT2<:CliffordNumber} 
   mapreduce((b,c)->isapprox(b,c; kwargs...), (acc,e)->acc && e, Multivector(M), Multivector(N))
 end

test/KVectors_runtests.jl

@@ -76,12 +76,12 @@ using .KG3
   e₁, e₂, e₃ = alle(KG3, 3)[1:3]
   𝐼 = alle(KG3,3)[end]
 
-  a = sortbasis(1.0e₁ + 3.0e₃)
+  a = sort_basis(1.0e₁ + 3.0e₃)
 
   @test eltype(a) <: Blade
   @test size(a) == (2,)
   @test a == KVector([1.0,0.0,3.0]) |> Multivectors.prune
-  @test sortbasis(1.0e₂) == 1.0e₂
+  @test sort_basis(1.0e₂) == 1.0e₂
   @test convert(KVector, 1.0e₂) == KVector(1.0e₂)
   @test first(a) == a[1]
   @test firstindex(a) == 1
@@ -108,7 +108,7 @@ using .KG3
   @test a/2.0 == a*0.5
   @test !(a == a∧e₂(1.0))
 
-  @test coords(a) == scalar.(sortbasis(a+0.0e₂))
+  @test coords(a) == scalar.(sort_basis(a+0.0e₂))
   @test coords(a[1]) == [scalar(a[1]), 0.0, 0.0]
   @test Multivectors.prune(KVector(coords(a) .* basis_1blades(a))) == a
   @test norm(basis_1vector(a)) == sqrt(3.0)

test/runtests.jl

@@ -149,6 +149,10 @@ using .CGA
   @test i+j == Multivector(i) + Multivector(j) == j+i == Multivector(i) + j == i + Multivector(j)
   @test i+j+k == i+(j+k) == (i+k)+j 
 
+  v = i+j+k+1.0
+  @test Multivectors.Δⱼ(v, 1) == involute(v)
+  @test Multivectors.Δⱼ(v, 2) == ~v
+  @test prune(shirokov_inv(i+j+k)*(i+j+k)) ≈ 1.0
 end
 
 module PG3
@@ -421,7 +425,7 @@ k	j	−i	−1
   # Transform a 1-vector with the sandwich product.
   v = reverse(q)*(1.0e₁+1.0e₂+1.0e₃)*q
 
-  v′ = grade(v, 1) |> prune∘sortbasis
+  v′ = grade(v, 1) |> prune∘sort_basis
   @test v′⋅1.0e₃ == 1.0
   @test v′⋅1.0e₁ ≈ sqrt(2.0)