From 62528e6c092da0173cb9483ccfd15c979a6fd8d0 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 02:52:28 +0100 Subject: [PATCH 01/32] Fix and test broken conversion --- src/Quaternion.jl | 2 +- test/test_Quaternion.jl | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 91d06c07..a270ef26 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -19,7 +19,7 @@ Quaternion(a::Vector) = Quaternion(0, a[1], a[2], a[3]) convert(::Type{Quaternion{T}}, x::Real) where {T} = Quaternion(convert(T, x), convert(T, 0), convert(T, 0), convert(T, 0)) convert(::Type{Quaternion{T}}, z::Complex) where {T} = - Quaternion(convert(T, real(z)), convert(T, imag(z)), convert(T, 0), convert(T, 0)) + Quaternion(convert(T, real(z)), convert(T, Base.imag(z)), convert(T, 0), convert(T, 0)) convert(::Type{Quaternion{T}}, q::Quaternion{T}) where {T <: Real} = q convert(::Type{Quaternion{T}}, q::Quaternion) where {T} = Quaternion(convert(T, q.s), convert(T, q.v1), convert(T, q.v2), convert(T, q.v3), q.norm) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index dea5da7c..dff1c0f8 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -67,6 +67,15 @@ let # test rotations end end +@testset "conversions" begin + z = rand(ComplexF64) + q = convert(Quaternion{Float64}, z) + @test q.s == real(z) + @test q.v1 == imag(z) + @test q.v2 == 0 + @test q.v3 == 0 +end + for _ in 1:100 let # test specialfunctions c = Complex(randn(2)...) From 72349b83f77f464df2cef0d6e493313779ceeaf5 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 02:53:21 +0100 Subject: [PATCH 02/32] Add tests for all complex analytic functions --- test/test_Quaternion.jl | 61 +++++++++++++++++++++++++++++++++-------- 1 file changed, 50 insertions(+), 11 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index dff1c0f8..11179099 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -76,22 +76,61 @@ end @test q.v3 == 0 end -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 i 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 "complex analytic functions" begin + # all complex analytic functions can be extended to the quaternions + q, q2 = randn(Quaternion{Float64}, 2) + 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, + ] + # since every quaternion is conjugate to a complex number, + # one can establish correctness as follows: + @testset for fun in unary_funs + for i 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 - @test exp(log(q)) ≈ q - @test exp(zero(q)) ≈ one(q) + q = randn(Quaternion{Float64}) + @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 all(sincos(q) .≈ (sin(q), cos(q))) + @test all(sincospi(q) .≈ (sinpi(q), cospi(q))) + @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 +end +for _ in 1:100 let # test qrotation and angleaxis inverse ax = randn(3); ax = ax / norm(ax) Θ = π * rand() From 10182f8c03bfaaaa25caba9b8e86401aa94970c6 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 02:54:56 +0100 Subject: [PATCH 03/32] Add extensions of complex analytic functions --- src/Quaternion.jl | 48 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index a270ef26..f3c80bb1 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -110,6 +110,54 @@ end argq(q::Quaternion) = normalizeq(Quaternion(0, q.v1, q.v2, q.v3)) +# extend complex analytic function f to the quaternion q +# adapted from Theorem 5 of doi.org/10.1017/S0305004100055638 +function extend_analytic(f, q::Quaternion) + a = abs_imag(q) + z = complex(q.s, a) + w = f(z) + wr, wi = reim(w) + scale = wi / a + if a > 0 + return Quaternion(wr, scale * q.v1, scale * q.v2, scale * q.v3) + else # quaternion may be real or complex + return Quaternion(wr, oftype(scale, wi), zero(scale), zero(scale)) + end +end + +for f in ( + :exp2, :exp10, :expm1, :log2, :log10, :log1p, + :tan, :asin, :acos, :atan, :sinh, :cosh, :tanh, :asinh, :acosh, :atanh, + :sinpi, :cospi, +) + @eval Base.$f(q::Quaternion) = extend_analytic($f, q) +end + +for f in (:sincos, :sincospi) + @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(scale), zero(scale)), + Quaternion(cr, oftype(cscale, ci), zero(scale), zero(scale)), + ) + end + end + end +end + function exp(q::Quaternion) s = q.s se = exp(s) From b253d304861b7e17ca9792b6becd6ce8947a3233 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 03:14:19 +0100 Subject: [PATCH 04/32] Remove slower methods --- src/Quaternion.jl | 39 +++++++++++---------------------------- 1 file changed, 11 insertions(+), 28 deletions(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index f3c80bb1..771c055a 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -114,20 +114,25 @@ argq(q::Quaternion) = normalizeq(Quaternion(0, q.v1, q.v2, q.v3)) # adapted from Theorem 5 of doi.org/10.1017/S0305004100055638 function extend_analytic(f, q::Quaternion) a = abs_imag(q) - z = complex(q.s, a) + s = q.s + z = complex(s, a) w = f(z) wr, wi = reim(w) scale = wi / a + norm = _isexpfun(f) && abs(s) < eps(float(typeof(s))) if a > 0 - return Quaternion(wr, scale * q.v1, scale * q.v2, scale * q.v3) + return Quaternion(wr, scale * q.v1, scale * q.v2, scale * q.v3, norm) else # quaternion may be real or complex - return Quaternion(wr, oftype(scale, wi), zero(scale), zero(scale)) + return Quaternion(wr, oftype(scale, wi), zero(scale), zero(scale), norm) end end +_isexpfun(::Union{typeof(exp),typeof(exp2),typeof(exp10)}) = true +_isexpfun(::Any) = false + for f in ( - :exp2, :exp10, :expm1, :log2, :log10, :log1p, - :tan, :asin, :acos, :atan, :sinh, :cosh, :tanh, :asinh, :acosh, :atanh, + :sqrt, :exp, :exp2, :exp10, :expm1, :log2, :log10, :log1p, + :sin, :cos, :tan, :asin, :acos, :atan, :sinh, :cosh, :tanh, :asinh, :acosh, :atanh, :sinpi, :cospi, ) @eval Base.$f(q::Quaternion) = extend_analytic($f, q) @@ -158,17 +163,6 @@ for f in (:sincos, :sincospi) end end -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, abs(s) < eps(typeof(s))) -end - function log(q::Quaternion) q, a = normalizea(q) s = q.s @@ -182,20 +176,9 @@ function log(q::Quaternion) end end -function sin(q::Quaternion) - L = argq(q) - return (exp(L * q) - exp(-L * q)) / (2 * L) -end - -function cos(q::Quaternion) - L = argq(q) - return (exp(L * q) + exp(-L * q)) / 2 -end - +(^)(q::Quaternion, p::Real) = extend_analytic(Base.Fix2(^, p), q) (^)(q::Quaternion, w::Quaternion) = exp(w * log(q)) -sqrt(q::Quaternion) = exp(0.5 * log(q)) - function linpol(p::Quaternion, q::Quaternion, t::Real) p = normalize(p) q = normalize(q) From c5cde18beda911781e936741660f1cf8bfdf6df1 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 03:18:07 +0100 Subject: [PATCH 05/32] Avoid invalidations --- src/Quaternion.jl | 3 ++- test/test_Quaternion.jl | 4 ++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 771c055a..ee5cbaa4 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -176,7 +176,8 @@ function log(q::Quaternion) end end -(^)(q::Quaternion, p::Real) = extend_analytic(Base.Fix2(^, p), q) +(^)(q::Quaternion{T}, p::T) where {T<:Real} = extend_analytic(Base.Fix2(^, p), q) +(^)(q::Quaternion{T}, p::S) where {T<:Real,S<:Real} = extend_analytic(Base.Fix2(^, p), q) (^)(q::Quaternion, w::Quaternion) = exp(w * log(q)) function linpol(p::Quaternion, q::Quaternion, t::Real) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 11179099..4b1a477a 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -128,6 +128,10 @@ 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 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 for _ in 1:100 From 98df88a495cc6f96e5b3446f55ae9ec29fdd4cbb Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 03:18:29 +0100 Subject: [PATCH 06/32] Increment version number --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8342884d..d9c71c9d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Quaternions" uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" -version = "0.4.3" +version = "0.4.4" [deps] DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" From 22847afdf35f475c2a99a610bb2558c6d1e0266e Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 03:22:25 +0100 Subject: [PATCH 07/32] Add test for resolved issue --- test/test_Quaternion.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 4b1a477a..292db537 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -134,6 +134,11 @@ end @test q^7.8f0 ≈ exp(7.8f0 * log(q)) end +let + # https://github.com/JuliaGeometry/Quaternions.jl/issues/39 + @test exp(Quaternion(1,1,1,1)) ≈ exp(Quaternion(1.0,1.0,1.0,1.0)) +end + for _ in 1:100 let # test qrotation and angleaxis inverse ax = randn(3); ax = ax / norm(ax) From 06b472d11f1b6fa608841109cbcea2b28e834770 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 03:25:33 +0100 Subject: [PATCH 08/32] Add inverse-specific overloads to avoid computing inverse --- src/Quaternion.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index ee5cbaa4..8633185e 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -133,6 +133,7 @@ _isexpfun(::Any) = false 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(q::Quaternion) = extend_analytic($f, q) From b8932bce0b88bd678fcfd0140846a94706953682 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 03:32:46 +0100 Subject: [PATCH 09/32] Remove ambiguous case --- src/Quaternion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 8633185e..8f59293e 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -178,7 +178,7 @@ function log(q::Quaternion) end (^)(q::Quaternion{T}, p::T) where {T<:Real} = extend_analytic(Base.Fix2(^, p), q) -(^)(q::Quaternion{T}, p::S) where {T<:Real,S<:Real} = extend_analytic(Base.Fix2(^, p), q) +(^)(q::Quaternion{T}, p::AbstractFloat) where {T<:Real} = extend_analytic(Base.Fix2(^, p), q) (^)(q::Quaternion, w::Quaternion) = exp(w * log(q)) function linpol(p::Quaternion, q::Quaternion, t::Real) From 9de628721ea40483c20d2f16c0d2b82c015b8a01 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sat, 19 Feb 2022 03:39:43 +0100 Subject: [PATCH 10/32] Don't overload sincospi if unavailable --- src/Quaternion.jl | 2 +- test/test_Quaternion.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 8f59293e..f42c0ad6 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -139,7 +139,7 @@ for f in ( @eval Base.$f(q::Quaternion) = extend_analytic($f, q) end -for f in (:sincos, :sincospi) +for f in (@static(VERSION ≥ v"1.6" ? (:sincos, :sincospi) : (:sincos,))) @eval begin function Base.$f(q::Quaternion) a = abs_imag(q) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 292db537..14d60569 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -119,7 +119,7 @@ end @test exp10(log10(q)) ≈ q @test expm1(log1p(q)) ≈ q @test all(sincos(q) .≈ (sin(q), cos(q))) - @test all(sincospi(q) .≈ (sinpi(q), cospi(q))) + VERSION ≥ v"1.6" && @test all(sincospi(q) .≈ (sinpi(q), cospi(q))) @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)] From ed7fc54d0a026d40fdb9dc32b27e378e62c7eb9e Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sun, 20 Feb 2022 13:40:56 +0100 Subject: [PATCH 11/32] Simplify log implementation --- src/Quaternion.jl | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index ab800f07..b4c49629 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -166,15 +166,10 @@ end function log(q::Quaternion) q, a = normalizea(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) - else - return Quaternion(log(a), th, 0.0, 0.0) - end + 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{T}, p::T) where {T<:Real} = extend_analytic(Base.Fix2(^, p), q) From 62426d000b0b5b9b9569d08120a5067b4fb3e24f Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sun, 20 Feb 2022 17:25:30 +0100 Subject: [PATCH 12/32] Test cis and cispi --- test/test_Quaternion.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 96a114b8..acac8a18 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -106,7 +106,7 @@ end # all complex analytic functions can be extended to the quaternions q, q2 = randn(Quaternion{Float64}, 2) unary_funs = [ - sqrt, inv, exp, exp2, exp10, expm1, log, log2, log10, log1p, + 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, @@ -116,7 +116,7 @@ end @testset for fun in unary_funs for i in 1:100 c = randn(ComplexF64) - @test fun(Quaternion(c)) ≈ fun(c) + fun !== cis && @test fun(Quaternion(c)) ≈ fun(c) @test q2 * fun(q) * inv(q2) ≈ fun(q2 * q * inv(q2)) end end @@ -144,6 +144,8 @@ end @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)) + @test cis(q) ≈ exp(normalize(q - real(q)) * q) + VERSION ≥ v"1.6" && @test cispi(q) ≈ cis(π * q) end let From 0ddb9b97a0537179747e5ef57ce38ae054d8c656 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sun, 20 Feb 2022 17:25:38 +0100 Subject: [PATCH 13/32] Add sinpi and cospi tests --- test/test_Quaternion.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index acac8a18..7583d156 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -130,6 +130,8 @@ end @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))) VERSION ≥ v"1.6" && @test all(sincospi(q) .≈ (sinpi(q), cospi(q))) @test tan(q) ≈ cos(q) \ sin(q) ≈ sin(q) / cos(q) From 386cc78292e27f519eb8a84e9e58d224dfbfbb01 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sun, 20 Feb 2022 17:26:17 +0100 Subject: [PATCH 14/32] Add cis and cispi docstrings --- src/Quaternion.jl | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index b4c49629..2e8d7d3d 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -130,8 +130,35 @@ end _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. + """ + cispi(q::Quaternion) = extend_analytic(cispi, q) +end + for f in ( - :sqrt, :exp, :exp2, :exp10, :expm1, :log2, :log10, :log1p, + :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, From 0a9f1fe540dd41c8e5f0fa5de315b2477220087e Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sun, 20 Feb 2022 17:29:23 +0100 Subject: [PATCH 15/32] Actually overload method --- src/Quaternion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 2e8d7d3d..355dee72 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -154,7 +154,7 @@ if VERSION ≥ v"1.6" !!! Note This is not equivalent to `exp(π*im*q)`. See [cis(::Quaternion)](@ref) for details. """ - cispi(q::Quaternion) = extend_analytic(cispi, q) + Base.cispi(q::Quaternion) = extend_analytic(cispi, q) end for f in ( From 27634cd37670f2119be91c137bbcffa06c5c4caf Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sun, 20 Feb 2022 17:38:26 +0100 Subject: [PATCH 16/32] Use testset --- test/test_Quaternion.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 7583d156..a2653e65 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -150,8 +150,8 @@ end VERSION ≥ v"1.6" && @test cispi(q) ≈ cis(π * q) end -let - # https://github.com/JuliaGeometry/Quaternions.jl/issues/39 +# 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 From 2d2a945d9de90651d316a12aae3da666370eb092 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sun, 20 Feb 2022 18:11:31 +0100 Subject: [PATCH 17/32] Remove real exponent overloads It takes too many subsequent overloads to resolve ambiguities. Can be handled in a future PR --- src/Quaternion.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 355dee72..7b0c988d 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -199,8 +199,6 @@ function log(q::Quaternion) return Quaternion(log(a), q.v1 * scale, q.v2 * scale, q.v3 * scale) end -(^)(q::Quaternion{T}, p::T) where {T<:Real} = extend_analytic(Base.Fix2(^, p), q) -(^)(q::Quaternion{T}, p::AbstractFloat) where {T<:Real} = extend_analytic(Base.Fix2(^, p), q) (^)(q::Quaternion, w::Quaternion) = exp(w * log(q)) function linpol(p::Quaternion, q::Quaternion, t::Real) From 21390e91829d47e42ed446b78ab217c9a77a3d31 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sun, 20 Feb 2022 18:16:44 +0100 Subject: [PATCH 18/32] Fix and test zero branch of sincos and sincospi --- src/Quaternion.jl | 4 ++-- test/test_Quaternion.jl | 6 +++++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 7b0c988d..ee7f7742 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -183,8 +183,8 @@ for f in (@static(VERSION ≥ v"1.6" ? (:sincos, :sincospi) : (:sincos,))) ) else return ( - Quaternion(sr, oftype(sscale, si), zero(scale), zero(scale)), - Quaternion(cr, oftype(cscale, ci), zero(scale), zero(scale)), + Quaternion(sr, oftype(sscale, si), zero(sscale), zero(sscale)), + Quaternion(cr, oftype(cscale, ci), zero(cscale), zero(cscale)), ) end end diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index a2653e65..e876519c 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -133,7 +133,11 @@ end @test sinpi(q) ≈ sin(π * q) @test cospi(q) ≈ cos(π * q) @test all(sincos(q) .≈ (sin(q), cos(q))) - VERSION ≥ v"1.6" && @test all(sincospi(q) .≈ (sinpi(q), cospi(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)] From bf4085734e9825570a5aba10e3d3995f23bf9df1 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Sun, 20 Feb 2022 19:09:25 +0100 Subject: [PATCH 19/32] Add docstring for extend_analytic --- src/Quaternion.jl | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index b00e9b48..0ca6a00f 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -110,8 +110,27 @@ end argq(q::Quaternion) = normalizeq(Quaternion(0, q.v1, q.v2, q.v3)) -# extend complex analytic function f to the quaternion q -# adapted from Theorem 5 of doi.org/10.1017/S0305004100055638 +""" + extend_analytic(f, q::Quaternion) + +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 From bb7a599731599dbcb711fae6c67ef95fe10ff6be Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 22 Feb 2022 20:56:30 +0100 Subject: [PATCH 20/32] Move randn call into testset --- test/test_Quaternion.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 14411e97..cbe40ccf 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -172,7 +172,6 @@ end @testset "complex analytic functions" begin # all complex analytic functions can be extended to the quaternions - q, q2 = randn(Quaternion{Float64}, 2) 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, @@ -182,6 +181,7 @@ end # 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 i in 1:100 c = randn(ComplexF64) fun !== cis && @test fun(Quaternion(c)) ≈ fun(c) From 91f62d1ae30a71d4de421e2d88387abbede0a07e Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 22 Feb 2022 20:57:45 +0100 Subject: [PATCH 21/32] Move identity tests into a testset --- test/test_Quaternion.jl | 62 +++++++++++++++++++++-------------------- 1 file changed, 32 insertions(+), 30 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index cbe40ccf..f4a295fe 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -189,37 +189,39 @@ end end end - q = randn(Quaternion{Float64}) - @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)) + @testset "identities for analytic functions" + 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 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)) + @test cis(q) ≈ exp(normalize(q - real(q)) * q) + VERSION ≥ v"1.6" && @test cispi(q) ≈ cis(π * q) end - @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)) - @test cis(q) ≈ exp(normalize(q - real(q)) * q) - VERSION ≥ v"1.6" && @test cispi(q) ≈ cis(π * q) end # https://github.com/JuliaGeometry/Quaternions.jl/issues/39 From 3b5308816830383ff99bd271ea5444fe168ccb3c Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 22 Feb 2022 20:58:38 +0100 Subject: [PATCH 22/32] Run tests many times --- test/test_Quaternion.jl | 64 +++++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 31 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index f4a295fe..1cd0cf82 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -189,38 +189,40 @@ end end end - @testset "identities for analytic functions" - 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)) + @testset "identities for analytic functions" 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 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)) + @test cis(q) ≈ exp(normalize(q - real(q)) * q) + VERSION ≥ v"1.6" && @test cispi(q) ≈ cis(π * q) end - @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)) - @test cis(q) ≈ exp(normalize(q - real(q)) * q) - VERSION ≥ v"1.6" && @test cispi(q) ≈ cis(π * q) end end From b3d291ebdb1ce0624cb0b78c3ab16f8e06c1e512 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 22 Feb 2022 20:59:20 +0100 Subject: [PATCH 23/32] Move exp test with others --- test/test_Quaternion.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 1cd0cf82..ad9d9add 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -226,11 +226,6 @@ end end end -# 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 - for _ in 1:100 let # test qrotation and angleaxis inverse ax = randn(3); ax = ax / norm(ax) @@ -319,6 +314,11 @@ end @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 @testset "random quaternions" begin From a457fa3ecde1bf27df8cdbfc25213bd441f14cd9 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 22 Feb 2022 21:00:03 +0100 Subject: [PATCH 24/32] Move exp tests near other funs --- test/test_Quaternion.jl | 72 ++++++++++++++++++++--------------------- 1 file changed, 36 insertions(+), 36 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index ad9d9add..f9807039 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -226,6 +226,42 @@ end 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} + + # 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 + for _ in 1:100 let # test qrotation and angleaxis inverse ax = randn(3); ax = ax / norm(ax) @@ -285,42 +321,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} - - # 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 - @testset "random quaternions" begin @testset "quatrand" begin rng = Random.MersenneTwister(42) From 2bdb4fb77ad73cdeb6dfb49148422d131a1c6626 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 22 Feb 2022 21:08:28 +0100 Subject: [PATCH 25/32] Add additional comment --- src/Quaternion.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 3b005f9c..75d31e6f 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -146,7 +146,10 @@ function extend_analytic(f, q::Quaternion) norm = _isexpfun(f) && iszero(s) if a > 0 return Quaternion(wr, scale * q.v1, scale * q.v2, scale * q.v3, norm) - else # quaternion may be real or complex + else + # 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 From 1f620bb0284e335fc5e9cd340204852abde8e029 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 22 Feb 2022 21:25:27 +0100 Subject: [PATCH 26/32] Separate exponentiation tests --- test/test_Quaternion.jl | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index f9807039..91c794ba 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -216,10 +216,6 @@ 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 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)) @test cis(q) ≈ exp(normalize(q - real(q)) * q) VERSION ≥ v"1.6" && @test cispi(q) ≈ cis(π * q) end @@ -262,6 +258,30 @@ 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) From 84d7989d3ffdc5a68b9ac4d063a572a6af60b832 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 22 Feb 2022 21:41:27 +0100 Subject: [PATCH 27/32] Make log more efficient and accurate --- src/Quaternion.jl | 2 +- test/test_Quaternion.jl | 10 ++++++++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 75d31e6f..39bd2305 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -219,7 +219,7 @@ for f in (@static(VERSION ≥ v"1.6" ? (:sincos, :sincospi) : (:sincos,))) end function log(q::Quaternion) - q, a = normalizea(q) + a = abs(q) M = abs_imag(q) theta = atan(M, q.s) scale = theta / ifelse(iszero(M), oneunit(M), M) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 91c794ba..8b83c3c6 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -170,7 +170,7 @@ Base.:(/)(a::MyReal, b::Real) = a.val / b end end -@testset "complex analytic functions" begin +@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, @@ -189,7 +189,7 @@ end end end - @testset "identities for analytic functions" begin + @testset "identities" begin for _ in 1:100 q = randn(QuaternionF64) @test inv(q) * q ≈ q * inv(q) ≈ one(q) @@ -220,6 +220,12 @@ end 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 + end end @testset "exp" begin From 9ca83ab8b91c88a3fea29c1905fa896251d25673 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Tue, 22 Feb 2022 21:42:16 +0100 Subject: [PATCH 28/32] Move exp tests with others --- test/test_Quaternion.jl | 70 ++++++++++++++++++++--------------------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 8b83c3c6..d6d0e68b 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -225,42 +225,42 @@ end @testset "log" begin @test log(zero(QuaternionF64)) === Quaternion(-Inf, 0, 0, 0) end - 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} - - # 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)) + @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 From 3d20d43233d4b7639298c6f5139a08b740bf5142 Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 23 Feb 2022 10:26:46 +0100 Subject: [PATCH 29/32] Apply suggestions from code review Co-authored-by: Yuto Horikawa --- test/test_Quaternion.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 44a6c48b..9ce1742d 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -162,7 +162,7 @@ Base.:(/)(a::MyReal, b::Real) = a.val / b # since every quaternion is conjugate to a complex number, # one can establish correctness as follows: @testset for fun in unary_funs - for i in 1:100 + for _ in 1:100 c = randn(ComplexF64) @test fun(Quaternion(c)) ≈ fun(c) @test q2 * fun(q) * inv(q2) ≈ fun(q2 * q * inv(q2)) @@ -182,7 +182,7 @@ end # one can establish correctness as follows: @testset for fun in unary_funs q, q2 = randn(QuaternionF64, 2) - for i in 1:100 + 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)) From 6de82040d3bfe4297cee4b4176919803102cb08d Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 23 Feb 2022 10:32:48 +0100 Subject: [PATCH 30/32] Add function overloads to ReadMe --- README.md | 39 +++++++++++++++++++++++++++++++++++---- 1 file changed, 35 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 3f582d3d..7a9a5cd8 100644 --- a/README.md +++ b/README.md @@ -15,18 +15,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 From 540186de81b7e8c290e9882b82488e2e05207dca Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 23 Feb 2022 15:33:30 +0100 Subject: [PATCH 31/32] Remove useless spaces --- test/test_Quaternion.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_Quaternion.jl b/test/test_Quaternion.jl index 9ce1742d..fba6ff3c 100644 --- a/test/test_Quaternion.jl +++ b/test/test_Quaternion.jl @@ -232,30 +232,30 @@ end @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)) From 596749c22d023ef13c3a03df78f1799682786e8b Mon Sep 17 00:00:00 2001 From: Seth Axen Date: Wed, 23 Feb 2022 21:11:51 +0100 Subject: [PATCH 32/32] Increment version number --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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"