Extend all analytic functions in base for octonions

Project.toml

@@ -1,6 +1,6 @@
 name = "Octonions"
 uuid = "d00ba074-1e29-4f5e-9fd4-d67071d6a14d"
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/octonion.jl

@@ -126,52 +126,85 @@ function Base.isequal(q::Octonion, w::Octonion)
             isequal(q.v6, w.v6) & isequal(q.v7, w.v7))
 end
 
-function Base.exp(o::Octonion)
-  s = o.s
-  se = exp(s)
-  scale = se
-  th = abs_imag(o)
-  if th > 0
-    scale *= sin(th) / th
-  end
-  Octonion(se * cos(th),
-            scale * o.v1,
-            scale * o.v2,
-            scale * o.v3,
-            scale * o.v4,
-            scale * o.v5,
-            scale * o.v6,
-            scale * o.v7)
+"""
+    extend_analytic(f, o::Octonion)
+
+Evaluate the extension of the complex analytic function `f` to the octonions at `o`.
+
+Given ``o = s + a u``, where ``s`` is the real part, ``u`` is a pure unit octonion,
+and ``a \\ge 0`` is the magnitude of the imaginary part of ``o``,
+```math
+f(o) = \\Re(f(z)) + \\Im(f(z)) u,
+```
+is the extension of `f` to the octonions, where ``z = s + a i`` is a complex analog to
+``o``.
+
+See [^DentoniSce1973] and [^ColomboSabadini2020] for details.
+
+[^DentoniSce1973]: Dentoni, P. and Sce M. "Funzioni regolari nell'algebra di Cayley."
+                   Rendiconti del Seminario matematico della Università di Padova 50 (1973): 251-267.
+                   Translation: [^ColomboSabadini2020]
+[^ColomboSabadini2020]: Colombo, F., Sabadini, I., Struppa, D.C. (2020).
+                        Regular Functions in the Cayley Algebra.
+                        In: Michele Sce's Works in Hypercomplex Analysis.
+                        doi: [10.1007/978-3-030-50216-4_6](https://doi.org/10.1007/978-3-030-50216-4_6)
+"""
+function extend_analytic(f, o::Octonion)
+    # Adapted from Quaternions.jl
+    a = abs_imag(o)
+    s = o.s
+    z = complex(s, a)
+    w = f(z)
+    wr, wi = reim(w)
+    scale = wi / a
+    # o == real(o), so f(real(o)) may be real or complex, i.e. wi may be nonzero.
+    # we choose to embed complex numbers in the octonions by identifying the first
+    # imaginary octonion basis with the complex imaginary basis.
+    wi_octo = a > 0 ? map(x -> x * scale, imag_part(o)) : ntuple(_ -> zero(scale), Val(7))
+    return octo(wr, wi_octo...)
 end
 
-function Base.log(o::Octonion)
-  a = abs(o)
-  o = o / a
-  s = o.s
-  M = abs_imag(o)
-  th = atan(M, s)
-  if M > 0
-    M = th / M
-    return Octonion(log(a),
-                     o.v1 * M,
-                     o.v2 * M,
-                     o.v3 * M,
-                     o.v4 * M,
-                     o.v5 * M,
-                     o.v6 * M,
-                     o.v7 * M)
-  else
-    z = zero(th)
-    return Octonion(log(a), ifelse(iszero(a), z, th), z, z, z, z, z, z)
-  end
+for f in (:sqrt, :exp, :exp2, :exp10, :expm1, :log2, :log10, :log1p,
+    :sin, :cos, :tan, :asin, :acos, :atan, :sinh, :cosh, :tanh, :asinh, :acosh, :atanh,
+    :csc, :sec, :cot, :acsc, :asec, :acot, :csch, :sech, :coth, :acsch, :asech, :acoth,
+    :sinpi, :cospi,
+)
+    @eval Base.$f(o::Octonion) = extend_analytic($f, o)
 end
 
-Base.:^(o::Octonion, w::Octonion) = exp(w * log(o))
+for f in (@static(VERSION ≥ v"1.6" ? (:sincos, :sincospi) : (:sincos,)))
+    @eval begin
+        function Base.$f(o::Octonion)
+            a = abs_imag(o)
+            z = complex(o.s, a)
+            s, c = $f(z)
+            sr, si = reim(s)
+            cr, ci = reim(c)
+            sscale = si / a
+            cscale = ci / a
+            ov = imag_part(o)
+            si_octo = a > 0 ? map(x -> x * sscale, ov) : ntuple(_ -> zero(sscale), Val(7))
+            ci_octo = a > 0 ? map(x -> x * cscale, ov) : si_octo
+            return octo(sr, si_octo...), octo(cr, ci_octo...)
+        end
+    end
+end
 
