Ga3.jl

Ga3.jl

# Geometric Algebra G(3) Library in Julia

# Types
struct G3Empty end

struct G3Vector{T<:Number}
    e1::T
    e2::T
    e3::T
end

struct G3Bivector{T<:Number}
    e23::T
    e31::T
    e12::T
end

struct G3Spinor{T<:Number}
    s::T
    B::G3Bivector{T}
end

struct G3Pseudoscalar{T<:Number}
    e123::T
end

struct G3Oddspinor{T<:Number}
    v::G3Vector{T}
    P::G3Pseudoscalar{T}
end

struct G3Multivector{T<:Number}
    s::T
    v::G3Vector{T}
    B::G3Bivector{T}
    P::G3Pseudoscalar{T}
end

# Constructors
G3Vector{T}() where {T<:Number} = G3Vector{T}(zero(T), zero(T), zero(T))
G3Bivector{T}() where {T<:Number} = G3Bivector{T}(zero(T), zero(T), zero(T))
G3Pseudoscalar{T}() where {T<:Number} = G3Pseudoscalar{T}(zero(T))
G3Spinor{T}() where {T<:Number} = G3Spinor{T}(zero(T), G3Bivector{T}())
G3Oddspinor{T}() where {T<:Number} = G3Oddspinor{T}(G3Vector{T}(), G3Pseudoscalar{T}())
G3Multivector{T}() where {T<:Number} = G3Multivector{T}(zero(T), G3Vector{T}(), G3Bivector{T}(), G3Pseudoscalar{T}())

G3Vector{T}(x::Number, y::Number, z::Number) where {T<:Number} = G3Vector{T}(T(x), T(y), T(z))
G3Bivector{T}(x::Number, y::Number, z::Number) where {T<:Number} = G3Bivector{T}(T(x), T(y), T(z))
G3Spinor{T}(s::Number, x::Number, y::Number, z::Number) where {T<:Number} = G3Spinor{T}(T(s), G3Bivector{T}(x, y, z))

# Operations
import Base: +, *

+(a::T, b::T) where {T<:Number} = a + b
+(a::G3Vector{T}, b::G3Vector{T}) where {T<:Number} = G3Vector{T}(a.e1 + b.e1, a.e2 + b.e2, a.e3 + b.e3)
+(a::G3Bivector{T}, b::G3Bivector{T}) where {T<:Number} = G3Bivector{T}(a.e23 + b.e23, a.e31 + b.e31, a.e12 + b.e12)
+(a::G3Pseudoscalar{T}, b::G3Pseudoscalar{T}) where {T<:Number} = G3Pseudoscalar{T}(a.e123 + b.e123)
+(a::G3Spinor{T}, b::G3Spinor{T}) where {T<:Number} = G3Spinor{T}(a.s + b.s, a.B + b.B)
+(a::G3Oddspinor{T}, b::G3Oddspinor{T}) where {T<:Number} = G3Oddspinor{T}(a.v + b.v, a.P + b.P)
+(a::G3Multivector{T}, b::G3Multivector{T}) where {T<:Number} = G3Multivector{T}(a.s + b.s, a.v + b.v, a.B + b.B, a.P + b.P)

