Skip to content

Fix NaN arithmetic #81

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

Open
PatrickHaecker wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Fix NaN arithmetic #81

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a6caacc
Keep `infpromote` from converting an infinity to `BigFloat`
Sep 1, 2026
ad04aaf
Propagate `NaN` through the arithmetic
Sep 1, 2026
2d4604f
Return a `Bool` from `signbit(::ComplexInfinity)`
Sep 1, 2026
2b74b78
Compare every ordering against the float result
Sep 1, 2026
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
6 changes: 2 additions & 4 deletions src/Infinities.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,11 +107,9 @@ ComplexInfinity{T}(::Infinity) where T<:Real = ComplexInfinity{T}()
ComplexInfinity(::Infinity) = ComplexInfinity()
ComplexInfinity{T}(x::RealInfinity) where T<:Real = ComplexInfinity{T}(signbit(x))
ComplexInfinity(x::RealInfinity) = ComplexInfinity(signbit(x))
ComplexInfinity{T}(x::ComplexInfinity) where T<:Real = ComplexInfinity(T(signbit(x))) # ambiguity fix
ComplexInfinity{T}(x::ComplexInfinity) where T<:Real = ComplexInfinity(T(x.signbit)) # ambiguity fix

signbit(y::ComplexInfinity{Bool}) = y.signbit
signbit(y::ComplexInfinity{<:Integer}) = !(mod(y.signbit,2) == 0)
signbit(y::ComplexInfinity) = y.signbit
signbit(y::ComplexInfinity) = mod(y.signbit, 2) == 1

convert(::Type{ComplexInfinity{T}}, ::Infinity) where T = ComplexInfinity{T}()
convert(::Type{ComplexInfinity}, ::Infinity) = ComplexInfinity()
Expand Down
29 changes: 20 additions & 9 deletions src/algebra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,9 @@
@inline infpromote(x::RealInfinity, y::Union{Integer, Rational}) = (x, float(y))
@inline infpromote(x::Union{Integer, Rational}, y::RealInfinity) = (float(x), y)
@inline infpromote(x::RealInfinity, ::InfiniteCardinal) = (x, ∞)
# `Base` promotes every `Real` to `BigFloat`, which would convert the infinity away.
@inline infpromote(x::BigFloat, y::Union{Infinity,RealInfinity}) = (x, y)
@inline infpromote(x::Union{Infinity,RealInfinity}, y::BigFloat) = (x, y)


# sign
Expand All @@ -27,7 +30,8 @@
@inline __add(x, y::AllInfinities) = isinf(x) ? _infadd(toinf(x), y) : y
@inline __add(x::Integer, y::InfiniteCardinal) = max(x, y)

@inline _add(x, y) = __add(infpromote(x, y)...)
# A `NaN` argument is returned unchanged, as it is over the floats. Types with no `NaN` fold the test away.
@inline _add(x, y) = isnan(x) ? x : __add(infpromote(x, y)...)

+(x::Number, y::AllInfinities) = _add(x, y)
+(x::AllInfinities, y::Number) = _add(y, x)
Expand All @@ -45,15 +49,20 @@
# multiplication

@inline _sb(x) = signbit(x)
@inline _sb(x::Complex) = angle(x)/π # overloading `signbit` causes type piracy
@inline _sb(x::Complex) = angle(x)/π # overloading `signbit` causes type piracy

@dlfivefifty dlfivefifty Sep 3, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hahaha this is exactly what I just suggested!

@dlfivefifty dlfivefifty Sep 3, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Btw For every θ that is "natural" in base 2 (1, -1, 1/2, 1/4, 3/4, etc.) it turns out angle(exp(im*θ)) == θ. This seems like a a very lucky coincidence that rounding goes in the right direction!

@PatrickHaecker PatrickHaecker Sep 3, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hahaha this is exactly what I just suggested!

Yes, but it only covers part of it, because we still have no defined interval the value is in and we still have a very strange format so store an angle when we have an integer. So this needs more thought.

