Bug 16122

Summary: gij - Incorrect result due to computations in extended precision on x86
Product: gcc Reporter: Debian GCC Maintainers <debian-gcc>
Component: libgcjAssignee: Not yet assigned to anyone <unassigned>
Status: RESOLVED WONTFIX    
Severity: normal CC: gcc-bugs, java-prs, vincent-gcc
Priority: P2    
Version: 3.4.1   
Target Milestone: ---   
Host: Target: i486-linux
Build: Known to work:
Known to fail: Last reconfirmed: 2005-09-24 15:56:22

Description Debian GCC Maintainers 2004-06-21 19:58:29 UTC 
[forwarded from http://bugs.debian.org/255525]

3.3.4 release and 3.4.1 branch, bug submitter writes:

Concerning the following Java source: 
 
------------------------------------------------------------------ 
// $Id: test.java 3734 2004-06-21 15:36:52Z lefevre $ 
 
public class test { 
    public static void main(String[] args) throws Exception { 
        test t = new test(); 
        t.doTest(); 
    } 
 
    volatile double x, y, z, d; 
 
    public void doTest() { 
        x = 9007199254740994.0; /* 2^53 + 2 */ 
        y = 1.0 - 1/65536.0; 
        z = x + y; 
        d = z - x; 
        System.out.println("z = " + z); 
        System.out.println("d = " + d); 
    } 
} 
------------------------------------------------------------------ 
I've compiled it with "gcj -C test.java" (GCC 3.3.4). 
 
Both IBM's and Sun's JVM give the correct result: 
 
greux:~/wd/src/fp> /global/greux/lefevre/IBMJava2-131/jre/bin/java test 
z = 9.007199254740994E15 
d = 0.0 
 
greux:~/wd/src/fp> /usr/local/j2re1.4.1/bin/java test 
z = 9.007199254740994E15 
d = 0.0 
 
but not gij: 
 
greux:~/wd/src/fp> /usr/bin/gij test 
z = 9.007199254740996E15 
d = 2.0 
 
gij should switch the FPU of the x86 processor 
(Pentium III in my case) to rounding in double precision to avoid the 
effect of the "double rounding" (you may find some information about 
this effect here: <http://www.srware.com/linux_numerics.txt>).
Comment 1 Andrew Pinski 2004-06-21 20:24:18 UTC 
Actually this is not a bug unless you use strictfp but IIRC strictfp does nothing to the problem anyways.
Comment 2 Vincent Lefèvre 2004-06-22 06:04:29 UTC 
I don't know very much about Java, but I've talked with someone who has some
knowledge in Java and floating point, and according to him, Java FP arithmetic
operations must round exactly in double precision (this is also confirmed by the
documentation I could see). strictfp concerns the exponent range of the double
precision vs extended precision, but this is not the problem here.
Comment 3 Bryce McKinlay 2004-06-22 16:51:42 UTC 
Confirmed with HEAD. The problem doesn't occur with native compilation, it is a
gij issue only.
Comment 4 Andrew Haley 2006-02-14 15:45:51 UTC 
A bit more explanation.

The problem is caused by the fact that 9007199254740994.0 + 0.9999847412109375 is carried out in extended precision, and the result is rounded to 9007199254740995.  In double precision, the result of this calculation is rounded to 9007199254740994.

When the extended-precision value is rounded to double for storing, it is rounded (again) to 9007199254740996.

So, double rounding leads gij to return a value of d that is 2.0, whereas it should be 0.0.

Note however, that the true accurate value for d, calculated at infinite precision, is 1-(2^-16).  So, the absolute error for gcj is 1+(2^-16) and the absolute error with correct rounding is 1-(2^-16).  (I'm not surprised this hasn't been reported as a problem with any real applications!)

It might be worth setting the floating-point precision of gcj to double, but that would only fix the double-precision case, and I presume we'd still have the same double rounding problem for floats.  

And in any case, I do not know if libcalls would be affected by being entered with the FPU in round-to-double mode.  We might end up breaking things.
Comment 5 Vincent Lefèvre 2006-02-14 17:03:30 UTC 
(In reply to comment #4)
> Note however, that the true accurate value for d, calculated at infinite
> precision, is 1-(2^-16).  So, the absolute error for gcj is 1+(2^-16) and the
> absolute error with correct rounding is 1-(2^-16).  (I'm not surprised this
> hasn't been reported as a problem with any real applications!)

Note that some algorithms may be sensitive to this difference. I give an example in <http://www.vinc17.org/research/publi.en.html#Lef2005b> (The Euclidean division implemented with a floating-point division and a floor); the effect of extended precision is dealt with in Section 5.

A second problem is the reproducilibity of the results on various architectures. Under probabilistic hypotheses, there should be something like 1 case over 2048 that is incorrectly rounded (in the real world, this is much less).

> It might be worth setting the floating-point precision of gcj to double, but
> that would only fix the double-precision case, and I presume we'd still have
> the same double rounding problem for floats.  

Yes, however doubles are nowadays used much more often than floats, IMHO. I think that fixing the problem for the doubles would be sufficient (as it is probably too difficult to do better), though not perfect.

> And in any case, I do not know if libcalls would be affected by being entered
> with the FPU in round-to-double mode.  We might end up breaking things.

The only glibc function for which problems have been noticed is pow in corner cases. See <http://sources.redhat.com/bugzilla/show_bug.cgi?id=706>. And it is also inaccurate when the processor is configured in extended precision; so in any case, users shouldn't rely on it. I'd be interested in other cases, if any.

More information here: <http://www.vinc17.org/research/extended.en.html>
Comment 6 Andrew Pinski 2016-09-30 22:51:02 UTC 
Closing as won't fix as libgcj (and the java front-end) has been removed from the trunk.