Floating Point Arithmetic and Reproducibility

Floating point arithmetic and reproducibility

There are issues with reproducibility of even elementary floating-point operations in high-level languages, but they are usually controllable with various platform-specific operations such as setting compiler switches, using custom code to set floating-point controls and modes, or, if necessary, writing essential operations in assembly. As developed in comments, the specific problem you encountered may be that different C implementations use different precision to evaluate intermediate floating-point expressions. Often this can be controlled by compiler switches or by including casts and assignments in the expressions to require rounding to the nominal type (thus discarding excess precision).

However, the more complicated functions, such as exp and cos, are not routinely reproducible on different platforms. Although the 2008 IEEE-754 standard recommends that these functions be implemented with correct rounding, this task has not been completed for any math library with a run-time with a known bound. Nobody in the world has done the mathematics to accomplish this.

The CRlibm project has implemented some of the functions with known run-time bounds, but the work is incomplete. (Per Pascal Cuoq’s comment, when CRlibm does not have a proven run-time bound for correct rounding, it falls back to a result highly likely to be correctly rounded due to being computed with very high precision.) Figuring out how to deliver a correctly-rounded result in a bounded time and proving it is difficult for many functions. (Consider how you might prove that no value of cos(x), where x is any double value, is closer than some small distance e from the midpoint between two representable values. The midpoint is important because it is where rounding must change from returning one result to returning another, and e tells you how accurately and precisely you must calculate an approximation in order to provide correct rounding.)

The current state of affairs is that many of the functions in the math library are approximated, some accuracy looser than correct rounding is delivered, and different vendors use different implementations with different approximations. I am supposing that R uses some of these functions in its rexp implementation, and that it uses the native libraries of its target platforms, so it gets different results on different platforms.

To remedy this, you might consider using a common math library on the platforms you target (possibly CRlibm).

Reproducibility of floating point operation result

In the gaming industry this is referred to as deterministic lockstep, and is very important for real-time networked games where the clients and server need to be in agreement about the state of physics objects (players, projectiles, deformable terrain etc).

According to Glenn Fiedler's article on Floating Point Determinism, the answer is "a resoundingly limp maybe"; if you run the same binary on the same architecture and restrict the use of features that are less well specified than basic floating-point, then you can get the same results. Otherwise, if you use different compilers, or allow your code to use SSE or 80-bit floating point, then results will vary between different executables and different machines.

Yosef Kreinin recommends:

• scanning assembler output for algebraic optimisations and applying them to your source code;
• suppressing fused multiply-add and other advanced instructions (e.g. the sin trigonometric function);
• and using SSE or SSE2, or otherwise setting the FPU CSR to 64-bit. (Yes, this conflicts with Glenn Fiedler's recommendation.)

And of course, test your code on multiple different machines; take hashes of intermediate outputs, so you can tell just where and when your simulations are diverging.

Is SSE floating-point arithmetic reproducible?

SSE is fully specified*. Muller is an expert in floating point arithmetic; who are you going to trust, him or some guy on a gamedev forum?

(*) there are actually a few exceptions for non-IEEE-754 operations like rsqrtss, where Intel never fully specified the behavior, but that doesn't effect the IEEE-754 basic operations, and more importantly their behavior can't actually change at this point because it would break binary compatibility for too many things, so they're as good as specified.

Floating point math in python / numpy not reproducible across machines

Floating point calculations are not always reproducible.

You may get reproducible results for floating calculations across different machines if you use the same executable image, inputs, libraries built with the same compiler and identical compiler settings (switches).

However if you use a dynamically linked library you may get different results, because of numerous reasons. First of all, as Veedrac pointed in comments it might use different algorithms for its routines on different architectures. Second, a compiler might produce different code depending on switches (various optimizations, control settings). Even a+b+c yields non-deterministic results across machines and compilers, because we can not be sure about order of evaluation, precision in intermediate calculations.

Read here why it is not guaranteed to get identical results on different IEEE 754-1985 implementations. New standard (IEEE 754-2008) tries to go further, but it still doesn't guarantee identical results among different implementations, because for example it allows implementers to choose when tinyness (underflow exception) is detected

More information about floating point determinism can be found in this article.

In R, how can I get PRNG to give identical floating point numbers between platforms?

(Comments from Oct. 29 and Nov. 2 moved here and edited.)

I should note that such subtle reproducibility issues with pseudorandom number generators (PRNGs) can occur when floating-point arithmetic is involved. For instance, Intel's instruction set architecture might make use of 80-bit extended precision for internal arithmetic. Extended precision, though, is only one way (among a host of others) that floating-point arithmetic might lead to non-reproducible pseudorandom numbers. Consider that Intel's and Arm's instruction set architectures are different enough to cause reproducibility issues. (If I understand, an Arm instruction set is what is used in Apple's M1 chip.)

By contrast, integer arithmetic has fewer reproducibility problems.

Thus, if bit-for-bit reproducibility matters to you, you should try to find an R language PRNG that uses only integer operations. (Indeed, computers generate pseudorandom floating-point numbers via integers, not the other way around, and most PRNGs produce integers, not floating-point numbers.)

