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

Base: correctly rounded floats constructed from rationals #49749

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

nsajko
Copy link
Contributor

@nsajko nsajko commented May 11, 2023

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
Copy link
Contributor Author

nsajko commented May 11, 2023

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.

Copy link
Contributor

@Seelengrab Seelengrab left a comment

Choose a reason for hiding this comment

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

We have lots of screen space - please, let's use it not just vertically. We are fine with 92 characters per line (sometimes more if it's just a few characters over) after all.

base/mpfr.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
test/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
@nsajko nsajko force-pushed the exact_rational_to_float branch 2 times, most recently from 2a10919 to 50d5aa8 Compare May 12, 2023 13:00
@nsajko
Copy link
Contributor Author

nsajko commented May 12, 2023

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.

Copy link
Contributor

@Seelengrab Seelengrab left a comment

Choose a reason for hiding this comment

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

Thank you for sticking with this! I haven't added the style points to every occurence, but please apply it uniformly.

We can worry about performance afterwards, since this PR got started with correctness. Talking about potential performance improvements without benchmarks is a bit of a moot point.

base/float.jl Outdated Show resolved Hide resolved
base/mpfr.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
test/rational.jl Outdated Show resolved Hide resolved
test/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/mpfr.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
base/rational.jl Outdated Show resolved Hide resolved
Copy link
Contributor

@Seelengrab Seelengrab left a comment

Choose a reason for hiding this comment

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

Looks good!

I think the gain in correctness is nice, now let's hope this doesn't cost too much performance 😬

@nsajko
Copy link
Contributor Author

nsajko commented May 17, 2023

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
Copy link
Contributor

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

base/rational.jl Outdated Show resolved Hide resolved
@nsajko
Copy link
Contributor Author

nsajko commented May 17, 2023

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 nsajko force-pushed the exact_rational_to_float branch 2 times, most recently from 424dc82 to dd48fdd Compare May 17, 2023 16:18
@nsajko
Copy link
Contributor Author

nsajko commented May 17, 2023

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

c = rational_to_float_components(num, den, prec, BigInt, rm)

││││┌ @ 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
Copy link
Contributor

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
Copy link
Contributor Author

nsajko commented May 18, 2023

Is there still something to do here from my side?

@nsajko
Copy link
Contributor Author

nsajko commented May 18, 2023

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.

Copy link
Contributor

@Seelengrab Seelengrab left a comment

Choose a reason for hiding this comment

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

Other than these, yeah, looks good from my POV 👍

base/float.jl Outdated Show resolved Hide resolved
test/rational.jl Outdated Show resolved Hide resolved
@nsajko nsajko force-pushed the exact_rational_to_float branch 2 times, most recently from 479e877 to 508cfdd Compare May 19, 2023 10:33
@giordano giordano added rationals The Rational type and values thereof maths Mathematical functions forget me not PRs that one wants to make sure aren't forgotten labels May 24, 2023
@nsajko
Copy link
Contributor Author

nsajko commented Jan 14, 2024

Added a test, apart from addressing the comments.

@nsajko
Copy link
Contributor Author

nsajko commented Jan 14, 2024

Expanded the / doc string with relevant information.

(; integral_significand, exponent)
end

function to_float_components(::Type{T}, num, den, precision, max_subnormal_exp, romo, sb) where {T}
Copy link
Member

Choose a reason for hiding this comment

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

does this method need to exist? if this is purely for internal use, surely the caller can be responsible for converting num to the type they want.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It's also used in mpfr.jl

base/rational.jl Outdated
(numerator(y), denominator(y))
end

function rational_to_floating_point(::Type{F}, x, rm = RoundNearest, prec = precision(F)) where {F}
Copy link
Member

Choose a reason for hiding this comment

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

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

We need the rounding modes for BigFloat. Also, having rounding modes other than the default one enables more thorough testing that applies to all rounding modes.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Most of the logic for the rounding modes is now in rounding.jl anyway.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

to_floating_point_impl(::Type{T}, ::Type{S}, num, den, rm, prec) where {T<:Base.IEEEFloat,S} could be simplified by requiring rm == RoundNearest, but then the tests would suffer.

Copy link
Contributor Author

@nsajko nsajko Jan 14, 2024

Choose a reason for hiding this comment

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

giving the user the ability

