diff --git a/Project.toml b/Project.toml index e325906..860b422 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Quaternions" uuid = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0" -version = "0.7.3" +version = "0.7.4" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/Quaternion.jl b/src/Quaternion.jl index 3d0f729..2fb3d12 100644 --- a/src/Quaternion.jl +++ b/src/Quaternion.jl @@ -136,11 +136,42 @@ Quaternion{Int64}(1, -2, -3, -4) ``` """ Base.conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3) -Base.abs(q::Quaternion) = sqrt(abs2(q)) +function Base.abs(q::Quaternion) + a = max(abs(q.s), abs(q.v1), abs(q.v2), abs(q.v3)) + if isnan(a) && isinf(q) + return typeof(a)(Inf) + elseif iszero(a) || isinf(a) + return a + else + return sqrt(abs2(q / a)) * a + end +end Base.float(q::Quaternion{T}) where T = convert(Quaternion{float(T)}, q) -abs_imag(q::Quaternion) = sqrt(q.v2 * q.v2 + (q.v1 * q.v1 + q.v3 * q.v3)) # ordered to match abs2 +function abs_imag(q::Quaternion) + a = max(abs(q.v1), abs(q.v2), abs(q.v3)) + if isnan(a) && (isinf(q.v1) | isinf(q.v2) | isinf(q.v3)) + return oftype(a, Inf) + elseif iszero(a) || isinf(a) + return a + else + return sqrt((q.v1 / a)^2 + (q.v2 / a)^2 + (q.v3 / a)^2) * a + end +end Base.abs2(q::Quaternion) = RealDot.realdot(q,q) -Base.inv(q::Quaternion) = conj(q) / abs2(q) +function Base.inv(q::Quaternion) + if isinf(q) + return quat( + copysign(zero(q.s), q.s), + flipsign(-zero(q.v1), q.v1), + flipsign(-zero(q.v2), q.v2), + flipsign(-zero(q.v3), q.v3), + ) + end + a = max(abs(q.s), abs(q.v1), abs(q.v2), abs(q.v3)) + p = q / a + iq = conj(p) / (a * abs2(p)) + return iq +end Base.isreal(q::Quaternion) = iszero(q.v1) & iszero(q.v2) & iszero(q.v3) Base.isfinite(q::Quaternion) = isfinite(q.s) & isfinite(q.v1) & isfinite(q.v2) & isfinite(q.v3) @@ -182,9 +213,28 @@ function Base.:*(q::Quaternion, w::Quaternion) return Quaternion(s, v1, v2, v3) end -Base.:/(q::Quaternion, w::Quaternion) = q * inv(w) +function Base.:/(q::Quaternion{T}, w::Quaternion{T}) where T + # handle over/underflow while matching the behavior of /(a::Complex, b::Complex) + a = max(abs(w.s), abs(w.v1), abs(w.v2), abs(w.v3)) + if isinf(w) + if isfinite(q) + return quat( + zero(T)*sign(q.s)*sign(w.s), + -zero(T)*sign(q.v1)*sign(w.v1), + -zero(T)*sign(q.v2)*sign(w.v2), + -zero(T)*sign(q.v3)*sign(w.v3), + ) + end + return quat(T(NaN), T(NaN), T(NaN), T(NaN)) + end + p = w / a + return (q * conj(p)) / RealDot.realdot(w, p) +end Base.:(==)(q::Quaternion, w::Quaternion) = (q.s == w.s) & (q.v1 == w.v1) & (q.v2 == w.v2) & (q.v3 == w.v3) +function Base.isequal(q::Quaternion, w::Quaternion) + isequal(q.s, w.s) & isequal(q.v1, w.v1) & isequal(q.v2, w.v2) & isequal(q.v3, w.v3) +end """ extend_analytic(f, q::Quaternion) diff --git a/test/Quaternion.jl b/test/Quaternion.jl index 6d9d335..b02c0ca 100644 --- a/test/Quaternion.jl +++ b/test/Quaternion.jl @@ -59,6 +59,13 @@ end @test Quaternion(1, 2, 3, 4) != Quaternion(1, 2, 3, 5) end + @testset "isequal" begin + @test isequal(Quaternion(1, 2, 3, 4), Quaternion(1.0, 2.0, 3.0, 4.0)) + @test !isequal(Quaternion(1, 2, 3, 4), Quaternion(5, 2, 3, 4)) + @test isequal(Quaternion(NaN, -0.0, Inf, -Inf), Quaternion(NaN, -0.0, Inf, -Inf)) + @test !isequal(Quaternion(NaN, 0.0, Inf, -Inf), Quaternion(NaN, -0.0, Inf, -Inf)) + end + @testset "convert" begin @test convert(Quaternion{Float64}, 1) === Quaternion(1.0) @test convert(Quaternion{Float64}, Quaternion(1, 2, 3, 4)) === @@ -153,7 +160,31 @@ end @test conj(conj(q)) === q @test conj(conj(qnorm)) === qnorm @test float(Quaternion(1, 2, 3, 4)) === float(Quaternion(1.0, 2.0, 3.0, 4.0)) - @test Quaternions.abs_imag(q) == abs(Quaternion(0, q.v1, q.v2, q.v3)) + @test Quaternions.abs_imag(q) ≈ abs(Quaternion(0, q.v1, q.v2, q.v3)) + end + + @testset "abs/abs_imag don't over/underflow" begin + for x in [1e-300, 1e300, -1e-300, -1e300] + @test abs(quat(x, 0, 0, 0)) == abs(x) + @test abs(quat(0, x, 0, 0)) == abs(x) + @test abs(quat(0, 0, x, 0)) == abs(x) + @test abs(quat(0, 0, 0, x)) == abs(x) + @test Quaternions.abs_imag(quat(0, x, 0, 0)) == abs(x) + @test Quaternions.abs_imag(quat(0, 0, x, 0)) == abs(x) + @test Quaternions.abs_imag(quat(0, 0, 0, x)) == abs(x) + end + @test isnan(abs(quat(NaN, NaN, NaN, NaN))) + @test abs(quat(NaN, Inf, NaN, NaN)) == Inf + @test abs(quat(-Inf, NaN, NaN, NaN)) == Inf + @test abs(quat(0.0)) == 0.0 + @test abs(quat(Inf)) == Inf + @test abs(quat(1, -Inf, 2, 3)) == Inf + @test isnan(Quaternions.abs_imag(quat(0, NaN, NaN, NaN))) + @test Quaternions.abs_imag(quat(0, Inf, NaN, NaN)) == Inf + @test Quaternions.abs_imag(quat(0, NaN, -Inf, NaN)) == Inf + @test Quaternions.abs_imag(quat(0.0)) == 0.0 + @test Quaternions.abs_imag(quat(0.0, 0.0, Inf, 0.0)) == Inf + @test Quaternions.abs_imag(quat(0, 1, -Inf, 2)) == Inf end @testset "algebraic properties" begin @@ -171,6 +202,21 @@ end end end + @testset "inv does not under/overflow" begin + x = 1e-300 + y = inv(x) + @test isequal(inv(quat(x, 0.0, 0.0, 0.0)), quat(y, -0.0, -0.0, -0.0)) + @test isequal(inv(quat(0.0, x, 0.0, 0.0)), quat(0.0, -y, -0.0, -0.0)) + @test isequal(inv(quat(0.0, 0.0, x, 0.0)), quat(0.0, -0.0, -y, -0.0)) + @test isequal(inv(quat(0.0, 0.0, 0.0, x)), quat(0.0, -0.0, -0.0, -y)) + @test isequal(inv(quat(y, 0.0, 0.0, 0.0)), quat(x, -0.0, -0.0, -0.0)) + @test isequal(inv(quat(0.0, y, 0.0, 0.0)), quat(0.0, -x, -0.0, -0.0)) + @test isequal(inv(quat(0.0, 0.0, y, 0.0)), quat(0.0, -0.0, -x, -0.0)) + @test isequal(inv(quat(0.0, 0.0, 0.0, y)), quat(0.0, -0.0, -0.0, -x)) + @test isequal(inv(quat(-Inf, 1, -2, 3)), quat(-0.0, -0.0, 0.0, -0.0)) + @test isequal(inv(quat(1, -2, Inf, 3)), quat(0.0, 0.0, -0.0, -0.0)) + end + @testset "isreal" begin @test isreal(Quaternion(1, 0, 0, 0)) @test !isreal(Quaternion(2, 1, 0, 0)) @@ -275,6 +321,32 @@ end @test q2 \ q ≈ inv(q2) * q @test q / x ≈ x \ q ≈ inv(x) * q end + @testset "no overflow/underflow" begin + @testset for x in [1e-300, 1e300, -1e-300, -1e300] + @test quat(x) / quat(x) == quat(1) + @test quat(x) / quat(0, x, 0, 0) == quat(0, -1, 0, 0) + @test quat(x) / quat(0, 0, x, 0) == quat(0, 0, -1, 0) + @test quat(x) / quat(0, 0, 0, x) == quat(0, 0, 0, -1) + @test quat(0, x, 0, 0) / quat(x, 0, 0, 0) == quat(0, 1, 0, 0) + @test quat(0, x, 0, 0) / quat(0, x, 0, 0) == quat(1, 0, 0, 0) + @test quat(0, x, 0, 0) / quat(0, 0, x, 0) == quat(0, 0, 0, -1) + @test quat(0, x, 0, 0) / quat(0, 0, 0, x) == quat(0, 0, 1, 0) + end + @testset for T in [Float32, Float64] + o = one(T) + z = zero(T) + inf = T(Inf) + nan = T(NaN) + @testset for s in [1, -1], t in [1, -1] + @test isequal(quat(o) / quat(s*inf), quat(s*z, -z, -z, -z)) + @test isequal(quat(o) / quat(s*inf, t*o, z, t*z), quat(s*z, -t*z, -z, -t*z)) + @test isequal(quat(o) / quat(s*inf, t*nan, t*z, z), quat(s*z, nan, -t*z, -z)) + @test isequal(quat(o) / quat(s*inf, t*inf, t*z, z), quat(s*z, -t*z, -t*z, -z)) + end + @test isequal(quat(inf) / quat(inf, 1, 2, 3), quat(nan, nan, nan, nan)) + @test isequal(quat(inf) / quat(inf, 1, 2, -inf), quat(nan, nan, nan, nan)) + end + end end @testset "^" begin