Skip to content

Commit

Permalink
imported Leibniz and improved expm1, cosh, sinh, isapprox
Browse files Browse the repository at this point in the history
  • Loading branch information
@chakravala
chakravala committed Jun 22, 2019
1 parent f76dd85 commit 7295818
Show file tree
Hide file tree
Showing 5 changed files with 75 additions and 48 deletions.
1 change: 1 addition & 0 deletions Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ AbstractTensors = "a8e43f4a-99b7-5565-8bf1-0165161caaea"
Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
ComputedFieldTypes = "459fdd68-db75-56b8-8c15-d717a790f88e"
DirectSum = "22fd7b30-a8c0-5bf2-aabe-97783860d07c"
Leibniz = "edad4870-8a01-11e9-2d75-8f02e448fc59"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Reduce = "93e0c654-6965-5f22-aba9-9c1ae6b3c259"
Requires = "ae029012-a4dd-5104-9daa-d747884805df"
Expand Down
14 changes: 13 additions & 1 deletion src/Grassmann.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -323,8 +323,20 @@ end

# ParaAlgebra

using Leibniz
import Leibniz: ∂, d

generate_product_algebra(:(Leibniz.Operator),:svec)

function (V::Signature{N})(d::Leibniz.Derivation{T,O}) where {N,T,O}
x=sum([(∂(V,k)^2) for k ∈ 1:N])^div(isodd(O) ? O-1 : O,2)
isodd(O) ? sum([(x*∂(V,k))*getbasis(V,1<<(k-1)) for k ∈ 1:N]) : x*getbasis(V,0)
end

∂(ω::T) where T<:TensorAlgebra{V} where V = ω⋅V(∇)
d(ω::T) where T<:TensorAlgebra{V} where V = V(∇)∧ω

function __init__()
@require Leibniz="edad4870-8a01-11e9-2d75-8f02e448fc59" generate_product_algebra(:(Leibniz.Operator),:svec)
@require Reduce="93e0c654-6965-5f22-aba9-9c1ae6b3c259" begin
import Reduce: RExpr
*(a::RExpr,b::Basis{V}) where V = SValue{V}(a,b)
Expand Down
34 changes: 27 additions & 7 deletions src/composite.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,9 @@ export exph, log_fast, logh_fast

## exponential & logarithm function

@inline Base.expm1(t::Basis{V,0}) where V = SValue{V}(ℯ-1)
@inline Base.expm1(t::T) where T<:TensorTerm{V,0} where V = SValue{V}(expm1(value(t)))

function Base.expm1(t::T) where T<:TensorAlgebra{V} where V
S,term,f = t,(t^2)/2,frobenius(t)
norms = MVector(f,frobenius(term),f)
Expand All @@ -22,13 +25,15 @@ function Base.expm1(t::T) where T<:TensorAlgebra{V} where V
end

@eval function Base.expm1(b::MultiVector{T,V}) where {T,V}
$(insert_expr((:N,:t,:out,:bs,:bn),:mvec,:T,Float64)...)
B = value(b)
sb,nb = scalar(b),norm(B)
sb ≈ nb && (return SValue{V}(expm1(sb)))
$(insert_expr((:N,:t,:out,:bs,:bn),:mvec,:T,Float64)...)
S = zeros(mvec(N,t))
term = zeros(mvec(N,t))
S .= value(b)
S .= B
out .= value(b^2)/2
norms = MVector(norm(S),norm(out),norm(term))
norms = MVector(nb,norm(out),norm(term))
k::Int = 3
@inbounds while (norms[2]<norms[1] || norms[2]>1) && k ≤ 10000
S += out
Expand Down Expand Up @@ -130,11 +135,20 @@ for base ∈ (2,10)
@eval Base.$fe(t::T) where T<:TensorAlgebra = exp(log($base)*t)
end