to be clear, that particular method is not meant to be public (I think the PR doesn't add any new public names)

Copy link
Member

Choose a reason for hiding this comment

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

Can't we do the BigFloat ones by just doing arbitrary precision BigFloat division? I'm not really sure if it's worth adding a bunch of code to make BigFloat(::Rational) fast.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I don't think that would be an improvement, by any metric. The relevant code was already added in #50691 (into rounding.jl), so now it doesn't cost anything to support the known rounding modes for BigFloat, although the IEEEFloat method could be simplified a bit if necessary, at the cost of some testing, as mentioned above.

Copy link
Contributor Author

@nsajko nsajko Jan 14, 2024

Choose a reason for hiding this comment

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

by any metric

Maybe the way you and @Joel-Dahne propose would be faster, though? But still, having all float types under a single code path seems better as any bugs should manifest (and be fixed) sooner, so it's more maintainable.

@nsajko
Copy link
Contributor Author

nsajko commented Jan 14, 2024

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
Copy link

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
Copy link
Contributor Author

nsajko commented Jan 15, 2024

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
Copy link

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.

@nsajko nsajko changed the title Base: make constructing a float from a rational exact Base: correctly rounded floats constructed from rationals Jan 15, 2024
Copy link

@Joel-Dahne Joel-Dahne left a comment

Choose a reason for hiding this comment

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

In general I think the change looks good! As far as I can tell the most important methods are to_float_components, to_floating_point_fallback and the to_floating_point_impl implementation for machine floats. The rest is more or less gluing code and tests (both of which there are a lot of).

In general the coding style is quite different from my personal one. That is of course fine, but I have some general comments:

  • Often times temporary names are introduced only to be used once. For example the exp variable in to_floating_point_fallback.
  • There are a lot of let-blocks that I don't really see the point of. For example in to_floating_point_fallback again.
  • 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?

Comment on lines +134 to +131
is_zero = num_is_zero & !den_is_zero
is_inf = !num_is_zero & den_is_zero

Choose a reason for hiding this comment

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

Here I think it would be natural to make an early return on these checks. So replace them with

num_is_zero & !den_is_zero && return zero(T)
!num_is_zero & den_is_zero && return sb * T(Inf)
num_is_zero & den_is_zero && return T(NaN)

That removes the need for the if-statement further down and , which, for me, simplifies the flow.

Similar things could be done in the other implementations, though it makes less of a difference there.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I prefer to avoid return when possible. It's a jump/goto statement, breaking nice control-flow properties known as structured programming. I do use jumps in control flow, but only when necessary.

Copy link
Contributor Author

@nsajko nsajko Jan 15, 2024

Choose a reason for hiding this comment

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

I used to be a Go nut (as in "Golang aficionado"), so I see where you're coming from, though.

Choose a reason for hiding this comment

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

It's only a matter of style of course. I'm not going to fight you over it 😄

base/rational.jl Outdated Show resolved Hide resolved
# `BigInt` is a safe default.
to_float_promote_type(::Type{F}, ::Type{S}) where {F,S} = BigInt

const BitIntegerOrBool = Union{Bool,Base.BitInteger}
Copy link

@Joel-Dahne Joel-Dahne Jan 15, 2024

Choose a reason for hiding this comment

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

This variable is only used once, I think it is more natural to just use the union directly. Though in this case I guess it does make the method declaration below slightly too long to fit on one line.

Copy link
Contributor Author

@nsajko nsajko Jan 15, 2024

Choose a reason for hiding this comment

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

My reasoning here was basically:

  1. I'm not exactly sure what's the preferred way to break lines (or format code in general) in Julia source in the Julia project, so I prefer to avoid it completely
  2. Adding a constant is basically free (as long as it's not huge, doesn't mess with precompilation somehow, etc.)

Choose a reason for hiding this comment

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

I understand, it also explains some of the style of code in other places!

base/mpfr.jl Outdated Show resolved Hide resolved
base/rational_to_float.jl Outdated Show resolved Hide resolved
test/rational.jl Show resolved Hide resolved
test/rational.jl Outdated Show resolved Hide resolved
test/rational.jl Outdated Show resolved Hide resolved
test/rational.jl Outdated Show resolved Hide resolved
base/rational_to_float.jl Show resolved Hide resolved
@Joel-Dahne
Copy link

Joel-Dahne commented Jan 15, 2024

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 nsajko added the bugfix This change fixes an existing bug label Jan 15, 2024
@nsajko
Copy link
Contributor Author

nsajko commented Jan 15, 2024

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
Copy link

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

@nsajko
Copy link
Contributor Author

nsajko commented Jan 16, 2024

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
Copy link
Member

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

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 JuliaLang#45213

Updates JuliaLang#50940

Updates JuliaLang#52507

Fixes JuliaLang#52394

Closes JuliaLang#52395

Fixes JuliaLang#52859
@nsajko
Copy link
Contributor Author

nsajko commented May 12, 2024

Fixed conflict. Could this be merged as per #53641 (comment) (caused by correctness fix, and the performance regression mainly applies to bignums)? The performance for bignums can be fixed later, this PR is already quite large

@nsajko
Copy link
Contributor Author

nsajko commented May 12, 2024

The test failures are real.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bugfix This change fixes an existing bug forget me not PRs that one wants to make sure aren't forgotten maths Mathematical functions rationals The Rational type and values thereof
Projects
None yet
7 participants