generalized exp and inv methods fixed #15 #13

src/algebra.jl

@@ -985,57 +985,64 @@ for Blade ∈ MSB
         div(a::$Blade{T,V,G},m) where {T,V,G} = $Blade{T,V,G}(div.(value(a),m))
     end
 end
-@pure /(a,b::T) where T<:TensorTerm = a*inv(b)
-for Blade ∈ MSB
-    @eval /(a,b::T) where T<:$Blade = a*inv(b)
+function inv(m::MultiVector{T,V}) where {T,V}
+    rm = ~m
+    d = m*rm
+    for k ∈ 0:ndims(V)
+        dk = d(k)
+        frobenius(dk) ≈ frobenius(d) && (return rm/dk)
+    end
+    throw(error("inv($m) is undefined"))
 end
+rem(a::MultiVector{T,V},m) where {T,V} = MultiVector{T,V}(rem.(value(a),m))
+div(a::MultiVector{T,V},m) where {T,V} = MultiVector{T,V}(div.(value(a),m))
+@pure /(a,b::T) where T<:TensorAlgebra = a*inv(b)
 for Term ∈ (:TensorTerm,MSB...,:MultiVector,:MultiGrade)
     @eval begin
         @pure /(a::S,b::T) where {S<:$Term,T<:Number} = a*inv(b)
         @pure /(a::S,b::UniformScaling) where S<:$Term = a*inv(vectorspace(a)(b))
     end
 end
-rem(a::MultiVector{T,V},m) where {T,V} = MultiVector{T,V}(rem.(value(a),m))
-div(a::MultiVector{T,V},m) where {T,V} = MultiVector{T,V}(div.(value(a),m))
 
 ## exponential & logarithm function
 
 import Base: exp, log
 
 function exp(t::T) where T<:TensorAlgebra{V} where V
     terms = TensorAlgebra{V}[t,(t^2)/2]
-    norms = abs.(terms)
+    norms = frobenius.(terms)
     k = 3
     while norms[end]<norms[end-1] || norms[end]>1
         push!(terms,(terms[end]*t)/k)
-        push!(norms,abs(terms[end]))
+        push!(norms,frobenius(terms[end]))
         k += 1
     end
     return 1+sum(terms[1:end-1])
 end
 
 function log(t::T) where T<:TensorAlgebra{V} where V
-    frob = abs(t)
+    frob = frobenius(t)
     k = 3
+    τ = t-1
     if frob ≤ 5/4
-        prods = TensorAlgebra{V}[t-1,(t-1)^2]
-        terms = TensorAlgebra{V}[t,prods[2]/2]
-        norms = [frob,abs(terms[2])]
+        prods = TensorAlgebra{V}[τ,τ^2]
+        terms = TensorAlgebra{V}[τ,prods[2]/2]
+        norms = [frob,frobenius(terms[2])]
         while (norms[end]<norms[end-1] || norms[end]>1) && k ≤ 300
-            push!(prods,prods[end]*(t-1))
+            push!(prods,prods[end]*τ)
             push!(terms,prods[end]/(k*(-1)^(k+1)))
-            push!(norms,abs(terms[end]))
+            push!(norms,frobenius(terms[end]))
             k += 1
         end
     else
-        s = inv(t*inv(t-1))
+        s = inv(t*inv(τ))
         prods = TensorAlgebra{V}[s,s^2]
         terms = TensorAlgebra{V}[s,2prods[2]]
-        norms = abs.(terms)
+        norms = frobenius.(terms)
         while (norms[end]<norms[end-1] || norms[end]>1) && k ≤ 300
             push!(prods,prods[end]*s)
             push!(terms,k*prods[end])
-            push!(norms,abs(terms[end]))
+            push!(norms,frobenius(terms[end]))
             k += 1
         end
     end

src/multivectors.jl

@@ -313,7 +313,7 @@ end
 ## Generic
 
 import Base: isinf, isapprox
-export basis, grade, hasinf, hasorigin, isorigin, scalar
+export basis, grade, hasinf, hasorigin, isorigin, scalar, frobenius, norm, norm2
 
 const VBV = Union{MValue,SValue,MBlade,SBlade,MultiVector}
 
@@ -337,9 +337,11 @@ const VBV = Union{MValue,SValue,MBlade,SBlade,MultiVector}
 @pure hasorigin(t::Union{MValue,SValue}) = hasorigin(basis(t))
 @pure hasorigin(m::TensorAlgebra) = hasorigin(vectorspace(m))
 
+@inline frobenius(t::T) where T<:TensorAlgebra = norm(value(t))
+
 function isapprox(a::S,b::T) where {S<:TensorMixed{T1},T<:TensorMixed{T2}} where {T1, T2}
     rtol = Base.rtoldefault(T1, T2, 0)    
-    return norm(value(a-b)) <= rtol * max(norm(value(a)), norm(value(b)))
+    return frobenius(a-b) <= rtol * max(frobenius(a), frobenius(b))
 end
 isapprox(a::S,b::T) where {S<:TensorMixed,T<:TensorTerm} = isapprox(a, MultiVector(b))
 isapprox(b::S,a::T) where {S<:TensorTerm,T<:TensorMixed} = isapprox(a, MultiVector(b))
@@ -378,6 +380,8 @@ for Blade ∈ MSB
     end
 end
 
+isscalar(t::MultiVector) = frobenius(t) == scalar(t)
+
 ## MultiGrade{N}
 
 struct MultiGrade{V} <: TensorAlgebra{V}

test/issuestests.jl

@@ -54,10 +54,10 @@ end
     i, j, k = hyperplanes(ℝ^3)
     alpha = 0.5π
 
-    @test exp(alpha/2*(i)) == 0.7071067811865476 - 0.7071067811865475v23
+    @test exp(alpha/2*(i)) ≈ sqrt(2)*(1+i)/2
 
     a, b, c = 1/sqrt(2) * [1, 1, 0]
-    @test_broken exp(alpha/2*(a*i + b*j + c*k))
+    @test exp(alpha/2*(a*i + b*j + c*k)) ≈ (sqrt(2)+j+i)/2
 end
 
 @testset "Issue #20: geometric product of null basis and negative origin" begin