@inline Base.sqrt(t::T) where T<:TensorAlgebra = exp(log(t)/2)
@inline Base.cbrt(t::T) where T<:TensorAlgebra = exp(log(t)/3)
for (qrt,n) ∈ ((:sqrt,2),(:cbrt,3))
@eval begin
@inline Base.$qrt(t::Basis{V,0} where V) = t
@inline Base.$qrt(t::T) where T<:TensorTerm{V,0} where V = SValue{V}($qrt(value(t)))
@inline function Base.$qrt(t::T) where T<:TensorAlgebra
isscalar(t) ? $qrt(scalar(t)) : exp(log(t)/$n)
end
end
end

## trigonometric

@inline Base.cosh(t::T) where T<:TensorTerm{V,0} where V = SValue{V}(cosh(value(t)))

function Base.cosh(t::T) where T<:TensorAlgebra{V} where V
τ = t^2
S,term = τ/2,(τ^2)/24
Expand All @@ -153,6 +167,8 @@ function Base.cosh(t::T) where T<:TensorAlgebra{V} where V
end

@eval function Base.cosh(b::MultiVector{T,V,E}) where {T,V,E}
sb,nb = scalar(b),frobenius(b)
sb ≈ nb && (return SValue{V}(cosh(sb)))
$(insert_expr((:N,:t,:out,:bs,:bn),:mvec,:T,Float64)...)
τ::MultiVector{T,V,E} = b^2
B = value(τ)
Expand Down Expand Up @@ -189,6 +205,8 @@ end
return MultiVector{t,V}(S)
end

@inline Base.sinh(t::T) where T<:TensorTerm{V,0} where V = SValue{V}(sinh(value(t)))

function Base.sinh(t::T) where T<:TensorAlgebra{V} where V
τ,f = t^2,frobenius(t)
S,term = t,(t*τ)/6
Expand All @@ -206,6 +224,8 @@ function Base.sinh(t::T) where T<:TensorAlgebra{V} where V
end

@eval function Base.sinh(b::MultiVector{T,V,E}) where {T,V,E}
sb,nb = scalar(b),frobenius(b)
sb ≈ nb && (return SValue{V}(sinh(sb)))
$(insert_expr((:N,:t,:out,:bs,:bn),:mvec,:T,Float64)...)
τ::MultiVector{T,V,E} = b^2
B = value(τ)
Expand Down Expand Up @@ -268,7 +288,7 @@ Base.cosc(t::T) where T<:TensorAlgebra{V} where V = iszero(t) ? zero(V) : (x=(1*

exph(t) = cosh(t)+sinh(t)

function log_fast(t)
function log_fast(t::T) where T<:TensorAlgebra{V} where V
term = zero(V)
norm = MVector(0.,0.)
while true
Expand All @@ -280,7 +300,7 @@ function log_fast(t)
return term
end

function logh_fast(t)
function logh_fast(t::T) where T<:TensorAlgebra{V} where V
term = zero(V)
norm = MVector(0.,0.)
while true
Expand Down
55 changes: 34 additions & 21 deletions src/multivectors.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ end

@inline indices(b::Basis{V}) where V = indices(bits(b),ndims(V))

@pure Basis{V}() where V = getbasis(V,0)
@pure Basis{V}(i::Bits) where V = getbasis(V,i)
Basis{V}(b::BitArray{1}) where V = getbasis(V,bit2int(b))

Expand Down Expand Up @@ -338,51 +339,63 @@ const VBV = Union{MValue,SValue,MBlade,SBlade,MultiVector}
@pure hasorigin(m::TensorAlgebra) = hasorigin(vectorspace(m))

@inline frobenius(t::T) where T<:TensorAlgebra = norm(value(t))
@inline Base.iszero(t::T) where T<:TensorAlgebra = frobenius(t) == 0
@inline Base.isone(t::T) where T<:TensorAlgebra = frobenius(t) == scalar(t) == 1
@inline Base.iszero(t::T) where T<:TensorAlgebra = frobenius(t) ≈ 0
@inline Base.isone(t::T) where T<:TensorAlgebra = frobenius(t) ≈ scalar(t) ≈ 1

function isapprox(a::S,b::T) where {S<:TensorMixed{T1},T<:TensorMixed{T2}} where {T1, T2}
rtol = Base.rtoldefault(T1, T2, 0)
for A ∈ (:TensorTerm,MSB...), B ∈ (:TensorTerm,MSB...)
@eval isapprox(a::S,b::T) where {S<:$A,T<:$B} = vectorspace(a)==vectorspace(b) && (grade(a)==grade(b) ? frobenius(a)≈frobenius(b) : (iszero(a) && iszero(b)))
end
isapprox(a::S,b::T) where {S<:MultiVector,T<:MultiVector} = vectorspace(a)==vectorspace(b) && value(a) ≈ value(b)
function isapprox(a::S,b::T) where {S<:TensorAlgebra,T<:TensorAlgebra}
rtol = Base.rtoldefault(valuetype(a), valuetype(b), 0)
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))
isapprox(a::S,b::T) where {S<:TensorTerm,T<:TensorTerm} = isapprox(MultiVector(a), MultiVector(b))
function isapprox(a::S,b::T) where {S<:TensorAlgebra,T<:Number}
rtol = Base.rtoldefault(valuetype(a), T, 0)
return frobenius(a-b) <= rtol * max(frobenius(a), b)
end
isapprox(a::S,b::T) where {S<:Number,T<:TensorAlgebra} = isapprox(b,a)

