Comment by codeflo | Hacker News Reader

Comment by codeflo | Hacker News Reader

There's a widespread misconception (not shared by the article author) that floating point arithmetic is "imprecise" in the sense that the result is off by some amount of random noise. On the contrary, IEEE floating point results are precisely specified to produce the closest representable value to the mathematically exact result.

For efficiency reasons, library functions (and sometimes, sadly, hardware implementations) don't always guarantee this (closest representable) exact result. That's fine if the function is then called "approximate_inverse_square_root" or something. Sometimes being off by some epsilon is a good tradeoff for 10x efficiency. I'm unhappy when such a tradeoff is smuggled into my math library without warning.

I checked Rust's source code to confirm the article's claim that it doesn't have a precise log function, and indeed, here's the implementation:

    pub fn log(self, base: f64) -> f64 {
        self.ln() / base.ln()
    }

I'm sure other math libraries aren't better. But this is not a correct implementation of arbitrary-base logarithm. A function like that perhaps shouldn't even be offered in the standard library at all (since it's so trivial to begin with), or at the very least, not with that name. If a programmer wants to opt-in to a fast but slightly wrong value, they should do so explicitly, in my opinion.

Thank you! This point has been driving me mad for the last couple of years.

While answering a sibling comment, I found a paper analyzing various libms for their precision: https://homepages.loria.fr/pzimmermann/papers/glibc238-20230... (2023). I only skimmed it, but it looks like only LLVM's libm guarantees 0.5 ulp for every supported single-precision operation, and everyone else is wildly off. That's a pretty good result, though -- it means there's at least one reasonably compliant libm :)

Another sibling comments says that guaranteeing 0.5 ulp for double-precision operations is nigh impossible (and then another says it is after all). I don't have the knowledge to confirm which is true, but it's possible that this is the best we can get.

Table--Maker's dilemma [0] says that we probably can't get correct rounding for some functions. That said, you can still calculate every function within 1 ULP by using extended precision (for example the double-double arithmetic, which uses two doubles to represent a single fp number) and rounding.

But I don't think anyone does this in practice. Simulations don't need perfect precision, anything involving real data has to deal with measurement errors which are usually much larger than f64 ULP scales, and then there's configurable arbitrary-precision for when it's really needed.

[0]: https://people.eecs.berkeley.edu/~wkahan/LOG10HAF.TXT, which coincidentally is also about the log function

> On the contrary, IEEE floating point results are precisely specified to produce the closest representable value to the mathematically exact result.

Not sure about Rust, but C++ only guarantees that for four arithmetic operations and `sqrt`. Not `pow`, not `log`. Not even when when all inputs/outputs are integers. Sometimes you can even get `(int)pow(10, 2) == 99`, which I believe is fully standard-compliant.

Yeah, those operations (+fma) are the only ones for which IEEE 754 mandates correct rounding.

Calculating correct rounding for other functions without hardware support is both challenging to implement and is often really slow. And that's before one gets into special functions. Or, even worse, inverses of special functions. Sometimes you might be lucky when those don't over/underflow around the edges of the domain

It's still challenging to implement with hardware support but you don't notice because it challenges Intel instead of you. There is a reason that no instruction set since x87 implements transcendentals - and x87 did them in microcode.

True!

Older numerical libraries such as cephes had to work with heterogeneous non-IEEE floating point implementations: you can still see remnants of VAX tests in the docs, for example. That, too, must've been quite a struggle

I am curious what the correct algorithm would be, if for no other reason than to find out how on god's green earth it could be slower than using a division as an erstwhile shortcut. 8I

In this particular case, there never exists any need to compute a logarithm of an arbitrary base, so such a function must not exist in a library.

If someone had the weird idea to write some formula with logarithms in a random base they should convert that formula to use binary logarithms, because there is no advantage in using another base.

Nonetheless, if one would want to include such a function, or any other function, one could reduce the argument to some small range using mathematical properties of the function and then some polynomial approximation should be determined, in order to use it to compute the function values.

If you determine a polynomial approximation for a function, whichever function it is, you can compute the entire function faster and more accurately than when you compute the function as a composition of other functions, one or more of which are computed by their own polynomial approximations.

In order to ensure that the function is rounded correctly, one must determine how many polynomial terms are needed and how big must be the numbers used for storing the coefficients.

For a function with a single argument, like "log2", it is much easier to do an exhaustive search to find all the argument values that would need more polynomial terms or more bits in the coefficients, to guarantee correct rounding.

For functions with 2 arguments, it can be difficult or impossible to find all problematic argument pairs, unless the function has some special mathematical properties that would allow a formal proof about the accuracy needed in polynomial evaluation to ensure correct rounding.

Reply on news.ycombinator.com