Hacker News Reader

> Things are much harder with 64 bits, unfortunately.

Solved and published over 15 years ago. The underlying methods were published over 20 years ago (V. Lef`evre, J.M. Muller, and A. Tisserand. Towards correctly rounded transcendentals. and Lefvre's associated PhD thesis).

https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804/file...

The Correctly Rounded Math Library (CRLIBM) is available for several languages.

CRLIBM is great, but it is significantly slower than alternatives since it requires very high internal precision. The cool thing with this paper is that it gives correct rounding without as much of a cost for 32 bit.

That's a bit glib. I don't believe that the bounds are fully proven for all the transcendentals (and I believe that exponentiation is only "proven" over restricted subdomains).

Single precision floating point is a small enough domain that it can be fully enumerated--you can prove that everything rounds correctly and find the problematic cases. 24-bits just isn't that large to modern computers.

Double precision floating point is still too large to fully enumerate. So, there could be an exception lurking in there that goes out a long way.

In reality, it seems that rounding requires about 2n + log2 n bits and never 3n bits (I don't think we've ever found an exception to this). However, if you can prove that, you should submit it as it would be a significant advance in the world of mathematics.

In general, there shouldn't be a rule like 3n bits. For random functions, 1 in 8^n values should require more than 3n bits. Since n bit inputs give 2^n possible values, this implies that for bit length n, there is a 4^-n chance that any given function will require more than 3n bits.

Except that this isn't the case. We know empirically because we have enumerated the single precision space completely.

So, the problem isn't bounded by simple probabilities. There are other properties involved that cause it to diverge from the "naive" probabilistic estimates. What those properties are we don't know yet.

The "rounding" problem is called the "Table Maker's Dilemma" and it descends into some pretty fundamental transcendental number theory.

This is not an "easy" problem. Proof that PI was transcendental is relatively modern (1885).

This doesn't contradict what I said. I said that the odds of a function needing more than 3n digits was 4^-n. For single precision, that means you would expect almost all functions to be OK with 3n bits (1 in 2^46 wouldn't be).

if it was a randomly selected function, which it is not.

In practice the low-order bits of transcendental functions are pretty much indistinguishable from random for _most_ inputs. You do have to be very careful close to especially nice values though, because this falls apart catastrophically there (e.g. logarithms close to 1)--it should be pretty obvious from the Taylor expansions why this happens. The good news is that when this does happen, almost by construction you usually don't end up very close to an exact halfway point for rounding.

> In practice the low-order bits of transcendental functions are pretty much indistinguishable from random for _most_ inputs.

In my experience, low order bits of transcendental functions make terrible rng sources (even discounting their runtime). Your example of the logarithm is great -- for large inputs, you need large changes to flip even low order bits.

If you're just looking for hard-to-round inputs, you can also enumerate the double precision space completely with some cleverness and modest patience.

Elkies [0] computes a piecewise linear approximation to his (number-theoretic) function of interest and does a lattice reduction on each piece to filter out the vast majority of the space. The same approach, works almost everywhere on all of the transcendental functions I've tried, and I believe the CRlibm group uses it.

[0] https://arxiv.org/pdf/math/0005139.pdf

Oooh, that's an interesting take.

Apple had a library that did that for 80-bit floating point, for both the 6502 (!) and 68000: SANE (https://en.wikipedia.org/wiki/Standard_Apple_Numerics_Enviro...)

And yes, performance wasn’t good, especially compared to the x87 that IBM PCs used. It did claim to round the last bit correctly, though, where early x87s could have large errors (https://randomascii.wordpress.com/2014/10/09/intel-underesti...)

> No more getting slightly different result when changing platform, OS, or even just updating libm

Optimizing compilers are usually used with flags that cause them to break IEEE754.

Still, it would be nice if results from low-optimization runs might be bit-exact.

Yes, indeed it is the case with fused multiply-add operations, as I wrote below in another thread. By default at -O3, unfortunately, gcc and icc both take responsibility for deciding whether or not to use FMA instructions when they see (a x b + c) expressions. As a result we can even see slightly different result across compiles!!

However, as far as I understand, specifying a proper language standard (for example -std=c99, -std=c++14, etc.) restores sane behavior (see -ffast-math, and more specifically -fexcess-precision).

Are there other cases where -O3 breaks IEEE754 on modern setups?

Fast math flags are actually far more complicated in C/C++ compilers than just a simple -ffast-math flag, although I suppose most users may not be acquainted with the degree of control possible. In LLVM, for example, fast math flags are actually a set of 7 independently-toggled flags (no-signed-zero, no-infinity, no-NaN, allow FMAs, allow reassociation, allow approximate reciprocal, and allow approximate functions), that at the C/C++ level can be toggled either by command-line options or on a per-scope basis. ICC actually has two different levels of fast math [1].

The dark underside of IEEE 754 is that full support for it also requires features that most languages don't provide access for: rounding mode and exception handling (aka sticky bits). In C/C++, using these features correctly requires STDC FENV_ACCESS support which is only in ICC and recent clang (neither GCC nor MSVC support it). And there's usually two hardware bits for handling denormals (flush-to-zero and denormals-are-zero) that aren't in IEEE 754, which tend to be on by default when you're compiling for some targets (most notably GPUs).

[1] https://software.intel.com/content/www/us/en/develop/documen... -- note fp-model=fast=[1|2]

What is a "per scope basis"?

You can change it between any pair of curly braces.

E.g.:

  void foo(int N) {
    for (int i = 0; i < N; i++) {
      #pragma STDC FENV_ACCESS ON
      // This applies only to the for loop...
    }
    // But not here, for example!
  }

> Yes, indeed it is the case with fused multiply-add operations, as I wrote below in another thread

This is not always true. Some FMA operations round after the multiply, and some do not.

> Optimizing compilers are usually used with flags that cause them to break IEEE754.

For languages that don't mandate IEEE-754. There are plenty of languages that do, like Java, JavaScript, C#, and WebAssembly.

It's worth noting that usually such optimizations are things like enabling algebraic reassociation and commutation (of primitives like +, -, /, ==, etc) operations, fused-multiply add, and others. It's really not clear to me how much disabling all those optimization costs in the real world.

Things get crazy when compilers want to do major reorganization of loops, optimizing stencils, polyhedral optimizations. I am not an expert here, but I think these rely on properties of float arithmetic that don't hold in all cases, i.e. they rely on UB.

> It's really not clear to me how much disabling all those optimization costs in the real world.

I assure you that SPECcpu scores are a lot lower, if that's what you care about. And yes, it's extremely helpful to do major reorganization of loops if the program has multi-dimensional loops in the wrong order.