"""
scalar(multivector)
Return the scalar (grade 0) part of any multivector.
"""
scalar(t::T) where T<:TensorTerm{V,0} where V = t
scalar(t::T) where T<:TensorTerm{V} where V = zero(V)
scalar(t::MultiVector) = t(0,1)
@inline scalar(t::T) where T<:TensorTerm{V,0} where V = t
@inline scalar(t::T) where T<:TensorTerm{V} where V = zero(V)
@inline scalar(t::MultiVector) = t(0,1)
for Blade ∈ MSB
@eval begin
scalar(t::$Blade{T,V,0} where {T,V}) = a(1)
scalar(t::$Blade{T,V} where T) where V = zero(V)
@inline scalar(t::$Blade{T,V,0} where {T,V}) = t(1)
@inline scalar(t::$Blade{T,V} where T) where V = zero(V)
end
end

vector(t::T) where T<:TensorTerm{V,1} where V = t
vector(t::T) where T<:TensorTerm{V} where V = zero(V)
vector(t::MultiVector) = t(1)
@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(1)
for Blade ∈ MSB
@eval begin
vector(t::$Blade{T,V,1} where {T,V}) = a
vector(t::$Blade{T,V} where T) where V = zero(V)
@inline vector(t::$Blade{T,V,1} where {T,V}) = t
@inline vector(t::$Blade{T,V} where T) where V = zero(V)
end
end

volume(t::T) where T<:TensorTerm{V,G} where {V,G} = G == ndims(V) ? t : zero(V)
volume(t::MultiVector{T,V} where T) where V = t(ndims(V),1)
@inline volume(t::T) where T<:TensorTerm{V,G} where {V,G} = G == ndims(V) ? t : zero(V)
@inline volume(t::MultiVector{T,V} where T) where V = t(ndims(V),1)
for Blade ∈ MSB
@eval begin
volume(t::$Blade{T,V,G} where T) where {V,G} = G == ndims(V) ? t : zero(V)
@inline volume(t::$Blade{T,V,G} where T) where {V,G} = G == ndims(V) ? t : zero(V)
end
end

isscalar(t::MultiVector) = frobenius(t) == scalar(t)
@inline isscalar(t::T) where T<:TensorTerm = grade(t) == 0 || iszero(t)
@inline isscalar(t::T) where T<:SBlade = grade(t) == 0 || iszero(t)
@inline isscalar(t::T) where T<:MBlade = grade(t) == 0 || iszero(t)
@inline isscalar(t::MultiVector) = frobenius(t) ≈ scalar(t)

@inline Base.real(t::T) where T<:TensorAlgebra = scalar(t)
@inline Base.imag(t::T) where T<:TensorAlgebra = t-real(t)

## MultiGrade{N}

Expand Down
19 changes: 0 additions & 19 deletions src/symbolic.jl
View file

This file was deleted.

0 comments on commit 7295818

Please sign in to comment.