The problematic equation has been re-formulated as a test case for the test suite of the GNU Scientific Library. It fails there, too. See the attached e-mail to Brian Gough <bjg@network-theory.co.uk>, author of the GNU Scientific Library.
Brian Gough (bjg@network-theory.co.uk) wrote on 2004-01-08: The coefficient 0.1 does not have an exact representation in binary, so the number used will depend on the precision and the rounding mode. This will cause different results. I suspect that GCJ on Intel uses the 80-bit extended-precision floating point registers (which are the gcc default, in C also) while the other virtual machines probably use strict double-precision. I don't know much about GCJ but I'm sure the developers are aware of the issue, so there may be an option to control it. > Could it be that the test for (disc == 0) should actually be replaced by > something like (fabs(disc - 0.0) < EPSILON)? If so, which value is > appropriate for EPSILON? The best value is 0, as there is no natural choice for epsilon (it introduces an arbitrary constant). My recommendation: choose test coefficients which are exact in binary (e.g. in this case multiply the polynomial by 5) to avoid the problem.
This nice article gives lots of insight into floating-point-things. http://docs.sun.com/source/806-3568/ncg_goldberg.html hth dvholten
here is another link to some fp-insights: http://en.wikipedia.org/wiki/Loss_of_significance IMO, the way to compute and use the discriminant is suspicious at least - in any implementation: to perform some fp-arithmetic, subtract the results and compare to (exact) 0 is a guess, nothing more. I think, a systematic search for single solutions will produce a large number of slightly wrong results. The current way of computing d = b * b - 4 * a * c if ( d == 0 ) ... leaves too much room for different implementations: an aggressive optimizer might use the (extended?) precision value of d from top-of-stack, instead of reloading (rounded) double-precision d. in the case, where b*b and 4*a*c is almost equal, the further use of d ( or sqrt(d) ) is even more inaccurate. so, it seems, a 'real good' algorithm is somewhat sophisticated - or just use this simpler one with a defined epsilon-constant. dvholten
I'm currently hacking this code (see mailing list), and I don't think this is much of a problem. "Floating point operations are like moving around piles of sand, every time you do it, you lose some sand and picke up some dirt." An exact result should never be expected. Anyway, this routine is likely not intended as the end-all of quad solvers, but rather as a utility for calculating bezier intersections, in which case it fits the bill.
A reply to dvholten's comment that "an aggressive optimizer might use the (extended?) precision value": I think this optimization would be a bit too aggressive, as it would be in conflict to the JVM specification. The dsub instruction mandates value set conversion to IEEE 64-bit. See http://java.sun.com/docs/books/vmspec/2nd-edition/html/Instructions2.doc3.html#dsub
Actually, the instruction in question is dcmpl or dcmpg, not dsub. But the dcmpX arguments must equally be converted to the closest value representable in 64-bit IEEE 754: http://java.sun.com/docs/books/vmspec/2nd-edition/html/Instructions2.doc3.html#dcmpop
some more on this... my understanding of (non-) strictfp is, that intermediate results and local variables can be of any precision the implementation feels right. So, it would be completly legal (and not too agressive) to use the intermediate value of d for comparison. With 'agressive optimizer' i originally meant something like java -> native compilers (like gnu java), which may introduce this kind of optimizing in a late code emitting step, when fine-print semantics of the source language are not available. there is another problem in this innocent function: assuming we find a (single) root. what, if some suspicous mind gets the idea to stick this into the original polynom to get a zero-result? how to compute this? y = a*x*x + b*x +c ? or y = (a*x)+b)*x +c ? or the same in a loop? this would give (upto) 3 different results - some of which may even be zero... maybe the whole function is even useless altogether - lets assume we let the values for a, b and c run from -MAX_DOUBLE to +MAX_DOUBLE and and look at all the results. The result-set is large, but it is finite. A large part of it gives the trivial 'no root'-cases - that's simple. Then we have a large part of trivial 'double-root' cases. How large is the remaining part? How much of this is numerically correct? 10^-3 ? 10^-6? 10^-9 ? less ?? Other functions, like sin() or sqrt() are much much better in their outcome. dvholten
Could this bring some clarity to the issue? http://www.srware.com/linux_numerics.txt Is the FPU in 80-bit precision mode?
The method java.awt.geom.QuadCurve2D.solveQuadratic sometimes does not return the right results for the following equation: (x^2)/10 + 20x + 1000 = 0 Expected result: -100 This equation is tested by Mauve; see gnu.testlet.java.awt.geom.QuadCurve2D. The problem occurs with various versions of gcj, compiling for IA-32. However, it does NOT occur when executing the same Mauve testlet on jamvm, using the same implementation for java.awt.geom.QuadCurve2D. It turned out that a small Java test program for computing the equation's discriminant gives a very small negative number when run on gcj/IA-32, but exactly zero when run on jamvm or the Sun J2SE 1.4.1_01. This bug report comes from GNU Classpath: https://savannah.gnu.org/bugs/index.php?func=detailitem&item_id=7123 (Please see the very large audit trail there for more information.)
This is now also tracked in the GCC AWT bugzilla as http://gcc.gnu.org/bugzilla/show_bug.cgi?id=16825
Confirmed, I have seen the discussion of this bug before somewhere but I had forgot where, I think the java-patches or the java list.
Created attachment 7507 [details] Java testcase
Created attachment 7508 [details] C test case
Created attachment 7509 [details] C test case with FPU set to double-precision rounding
This isn't a Awt or libgcj bug itself, but rather an effect of the extended-mode of FPU being the default leading to inconsistent results with 64-bit doubles. A suggestion for a 'fix' is to set the FPU to 64-bit double rounding for Java programs, since Java does not have a 'long double' type. This would both stricter and give more consistent results with other platforms.
Note the java "standard" does not require the FPU set to double-precision except when working in strict mode (this was changed after 1.0 IIRC).
Having done some testing (see comments in GCC bugzilla), I've determined this is not a bug in classpath, but rather a rounding inaccuracy because the processor is in 80-bit precision mode and not using 64-bit rounding. I'm not sure as to whether this is allowed or not, since the API does not specifiy if the method should use strictfp or not.
Regardless whether there's strictfp used or not having FPU in 80-bit mode is not allowed by java specs. There used to be a proposal to introduce a special keyword "whatever math platform has"-like, but it didn't get thru. Instead the "default" requirements have been relaxed (to pretty much exactly match x86 with FPU in 64bit presision mode, at least for doubles) and stricfp keywork added (but none of the free JVMs implements it anyway AFAIK) Can you please try it with SableVM (on x86)? I am 99% sure it'd work the same as Sun's VM. I'll fix it if it didn't.
Gadek: Yes, it works with SableVM, as well as JamVM and gij . (which saves the intermediate results to memory, and thus does the 64-bit rounding) Kaffe and gcj do display this bug though. So if you're certain this is wrong, then this bug should be delegated to them.
AFAIK, the gcj people already know that their fp implementation is violating the JVM spec. Nevertheless, I think we should use an algorithm with better numerical properties, e.g. the one mentioned by dvholten (see the end of http://en.wikipedia.org/wiki/Loss_of_significance). Too bad I don't have the time myself right now...
*** Bug 22724 has been marked as a duplicate of this bug. ***
We need to reimplement QuadCurve2D.solveQuadratic, possibly using the algorithm described here: http://www.library.cornell.edu/nr/bookcpdf/c5-6.pdf
QuadCurve2D.solveQuadratic already uses the algorithm in Numerical Recipes. *** This bug has been marked as a duplicate of 323 ***