Base: correctly rounded floats constructed from rationals

[nsajko]

Constructing a floating-point number from a Rational should now be correctly rounded.

Implementation approach:

1. Convert the (numerator, denominator) pair to a (sign bit, integral significand, exponent) triplet using integer arithmetic. The integer type in question must be wide enough.

2. Convert the above triplet into an instance of the chosen FP type. There is special support for IEEE 754 floating-point and for BigFloat, otherwise a fallback using ldexp is used.

As a bonus, constructing a BigFloat from a Rational should now be thread-safe when the rounding mode and precision are provided to the constructor, because there is no access to the global precision or rounding mode settings.

Updates #45213

Updates #50940

Updates #52507

Fixes #52394

Closes #52395

Fixes #52859

[nsajko]

Regarding performance:

1. the ldexp calls are not essential in the cases of the IEEE floating-point formats, ldexp could be replaced with a bit of bit twiddling for those common cases. I can do that now or in a future PR.
2. BigInt should have it's own implementation for performance reasons, although I'm not sure whether it'd be better to put that in Base or in MutableArithmetics.

[nsajko]

The second commit is the mentioned optimization of using bit twiddling instead of ldexp for IEEE floats. If you're fine with having it in this PR, I can squash them, otherwise I'll move it to a subsequent PR.

Regarding the 92 character line length limit, does it apply to comments equally? I usually try to have comments be narrower than code, but if you want I'll rewrap them to 92 columns.

[nsajko]

performance

Profiling reveals that a considerable amount of time is spent in the line

y = Rational(promote(numerator(x), denominator(x))...)

The Rational call spends almost all of its time in divgcd, however this is unnecessary and could be avoided if we called unsafe_rational instead. I think this will be the single biggest performance win, but there are other possible optimization opportunities that should be looked into, such as optimizing rational_to_float_components_impl to do just one division and then round manually (like it was before), and the optimization of replacing ldexp with bit-level operations for IEEE floats.

[Seelengrab]

Which cases spend their time in promotion in particular? IIRC this should be a noop if there's nothing to pomote/convert :eyes:

[nsajko]

Which cases spend their time in promotion in particular? IIRC this should be a noop if there's nothing to pomote/convert

The problem is not promote, it's the Rational call in the same line. I already replaced it with unsafe_rational locally, to avoid the costly but unnecessary divgcd.

[nsajko]

JET shows that this line causes run time dispatch, I think because the type of rm is unknown:

https://github.com/JuliaLang/julia/blob/dd48fddc2b9a25defac1043d64fb088b9d3e8a87/base/mpfr.jl#L329

││││┌ @ mpfr.jl:329 Base.MPFR.rational_to_float_components(%18, %27, prec, Base.MPFR.BigInt, %102)
│││││ runtime dispatch detected: Base.MPFR.rational_to_float_components(%18::Int64, %27::Int64, prec::Int64, Base.MPFR.BigInt, %102::RoundingMode)::Base.FloatComponentsResult{BigInt, Int64}

There doesn't seem to be a way around this without introducing a new type for roundings and using it in rational_to_float_components_impl instead of RoundingMode. This doesn't seem too difficult but it may be overkill?

[Seelengrab]

This is just for BigFloat, right? I don't think this dispatch will be particularly problematic then - any computation involving BigFloat of sizes/precisions where it's worth using BigFloat instead of something like Float128 will probably be bottlenecked by MPFR/GMP operations under the hood instead of this dispatch (not to mention the other allocating BigFloat/BigInt operations in here).

[nsajko]

Is there still something to do here from my side?

[nsajko]

Another small optimization would be using ndigits0zpb instead of ndigits. Not sure if I should do this now or wait until this PR is merged. EDIT: probably top_set_bit would be even better than ndigits.

[Seelengrab]

The failures seem to be related to package extensions, not this PR:

https://buildkite.com/julialang/julia-master/builds/24099#01883392-fcbf-4f7f-b20c-6871fbe4a277/748-1732

[nsajko]

This seems like a pretty complex (and slow) implementation of something that ought to be pretty simple.

