Skip to content
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.

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

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

Merged
merged 5 commits into from
Jan 24, 2023
Merged
Changes from 1 commit
Commits
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
Prev Previous commit
Next Next commit
Improve numeric stability
sethaxen committed Jan 23, 2023
commit 71c77fea57ec4d0137cdbe94624dbade0a3f4c34
50 changes: 46 additions & 4 deletions src/octonion.jl
Original file line number Diff line number Diff line change
@@ -35,11 +35,28 @@ imag_part(o::Octonion) = (o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7)
(*)(x::Real, o::Octonion) = o * x

conj(o::Octonion) = Octonion(o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7)
abs(o::Octonion) = sqrt(abs2(o))
abs(o::Octonion) = hypot(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7)
float(q::Octonion{T}) where T = convert(Octonion{float(T)}, q)
abs_imag(o::Octonion) = sqrt((o.v4 * o.v4 + (o.v2 * o.v2 + o.v6 * o.v6)) + ((o.v1 * o.v1 + o.v5 * o.v5) + (o.v3 * o.v3 + o.v7 * o.v7))) # ordered to match abs2
inv(o::Octonion) = conj(o) / abs2(o)
abs_imag(o::Octonion) = hypot(o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7)
abs2(o::Octonion) = RealDot.realdot(o, o)
function inv(o::Octonion)
if isinf(o)
return octo(
copysign(zero(o.s), o.s),
flipsign(-zero(o.v1), o.v1),
flipsign(-zero(o.v2), o.v2),
flipsign(-zero(o.v3), o.v3),
flipsign(-zero(o.v4), o.v4),
flipsign(-zero(o.v5), o.v5),
flipsign(-zero(o.v6), o.v6),
flipsign(-zero(o.v7), o.v7),
)
end
a = max(abs(o.s), maximum(abs, imag_part(o)))
p = o / a
io = conj(p) / (a * abs2(p))
return io
end

isreal(o::Octonion) = iszero(o.v1) & iszero(o.v2) & iszero(o.v3) & iszero(o.v4) & iszero(o.v5) & iszero(o.v6) & iszero(o.v7)
isfinite(o::Octonion) = isfinite(real(o)) & isfinite(o.v1) & isfinite(o.v2) & isfinite(o.v3) & isfinite(o.v4) & isfinite(o.v5) & isfinite(o.v6) & isfinite(o.v7)
@@ -79,10 +96,35 @@ function (*)(o::Octonion, w::Octonion)
return Octonion(s, v1, v2, v3, v4, v5, v6, v7)
end

(/)(o::Octonion, w::Octonion) = o * inv(w)
function Base.:/(o::Octonion{T}, w::Octonion{T}) where T
# handle over/underflow while matching the behavior of /(a::Complex, b::Complex)
a = max(abs(real(w)), maximum(abs, imag_part(w)))
if isinf(w)
if isfinite(o)
return octo(
zero(T)*sign(o.s)*sign(w.s),
-zero(T)*sign(o.v1)*sign(w.v1),
-zero(T)*sign(o.v2)*sign(w.v2),
-zero(T)*sign(o.v3)*sign(w.v3),
-zero(T)*sign(o.v4)*sign(w.v4),
-zero(T)*sign(o.v5)*sign(w.v5),
-zero(T)*sign(o.v6)*sign(w.v6),
-zero(T)*sign(o.v7)*sign(w.v7),
)
end
return octo(T(NaN), T(NaN), T(NaN), T(NaN), T(NaN), T(NaN), T(NaN), T(NaN))
end
p = w / a
return (o * conj(p)) / RealDot.realdot(w, p)
end

(==)(q::Octonion, w::Octonion) = (q.s == w.s) & (q.v1 == w.v1) & (q.v2 == w.v2) & (q.v3 == w.v3) &
(q.v4 == w.v4) & (q.v5 == w.v5) & (q.v6 == w.v6) & (q.v7 == w.v7)
function Base.isequal(q::Octonion, w::Octonion)
return (isequal(q.s, w.s) & isequal(q.v1, w.v1) & isequal(q.v2, w.v2) &
isequal(q.v3, w.v3) & isequal(q.v4, w.v4) & isequal(q.v5, w.v5) &
isequal(q.v6, w.v6) & isequal(q.v7, w.v7))
end

function exp(o::Octonion)
s = o.s
88 changes: 87 additions & 1 deletion test/octonion.jl
Original file line number Diff line number Diff line change
@@ -57,6 +57,21 @@ end
@test Octonion(1.0, 2, 3, 4, 5, 6, 7, 8) != Octonion(1, 2, 3, 4, 1, 2, 3, 4)
end

@testset "isequal" begin
@test isequal(Octonion(1:8...), Octonion(1.0:8.0...))
x = ntuple(identity, 8)
o = Octonion(x...)
for i in 1:8
x2 = Base.setindex(x, 0, i)
o2 = Octonion(x2...)
@test !isequal(o, o2)
end
o = Octonion(NaN, -0.0, Inf, -Inf, 0.0, NaN, Inf, -Inf)
@test isequal(o, o)
@test !isequal(o, Octonion(NaN, 0.0, Inf, -Inf, 0.0, NaN, Inf, -Inf))
@test !isequal(o, Octonion(NaN, -0.0, Inf, -Inf, -0.0, NaN, Inf, -Inf))
end

