From a03397a0ebdb04094d631fcc535639019e34690f Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 25 Jan 2023 00:54:43 +0100 Subject: [PATCH 1/5] Add extend_analytic --- src/octonion.jl | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/src/octonion.jl b/src/octonion.jl index b237e53..bce00bd 100644 --- a/src/octonion.jl +++ b/src/octonion.jl @@ -126,6 +126,44 @@ function Base.isequal(q::Octonion, w::Octonion) isequal(q.v6, w.v6) & isequal(q.v7, w.v7)) end +""" + 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.exp(o::Octonion) s = o.s se = exp(s) From 47d76915dc39b4894860443702cb6635528185d1 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 25 Jan 2023 00:56:28 +0100 Subject: [PATCH 2/5] Use extend_analytic for exp and sqrt --- src/octonion.jl | 22 ++------------ test/octonion.jl | 75 ++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 75 insertions(+), 22 deletions(-) diff --git a/src/octonion.jl b/src/octonion.jl index bce00bd..3017b4d 100644 --- a/src/octonion.jl +++ b/src/octonion.jl @@ -164,22 +164,8 @@ function extend_analytic(f, o::Octonion) return octo(wr, wi_octo...) 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) +for f in (:sqrt, :exp) + @eval Base.$f(o::Octonion) = extend_analytic($f, o) end function Base.log(o::Octonion) @@ -206,10 +192,6 @@ end Base.:^(o::Octonion, w::Octonion) = exp(w * log(o)) -function Base.sqrt(o::Octonion) - exp(0.5 * log(o)) -end - 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} diff --git a/test/octonion.jl b/test/octonion.jl index f84ec2e..7e5a28c 100644 --- a/test/octonion.jl +++ b/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 From 451a3ac00810cc53b3d35fc8c7231198810f7c5a Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 25 Jan 2023 00:57:16 +0100 Subject: [PATCH 3/5] Use extend_analytic approach for more functions --- src/octonion.jl | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/octonion.jl b/src/octonion.jl index 3017b4d..0962c94 100644 --- a/src/octonion.jl +++ b/src/octonion.jl @@ -164,10 +164,32 @@ function extend_analytic(f, o::Octonion) return octo(wr, wi_octo...) end -for f in (:sqrt, :exp) +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 +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.log(o::Octonion) a = abs(o) o = o / a From 2e45efa504bb0b314204811de31fe9f923970b23 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 25 Jan 2023 00:58:01 +0100 Subject: [PATCH 4/5] Update log --- src/octonion.jl | 29 ++++++++++------------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/src/octonion.jl b/src/octonion.jl index 0962c94..0a59681 100644 --- a/src/octonion.jl +++ b/src/octonion.jl @@ -190,26 +190,17 @@ for f in (@static(VERSION ≥ v"1.6" ? (:sincos, :sincospi) : (:sincos,))) end end +# ~2x faster than extend_analytic(log, o) 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 + 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)) From 8c0525c409ef6781d715361285640e762fdec304 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 25 Jan 2023 01:01:15 +0100 Subject: [PATCH 5/5] Increment version number --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 25abb17..ce3e81f 100644 --- a/Project.toml +++ b/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"