"Make everything as simple as possible, but not simpler." Note the "as possible" part.

simonbyrne's draft implementations, though imperfect and incomplete, seem promising, efficient, and simple

Irrelevant. Obviously the goal of this PR is exactness.

I think that it is important to try to keep complexity down for maintainability and to reduce bugs.

Maintainability is a worthy goal, but not at the expense of correctness.

This PR can probably achieve its goals with a much simpler implementation

I assume this was revealed to you in a dream? Where's the PR, though?

Most complex implementations introduce bugs that are regressions. This is no exception

Thanks for pointing out the edge case, it seems easy enough to fix.

[LilithHafner]

I will be able to join the triage call this week briefly, if at all. My order of preference for conversion of rationals to floating point is

1. Accurate to .5 ulp and roughly as fast as or faster than floating point division
2. Accurate to 1 ulp and roughly as fast as floating point division (current, maybe actually 1.5 ulps off?)
3. Accurate to .5 ulp and substantially slower than floating point division (pr, with a little work)
4. Usually accurate to .5 ulp, but occasionally throws and substantially slower than floating point division (pr, now)

Code complexity (inc. maintainability, sysimg size, bug surface, and build time) is another important factor. I don't like added complexity, but if it is necessary to get us to option 1, I won't object too much.

[oscardssmith]

Triage thinks the performance tradeoff is worth it (once the bugs are fixed). Rationals are already pretty slow.

[Joel-Dahne]

Sorry for joining the conversation late, I was pointed here from #52859.

Regarding this part in the initial comment

However, the requirement for exactness means that the BigFloat precision would need to be manipulated, something that seems to be problematic in Julia. So implement the division using integer arithmetic and ldexp.

It's not a problem to manipulate the precision if you call MPFR directly. I believe the rational_to_bigfloat method could be replaced with the following code

function BigFloat(x::Rational{T}, r::MPFRRoundingMode=ROUNDING_MODE[]; precision::Integer=DEFAULT_PRECISION[]) where {T}
    res = BigFloat(precision = precision)
    if T <: Union{UInt8,UInt16,UInt32,Culong}
        # 64 bits is enough to represent the numerator exactly
        num = BigFloat(numerator(x), precision = 64)::BigFloat
        @ccall Base.MPFR.libmpfr.mpfr_div_ui(
            res::Ref{BigFloat},
            num::Ref{BigFloat},
            denominator(x)::Culong,
            r::Base.MPFR.MPFRRoundingMode,
        )::Int32
    elseif T <: Union{Bool,Int8,Int16,Int32,Clong}
        # 64 bits is enough to represent the numerator exactly
        num = BigFloat(numerator(x), precision = 64)::BigFloat
        @ccall Base.MPFR.libmpfr.mpfr_div_si(
            res::Ref{BigFloat},
            num::Ref{BigFloat},
            denominator(x)::Clong,
            r::Base.MPFR.MPFRRoundingMode,
        )::Int32
    else
        num_big = convert(BigInt, numerator(x))
        # FIXME: How to efficiently determine what precision is needed
        # to represent the numerator exactly
        num_big_prec = ndigits(num_big, base = UInt8(2), pad = false)
        num = BigFloat(num_big, precision = num_big_prec)::BigFloat
        @assert BigInt(num) == num_big # FIXME: Remove this
        @ccall Base.MPFR.libmpfr.mpfr_div_z(
            res::Ref{BigFloat},
            num::Ref{BigFloat},
            convert(BigInt, denominator(x))::Ref{BigInt},
            r::Base.MPFR.MPFRRoundingMode,
        )::Int32
    end
    return res
end

That is, you convert the numerator to an exact BigFloat and then use the division by integer functions from MPFR. The only non-trivial part is to determine what precision you need to store the numerator exactly. This offloads all the heavy lifting to MPFR, which means we don't have to maintain that code :smile:

[Joel-Dahne]

After some digging I found that ndigits for BigInt uses Base.GMP.MPZ.sizeinbase internally, so that should be plenty fast!

