Minor correctness improvements

Project.toml

@@ -1,7 +1,7 @@
 name = "SimpleGA"
 uuid = "6f159fe9-4850-48a0-97c4-7cf487b3e6ea"
 authors = ["Monumo Ltd.", "Tom Gillam <tpgillam@googlemail.com>"]
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/index.md

@@ -212,7 +212,7 @@ M -> M * e4  # Map from odd to even.
     * `STA.I4 = g0 * g1 * g2 * g3`
 * Extra functions: `bar(M) = g0 * M * g0 `
 
-The algebra G(1,3) is the geometric algebra for Minkowski spacetime, better known as the Spacetime Algebra (STA). STA is the standard abbreviationn for this algebra so is used for the namespace.
+The algebra G(1,3) is the geometric algebra for Minkowski spacetime, better known as the Spacetime Algebra (STA). STA is the standard abbreviation for this algebra so is used for the namespace.
 
 Even elements in the STA are isomorphic to 2X2 complex matrices (the full Pauli algebra), which is essentially how they are implemented here. (We acually unwrap the matrices and store the entries explicitly in a struct, which turns out to be slightly more efficient.)
 
@@ -228,12 +228,12 @@ Note that the STA is also perfect for the conformal model of 2D Euclidean space,
 ## G(3,1) and 2D conformal geometry
 * Name: `GA31`
 * Basis: `GA31.basis`
-* Basis vectors: `[GA31.e1, GA31.e2, GA31.e3. GA31.f3]`
+* Basis vectors: `[GA31.e1, GA31.e2, GA31.e3, GA31.f3]`
 * Other basis elements: `GA31.I4 = e1 * e2 * e3 * f3`
 
 The algebra G(3,1) is the close relative of the STA, just with signatures flipped. Under the hood the implementation is almost identical, based on 2X2 complex matrices. Only a handful of signs are changed.
 
-The main reason for including G(3,1) as a separate algebra is its use as the conformal algebra for the 2D plane. It is more convenient for 2D work to have an algebra where the 2D basis vectors have positive norm. As this is the main use case, we have not copied the STA convention of displaying the basis vectors with Greek letters. Instead we have adopted the standard convention of ```e1, e2, e3```s have positive square and ```f3``` having negative square.
+The main reason for including G(3,1) as a separate algebra is its use as the conformal algebra for the 2D plane. It is more convenient for 2D work to have an algebra where the 2D basis vectors have positive norm. As this is the main use case, we have not copied the STA convention of displaying the basis vectors with Greek letters. Instead we have adopted the standard convention of `e1, e2, e3`s have positive square and `f3` having negative square.
 
 
 ## G(3,0,1), the PGA

src/cga.jl