*(a::T, b::T) where {T<:Number} = a * b
*(s::T, v::G3Vector{T}) where {T<:Number} = G3Vector{T}(s * v.e1, s * v.e2, s * v.e3)
*(v::G3Vector{T}, s::T) where {T<:Number} = G3Vector{T}(s * v.e1, s * v.e2, s * v.e3)
*(s::T, B::G3Bivector{T}) where {T<:Number} = G3Bivector{T}(s * B.e23, s * B.e31, s * B.e12)
*(B::G3Bivector{T}, s::T) where {T<:Number} = G3Bivector{T}(s * B.e23, s * B.e31, s * B.e12)
*(s::T, q::G3Spinor{T}) where {T<:Number} = G3Spinor{T}(s * q.s, s * q.B)
*(q::G3Spinor{T}, s::T) where {T<:Number} = G3Spinor{T}(s * q.s, s * q.B)
*(s::T, P::G3Pseudoscalar{T}) where {T<:Number} = G3Pseudoscalar{T}(s * P.e123)
*(P::G3Pseudoscalar{T}, s::T) where {T<:Number} = G3Pseudoscalar{T}(s * P.e123)
*(s::T, o::G3Oddspinor{T}) where {T<:Number} = G3Oddspinor{T}(s * o.v, s * o.P)
*(o::G3Oddspinor{T}, s::T) where {T<:Number} = G3Oddspinor{T}(s * o.v, s * o.P)
*(a::G3Vector{T}, b::G3Vector{T}) where {T<:Number} = G3Spinor{T}(
    a.e1*b.e1 + a.e2*b.e2 + a.e3*b.e3,
    G3Bivector{T}(a.e2*b.e3 - b.e2*a.e3, a.e3*b.e1 - b.e3*a.e1, a.e1*b.e2 - b.e1*a.e2)
)
*(a::G3Bivector{T}, b::G3Bivector{T}) where {T<:Number} = G3Spinor{T}(
    -a.e23*b.e23 - a.e31*b.e31 - a.e12*b.e12,
    G3Bivector{T}(a.e12*b.e31 - a.e31*b.e12, a.e23*b.e12 - a.e12*b.e23, a.e31*b.e23 - a.e23*b.e31)
)
*(v::G3Vector{T}, B::G3Bivector{T}) where {T<:Number} = G3Oddspinor{T}(
    G3Vector{T}(v.e3*B.e31 - v.e2*B.e12, v.e1*B.e12 - v.e3*B.e23, v.e2*B.e23 - v.e1*B.e31),
    G3Pseudoscalar{T}(v.e1*B.e23 + v.e2*B.e31 + v.e3*B.e12)
)
*(B::G3Bivector{T}, v::G3Vector{T}) where {T<:Number} = G3Oddspinor{T}(
    G3Vector{T}(-v.e3*B.e31 + v.e2*B.e12, -v.e1*B.e12 + v.e3*B.e23, -v.e2*B.e23 + v.e1*B.e31),
    G3Pseudoscalar{T}(v.e1*B.e23 + v.e2*B.e31 + v.e3*B.e12)
)
*(v::G3Vector{T}, P::G3Pseudoscalar{T}) where {T<:Number} = G3Bivector{T}(v.e1*P.e123, v.e2*P.e123, v.e3*P.e123)
*(P::G3Pseudoscalar{T}, v::G3Vector{T}) where {T<:Number} = G3Bivector{T}(v.e1*P.e123, v.e2*P.e123, v.e3*P.e123)
*(B::G3Bivector{T}, P::G3Pseudoscalar{T}) where {T<:Number} = G3Vector{T}(-B.e23*P.e123, -B.e31*P.e123, -B.e12*P.e123)
*(P::G3Pseudoscalar{T}, B::G3Bivector{T}) where {T<:Number} = G3Vector{T}(-B.e23*P.e123, -B.e31*P.e123, -B.e12*P.e123)
*(a::G3Spinor{T}, b::G3Spinor{T}) where {T<:Number} = let AB = a.B * b.B
    G3Spinor{T}(a.s*b.s + AB.s, AB.B + (a.s * b.B) + (b.s * a.B))
end
*(a::G3Pseudoscalar{T}, b::G3Pseudoscalar{T}) where {T<:Number} = -a.e123 * b.e123

# Wedge Product (∧)
∧(s::T, t::T) where {T<:Number} = s * t
∧(s::T, v::G3Vector{T}) where {T<:Number} = s * v
∧(v::G3Vector{T}, s::T) where {T<:Number} = s * v
∧(a::G3Vector{T}, b::G3Vector{T}) where {T<:Number} = G3Bivector{T}(a.e2*b.e3 - b.e2*a.e3, a.e3*b.e1 - b.e3*a.e1, a.e1*b.e2 - b.e1*a.e2)
∧(v::G3Vector{T}, B::G3Bivector{T}) where {T<:Number} = G3Pseudoscalar{T}(v.e1*B.e23 + v.e2*B.e31 + v.e3*B.e12)
∧(B::G3Bivector{T}, v::G3Vector{T}) where {T<:Number} = G3Pseudoscalar{T}(v.e1*B.e23 + v.e2*B.e31 + v.e3*B.e12)