@PatrickHaecker PatrickHaecker Sep 3, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Btw For every θ that is "natural" in base 2 (1, -1, 1/2, 1/4, 3/4, etc.) it turns out angle(exp(im*θ)) == θ. This seems like a a very lucky coincidence that rounding goes in the right direction!

You are definitely onto something here. While it is true for all your numbers, it interestingly is not true, e.g., for 1/8:

julia> angle(exp(im*π/8))/π
0.12500000000000003

So although 1/8 is exact in floating points, the calculation is not exact. That not only means that our tests are incomplete here (I already added some, but have not yet pushed). It also means that we need to check whether this loss of precision is a problem anywhere. And there are quite some locations where we do the exact comparison.

But I need to do this with a fresh mind tomorrow.

Especially as this lead me to toinf and I currently think we have a bug in toinf:

julia> angle(toinf(complex(0,Inf)))
4.934802200544679

julia> angle(complex(0,Inf))
1.5707963267948966

And I don't understand how we can have such a large bug here without the test suite detecting it. So I probably need to have a fresh look tomorrow to see why this is not a bug and only my current understanding is incomplete.

@PatrickHaecker PatrickHaecker Sep 4, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It was indeed a bug and there is an explanation why the tests haven't caught it. Fix, details and off-axis tests in #85. So for all I know this conversation is resolved.

@dlfivefifty dlfivefifty Sep 4, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The really important case is that for r > 0 where r is any real we always have _sb(r+0im) = zero(r), _sb(-r±0im) = one(r), _sb(r*im) = one(r)/2, _sb(-r*im) = -one(r)/2

so we want this to be correct for BigFloat and other types too. Actually, I'm not sure what the types should be (I guess for <: Real we want them to be Bool)

@PatrickHaecker PatrickHaecker Sep 4, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I assume you mean == as comparison. What do you mean with "types"? The storage type or the return type?

I am currently trying a combination of UInt64 as a storage type and using sign with a Float64 as a possible user interface in addition to other methods for Base functions.
This seems to work well up to now (unique values, some precision everywhere, wrapping for free, usage without rounding possible), but I can't yet see whether this will work in the whole package, which is the relevant question in the end.

@dlfivefifty dlfivefifty Sep 4, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a bit of a fringe question so I don't think any choice makes a big practical choice. But what I mean is if _sb(x::Real) = signbit(x) then the return type is a Bool, but for _sb(::Complex{<:Union{Int,Float64}) it should probably be a Float64. But then what about _sb(::Complex32) and _sb(::Complex{BigFloat})? Probably Float32 and BigFloat. (which I believe your implementation does)

@dlfivefifty dlfivefifty Sep 4, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess the concern is: do we actually know that rounding works correctly for all BigFloat? Eg, we have

julia> setprecision(10) do
       angle(one(BigFloat)im)
       end
1.5703

julia> setprecision(20) do
       angle(one(BigFloat)im)
       end
1.570797

julia> setprecision(40) do
       angle(one(BigFloat)im)
       end
1.5707963267942

There is no obvious guarantee that when we divide this by π we get exactly 1/2.

Except.....in tests we do up to precision 10_000:

julia> for p in 1:10_000
       setprecision(p) do
       @test angle(one(BigFloat)im)/π == 1/2
       end
       end

Which is to say: I think by luck we are fine.

@PatrickHaecker PatrickHaecker Sep 5, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I assume you talking about #81 and not about the possible future solution.

I think in general we can't be exact with floats. The user can always have any arbitrary transformations between two points in time and when they compare then according to what should be identical with two mathematical reals it can be different in our implementation.

That being said, I think in this case it is probably more than luck. I am speculating, because I haven't followed into MPFR's implementation, but your example boils down to:

atan(one(BigFloat), zero(BigFloat))/π == 1/2

which is

atan(one(BigFloat), zero(BigFloat)) == π/2

