How to make gcc on SUN calculate floating points the same way as in Linux

Question

I have a project where I have to perform some mathematics calculations with double variables. The problem is that I get different results on SUN Solaris 9 and Linux. There are a lot of ways (explained here and other forums) how to make Linux work as Sun, but not the other way around. I cannot touch the Linux code, so it is only SUN I can change. Is there any way to make SUN to behave as Linux?

The code I run(compile with gcc on both systems):

int hash_func(char *long_id)          
{                                 
    double      product, lnum, gold;
    while (*long_id)
        lnum = lnum * 10.0 + (*long_id++ - '0');
    printf("lnum  => %20.20f\n", lnum);
    lnum = lnum * 10.0E-8;
    printf("lnum  => %20.20f\n", lnum);
    gold = 0.6125423371582974;
    product = lnum * gold;
    printf("product => %20.20f\n", product);
    ...
}

if the input is 339886769243483

the output in Linux:

lnum  => 339886769243**483**.00000000000000000000

lnum  => 33988676.9243**4829473495483398**

product => 20819503.600158**59827399253845**

When on SUN:

lnum => 339886769243483.00000000000000000000

lnum => 33988676.92434830218553543091

product = 20819503.600158**60199928283691**

Note: The result is not always different, moreover most of the times it is the same. Just 10 15-digit numbers out of 60000 have this problem.

Please help!!!

Answer 1

Two things:

• you need to read this: http://docs.sun.com/source/806-3568/ncg_goldberg.html

• you need to decide what your requirements are for numerical precision in your application

Answer 2

These values differ by less than one part in 2⁵². But a double-precision number has just 52 bits behind the radix point. So they differ by just the last bit. It may be that one machine isn't always rounding correctly, but you're doing two multiplications, so the answer is going to have that much error anyway.

Answer 3

The real answer here is another question: why do you think you need this? There may be a better way to accomplish what you're trying to do that doesn't depend on intricate details of platform floating-point. Having said that...

It's unfortunate that you can't change the Linux code, since it's really the Linux results that are deficient here. The SUN results are as good as they could possibly be: they're correctly rounded; each multiplication gives the unique (in this case) C double that's closest to the result. In contrast, the first Linux multiplication does not give a correctly rounded result.

Your Linux results come from a 32-bit system on x86 hardware, right? The results you show are consistent with, and likely caused by, the phenomenon of 'double rounding': the result of the first multiplication is first rounded to 64-bit precision (the precision used internally by the Intel x87 FPU), and then re-rounded to the usual 53-bit precision of a double. Most of the time (around 1999 times out of 2000 or so on average) this double round has the same effect as a single round to 53-bit precision would have had, but occasionally it can produce a different result, and that's what you're seeing here.

As you say, there are ways to fix the Linux results to match the Solaris ones: one of these is to use appropriate compiler flags to force the use of SSE2 instructions for floating-point operations if possible. The recent 4.5 release of gcc also fixes the difference by means of a new -fexcess-precision flag, though the fix may impact performance when not using SSE2.

[Edit: after several rereads of the gcc manuals, the gcc-patches mailing list thread at http://gcc.gnu.org/ml/gcc-patches/2008-11/msg00105.html, and the related gcc bug report, it's still not clear to me whether use of -fexcess-precision=standard does in fact eliminate double rounding on x87 systems; I think the answer depends on the value of FLT_EVAL_METHOD. I don't have a 32-bit Linux/x86 machine handy to test this on.]

But I don't know how you'd fix the Solaris results to match the Linux ones, and I'm not sure why you'd want to: you'd be making the Solaris results less accurate instead of making the Linux results more accurate.

[Edit: caf has a good suggestion here. On Solaris, try deliberately using long double for intermediate results, then forcing back to double. If done right, this should reproduce the double rounding effect that you're seeing in Linux.]

See David Monniaux's excellent paper The pitfalls of verifying floating-point computations for a good explanation of double rounding. It's essential reading after the Goldberg article mentioned in an earlier answer.