C library for multiple-precision floating-point arithmetic with correct rounding
mpfr.org
mpfr.org
SLEEF itself is used by PyTorch.
What exactly is correct rounding? Sounds kinda like perfect approximation.
I search/skimmed all references to “round” on the linked page, but didn’t find a definition of what “correct” was. Wikipedia page on rounding doesn’t use the term correct anywhere. Is this a mathematical adjective in this case? Or just the opinion of the authors?
I’m well aware of banker’s rounding, and rounding up, and round down, and rounding towards zero, and rounding towards signed infinity, etc. But “correct” as a formal modifier is new to me. What did I miss?
> I’m well aware of banker’s rounding, and rounding up, and round down, and rounding towards zero, and rounding towards signed infinity, etc.
MPFR supports most of those rounding modes you mentioned:
https://www.mpfr.org/mpfr-current/mpfr.html#Rounding
> But “correct” as a formal modifier is new to me
It just means that once you select a rounding mode, it is guaranteed (modulo bugs) to give you the correct answer rounded according to your selected rounding mode.
At first one can think: "obviously, who wants incorrect rounding?". And indeed all recent floating-point hardware implements correct rounding... But only for 32- or 64-bit floats, and only for addition, subtraction, multiplication and division.
By contrast, most libc/libm implementations do not have correct rounding for transcendental functions (cos, sin, tan, exp, log, etc.). There are some exceptions (Rlibm, CRlibm).
MPFR implements correct rounding
- at any bit width and
- for all transcendental functions.
For another example, interval arithmetic is a really cool approach to keeping an eye on your accuracy. There, you want to round your upper bounds up, and round your lower bounds down.
For the ordinary arithmetic operations, including for square root and for reciprocal square root, there are algorithms to do correct rounding without excessive computation. Because of that, the standard for floating-point arithmetic specifies that a conforming implementation must use correct rounding for those operations.
On the other hand, for transcedental functions, like exponential, logarithm, sine and so on, there is no general method to guarantee that the result will be correctly rounded without having to compute an unlimited number of result bits before being able to round correctly.
Because of the difficulty, many mathematical libraries and also all hardware implementations of transcedental functions round correctly many of the possible results, but not all of them.
Nevertheless, for each usual transcedental functions and for the numbers of bits used frequently, i.e. single-, double- and quadruple-precision, there are known algorithms to round correctly, which have usually been found by exhaustive searches for corner cases, when more result bits must be computed than for most other results.
That may not be what it does, but that's certainly what I'd expect "correct rounding" to do.
There is a similar concept which is supported experimentally by MPFR, which is faithful rounding. With faithful rounding, a number is rounded either to the closest value up or to the closest value down, but it is not necessarily the nearest of them. That is, if a and b are representable and no value between them is representable, and if a < x < b, then with faithful rounding x will be rounded to any of a or b. This can sometimes be faster, but is not reproducible any more.
In fact I believe MPFR depends on GMP.
Includes not only infinite-precision [correction: arbitrary-precision] floats (not sure about the correct rounding part) but also integers (so subsumes GMP), rationals, complex numbers, and polynomials.
The original library claim for correct rounding is very important for many tasks, and is the entire reason for that library. I doubt your library does it, it it would claim it, since it's extremely hard to do and would be a great selling point.
Here's [1] a slightly older paper with lots of details on the problems in getting such a library to work. There's massive literature on this
[1] https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804/file...
When reproducibility doesn't matter, it's basically always better to use a not-necessarily-correctly-rounded implementation at somewhat higher precision, as this gets you smaller overall errors for less effort. But there are a lot of domains where reproducibility matters, and that's where correct rounding comes in.
If it doesn't matter, then high precision rarely matters either. Once you have results with unknown noise in them, and you do a few operations on those, the noise compounds, at approximately sqrt(N) bits for N steps of computation (rough estimate, it certainly depends on algorithm selection, compiler nuances, stability of problem, etc.). So high precision degrades quickly when doing any work without using a reproducible and understood rounding mode.
>correct rounding is that it is the easiest mechanism to unlock reproducibility
Any rounding method works works for reproducibility - truncate, round to zero, Bankers rounding, etc., all (roughly) equivalently hard/easy to implement. However, pretty much none of these are possible without the ability to do correcting rounding first as a result of the Table Makers' Dilemma.
For example, Banker's rounding is generally better for doing computation than the usual correctly rounded. Knuth has some papers on this somewhere which are fascinating.
Also note most of the non-reproducible issues in numerical computation are compiler added - compile code with GCC version X, get one result, compile with later version, and results may change. Or across compilers. Negative 0, denormals, rounding modes, all need carefully dealt with. It can be done, but is not trivial. Keeping floating point solid is very, very tricky over compilers and time.
This isn't really true for any of the numerous numerical algorithms that satisfy backwards error bounds; there you would always prefer to have a well-rounded arithmetic at higher precision, if reproducibility is not required.
> Any rounding method works works for reproducibility
Only if you round _correctly_ in that direction.
> Banker's rounding is generally better for doing computation than the usual correctly rounded.
Banker's rounding _is_ the default IEEE 754 rounding mode.
This, FWIW, is what GCC uses MPFR for.
"If you input 0.1 into atan, the result should round up when outputting 5 digits, down with 8 digits and up with 12 digits"?
Also, "round up" is not any easier than "round correctly" if you require it to be consistent. Say you compute a result like 9.004 and want to round to two digits - if the exact result would have been 8.991, then rounding your computation as 9.01 is different from the standard-specified result of 9.00.
So your standard would literally have to specify the desired output for every combination of operation, input, and output precision - obviously impossible.
True. I've edited my comment to reflect the error.
> I doubt your library does it
Just for the record, CLN is not mine. I'm just a (casual) user. And you're right, it probably does not do correct rounding. I was unaware of the significance of this.
The above notwithstanding, I think CLN is pretty cool and deserves more attention than it gets.
Agreed - I think all such libraries are interesting. Knowing the pitfalls of them is also useful and not obvious.
And in case its relevant, in this setting, students are only allowed to use single (32-bit) precision floats and operations.
That's actually a good method that lots of numerical analysts use. I think C99 has several rounding modes - I often have code that runs the same computation under each mode and look at what your doing.
>not all crappy code generates a wide window
True, often you need to find the values that make it puke. I think there was a paper recently on doing this automatically.... If I can refind it I'll post - it was cool stuff.
>students are only allowed to use single (32-bit) precision floats and operations
For stuff this small you can often brute force all values or at least all edge case values. I sometimes do the trick of converting 32 bit integers (watch out for weird C aliasing rules) to hit lots of numbers near precision edges and run them through, denormals, other edge cases, and pseudo-fuzz the results. This will sometimes turn up more.
If you teach, and have not worked through some numerical analysis books that provide good background for this stuff, read Goldberg's paper "What Every Computer Scientist Should Know About Floating-Point Arithmetic" [1]. The next good thing to dig through is Higham's book "Accuracy and Stability of Numerical Algorithms" [2]
Those will give you lots of places to look for issues, and the ability to analyze algorithms more concretely for floating point issues.
[1] https://docs.oracle.com/cd/E19957-01/800-7895/800-7895.pdf [2] http://ftp.demec.ufpr.br/CFD/bibliografia/Higham_2002_Accura...
In my class, FP computation is something I try to teach the students to be mindful of, but it is hardly the focus of the class (about certain kinds of data visualization). FP computation arises in lots of different places in many different functions, often in the context of a larger data processing pipeline that actually starts and ends with 16-bit integers (e.g. medical imaging data). So I'm not sure how I'd thoroughly exercise the FP-related code in the way that you describe. Hence my current use of something I can easily control at the outset of the whole pipeline (the direction of rounding).
Do you have any ideas on injecting randomized rounding mode changes prior to FP ops in the assembly, or something functionally equivalent for characterizing C code (ideally without code changes)? I just found [1] which seems extremely related, but I'll have to look into whether the code still works.
[2] https://tel.archives-ouvertes.fr/tel-03110553/document
https://www.mpfr.org/mpfr-4.0.1/timings.html
Also, MPFR depends on GMP which has integers and rationals. Then there is MPC which in turn depends on MPFR and provides complex numbers.
Also another amazing piece of software behind quickjs:
I personally wish more of HN were actually this kind of thing; less about another web tech piece that will be absolute in 3 months anyway.
But hey variety is good and it takes all kinds of us and all kinds of what we all do.
Also coding a simple emulator, like for the chip-8 for example, and learning what a carry flag does would help with basic understanding of the problem.
These days I'm super lazy and just Python the big int problems on things like Advent of Code. Much more boring. And because of floating point issues, I haven't seen them give something in the real domain, as compared to whole numbers.
The mathematical setup work is not rigourous enough. I allow myself to say that, because I know my maths teachers would have rejected such work.
Solution, cleanup and redo, this time very rigourously, the maths setup work of this document. Yes, it will probably add a few pages.
I had a look at SLEEF, but I could not find the mathematical proof document of the accuracy of its approximants.
If so, "because my maths teachers" is not a very good rebuttal to what is obviously a form of program documentation and not peer reviewed publication.
The document could use some editorial cleanup, but so don't see anything systematically wrong with the math. It's understandable by people versed in the language of numerical analysis, which is certainly the target audience.
I don't think there is an accuracy proof document of SLEEF approximants.