@@ -35,125 +35,120 @@ const f4 = Odd{Int8}(
 const I5 = e1 * e2 * e3 * e4 * f4
 const basis = SA[e1, e2, e3, e4, f4]
 
-#Sets tolerance for not displaying results. Adding 1 to comparison seems to work well.
-approxzero(x::Real) = isapprox(1 + x, 1.0)
-
-function mv_to_text(a::Even)
+function Base.show(io::IO, a::Even)
     res = ""
     scl = tr(a)
-    tp = approxzero(scl) ? "" : " + " * string(scl)
+    tp = iszero(scl) ? "" : " + " * string(scl)
     res *= tp
     scl = dot(a, -e1 * e2)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1e2"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1e2"
     res *= tp
     scl = dot(a, -e2 * e3)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e2e3"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e2e3"
     res *= tp
     scl = dot(a, -e3 * e1)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e3e1"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e3e1"
     res *= tp
     scl = dot(a, -e1 * e4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1e4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1e4"
     res *= tp
     scl = dot(a, -e2 * e4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e2e4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e2e4"
     res *= tp
     scl = dot(a, -e3 * e4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e3e4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e3e4"
     res *= tp
     scl = dot(a, e1 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1f4"
     res *= tp
     scl = dot(a, e2 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e2f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e2f4"
     res *= tp
     scl = dot(a, e3 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e3f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e3f4"
     res *= tp
     scl = dot(a, e4 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e4f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e4f4"
     res *= tp
     scl = dot(a, -e1 * I5)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "I5e1"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "I5e1"
     res *= tp
     scl = dot(a, -e2 * I5)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "I5e2"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "I5e2"
     res *= tp
     scl = dot(a, -e3 * I5)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "I5e3"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "I5e3"
     res *= tp
     scl = dot(a, -e4 * I5)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "I5e4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "I5e4"
     res *= tp
     scl = dot(a, f4 * I5)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "I5f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "I5f4"
     res *= tp
     if (length(res) == 0)
         res = "0.0"
     else
         res = chop(res; head=3, tail=0)
     end
-    return res
+    return print(io, res)
 end
 
-function mv_to_text(a::Odd)
+function Base.show(io::IO, a::Odd)
     res = ""
     scl = dot(a, e1)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1"
     res *= tp
     scl = dot(a, e2)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e2"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e2"
     res *= tp
     scl = dot(a, e3)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e3"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e3"
     res *= tp
     scl = dot(a, e4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e4"
     res *= tp
     scl = dot(a, -f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "f4"
     res *= tp
     scl = dot(a, -e1 * e2 * e3)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1e2e3"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1e2e3"
     res *= tp
     scl = dot(a, -e1 * e2 * e4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1e2e4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1e2e4"
     res *= tp
     scl = dot(a, -e1 * e3 * e4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1e3e4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1e3e4"
     res *= tp
     scl = dot(a, -e2 * e3 * e4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e2e3e4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e2e3e4"
     res *= tp
     scl = dot(a, e1 * e2 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1e2f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1e2f4"
     res *= tp
     scl = dot(a, e1 * e3 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1e3f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1e3f4"
     res *= tp
     scl = dot(a, e2 * e3 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e2e3f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e2e3f4"
     res *= tp
     scl = dot(a, e1 * e4 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e1e4f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e1e4f4"
     res *= tp
     scl = dot(a, e2 * e4 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e2e4f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e2e4f4"
     res *= tp
     scl = dot(a, e3 * e4 * f4)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "e3e4f4"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "e3e4f4"
     res *= tp
     scl = dot(a, -I5)
-    tp = approxzero(scl) ? "" : " + " * string(scl) * "I5"
+    tp = iszero(scl) ? "" : " + " * string(scl) * "I5"
     res *= tp
     if (length(res) == 0)
         res = "0.0"
     else
         res = chop(res; head=3, tail=0)
     end
-    return res
+    return print(io, res)
 end
 
-include("show.jl")
-
 end #Module

src/core20.jl

@@ -38,12 +38,12 @@ Base.:(-)(a::Odd, b::Odd) = Odd(a.c1 - b.c1)
 
 #Scalar addition / subtraction. Other cases are in GAcommon
 #Relies on Julia's promotion rules to do the sensible thing.
-Base.:(+)(num::Number, a::Even) = Even(a.c1 + num)
-Base.:(-)(num::Number, a::Even) = Even(-a.c1 + num)
+Base.:(+)(num::Real, a::Even) = Even(a.c1 + num)
+Base.:(-)(num::Real, a::Even) = Even(-a.c1 + num)
 
 #Multiplication
-Base.:(*)(num::Number, a::Even) = Even(num * a.c1)
-Base.:(*)(num::Number, a::Odd) = Odd(num * a.c1)
+Base.:(*)(num::Real, a::Even) = Even(num * a.c1)
+Base.:(*)(num::Real, a::Odd) = Odd(num * a.c1)
 Base.:(*)(a::Even, b::Even) = Even(a.c1 * b.c1)
 Base.:(*)(a::Even, b::Odd) = Odd(conj(a.c1) * b.c1)
 Base.:(*)(a::Odd, b::Even) = Odd(a.c1 * b.c1)

src/core24.jl

@@ -48,12 +48,12 @@ Base.:(-)(a::Odd, b::Odd) = Odd(a.m - b.m)
 
 #Scalar addition / subtraction. Other cases are in GAcommon
 #Relies on Julia's promotion rules to do the sensible thing.
-Base.:(+)(num::Number, a::Even) = Even(a.m + num * one(a.m))
-Base.:(-)(num::Number, a::Even) = Even(-a.m + num * one(a.m))
+Base.:(+)(num::Real, a::Even) = Even(a.m + num * one(a.m))
+Base.:(-)(num::Real, a::Even) = Even(-a.m + num * one(a.m))
 
 #Multiplication
-Base.:(*)(num::Number, a::Even) = Even(num * a.m)
-Base.:(*)(num::Number, a::Odd) = Odd(num * a.m)
+Base.:(*)(num::Real, a::Even) = Even(num * a.m)
+Base.:(*)(num::Real, a::Odd) = Odd(num * a.m)
 Base.:(*)(a::Even, b::Even) = Even(a.m * b.m)
 Base.:(*)(a::Even, b::Odd) = Odd(a.m * b.m)
 Base.:(*)(a::Odd, b::Even) = Odd(a.m * g2.m * conj(b.m) * g2.m)

src/core30.jl

@@ -50,12 +50,12 @@ Base.:(-)(a::Odd, b::Odd) = Odd(a.w - b.w, a.x - b.x, a.y - b.y, a.z - b.z)
 
 # Scalar addition / subtraction. Other cases are in GAcommon
 # Relies on Julia's promotion rules to do the sensible thing.
-Base.:(+)(num::Number, a::Even) = Even(a.w + num, a.x, a.y, a.z)
-Base.:(-)(num::Number, a::Even) = Even(-a.w + num, -a.x, -a.y, -a.z)
+Base.:(+)(num::Real, a::Even) = Even(a.w + num, a.x, a.y, a.z)
+Base.:(-)(num::Real, a::Even) = Even(-a.w + num, -a.x, -a.y, -a.z)
 
 # Multiplication
-Base.:(*)(num::Number, a::Even) = Even(num * a.w, num * a.x, num * a.y, num * a.z)
-Base.:(*)(num::Number, a::Odd) = Odd(num * a.w, num * a.x, num * a.y, num * a.z)
+Base.:(*)(num::Real, a::Even) = Even(num * a.w, num * a.x, num * a.y, num * a.z)
+Base.:(*)(num::Real, a::Odd) = Odd(num * a.w, num * a.x, num * a.y, num * a.z)
 
 function Base.:(*)(a::Even, b::Even)
     return Even(

src/core31.jl

@@ -72,12 +72,12 @@ Base.:(-)(a::Odd, b::Odd) = Odd(a.c1 - b.c1, a.c2 - b.c2, a.c3 - b.c3, a.c4 - b.
 
 #Scalar addition / subtraction. Other cases are in GAcommon
 #Relies on Julia's promotion rules to do the sensible thing.
-Base.:(+)(num::Number, a::Even) = Even(a.c1 + num, a.c2, a.c3, a.c4 + num)
-Base.:(-)(num::Number, a::Even) = Even(-a.c1 + num, -a.c2, -a.c3, -a.c4 + num)
+Base.:(+)(num::Real, a::Even) = Even(a.c1 + num, a.c2, a.c3, a.c4 + num)
+Base.:(-)(num::Real, a::Even) = Even(-a.c1 + num, -a.c2, -a.c3, -a.c4 + num)
 
 #Multiplication
-Base.:(*)(num::Number, a::Even) = Even(num * a.c1, num * a.c2, num * a.c3, num * a.c4)
-Base.:(*)(num::Number, a::Odd) = Odd(num * a.c1, num * a.c2, num * a.c3, num * a.c4)
+Base.:(*)(num::Real, a::Even) = Even(num * a.c1, num * a.c2, num * a.c3, num * a.c4)
+Base.:(*)(num::Real, a::Odd) = Odd(num * a.c1, num * a.c2, num * a.c3, num * a.c4)
 
 function Base.:(*)(a::Even, b::Even)
     return Even(
@@ -162,6 +162,16 @@ end
 LinearAlgebra.norm(a::Even) = sqrt(abs(dot(a, a)))
 LinearAlgebra.norm(a::Odd) = sqrt(abs(dot(a, a)))
 
+"""
+    _as_even(x::Complex) -> Even
+
+Wrap a complex number `x` as an Even type.
+
+We don't expose this so as to avoid leaking our representation abstraction in the definition
+of +, - and *.
+"""
+_as_even(x::Complex) = Even(x, zero(x), zero(x), x)
+
 #Exponentiation
 function SimpleGA.bivector_exp(a::Even)
     a = project(a, 2)
@@ -170,7 +180,7 @@ function SimpleGA.bivector_exp(a::Even)
     return if iszero(fct)
         1 + a
     else
-        cosh(fct) + sinh(fct) / fct * a
+        _as_even(cosh(fct)) + _as_even(sinh(fct) / fct) * a
     end
 end
 
@@ -180,7 +190,7 @@ function Base.exp(a::Even)
     return if iszero(fct)
         R
     else
-        exp(fct) * R
+        _as_even(exp(fct)) * R
     end
 end
 

src/core33.jl

@@ -57,12 +57,12 @@ Base.:(-)(a::Odd, b::Odd) = Odd(a.p - b.p, a.m - b.m)
 
 #Scalar addition / subtraction. Other cases are in GAcommon
 #Relies on Julia's promotion rules to do the sensible thing.
-Base.:(+)(num::Number, a::Even) = Even(a.p + num * one(a.p), a.m + num * one(a.m))
-Base.:(-)(num::Number, a::Even) = Even(-a.p + num * one(a.p), -a.m + num * one(a.m))
+Base.:(+)(num::Real, a::Even) = Even(a.p + num * one(a.p), a.m + num * one(a.m))
+Base.:(-)(num::Real, a::Even) = Even(-a.p + num * one(a.p), -a.m + num * one(a.m))
 
 #Multiplication
-Base.:(*)(num::Number, a::Even) = Even(num * a.p, num * a.m)
-Base.:(*)(num::Number, a::Odd) = Odd(num * a.p, num * a.m)
+Base.:(*)(num::Real, a::Even) = Even(num * a.p, num * a.m)
+Base.:(*)(num::Real, a::Odd) = Odd(num * a.p, num * a.m)
 Base.:(*)(a::Even, b::Even) = Even(a.p * b.p, a.m * b.m)
 Base.:(*)(a::Even, b::Odd) = Odd(a.p * b.p, a.m * b.m)
 Base.:(*)(a::Odd, b::Even) = Odd(a.p * b.m, a.m * b.p)

src/core40.jl

@@ -44,12 +44,12 @@ Base.:(-)(a::Even, b::Even) = Even(a.qp - b.qp, a.qm - b.qm)
 Base.:(-)(a::Odd, b::Odd) = Odd(a.qp - b.qp, a.qm - b.qm)
 
 #Scalar addition / subtraction. Other cases are in GAcommon
-Base.:(+)(num::Number, a::Even) = Even(a.qp + num, a.qm + num)
-Base.:(-)(num::Number, a::Even) = Even(-a.qp + num, -a.qm + num)
+Base.:(+)(num::Real, a::Even) = Even(a.qp + num, a.qm + num)
+Base.:(-)(num::Real, a::Even) = Even(-a.qp + num, -a.qm + num)
 
 #Multiplication
-Base.:(*)(num::Number, a::Even) = Even(num * a.qp, num * a.qm)
-Base.:(*)(num::Number, a::Odd) = Odd(num * a.qp, num * a.qm)
+Base.:(*)(num::Real, a::Even) = Even(num * a.qp, num * a.qm)
+Base.:(*)(num::Real, a::Odd) = Odd(num * a.qp, num * a.qm)
 Base.:(*)(a::Even, b::Even) = Even(a.qp * b.qp, a.qm * b.qm)
 Base.:(*)(a::Even, b::Odd) = Odd(a.qp * b.qp, a.qm * b.qm)
 Base.:(*)(a::Odd, b::Even) = Odd(a.qp * b.qm, a.qm * b.qp)

src/corecga.jl

@@ -70,12 +70,12 @@ Base.:(-)(a::Even, b::Even) = Even(a.q1 - b.q1, a.q2 - b.q2, a.q3 - b.q3, a.q4 -
 Base.:(-)(a::Odd, b::Odd) = Odd(a.q1 - b.q1, a.q2 - b.q2, a.q3 - b.q3, a.q4 - b.q4)
 
 #Scalar addition / subtraction. Other cases are in GAcommon
-Base.:(+)(num::Number, a::Even) = Even(a.q1 + num, a.q2, a.q3, a.q4 + num)
-Base.:(-)(num::Number, a::Even) = Even(-a.q1 + num, -a.q2, -a.q3, -a.q4 + num)
+Base.:(+)(num::Real, a::Even) = Even(a.q1 + num, a.q2, a.q3, a.q4 + num)
+Base.:(-)(num::Real, a::Even) = Even(-a.q1 + num, -a.q2, -a.q3, -a.q4 + num)
 
 #Multiplication
-Base.:(*)(num::Number, a::Even) = Even(num * a.q1, num * a.q2, num * a.q3, num * a.q4)
-Base.:(*)(num::Number, a::Odd) = Odd(num * a.q1, num * a.q2, num * a.q3, num * a.q4)
+Base.:(*)(num::Real, a::Even) = Even(num * a.q1, num * a.q2, num * a.q3, num * a.q4)
+Base.:(*)(num::Real, a::Odd) = Odd(num * a.q1, num * a.q2, num * a.q3, num * a.q4)
 
 function Base.:(*)(a::Even, b::Even)
     return Even(

src/corepga.jl

@@ -48,12 +48,12 @@ Base.:(-)(a::Even, b::Even) = Even(a.q - b.q, a.n - b.n)
 Base.:(-)(a::Odd, b::Odd) = Odd(a.q - b.q, a.n - b.n)
 
 #Scalar addition / subtraction. Other cases are in GAcommon
-Base.:(+)(num::Number, a::Even) = Even(a.q + num, a.n)
-Base.:(-)(num::Number, a::Even) = Even(-a.q + num, -a.n)
+Base.:(+)(num::Real, a::Even) = Even(a.q + num, a.n)
+Base.:(-)(num::Real, a::Even) = Even(-a.q + num, -a.n)
 
 #Multiplication
-Base.:(*)(num::Number, a::Even) = Even(num * a.q, num * a.n)
-Base.:(*)(num::Number, a::Odd) = Odd(num * a.q, num * a.n)
+Base.:(*)(num::Real, a::Even) = Even(num * a.q, num * a.n)
+Base.:(*)(num::Real, a::Odd) = Odd(num * a.q, num * a.n)
 Base.:(*)(a::Even, b::Even) = Even(a.q * b.q, a.q * b.n + a.n * b.q)
 Base.:(*)(a::Even, b::Odd) = Odd(a.q * b.q, a.q * b.n - a.n * b.q)
 Base.:(*)(a::Odd, b::Even) = Odd(a.q * b.q, a.q * b.n + a.n * b.q)