From 37ffd8aa7da6ad2c95bd053889de69011464005a Mon Sep 17 00:00:00 2001 From: Michael Reed <18372368+chakravala@users.noreply.github.com> Date: Mon, 17 Jun 2019 20:07:29 -0400 Subject: [PATCH] generalized exp and inv methods fixed #15 #13 --- src/algebra.jl | 39 +++++++++++++++++++++++---------------- src/multivectors.jl | 8 ++++++-- test/issuestests.jl | 4 ++-- 3 files changed, 31 insertions(+), 20 deletions(-) diff --git a/src/algebra.jl b/src/algebra.jl index 14eb3f8..2e11ebf 100644 --- a/src/algebra.jl +++ b/src/algebra.jl @@ -985,18 +985,24 @@ 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 @@ -1004,38 +1010,39 @@ 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]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]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]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 diff --git a/src/multivectors.jl b/src/multivectors.jl index e464d13..1920c18 100644 --- a/src/multivectors.jl +++ b/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} diff --git a/test/issuestests.jl b/test/issuestests.jl index fd5acf1..1b37ed4 100644 --- a/test/issuestests.jl +++ b/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