and I strongly assume that this special cases such that atan will return the best rounding of π/2 according to the current precision if the first argument is positive and the second is 0. This triggers reliably because 1 and 0 are exact. π and π/2 will then have the same number in the significant and only differ in the mantissa by exactly 1 leading to 0.5 exactly in the division.

So while it would be possible to implement that in a different way, I think an implementation which cares about precision will always return the exact value here, despite being a floating point operation. That's just a special case which is exact when using binary floating points. It wouldn't for some other base in general, but it is with base 2.

@inline _sb(x::ComplexInfinity) = x.signbit # the whole angle, not just its sign

@inline __mul(x, y::AllInfinities) = RealInfinity(_sb(x) ⊻ _sb(y))
@inline __mul(x, y::ComplexInfinity) = ComplexInfinity(_sb(x) + _sb(y))
@inline __mul(x, y::ComplexInfinity{Bool}) = ComplexInfinity(_sb(x) ⊻ _sb(y))
@inline __mul(x::Complex, y::ComplexInfinity{Bool}) = ComplexInfinity(_sb(x) + _sb(y))
@inline __mul(x::Integer, y::InfiniteCardinal) = x > 0 ? y : throw(ArgumentError("Cannot multiply $x * $y"))

@inline _mul(x, y) = iszero(x) ? throw(ArgumentError("Cannot multiply $x * $y")) : __mul(infpromote(x, y)...)
@inline function _mul(x, y)
isnan(x) && return x
iszero(x) && throw(ArgumentError("Cannot multiply $x * $y"))
__mul(infpromote(x, y)...)
end

*(x::Number, y::AllInfinities) = _mul(x, y)
*(x::AllInfinities, y::Number) = _mul(y, x)
Expand All @@ -68,6 +77,7 @@

# mod
@inline function _mod(x::Real, y::IntegerInfinities)
isnan(x) && return x
signbit(x) == signbit(y) || throw(ArgumentError("mod($x,$y) is unbounded"))
x
end
Expand All @@ -76,14 +86,14 @@ mod(::IntegerInfinities, ::Real) = NotANumber()
mod(::IntegerInfinities, ::IntegerInfinities) = NotANumber()

# fld, cld, div
_divinf(T) = zero(T)
_fldinf(x) = signbit(x) ? -one(x) : zero(x)
_cldinf(x) = signbit(x) ? zero(x) : one(x)
div(::T, ::IntegerInfinities) where T <: Real = _divinf(T)
_divinf(x) = isnan(x) ? x : zero(x)
_fldinf(x) = isnan(x) ? x : signbit(x) ? -one(x) : zero(x)
_cldinf(x) = isnan(x) ? x : signbit(x) ? zero(x) : one(x)
div(x::Real, ::IntegerInfinities) = _divinf(x)
fld(x::Real, ::IntegerInfinities) = _fldinf(x)
cld(x::Real, ::IntegerInfinities) = _cldinf(x)

