Extend all analytic functions in base for quaternions (#56)

* Fix and test broken conversion

* Add tests for all complex analytic functions

* Add extensions of complex analytic functions

* Remove slower methods

* Avoid invalidations

* Increment version number

* Add test for resolved issue

* Add inverse-specific overloads to avoid computing inverse

* Remove ambiguous case

* Don't overload sincospi if unavailable

* Simplify log implementation

* Test cis and cispi

* Add sinpi and cospi tests

* Add cis and cispi docstrings

* Actually overload method

* Use testset

* Remove real exponent overloads

It takes too many subsequent overloads to resolve ambiguities. Can be handled in a future PR

* Fix and test zero branch of sincos and sincospi

* Add docstring for extend_analytic

* Move randn call into testset

* Move identity tests into a testset

* Run tests many times

* Move exp test with others

* Move exp tests near other funs

* Add additional comment

* Separate exponentiation tests

* Make log more efficient and accurate

* Move exp tests with others

* Apply suggestions from code review

Co-authored-by: Yuto Horikawa <hyrodium@gmail.com>

* Add function overloads to ReadMe

* Remove useless spaces

* Increment version number

Co-authored-by: Yuto Horikawa <hyrodium@gmail.com>

Project.toml

@@ -1,6 +1,6 @@
 name = "Quaternions"
 uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
-version = "0.4.9"
+version = "0.5.0"
 
 [deps]
 DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"

README.md

@@ -17,18 +17,49 @@ Implemented functions are:
     conj
     abs
     abs2
-    exp
-    log
     normalize
     normalizea  (return normalized quaternion and absolute value as a pair)
     angleaxis  (taken as an orientation, return the angle and axis (3 vector) as a tuple)
     angle
     axis
+    sqrt
     exp
-    log
+    exp2
+    exp10
+    expm1
+    log2
+    log10
+    log1p
+    cis
+    cispi
     sin
     cos
-    sqrt
+    tan
+    asin
+    acos
+    atan
+    sinh
+    cosh
+    tanh
+    asinh
+    acosh
+    atanh
+    csc
+    sec
+    cot
+    acsc
+    asec
+    acot
+    csch
+    sech
+    coth
+    acsch
+    asech
+    acoth
+    sinpi
+    cospi
+    sincos
+    sincospi
     linpol  (interpolate between 2 normalized quaternions)
     rand
     randn

src/Quaternion.jl

@@ -115,43 +115,118 @@ end
 
 argq(q::Quaternion) = normalizeq(Quaternion(0, q.v1, q.v2, q.v3))
 
-function exp(q::Quaternion)
-    s = q.s
-    se = exp(s)
-    scale = se
-    th = abs_imag(q)
-    if th > 0
-        scale *= sin(th) / th
-    end
-    Quaternion(se * cos(th), scale * q.v1, scale * q.v2, scale * q.v3, iszero(s))
-end
+"""
+    extend_analytic(f, q::Quaternion)
 
-function log(q::Quaternion)
-    q, a = normalizea(q)
+Evaluate the extension of the complex analytic function `f` to the quaternions at `q`.
+
+Given ``q = s + a u``, where ``s`` is the real part, ``u`` is a pure unit quaternion,
+and ``a \\ge 0`` is the magnitude of the imaginary part of ``q``,
+
+```math
+f(q) = \\Re(f(z)) + \\Im(f(z)) u,
+```
+is the extension of `f` to the quaternions, where ``z = a + s i`` is a complex analog to
+``q``.
+
+See Theorem 5 of [^Sudbery1970] for details.
+
+[^Sudbery1970]
+    Sudbery (1979). Quaternionic analysis. Mathematical Proceedings of the Cambridge 
+    Philosophical Society,85, pp 199­225
+    doi:[10.1017/S030500410005563](https://doi.org/10.1017/S0305004100055638)
+"""
+function extend_analytic(f, q::Quaternion)
+    a = abs_imag(q)
     s = q.s
-    M = abs_imag(q)
-    th = atan(M, s)
-    if M > 0
-        M = th / M
-        return Quaternion(log(a), q.v1 * M, q.v2 * M, q.v3 * M)
+    z = complex(s, a)
+    w = f(z)
+    wr, wi = reim(w)
+    scale = wi / a
+    norm = _isexpfun(f) && iszero(s)
+    if a > 0
+        return Quaternion(wr, scale * q.v1, scale * q.v2, scale * q.v3, norm)
     else
-        return Quaternion(log(a), th, 0.0, 0.0)
+        # q == real(q), so f(real(q)) may be real or complex, i.e. wi may be nonzero.
+        # we choose to embed complex numbers in the quaternions by identifying the first
+        # imaginary quaternion basis with the complex imaginary basis.
+        return Quaternion(wr, oftype(scale, wi), zero(scale), zero(scale), norm)
     end
 end
 
-function sin(q::Quaternion)
-    L = argq(q)
-    return (exp(L * q) - exp(-L * q)) / (2 * L)
+_isexpfun(::Union{typeof(exp),typeof(exp2),typeof(exp10)}) = true
+_isexpfun(::Any) = false
+
+"""
+    cis(q::Quaternion)
+
+Return ``\\exp(u * q)``, where ``u`` is the normalized imaginary part of `q`.
+
+Let ``q = s + a u``, where ``s`` is the real part, ``u`` is a pure unit quaternion,
+and ``a \\ge 0`` is the magnitude of the imaginary part of ``q``.
+
+!!! Note
+    This is the extension of `cis(z)` for complex `z` to the quaternions and is not
+    equivalent to `exp(im * q)`. As a result, `cis(Quaternion(z)) ≠ cis(z)` when
+    `imag(z) < 0`.
+"""
+cis(q::Quaternion)
+
+if VERSION ≥ v"1.6"
+    """
+        cispi(q::Quaternion)
+
+    Compute `cis(π * q)` more accurately.
+
+    !!! Note
+        This is not equivalent to `exp(π*im*q)`. See [cis(::Quaternion)](@ref) for details.
+    """
+    Base.cispi(q::Quaternion) = extend_analytic(cispi, q)
 end
 
-function cos(q::Quaternion)
-    L = argq(q)
-    return (exp(L * q) + exp(-L * q)) / 2
+for f in (
+    :sqrt, :exp, :exp2, :exp10, :expm1, :log2, :log10, :log1p, :cis,
+    :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(q::Quaternion) = extend_analytic($f, q)
 end
 
-(^)(q::Quaternion, w::Quaternion) = exp(w * log(q))
+for f in (@static(VERSION ≥ v"1.6" ? (:sincos, :sincospi) : (:sincos,)))
+    @eval begin
+        function Base.$f(q::Quaternion)
+            a = abs_imag(q)
+            z = complex(q.s, a)
+            s, c = $f(z)
+            sr, si = reim(s)
+            cr, ci = reim(c)
+            sscale = si / a
+            cscale = ci / a
+            if a > 0
+                return (
+                    Quaternion(sr, sscale * q.v1, sscale * q.v2, sscale * q.v3),
+                    Quaternion(cr, cscale * q.v1, cscale * q.v2, cscale * q.v3),
+                )
+            else
+                return (
+                    Quaternion(sr, oftype(sscale, si), zero(sscale), zero(sscale)),
+                    Quaternion(cr, oftype(cscale, ci), zero(cscale), zero(cscale)),
+                )
+            end
+        end
+    end
+end
 
-sqrt(q::Quaternion) = exp(0.5 * log(q))
+function log(q::Quaternion)
+    a = abs(q)
+    M = abs_imag(q)
+    theta = atan(M, q.s)
+    scale = theta / ifelse(iszero(M), oneunit(M), M)
+    return Quaternion(log(a), q.v1 * scale, q.v2 * scale, q.v3 * scale)
+end
+
+(^)(q::Quaternion, w::Quaternion) = exp(w * log(q))
 
 function linpol(p::Quaternion, q::Quaternion, t::Real)
     p = normalize(p)

test/test_Quaternion.jl

@@ -156,22 +156,139 @@ Base.:(/)(a::MyReal, b::Real) = a.val / b
 # this used to throw an error
 @test qrotation([1, 0, 0], MyReal(1.5)) == qrotation([1, 0, 0], 1.5)
 
-for _ in 1:100
-    let # test specialfunctions
-        c = Complex(randn(2)...)
-        q, q2 = sample(Quaternion{Float64}, 4)
-        unary_funs = [exp, log, sin, cos, sqrt, inv, conj, abs2, norm]
-        # since every quaternion is conjugate to a complex number,
-        # one can establish correctness as follows:
-        for fun in unary_funs
-            @test fun(Quaternion(c)) ≈ Quaternion(fun(c))
+@testset "non-analytic functions" begin
+    q, q2 = randn(Quaternion{Float64}, 2)
+    unary_funs = [conj, abs, abs2, norm, sign]
+    # since every quaternion is conjugate to a complex number,
+    # one can establish correctness as follows:
+    @testset for fun in unary_funs
+        for _ in 1:100
+            c = randn(ComplexF64)
+            @test fun(Quaternion(c)) ≈ fun(c)
             @test q2 * fun(q) * inv(q2) ≈ fun(q2 * q * inv(q2))
         end
+    end
+end
+
+@testset "extended complex analytic functions" begin
+    # all complex analytic functions can be extended to the quaternions
+    unary_funs = [
+        sqrt, inv, exp, exp2, exp10, expm1, log, log2, log10, log1p, cis,
+        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,
+    ]
+    # since every quaternion is conjugate to a complex number,
+    # one can establish correctness as follows:
+    @testset for fun in unary_funs
+        q, q2 = randn(QuaternionF64, 2)
+        for _ in 1:100
+            c = randn(ComplexF64)
+            fun !== cis && @test fun(Quaternion(c)) ≈ fun(c)
+            @test q2 * fun(q) * inv(q2) ≈ fun(q2 * q * inv(q2))
+        end
+    end
 
-        @test exp(log(q)) ≈ q
-        @test exp(zero(q)) ≈ one(q)
+    @testset "identities" begin
+        for _ in 1:100
+            q = randn(QuaternionF64)
+            @test inv(q) * q ≈ q * inv(q) ≈ one(q)
+            @test sqrt(q) * sqrt(q) ≈ q
+            @test exp(log(q)) ≈ q
+            @test exp(zero(q)) ≈ one(q)
+            @test log(one(q)) ≈ zero(q)
+            @test exp2(log2(q)) ≈ q
+            @test exp10(log10(q)) ≈ q
+            @test expm1(log1p(q)) ≈ q
+            @test sinpi(q) ≈ sin(π * q)
+            @test cospi(q) ≈ cos(π * q)
+            @test all(sincos(q) .≈ (sin(q), cos(q)))
+            @test all(sincos(zero(q)) .≈ (sin(zero(q)), cos(zero(q))))
+            if VERSION ≥ v"1.6"
+                @test all(sincospi(q) .≈ (sinpi(q), cospi(q)))
+                @test all(sincospi(zero(q)) .≈ (sinpi(zero(q)), cospi(zero(q))))
+            end
+            @test tan(q) ≈ cos(q) \ sin(q) ≈ sin(q) / cos(q)
+            @test tanh(q) ≈ cosh(q) \ sinh(q) ≈ sinh(q) / cosh(q)
+            @testset for (f, finv) in [(sin, csc), (cos, sec), (tan, cot), (sinh, csch), (cosh, sech), (tanh, coth)]
+                @test f(q) ≈ inv(finv(q))
+            end
+            @testset for (f, finv) in [(asin, acsc), (acos, asec), (atan, acot), (asinh, acsch), (acosh, asech), (atanh, acoth)]
+                @test f(q) ≈ finv(inv(q))
+            end
+            @test cis(q) ≈ exp(normalize(q - real(q)) * q)
+            VERSION ≥ v"1.6" && @test cispi(q) ≈ cis(π * q)
+        end
     end
 
+    @testset "additional properties" begin
+        @testset "log" begin
+            @test log(zero(QuaternionF64)) === Quaternion(-Inf, 0, 0, 0)
+        end
+
+        @testset "exp" begin
+            @test exp(Quaternion(0, 0, 0, 0)) === Quaternion(1., 0., 0., 0., true)
+            @test exp(Quaternion(2, 0, 0, 0)) === Quaternion(exp(2), 0, 0, 0, false)
+            @test exp(Quaternion(0, 2, 0, 0)) === Quaternion(cos(2), sin(2), 0, 0, true)
+            @test exp(Quaternion(0, 0, 2, 0)) === Quaternion(cos(2), 0, sin(2), 0, true)
+            @test exp(Quaternion(0, 0, 0, 2)) === Quaternion(cos(2), 0, 0, sin(2), true)
+
+            @test norm(exp(Quaternion(0, 0, 0, 0))) ≈ 1
+            @test norm(exp(Quaternion(2, 0, 0, 0))) ≠ 1
+            @test norm(exp(Quaternion(0, 2, 0, 0))) ≈ 1
+            @test norm(exp(Quaternion(0, 0, 2, 0))) ≈ 1
+            @test norm(exp(Quaternion(0, 0, 0, 2))) ≈ 1
+
+            @test exp(Quaternion(0., 0., 0., 0.)) === Quaternion(1., 0., 0., 0., true)
+            @test exp(Quaternion(2., 0., 0., 0.)) === Quaternion(exp(2), 0, 0, 0, false)
+            @test exp(Quaternion(0., 2., 0., 0.)) === Quaternion(cos(2), sin(2), 0, 0, true)
+            @test exp(Quaternion(0., 0., 2., 0.)) === Quaternion(cos(2), 0, sin(2), 0, true)
+            @test exp(Quaternion(0., 0., 0., 2.)) === Quaternion(cos(2), 0, 0, sin(2), true)
+
+            @test norm(exp(Quaternion(0., 0., 0., 0.))) ≈ 1
+            @test norm(exp(Quaternion(2., 0., 0., 0.))) ≠ 1
+            @test norm(exp(Quaternion(0., 2., 0., 0.))) ≈ 1
+            @test norm(exp(Quaternion(0., 0., 2., 0.))) ≈ 1
+            @test norm(exp(Quaternion(0., 0., 0., 2.))) ≈ 1
+
+            @test exp(Quaternion(0,0,0,0)) isa Quaternion{Float64}
+            @test exp(Quaternion(0.,0,0,0)) isa Quaternion{Float64}
+            @test exp(Quaternion(0//1,0,0,0)) isa Quaternion{Float64}
+            @test exp(Quaternion(BigFloat(0),0,0,0)) isa Quaternion{BigFloat}
+
+            # https://github.com/JuliaGeometry/Quaternions.jl/issues/39
+            @testset "exp(::Quaternion{Int})" begin
+                @test exp(Quaternion(1,1,1,1)) ≈ exp(Quaternion(1.0,1.0,1.0,1.0))
+            end
+        end
+    end
+end
+
+@testset "^" begin
+    @testset "^(::Quaternion, ::Real)" begin
+        for _ in 1:100
+            q = randn(QuaternionF64)
+            @test q^2.0 ≈ q * q
+            @test q^1.0 ≈ q
+            @test q^-1.0 ≈ inv(q)
+            @test q^1.3 ≈ exp(1.3 * log(q))
+            @test q^7.8 ≈ exp(7.8 * log(q))
+            @test q^1.3f0 ≈ exp(1.3f0 * log(q))
+            @test q^7.8f0 ≈ exp(7.8f0 * log(q))
+        end
+    end
+    @testset "^(::Quaternion, ::Quaternion)" begin
+        @test Quaternion(ℯ,0,0,0)^Quaternion(0,0,π/2,0) ≈ Quaternion(0,0,1,0)
+        @test Quaternion(3.5,0,0,2.3)^Quaternion(0.2,0,0,1.7) ≈
+            Quaternion(real((3.5+2.3im)^(0.2+1.7im)),0,0,imag((3.5+2.3im)^(0.2+1.7im)))
+        for _ in 1:100
+            q, p = randn(QuaternionF64, 2)
+            @test q^p ≈ exp(p * log(q))
+        end
+    end
+end
+
+for _ in 1:100
     let # test qrotation and angleaxis inverse
         ax = randn(3); ax = ax / norm(ax)
         Θ = π * rand()
@@ -230,37 +347,6 @@ for _ in 1:100
     end
 end
 
-@testset "exp" begin
-    @test exp(Quaternion(0, 0, 0, 0)) === Quaternion(1., 0., 0., 0., true)
-    @test exp(Quaternion(2, 0, 0, 0)) === Quaternion(exp(2), 0, 0, 0, false)
-    @test exp(Quaternion(0, 2, 0, 0)) === Quaternion(cos(2), sin(2), 0, 0, true)
-    @test exp(Quaternion(0, 0, 2, 0)) === Quaternion(cos(2), 0, sin(2), 0, true)
-    @test exp(Quaternion(0, 0, 0, 2)) === Quaternion(cos(2), 0, 0, sin(2), true)
-
-    @test norm(exp(Quaternion(0, 0, 0, 0))) ≈ 1
-    @test norm(exp(Quaternion(2, 0, 0, 0))) ≠ 1
-    @test norm(exp(Quaternion(0, 2, 0, 0))) ≈ 1
-    @test norm(exp(Quaternion(0, 0, 2, 0))) ≈ 1
-    @test norm(exp(Quaternion(0, 0, 0, 2))) ≈ 1
-
-    @test exp(Quaternion(0., 0., 0., 0.)) === Quaternion(1., 0., 0., 0., true)
-    @test exp(Quaternion(2., 0., 0., 0.)) === Quaternion(exp(2), 0, 0, 0, false)
-    @test exp(Quaternion(0., 2., 0., 0.)) === Quaternion(cos(2), sin(2), 0, 0, true)
-    @test exp(Quaternion(0., 0., 2., 0.)) === Quaternion(cos(2), 0, sin(2), 0, true)
-    @test exp(Quaternion(0., 0., 0., 2.)) === Quaternion(cos(2), 0, 0, sin(2), true)
-
-    @test norm(exp(Quaternion(0., 0., 0., 0.))) ≈ 1
-    @test norm(exp(Quaternion(2., 0., 0., 0.))) ≠ 1
-    @test norm(exp(Quaternion(0., 2., 0., 0.))) ≈ 1
-    @test norm(exp(Quaternion(0., 0., 2., 0.))) ≈ 1
-    @test norm(exp(Quaternion(0., 0., 0., 2.))) ≈ 1
-
-    @test exp(Quaternion(0,0,0,0)) isa Quaternion{Float64}
-    @test exp(Quaternion(0.,0,0,0)) isa Quaternion{Float64}
-    @test exp(Quaternion(0//1,0,0,0)) isa Quaternion{Float64}
-    @test exp(Quaternion(BigFloat(0),0,0,0)) isa Quaternion{BigFloat}
-end
-
 @testset "random quaternions" begin
     @testset "quatrand" begin
         rng = Random.MersenneTwister(42)