Can I use rounding to ensure determinism of atomic floating point operations?

Question

I am developing a C application which needs floating-point determinism. I would also like the floating-point operations to be fairly fast. This includes standard transcendental functions not specified by IEEE754 like sine and logarithm. The software floating-point implementations I have considered are relatively slow, compared to hardware floating point, so I am considering simply rounding away one or two of the least significant bits from each each answer. The loss of precision is an adequate compromise for my application, but will this suffice to ensure deterministic results across platforms? All floating-point values will be doubles.

I realize order of operations is another potential source for variance in floating-point results. I have a way to address that already.

It would be terrific if there were software implementations of the major floating-point hardware implementations in use today, so I could test a hypothesis like this directly.

Answer 1

As I understand it, you have a software implementation of a transcendental function like sin(x), expressed in terms of IEEE standard operations such as floating point add and multiply, and you want to ensure that you get the same answer on all machines (or, at least, all the machines that you care about).

First, understand: this will not be portable to all machines. E.g. IBM mainframe hex floating point is not IEEE, and will not give the same answers. To get that exact, you would need to have a software implementation of the IEEEE standard operations like FP add and multiply.

I'm guessing that you only care about machines that implement IEEE standard floating point. And I am also guessing that you are not worried about NaNs, since NaNs were not completely standardized by IEEE 754-1985, and two opposite implementations arose: HP and MIPS, vedrsus almost everyone else.1

With those restrictions, how can you get variability in your calculations?

(1) If the code is being parallelized. Make sure that is not happening. (It's unlikely, but some machines might.) Parallelization is a major source of result variation in FP. At least one company I know, who cares about reproduceability and parallelism, refuses to use FP, and only uses integer.

(2) Ensure that the machine is set up appropriately.

E.g. most machines calculate in 32 or 64 bit precision (C original standard was 64 bit "double" everywhere. But Intel x86/x87 can calculate in 80 bit in registers, and round to 64 or 32 when spilling. 1 shows how to change the x86/x87 precision control from 80 bit to 64 bit, using inline assembly. Note that this code is assembly level and not portable - but most other machines already do computation in 32 or 64 bit precision, and you don't need to worry about the x87 80 bit.

(By the way, on x86 you you can only avoid all issues by using SSE FP; the old legacy Intel x87 FP can never give exactly the same answers (although if you set precision control (PC) to 64 bit rather than 80 bit, you will get the same results except if there was an intermediate overflow, since the exponent width is not affected, just the mantissa))

E.g. ensure that you are using the same underflow mode on all machines. I.e. ensure denorms or enabled, or oppositely that all machines are in flush to zero mode. Here it is a Dobson's choice: flush to zero modes are not standardized, but some machines, e.g. GPUs, simply have not had denormalized numbers. I.e. many machines have IEEE standard number FORMATS, but not actual IEEE standard arithmetic (with denorms). My druther is to require IEEE denorms, but if I were absolutely paranoid I would go with flush to zero, and force that flushing myself in the software.

(3) Ensure that you are using the same language ioptions. Older C programs do all calculations in "double" (64-bit), but it is now permissible to calculate in single precision. Whatever, you want to do it the same way on all machines.

(4) Some smaller items wrt your code:

Avoid big expressions that a compiler is likely to rearrange (if it doesn't implement strict FP switches properly)

Possible write all of your code in simple form like

double a = ...;
double b = ...;
double c = a *b;
double d = ...;
double e = a*d;
double f = c + e;

Rather than

f = (a*b) + (a*c);

which might be optimized to

f = a*(b+c);

I'll leave talking about compiler options for the last, because it is longer.

If you do all of these things, then your calculations should be absolutely repeatable. IEEE floating point is exact - it always gives the same answers. It is the rearranging of the calculations by the compiler on the way to the IEEE FP that introduces variability.

There should be no need for you to round off low order bits. But doing so also will not hurt, and may mask some issues. Remember: you may need to mask off at least one bit for every add....

(2) Compiler optimizations rearranging the code in different ways on different machines. As one commenter said, use whatever your compiler switches for strict FP are.

You might have to disable all optimization for the file containing your sin code.

You might have to use volatiles.

Hopefully there are compiler switches that are more specific. E.g. for gcc:

-ffp-contract=off --- disable fused multiply add, since not all of your target machines may have them.

-fexcess precision=standard --- disables stuff like Intel x86/x87 excess precision in internal registers

-std=c99 --- specifies fairly strict C language standard. Unfortunately not completely implemented, as I google it today

make sure you do not have optimizations enabled like -funsafe-math and -fassociativbe-math