-function Base.sqrt(o::Octonion)
-  exp(0.5 * log(o))
+# ~2x faster than extend_analytic(log, o)
+function Base.log(o::Octonion)
+    a = abs_imag(o)
+    theta = atan(a, o.s)
+    scale = theta / a
+    if a > 0
+        return octo(log(abs(o)), map(x -> x * scale, imag_part(o))...)
+    else
+        z = zero(scale)
+        return octo(log(abs(o.s)), oftype(scale, theta), z, z, z, z, z, z)
+    end
 end
 
+Base.:^(o::Octonion, w::Octonion) = exp(w * log(o))
+
 octorand(rng::AbstractRNG = Random.GLOBAL_RNG) = octo(randn(rng), randn(rng), randn(rng), randn(rng), randn(rng), randn(rng), randn(rng), randn(rng))
 
 function Base.rand(rng::AbstractRNG, ::Random.SamplerType{Octonion{T}}) where {T<:Real}

test/octonion.jl

@@ -491,7 +491,15 @@ end
     end
 
     @testset "analytic functions" begin
-        unary_funs = [sqrt, inv, exp, log]
+        # all complex analytic functions can be extended to the octonions
+        #! format: off
+        unary_funs = [
+            sqrt, inv, exp, exp2, exp10, expm1, log, log2, log10, log1p,
+            sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh,
+            csc, sec, cot, acsc, asec, acot, csch, sech, coth, acsch, asech, acoth,
+            sinpi, cospi,
+        ]
+        #! format: on
         # since every octonion is conjugate to a quaternion,
         # one can establish correctness as follows:
         @testset for fun in unary_funs
@@ -512,8 +520,71 @@ end
                 @test exp(log(o)) ≈ o
                 @test exp(zero(o)) === one(o)
                 @test log(one(o)) === zero(o)
+                @test exp2(log2(o)) ≈ o
+                @test exp10(log10(o)) ≈ o
+                @test expm1(log1p(o)) ≈ o
+                @test sinpi(o) ≈ sin(π * o)
+                @test cospi(o) ≈ cos(π * o)
+                @test all(sincos(o) .≈ (sin(o), cos(o)))
+                @test all(sincos(zero(o)) .≈ (sin(zero(o)), cos(zero(o))))
+                if VERSION ≥ v"1.6"
+                    @test all(sincospi(o) .≈ (sinpi(o), cospi(o)))
+                    @test all(sincospi(zero(o)) .≈ (sinpi(zero(o)), cospi(zero(o))))
+                end
+                @test tan(o) ≈ cos(o) \ sin(o) ≈ sin(o) / cos(o)
+                @test tanh(o) ≈ cosh(o) \ sinh(o) ≈ sinh(o) / cosh(o)
+                @testset for (f, finv) in [
+                    (sin, csc),
+                    (cos, sec),
+                    (tan, cot),
+                    (sinh, csch),
+                    (cosh, sech),
+                    (tanh, coth),
+                ]
+                    @test f(o) ≈ inv(finv(o))
+                end
+                @testset for (f, finv) in [
+                    (asin, acsc),
+                    (acos, asec),
+                    (atan, acot),
+                    (asinh, acsch),
+                    (acosh, asech),
+                    (atanh, acoth),
+                ]
+                    @test f(o) ≈ finv(inv(o))
+                end
+            end
+        end
+
+        @testset "additional properties" begin
+            @testset "log" begin
+                @test log(zero(OctonionF64)) === octo(-Inf)
+                @test log(one(OctonionF64)) === octo(0.0)
+                @test log(-one(OctonionF64)) ≈ _octo(log(complex(-1.0)))
+                x = rand()
+                @test log(octo(x)) ≈ octo(log(x))
+                @test log(octo(-x)) ≈ _octo(log(complex(-x)))
+            end
+
+            @testset "exp" begin
+                @test exp(octo(0)) === octo(1.0)
+                @test exp(octo(2)) === octo(exp(2))
+                @test norm(exp(octo(0))) ≈ 1
+                @test norm(exp(octo(2))) ≠ 1
+                @test exp(octo(0.0)) === octo(1.0)
+                for i in 2:8
+                    z = setindex!(zeros(8), 2, i)
+                    z2 = setindex!(zeros(8), sin(2), i)
+                    @test exp(octo(z...)) === octo(cos(2), z2[2:end]...)
+                    @test norm(exp(octo(z...))) ≈ 1
+                    @test exp(octo(2.0)) === octo(exp(2))
+                end
+                @test exp(octo(0)) isa OctonionF64
+                @test exp(octo(0.0)) isa OctonionF64
+                @test exp(octo(0//1)) isa OctonionF64
+                @test exp(octo(BigFloat(0))) isa Octonion{BigFloat}
+                @test exp(octo(fill(1, 8)...)) ≈ exp(octo(fill(1.0, 8)...))
             end
-            @test log(zero(OctonionF64)) === octo(-Inf)
         end
     end