_inffcd(x, y) = signbit(y) ? -x : x
_inffcd(x, y) = isnan(y) ? y : signbit(y) ? -x : x
for OP in (:fld,:cld,:div)
@eval begin
$OP(x::IntegerInfinities, y::Real) = _inffcd(x, y)
Expand All @@ -94,8 +104,9 @@ end
# power
# Although the base implementation can cover these cases, it can change overtime and yield inconsistent results.
# ref: https://github.com/JuliaMath/Infinities.jl/actions/runs/19993302836/
_infpow(::PositiveInfinity, p) = ifelse(iszero(p), one(p), ifelse(p > 0, +∞, +zero(p)))
_infpow(::PositiveInfinity, p) = isnan(p) ? p : ifelse(iszero(p), one(p), ifelse(p > 0, +∞, +zero(p)))
function _infpow(x::NegativeInfinity, p)
isnan(p) && return p
!isinteger(p) && throw(Base.Math.throw_exp_domainerror(x))
iszero(p) && return one(p)
isodd(p) && return ifelse(p > 0, -∞, -zero(p))
Expand Down
2 changes: 1 addition & 1 deletion src/ambiguities.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ for Typ in (Rational, )
for op in (:fld, :cld, :div)
@eval $op(x::InfiniteCardinal, y::$Typ) = _inffcd(x, y)
end
@eval div(::T, ::IntegerInfinities) where T <: $Typ = _divinf(T)
@eval div(x::$Typ, ::IntegerInfinities) = _divinf(x)
@eval fld(x::$Typ, ::IntegerInfinities) = _fldinf(x)
@eval cld(x::$Typ, ::IntegerInfinities) = _cldinf(x)
end
33 changes: 33 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -326,6 +326,9 @@ Base.iterate(s::CharString, i::Integer=1) = i ≤ length(s.chars) ? (s.chars[i],

@test signbit(ComplexInfinity(3))
@test !signbit(ComplexInfinity(100))
# `signbit` returns a `Bool` for every angle, as it does over the reals
@test signbit(ComplexInfinity(1.0)) === signbit(-ComplexInfinity()) === true
@test signbit(ComplexInfinity(0.5)) === signbit(ComplexInfinity()) === false
end

@testset "Set" begin
Expand Down Expand Up @@ -404,6 +407,8 @@ Base.iterate(s::CharString, i::Integer=1) = i ≤ length(s.chars) ? (s.chars[i],
@test T(-Inf) == inf == T(-Inf)
@test T(Inf) ≠ inf
end
@test T(2) + ∞ ≡ ∞ + T(2) ≡ ∞
@test T(2) * +∞ ≡ (+∞)^T(2) ≡ +∞
end
end

Expand Down Expand Up @@ -448,6 +453,20 @@ Base.iterate(s::CharString, i::Integer=1) = i ≤ length(s.chars) ? (s.chars[i],
@test sorted[1] === -∞ && sorted[2] === 1.0 && sorted[3] === ∞ && isnan(sorted[4])
end

@testset "NaN arithmetic" begin
for nan in (NaN, NaN32, NaN16, big(NaN)),
inf in (∞, +∞, -∞, ℵ₀, ComplexInfinity(), -ComplexInfinity())

for op in (+, -, *, div, fld, cld)
@test isnan(op(nan, inf)) && isnan(op(inf, nan))
end
@test isnan(mod(nan, inf)) # the other direction is `NotANumber` for every argument
end
for nan in (NaN, NaN32, NaN16, big(NaN)), inf in (+∞, -∞)
@test isnan(inf^nan)
end
end

@testset "ordinary values" begin
for inf in (∞, +∞, ℵ₀)
@test 1.0 < inf && !(inf < 1.0) && 1.0 ≤ inf && inf ≥ 1.0
Expand All @@ -460,6 +479,20 @@ Base.iterate(s::CharString, i::Integer=1) = i ≤ length(s.chars) ? (s.chars[i],
@test max(-∞, ∞) === ∞ && min(-∞, ∞) === -∞
end

@testset "against the floats" begin
values = (0, 1, -2, 1.5, -1.5, 0.0, -0.0, NaN, NaN32, Inf, -Inf,
prevfloat(Inf), nextfloat(-Inf), nextfloat(0.0))
for x in values, (inf, flt) in ((∞, Inf), (+∞, Inf), (-∞, -Inf), (ℵ₀, Inf))
for op in (<, ≤, >, ≥, ==, isless, isequal)
@test op(x, inf) == op(x, flt)
@test op(inf, x) == op(flt, x)
end
for op in (max, min)
@test isequal(op(x, inf), op(x, flt)) && isequal(op(inf, x), op(flt, x))
end
end
end

@testset "parsing" begin
@test tryparse(NegativeInfinity, "-∞") == NegativeInfinity()
@test tryparse(NegativeInfinity, " - ∞ ") == NegativeInfinity()
Expand Down
Loading