diff --git a/Project.toml b/Project.toml index f3690283..8a8ad4ee 100644 --- a/Project.toml +++ b/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" diff --git a/README.md b/README.md index 204e3158..3d9681a5 100644 --- a/README.md +++ b/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 diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 931d3372..39bd2305 100644 --- a/src/Quaternion.jl +++ b/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) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 8607e7ba..fba6ff3c 100644 --- a/test/test_Quaternion.jl +++ b/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)