For instance, for uniform variates, take the integer output of the Mersenne Twister algorithm without manipulating it. For Gaussian (and exponential) random variates, there is fortunately an algorithm by Karney that generates arbitrary-precision variates with only integer operations. Also, consider rational arithmetic built on underlying integer operations.

REFERENCES:

Karney, C.F.F., 2016. Sampling exactly from the normal distribution. ACM Transactions on Mathematical Software (TOMS), 42(1), pp.1-14.

Is specifying floating-point type sufficient to guarantee same results?

IEEE 754-2008 clause 11 describes what is necessary for reproducible floating-point results. This is largely:

• Bindings from the programming language to IEEE 754 operations (e.g., a*b performs the floating-point multiplication specified by IEEE 754).
• Ways to specify that reproducible results are desired. E.g., disable default language permissions to use wider precision than the nominal types of objects.
• Accurate conversions to and from external decimal character sequences.
• Avoiding some of the fancier features of IEEE 754.

These requirements are poorly supported in today’s compilers.

Adding .5 will not be a problem. All normal floating-point implementations represent .5 exactly and add it correctly. What will be a problem is that a language implementation may add .5 and keep the result precisely (more precisely than a usual double) while another implementation rounds the result to double. If math library routines are used (such as cos and log), that is another problem, because they are hard to compute well, and different implementations provide different approximations.

IEEE 754 is a good specification. Ideally, you would specify that implementations of your specification conform to IEEE 754.

C: Is the result of floating point arithmetic ALWAYS normalized?

Does the C language always store the floating-point numbers in the normalized form?

"It depends." As we'll see, it's more the hardware than the C language that determines this.

If the implementation uses something other than IEEE-754, there's not much we can say.

If the implementation does use IEEE-754, then all numbers are always stored normalized except the ones that aren't, namely the subnormals.

Does this also hold true for the results obtained after some arithmetic (addition, multiplication)?

Yes. (More on this below.)

Is it dependent on the language or the hardware - FPU?

It is typically dependent on the hardware. Most of the time, assuming the target processor supports floating-point at all, a C program is compiled straight to native floating-point instructions, without any language- or compiler-imposed extra processing. (This is by contrast to, for example, Java, which does have a language-imposed floating-point definition, which is implemented in part by the JVM.)

The C standard does have an optional section, "Annex F", which specifies a bunch of specific floating-point behavior, conformant with IEEE-754.

Now, if the C implementation adopts Annex F and is conformant with IEEE-754 (typically because the underlying hardware is, too), the answers to your first two questions become even easier. In IEEE-754 binary arithmetic, with one exception, there are no ambiguities of representation. Every number that can be represented in normalized form has exactly one normalized representation. Every number that cannot be represented in normalized form, but that can be represented as a subnormal, has exactly one subnormal representation. These constraints apply to every IEEE-754 floating-point number, including (naturally enough) the results of other operations.

(The exception, as Eric and Chux have reminded me in a comment, is zero, which IEEE-754 has two of, positive and negative.)

So the answer to "Is the result always normalized?" is "No" (because IEEE-754 most definitely has those subnormals, and of course zero), but if the question is "Does every number have a unique representation?", the answer is mostly "Yes." (Except, again, for zero. Or if you're one of the rare few who is doing something with the IEEE-754-2008 decimal formats, which are much less unique.) See also How to distinguish between 1 and zero floating-point values?

The last question, I suppose, is "How many C implementations adopt Annex F?", or, stated another way, "How many processors comply with IEEE-754?" For the CPU's on general-purpose computers (mainframes and personal computers), as far as I know the answer these days is "All of them". GPU's, on the other hand, are deliberately not quite compatible with IEEE-754 (because they can be even more insanely efficient that way). Microprocessors, for "embedded" work, I'm not so sure about. (Often they don't have viable floating-point at all.)

C/C++: Are IEEE 754 float addition/multiplication/... and int-to-float conversion standardized?

No, unless the macro __STD_IEC_559__ is defined.

Basically the standard does not require IEEE 754 compatible floating point, so most compilers will use whatever floating point support the hardware provides. If the hardware provides IEEE compatible floating point, most compilers for that target will use it and predefine the __STD_IEC_559__ macro.

If the macro is defined, then IEEE 754 guarantees the bit representation (but not the byte order) of float and double as 32-bit and 64-bit IEEE 754. This in turn guarantees bit-exact representation of double arithmetic (but note that the C standard allows float arithmetic to happen at either 32 bit or 64 bit precision).

The C standard requires that float to int conversion be the same as the trunc function if the result is in range for the resulting type, but unfortunately IEEE doesn't actually define the behavior of functions, just of basic arithmetic. The C spec also allows the compiler reorder operations in violation of IEEE754 (which might affect precision), but most that support IEEE754 will not do that wihout a command line option.

Anecdotal evidence also suggest that some compilers do not define the macro even though they should while other compilers define it when they should not (do not follow all the requirements of IEEE 754 strictly). These cases should probably be considered compiler bugs.

---

Related Topics