# Dual and Undual Operations
dual(v::G3Vector{T}) where {T<:Number} = G3Bivector{T}(v.e1, v.e2, v.e3)
dual(B::G3Bivector{T}) where {T<:Number} = G3Vector{T}(-B.e23, -B.e31, -B.e12)
dual(P::G3Pseudoscalar{T}) where {T<:Number} = P.e123
dual(s::T) where {T<:Number} = G3Pseudoscalar{T}(s)

undual(B::G3Bivector{T}) where {T<:Number} = G3Vector{T}(B.e23, B.e31, B.e12)
undual(v::G3Vector{T}) where {T<:Number} = G3Bivector{T}(-v.e1, -v.e2, -v.e3)
undual(s::T) where {T<:Number} = G3Pseudoscalar{T}(s)
undual(P::G3Pseudoscalar{T}) where {T<:Number} = P.e123

# Regressive Product (∨)
∨(a, b) = dual(undual(a) ∧ undual(b))

# Utilities
function Base.show(io::IO, v::G3Vector{T}) where {T<:Number}
    print(io, "$(v.e1)e₁ + $(v.e2)e₂ + $(v.e3)e₃")
end

function Base.show(io::IO, B::G3Bivector{T}) where {T<:Number}
    print(io, "$(B.e23)e₂₃ + $(B.e31)e₃₁ + $(B.e12)e₁₂")
end

function Base.show(io::IO, q::G3Spinor{T}) where {T<:Number}
    print(io, "$(q.s) + $(q.B)")
end

function Base.show(io::IO, P::G3Pseudoscalar{T}) where {T<:Number}
    print(io, "$(P.e123)e₁₂₃")
end

function Base.show(io::IO, o::G3Oddspinor{T}) where {T<:Number}
    print(io, "$(o.v) + $(o.P)")
end

function Base.show(io::IO, M::G3Multivector{T}) where {T<:Number}
    print(io, "$(M.s) + $(M.v) + $(M.B) + $(M.P)")
end

function Base.show(io::IO, ::G3Empty)
    print(io, "∅")
end

# Norm Squared Functions
norm_sqr(s::T) where {T<:Number} = s * s
norm_sqr(v::G3Vector{T}) where {T<:Number} = v.e1*v.e1 + v.e2*v.e2 + v.e3*v.e3
norm_sqr(B::G3Bivector{T}) where {T<:Number} = B.e23*B.e23 + B.e31*B.e31 + B.e12*B.e12
norm_sqr(P::G3Pseudoscalar{T}) where {T<:Number} = P.e123 * P.e123
norm_sqr(q::G3Spinor{T}) where {T<:Number} = q.s*q.s + norm_sqr(q.B)
norm_sqr(o::G3Oddspinor{T}) where {T<:Number} = norm_sqr(o.v) + norm_sqr(o.P)
norm_sqr(M::G3Multivector{T}) where {T<:Number} = M.s*M.s + norm_sqr(M.v) + norm_sqr(M.B) + norm_sqr(M.P)
norm_sqr_study(S::G3Multivector{T}) where {T<:Number} = S.s*S.s + S.P.e123*S.P.e123

# Norm Functions
norm_study(S::G3Multivector{T}) where {T<:Number} = sqrt(norm_sqr_study(S))
norm(x) = sqrt(norm_sqr(x))

# Normalization
normalize(x) = x * (1.0 / norm(x))

# Inverse Functions
inv(s::T) where {T<:Number} = 1.0 / s
inv(v::G3Vector{T}) where {T<:Number} = v * (1.0 / norm_sqr(v))
inv(B::G3Bivector{T}) where {T<:Number} = B * (-1.0 / norm_sqr(B))
inv(P::G3Pseudoscalar{T}) where {T<:Number} = G3Pseudoscalar{T}(-1.0 / P.e123)
inv(q::G3Spinor{T}) where {T<:Number} = conj(q) * (1.0 / norm_sqr(q))
inv(o::G3Oddspinor{T}) where {T<:Number} = conj(o) * (1.0 / norm_sqr(o))
study_conj(S::G3Multivector{T}) where {T<:Number} = G3Multivector{T}(S.s, zero(G3Vector{T}), zero(G3Bivector{T}), -S.P)
inv_study(S::G3Multivector{T}) where {T<:Number} = study_conj(S) * (1.0 / (S * study_conj(S)).s)
inv(M::G3Multivector{T}) where {T<:Number} = error("Inverse for general multivectors is not yet implemented.")
inv(b::G3Biparavector{T}) where {T<:Number} = inv(to_multivector(b))

