Incorrect rounding when converting irrationals to Float64 in extreme cases

[Joel-Dahne]

Conversion of irrationals to Float64 with a specified rounding mode is currently done by evaluating the irrational with BigFloat at 256 bits of precision and then rounding the result. This can lead to issues if 256 bits is not enough to correctly round the result. There are not examples like this in Base, but one could for example define the (fairly useless) irrational number $1 - \pi^{-300}$, which is clearly less than 1, for which we get (with Julia 1.10.0)

julia> Base.@irrational almost_one (1 - big(π)^-300)

julia> almost_one
almost_one = 1.0...

julia> Float64(almost_one, RoundDown)
1.0

julia> almost_one < 1
false

julia> BigFloat(almost_one)
1.0

julia> setprecision(() -> BigFloat(almost_one), 512)
0.999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999992837905

If one takes the exponent to be for example -100 then it works correctly, since 256 bits is then enough to get correct rounding.

This is clearly an extreme example and I don't think this would be a practical problem for anyone. But it does mean one has to be careful if the goal is to prove properties of an irrational, which would be relevant in for example IntervalArithmetic.jl.

[aravindh-krishnamoorthy]

I must preface this by saying that I'm working with Julia only since a few months. Indeed setting the precision to $512$ solves one of the problems. By using:

@assume_effects :total function (t::Type{T})(x::AbstractIrrational, r::RoundingMode) where T<:Union{Float32,Float64}
-    setprecision(BigFloat, 256) do
+    setprecision(BigFloat, 512) do
        T(BigFloat(x)::BigFloat, r)
    end
end

we get:

julia> using Revise

julia> using Base

julia> Revise.track(Base)

julia> Base.@irrational almost_one (1 - big(π)^-300)

julia> Float64(almost_one, RoundDown)
0.9999999999999999

julia> almost_one < 1
false

However, as we note, almost_one < 1 is still false.

julia> @enter almost_one < 1
In <(x, y) at promotion.jl:465
>465  <( x::Real, y::Real)     = (< )(promote(x,y)...)

About to run: (promote)(almost_one, 1)
1|debug> s
In promote(x, y) at promotion.jl:394
 394  function promote(x, y)
 395      @inline
>396      px, py = _promote(x, y)
 397      not_sametype((x,y), (px,py))
 398      px, py
 399  end

About to run: (Base._promote)(almost_one, 1)
1|debug> n
In promote(x, y) at promotion.jl:394
 394  function promote(x, y)
 395      @inline
 396      px, py = _promote(x, y)
>397      not_sametype((x,y), (px,py))
 398      px, py
 399  end

About to run: (tuple)(almost_one, 1)
1|debug> n
In promote(x, y) at promotion.jl:394
 394  function promote(x, y)
 395      @inline
 396      px, py = _promote(x, y)
 397      not_sametype((x,y), (px,py))
>398      px, py
 399  end

About to run: (tuple)(1.0, 1.0)
1|debug> n
In promote(x, y) at promotion.jl:394
 394  function promote(x, y)
 395      @inline
 396      px, py = _promote(x, y)
 397      not_sametype((x,y), (px,py))
>398      px, py
 399  end

About to run: return (1.0, 1.0)
1|debug> n
In <(x, y) at promotion.jl:465
>465  <( x::Real, y::Real)     = (< )(promote(x,y)...)

About to run: (Core._apply_iterate)(iterate, <, (1.0, 1.0))
1|debug> s
In <(x, y) at float.jl:587
>587  <( x::T, y::T) where {T<:IEEEFloat} = lt_float(x, y)

About to run: $(QuoteNode(Core.Intrinsics.lt_float))
1|debug> n
In <(x, y) at float.jl:587
>587  <( x::T, y::T) where {T<:IEEEFloat} = lt_float(x, y)

About to run: return false
1|debug> n
In <(x, y) at promotion.jl:465
>465  <( x::Real, y::Real)     = (< )(promote(x,y)...)

About to run: return false
1|debug> n
false

But, quirkily, even with setprecision(BigFloat, 256)

julia> almost_one < 1.0
true

julia> @enter almost_one < 1.0
In <(x, y) at irrationals.jl:93
>93  <(x::AbstractIrrational, y::Float64) = Float64(x,RoundUp) <= y

About to run: (Float64)(almost_one, RoundingMode{:Up}())
1|debug> s
In Type(x, r) at irrationals.jl:68
 65  end
 66  Rational{BigInt}(x::AbstractIrrational) = throw(ArgumentError("Cannot convert an AbstractIrrational to a Rational{BigInt}: use rationalize(BigInt, x) instead"))
 67
 68  @assume_effects :total function (t::Type{T})(x::AbstractIrrational, r::RoundingMode) where T<:Union{Float32,Float64}
>69      setprecision(BigFloat, 256) do
 70          T(BigFloat(x)::BigFloat, r)
 71      end
 72  end
 73
 74

About to run: $(QuoteNode(setprecision))
1|debug> n
In Type(x, r) at irrationals.jl:68
 65  end
 66  Rational{BigInt}(x::AbstractIrrational) = throw(ArgumentError("Cannot convert an AbstractIrrational to a Rational{BigInt}: use rationalize(BigInt, x) instead"))
 67
 68  @assume_effects :total function (t::Type{T})(x::AbstractIrrational, r::RoundingMode) where T<:Union{Float32,Float64}
>69      setprecision(BigFloat, 256) do
 70          T(BigFloat(x)::BigFloat, r)
 71      end
 72  end
 73
 74

About to run: return 1.0
1|debug> n
In <(x, y) at irrationals.jl:93
>93  <(x::AbstractIrrational, y::Float64) = Float64(x,RoundUp) <= y

About to run: (<=)(1.0, 1.0)
1|debug> n
In <(x, y) at irrationals.jl:93
>93  <(x::AbstractIrrational, y::Float64) = Float64(x,RoundUp) <= y

About to run: return true
1|debug> n
true

[Joel-Dahne]

Indeed increasing the precision fixes this specific case, but one could of course find other counter examples then! One could solve this by iteratively computing at higher precision, until one gets a result that is strictly between two Float64 values.

Thanks to your comment I also noticed that I gave the inequality in the wrong direction in my original comment. almost_one < 1 fails due to the conversion, not the rounding. As you mention almost_one < 1.0 works correctly. It is the below code that fails due to the rounding.

julia> almost_one > 1.0
true

[aravindh-krishnamoorthy]

Thank you; indeed for any exponent in setprecision(BigFloat, exponent), we can always find a number that will fail almost_one > 1.0 with the current method of handling irrational numbers. Perhaps the right way of handling this is to change the current method of comparison with floats?

On the other hand, I also feel that consistency is required between the results of almost_one < 1 and almost_one < 1.0. In both cases, one of the operands is irrational (a "higher" class), so whether 1 is interpreted as an Int64 or Float64 must not matter for the final result?