Make abs, abs_imag, inv, and / resistant to under/overflow (#122)

* Use hypot in abs and abs_imag

* Use approximate check

* Add tests against abs/exp under/overflow

* Increment patch number

* Ensure tests pass on 1.0

* Backport _hypot code for older Julia versions

* Mark test as non broken

* Implement abs without hypot

* Implement abs_imag efficiently

* Handle over/underflow for inv and div

* Make inv implementation more like Complex's

* Test abs_imag for all Julia versions

* Add tests for inv under/overflow

* Test more abs/abs_imag cases

* Simplify tests with isequal

* Use isequal

* Add isequal

* Add under/overflow tests for div

* Add missing abs_imag tests

* Apply suggestions from code review

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

* Update test/Quaternion.jl

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

* Update test/Quaternion.jl

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

* Increment patch number

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

Author: sethaxen

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"

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)

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