# Division
divide(x, y) = x * inv(y)

# Metric-Related Functions
## Right Complement
rc(s::T) where {T<:Number} = G3Pseudoscalar{T}(s)
rc(v::G3Vector{T}) where {T<:Number} = G3Bivector{T}(v.e1, v.e2, v.e3)
rc(B::G3Bivector{T}) where {T<:Number} = G3Vector{T}(B.e23, B.e31, B.e12)
rc(P::G3Pseudoscalar{T}) where {T<:Number} = P.e123
rc(q::G3Spinor{T}) where {T<:Number} = G3Oddspinor{T}(rc(q.B), rc(q.s))
rc(o::G3Oddspinor{T}) where {T<:Number} = G3Spinor{T}(rc(o.P), rc(o.v))
rc(M::G3Multivector{T}) where {T<:Number} = G3Multivector{T}(rc(M.P), rc(M.B), rc(M.v), rc(M.s))

## Aliases


const lc = rc
const star = rc

## Regressive Product
vee(a, b) = lc(wedge(rc(a), rc(b)))

## Dot Products
dot(x, y) = vee(x, star(y))
rdot(x, y) = vee(x, star(y))
ldot(x, y) = vee(star(x), y)

## Left Contraction (Simplified)
lcontract(a, b) = error("Full left contraction is not yet implemented.")

# Transcendental Functions
## Exponential Functions
exp(v::G3Vector{T}) where {T<:Number} = G3Paravector{T}(cosh(norm(v)), normalize(v) * sinh(norm(v)))
exp(B::G3Bivector{T}) where {T<:Number} = G3Spinor{T}(cos(norm(B)), normalize(B) * sin(norm(B)))
exp(q::G3Spinor{T}) where {T<:Number} = exp(q.s) * exp(q.B)
exp(p::G3Paravector{T}) where {T<:Number} = G3Paravector{T}(exp(p.s) * exp(p.v).s, exp(p.s) * exp(p.v).v)
exp(P::G3Pseudoscalar{T}) where {T<:Number} = G3Multivector{T}(cos(P.e123), zero(G3Vector{T}), zero(G3Bivector{T}), G3Pseudoscalar{T}(sin(P.e123)))
exp(o::G3Oddspinor{T}) where {T<:Number} = to_multivector(exp(o.v)) * exp(o.P)
exp(M::G3Multivector{T}) where {T<:Number} = error("Exponential for general multivectors is not yet implemented.")

## Square Root Functions
sqrt(S::G3Multivector{T}) where {T<:Number} = begin
    c = sqrt(0.5 * (S.s + norm_study(S)))
    G3Multivector{T}(c, zero(G3Vector{T}), zero(G3Bivector{T}), (0.5 / c) * S.P)
end
sqrt(q::G3Spinor{T}) where {T<:Number} = begin
    n = norm_sqr(q)
    if abs(n - 1.0) > 1e-4
        n = sqrt(n)
        q = q * (1.0 / n)
        s_sign = q.s < 0 ? -1.0 : 1.0
        q = G3Spinor{T}(q.s + s_sign, q.B)
        sqrt(n) * normalize(q)
    else
        s_sign = q.s < 0 ? -1.0 : 1.0
        q = G3Spinor{T}(q.s + s_sign, q.B)
        normalize(q)
    end
end
sqrt(B::G3Bivector{T}) where {T<:Number} = sqrt(G3Spinor{T}(0.0, B))
sqrt(P::G3Pseudoscalar{T}) where {T<:Number} = sqrt(G3Multivector{T}(0.0, zero(G3Vector{T}), zero(G3Bivector{T}), P))

## Logarithm Functions
log(q::G3Spinor{T}) where {T<:Number} = begin
    qn = norm_sqr(q)
    if abs(qn - 1.0) > 1e-4
        qn = sqrt(qn)
        q = q * (1.0 / qn)
        log(qn) + normalize(q.B) * acos(q.s)
    else
        zero(T) + normalize(q.B) * acos(q.s)
    end
end
log(B::G3Bivector{T}) where {T<:Number} = log(G3Spinor{T}(0.0, B))