@testset "convert" begin
@test convert(Octonion{Float64}, 1) === Octonion(1.0)
@test convert(Octonion{Float64}, Octonion(1, 2, 3, 4, 5, 6, 7, 8)) ===
@@ -142,10 +157,33 @@ end
@test conj(conj(q)) === q
@test conj(conj(qnorm)) === qnorm
@test float(Octonion(1:8...)) === Octonion(1.0:8.0...)
@test Octonions.abs_imag(q) ==
@test Octonions.abs_imag(q)
abs(Octonion(0, q.v1, q.v2, q.v3, q.v4, q.v5, q.v6, q.v7))
end

@testset "abs/abs_imag don't over/underflow" begin
for x in [1e-300, 1e300, -1e-300, -1e300]
for i in 1:8
z = Base.setindex(ntuple(zero, 8), x, i)
o = octo(z...)
@test abs(o) == abs(x)
i == 1 || @test Octonions.abs_imag(o) == abs(x)
end
end
@test isnan(abs(octo(fill(NaN, 8)...)))
@test abs(octo(NaN, Inf, fill(NaN, 6)...)) == Inf
@test abs(octo(-Inf, fill(NaN, 7)...)) == Inf
@test abs(octo(0.0)) == 0.0
@test abs(octo(Inf)) == Inf
@test abs(octo(1, -Inf, 3:8...)) == Inf
@test isnan(Octonions.abs_imag(octo(0, fill(NaN, 7)...)))
@test Octonions.abs_imag(octo(0, Inf, fill(NaN, 6)...)) == Inf
@test Octonions.abs_imag(octo(0, NaN, -Inf, fill(NaN, 5)...)) == Inf
@test Octonions.abs_imag(octo(0.0)) == 0.0
@test Octonions.abs_imag(octo(0.0, 0.0, Inf, fill(0, 5)...)) == Inf
@test Octonions.abs_imag(octo(0, 1, -Inf, 3:7...)) == Inf
end

@testset "algebraic properties" begin
for _ in 1:100, T in (Float32, Float64, Int32, Int64)
if T <: Integer
@@ -169,6 +207,25 @@ end
end
end

@testset "inv does not under/overflow" begin
x = 1e-300
y = inv(x)
for i in 1:8
z = zeros(8)
z[i] = x
z2 = vcat(0.0, fill(-0.0, 7))
z2[i] = i == 1 ? y : -y
o = octo(z...)
@test inv(octo(z...)) == octo(z2...)

z[i] = y
z2[i] = i == 1 ? x : -x
@test inv(octo(z...)) == octo(z2...)
end
@test isequal(inv(octo(-Inf, 1, -2, 3, -4, 5, -6, 7)), octo(-0.0, -0.0, 0.0, -0.0, 0.0, -0.0, 0.0, -0.0))
@test isequal(inv(octo(1, -2, Inf, 3, -4, 5, -6, 7)), octo(0.0, 0.0, -0.0, -0.0, 0.0, -0.0, 0.0, -0.0))
end

@testset "isreal" begin
@test isreal(octo(1))
@test !isreal(octo(1, 1, 0, 0, 0, 0, 0, 0))
@@ -361,6 +418,35 @@ end
@test o2 \ o ≈ inv(o2) * o
@test o / x ≈ x \ o ≈ inv(x) * o
end

@testset "no overflow/underflow" begin
@testset for x in [1e-300, 1e300, -1e-300, -1e300]
@test octo(x) / octo(x) == octo(1)
@testset for i in 2:8
z = Base.setindex(ntuple(zero, 8), x, i)
z2 = Base.setindex(ntuple(zero, 7), -1, i - 1)
@test octo(x) / octo(z...) == octo(0, z2...)
end
@test octo(0, x, zeros(6)...) / octo(x, 0, 0, 0, zeros(4)...) == octo(0, 1, 0, 0, zeros(4)...)
@test octo(0, x, zeros(6)...) / octo(0, x, 0, 0, zeros(4)...) == octo(1, 0, 0, 0, zeros(4)...)
@test octo(0, x, zeros(6)...) / octo(0, 0, x, 0, zeros(4)...) == octo(0, 0, 0, -1, zeros(4)...)
@test octo(0, x, zeros(6)...) / octo(0, 0, 0, x, zeros(4)...) == octo(0, 0, 1, 0, zeros(4)...)
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(octo(o) / octo(s*inf), octo(s*z, fill(-z, 7)...))
@test isequal(octo(o) / octo(s*inf, t*o, z, t*z, z, t*z, z, t*z), octo(s*z, -t*z, -z, -t*z, -z, -t*z, -z, -t*z))
@test isequal(octo(o) / octo(s*inf, t*nan, t*z, z, t*z, z, t*z, z), octo(s*z, nan, -t*z, -z, -t*z, -z, -t*z, -z))
@test isequal(octo(o) / octo(s*inf, t*inf, t*z, z, t*z, z, t*z, z), octo(s*z, -t*z, -t*z, -z, -t*z, -z, -t*z, -z))
end
@test isequal(octo(inf) / octo(inf, 1:7...), octo(fill(nan, 8)...))
@test isequal(octo(inf) / octo(inf, 1, 2, -inf, 4:7...), octo(fill(nan, 8)...))
end
end
end

@testset "^" begin