Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
sethaxen merged 24 commits into from
Jan 9, 2023
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
89463d7
Use hypot in abs and abs_imag
sethaxen Dec 10, 2022
9901662
Use approximate check
sethaxen Dec 10, 2022
65ed170
Add tests against abs/exp under/overflow
sethaxen Dec 10, 2022
1a25d87
Increment patch number
sethaxen Dec 10, 2022
ab0a01f
Ensure tests pass on 1.0
sethaxen Dec 10, 2022
4bfe65a
Merge branch 'main' into hypot
sethaxen Jan 7, 2023
8675e50
Backport _hypot code for older Julia versions
sethaxen Jan 7, 2023
b154432
Mark test as non broken
sethaxen Jan 7, 2023
f0e840f
Implement abs without hypot
sethaxen Jan 8, 2023
a2e138b
Implement abs_imag efficiently
sethaxen Jan 8, 2023
174c934
Handle over/underflow for inv and div
sethaxen Jan 8, 2023
690a08b
Make inv implementation more like Complex's
sethaxen Jan 8, 2023
f158626
Test abs_imag for all Julia versions
sethaxen Jan 8, 2023
4f38866
Add tests for inv under/overflow
sethaxen Jan 8, 2023
42219ea
Test more abs/abs_imag cases
sethaxen Jan 8, 2023
7c7f851
Simplify tests with isequal
sethaxen Jan 8, 2023
51296bc
Use isequal
sethaxen Jan 8, 2023
32ce61b
Add isequal
sethaxen Jan 8, 2023
a45cdf9
Add under/overflow tests for div
sethaxen Jan 8, 2023
083c86a
Add missing abs_imag tests
sethaxen Jan 8, 2023
481f73b
Apply suggestions from code review
sethaxen Jan 9, 2023
1b8a88f
Update test/Quaternion.jl
sethaxen Jan 9, 2023
bfa5e38
Update test/Quaternion.jl
sethaxen Jan 9, 2023
347c880
Increment patch number
sethaxen Jan 9, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
58 changes: 54 additions & 4 deletions src/Quaternion.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
74 changes: 73 additions & 1 deletion test/Quaternion.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)) ===
Expand Down Expand Up @@ -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
Expand All @@ -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))
Expand Down Expand Up @@ -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
Expand Down