[nsajko]

@Joel-Dahne I've been working on an improved version of this PR for some time, with significant changes (a rewrite?). I'll push the changes now, so take a look if you're still interested.

[nsajko]

The code is now much more modular, and it should be faster than before in the special case of IEEE 754 floats and BigFloats (although I haven't done any performance measurements yet). Apart from the extra FP-type specialization (ldexp is still used, but only as a fallback) mentioned in the previous sentence, the code is also much simpler now. TBH the previous version of this PR was kind of monstrous.

[nsajko]

@LilithHafner I now regret my tone in my response to you above. Sorry. While I suppose we still disagree regarding the exactness trade off, you we're definitely in the right when questioning the maintainability of the previous version of this PR.

[oscardssmith]

This is still a lot of code, but it is majority tests, and the implementation does look pretty reasonable for what it is.

[nsajko]

Added a test, apart from addressing the comments.

[nsajko]

Expanded the / doc string with relevant information.

[nsajko]

Oh, some of the tests are failing on 32-bit Linux and 32-bit Windows. EDIT: it's because div(::Int128, ::Int128) allocates on 32-bit Julia, which I didn't expect

[Joel-Dahne]

Let me try to take a closer look at this today!

Regarding the comment

do we actually care about giving the user the ability to specify rounding modes here? Most other functions in Base don't, and giving up that seems like it would simplify the code a lot.

Indeed most other functions in Base don't allow you to control the rounding, and I don't believe they should! Getting rounding correct is difficult (hence this PR) and it would probably fit better in a separate package. However, Base current does support rounding when converting to Float64, e.g. Float64(1 // 3, RoundUp). Since it is already supported I think it is a valid goal to also make it correct! For example this is used in IntervalArithmetic.jl when converting a rational interval to a Float64 interval.

[nsajko]

Wow! I had no idea the rounding mode argument was already supported for converting rationals to floats! I always assumed the rounding mode argument was only accepted when converting to BigFloat. Then clearly:

1. We do need the rounding modes to be supported explicitly in this PR
2. The PR needs to actually override the rounding mode methods

[Joel-Dahne]

I would also propose to say that this change makes constructing a float from a rational "correctly rounded", instead of "exact". It seems more clear to me at least, since most rationals are not exactly representable as floats.

[Joel-Dahne]

I also made a quick benchmark comparing this implementation for BigFloat with the one I proposed above. With the version in this PR I get

julia> @benchmark BigFloat($(5 // 7))
BenchmarkTools.Trial: 10000 samples with 130 evaluations.
 Range (min … max):  743.838 ns …  3.561 ms  ┊ GC (min … max):  0.00% … 31.21%
 Time  (median):       1.072 μs              ┊ GC (median):     0.00%
 Time  (mean ± σ):     1.920 μs ± 46.631 μs  ┊ GC (mean ± σ):  18.36% ±  0.86%

   ▄▆▄▂ ▁    ▃▆▇▇█▅▃▁▁  ▁▁    ▁▁ ▁▁  ▂▁        ▂▂    ▁▁▁▁      ▂
  ███████▇▆▅███████████████▆▆██████▇████▇▇▇▇▆▅▇███▆▆▇█████▇▇▅▇ █
  744 ns        Histogram: log(frequency) by time      2.09 μs <

 Memory estimate: 1.09 KiB, allocs estimate: 40.

julia> @benchmark BigFloat($(big(2)^200 // big(3)^200))
BenchmarkTools.Trial: 10000 samples with 63 evaluations.
 Range (min … max):  792.905 ns …  4.627 ms  ┊ GC (min … max):  0.00% … 49.30%
 Time  (median):       1.005 μs              ┊ GC (median):     0.00%
 Time  (mean ± σ):     1.959 μs ± 57.511 μs  ┊ GC (mean ± σ):  20.56% ±  0.71%

  ▁▄▄▄▆▇▆▃▃▆█▇▆▄▂ ▂▁▁▄▄   ▁▁▂▂▂▂▂▃▂▂▂        ▁▁▂▁▂▂▂▂▂▂▂▂▁▁    ▂
  ██████████████████████▇▇████████████▇▅▅▄▅▅▇███████████████▇▇ █
  793 ns        Histogram: log(frequency) by time      1.97 μs <

 Memory estimate: 1.16 KiB, allocs estimate: 38.

With my version I get

julia> @benchmark Base.MPFR.__BigFloat($(5 // 7))
BenchmarkTools.Trial: 10000 samples with 968 evaluations.
 Range (min … max):   78.898 ns …  27.099 μs  ┊ GC (min … max):  0.00% … 99.52%
 Time  (median):      88.857 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   102.076 ns ± 454.203 ns  ┊ GC (mean ± σ):  12.20% ±  2.97%

            ▁▂▂▂▁       ▃█▁▂▆▄                                   
  ▂▃▅▆▆▇▇▇████████▇▅▄▄▃▄██████▇▇█▅▄▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂ ▄
  78.9 ns          Histogram: frequency by time          111 ns <

 Memory estimate: 184 bytes, allocs estimate: 4.

julia> @benchmark Base.MPFR.__BigFloat($(big(2)^200 // big(3)^200))
BenchmarkTools.Trial: 10000 samples with 204 evaluations.
 Range (min … max):  376.975 ns … 498.437 μs  ┊ GC (min … max): 0.00% … 44.14%
 Time  (median):     410.777 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   551.419 ns ±   6.129 μs  ┊ GC (mean ± σ):  6.40% ±  0.58%

    ▁▆█▄                                                         
  ▃▅████▅▃▂▂▂▃▆█▆▄▃▅▃▃▂▂▂▃▂▂▂▂▂▃▃▃▃▃▂▂▂▂▂▂▂▁▂▂▂▁▂▂▂▂▂▂▁▂▂▁▂▁▁▂▂ ▃
  377 ns           Histogram: frequency by time          764 ns <

 Memory estimate: 368 bytes, allocs estimate: 10.

So 10 times faster for machine size integers and 2 times faster for large integers. Even more if one counts the GC time, which one probably should for BigFloat.

I don't think these performance gains are essential though, either version would be fine I think.

[nsajko]

Often times temporary names are introduced only to be used once. For example the exp variable in to_floating_point_fallback.

Got rid of exp, thanks.

There are a lot of let-blocks that I don't really see the point of. For example in to_floating_point_fallback again.

Yeah, some of the lets were definitely redundant. There's less of them now. Some of the remaining ones are because I wanted to give a local scope to if blocks that have their own variables, as the lack of a new scope can be surprising there.

There are many type annotations that I don't really see the point of. To take to_floating_point_fallback as an example again, surely the compiler can verify zero(T)::T without the type annotation?

There is no way to declare a return type for a function (as in, all the function's methods) in Julia. We don't even have a guarantee for the return type of constructors (#42372)! Thus I use a type assertion after calling some function whenever I can to catch bugs earlier, and to help contain possible type-instability (which can sometimes be introduced by the caller). Note that Julia is good with eliminating type assertions based on type inference, so they're basically free when they're not helpful.

So 10 times faster for machine size integers and 2 times faster for large integers.

Thanks, I'll profile this and potentially follow up.

[Joel-Dahne]

Sounds good! I don't have any further comments at this point!

[nsajko]

Part of the reason for this PR being slow for bignums is that the current div and divrem (not touched by this PR) are unnecessarily slow. Looking at the allocation profiler, divrem(::BigInt, ::Int, ::RoundingMode{:ToZero}) is really slow and allocates too much. It calls div, which does an unnecessary promotion on the denominator, then divrem does unnecessary copying arithmetic. And GMP actually provides mpz_tdiv_q_ui that directly calculates both the quotient and the remainder. There's a bunch more examples like this, but that's clearly a matter for a separate PR. To make matters worse, touching base/div.jl seems fraught with peril, as the dispatch situation between all the different function variants (div, rem, divrem, fld, cld, fldmod) is horrifying...

[oscardssmith]

I think the base/div.jl issues are significantly improved recently.