Bug 16825 - need to reimplement QuadCurve2D.solveQuadratic
Summary: need to reimplement QuadCurve2D.solveQuadratic
Status: RESOLVED DUPLICATE of bug 323
Alias: None
Product: classpath
Classification: Unclassified
Component: awt (show other bugs)
Version: unspecified
: P2 normal
Target Milestone: ---
Assignee: Thomas Fitzsimmons
: 22724 (view as bug list)
Depends on:
Reported: 2004-07-29 21:10 UTC by Mark Wielaard
Modified: 2006-08-02 10:56 UTC (History)
46 users (show)

See Also:
Known to work:
Known to fail:
Last reconfirmed: 2005-05-14 21:08:51

Java testcase (193 bytes, text/plain)
2004-11-09 23:17 UTC, Sven de Marothy
C test case (184 bytes, text/plain)
2004-11-09 23:17 UTC, Sven de Marothy
C test case with FPU set to double-precision rounding (265 bytes, text/plain)
2004-11-09 23:18 UTC, Sven de Marothy

Note You need to log in before you can comment on or make changes to this bug.
Description from-classpath 2004-01-07 14:15:10 UTC
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.
Comment 1 from-classpath 2004-03-16 00:08:40 UTC
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.
Comment 2 from-classpath 2004-05-02 20:18:17 UTC
This nice article gives lots of insight into floating-point-things.


Comment 3 from-classpath 2004-05-08 11:59:50 UTC
here is another link to some fp-insights:

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.

Comment 4 from-classpath 2004-05-16 22:17:20 UTC
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. 

Comment 5 from-classpath 2004-05-21 21:41:08 UTC
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
Comment 6 from-classpath 2004-05-21 22:08:51 UTC
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
Comment 7 from-classpath 2004-05-23 09:36:19 UTC
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.

Comment 8 from-classpath 2004-07-07 11:50:47 UTC
Could this bring some clarity to the issue?


Is the FPU in 80-bit precision mode?
Comment 9 Mark Wielaard 2004-07-29 21:10:35 UTC
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:

(Please see the very large audit trail there for more information.)
Comment 10 from-classpath 2004-07-29 21:16:22 UTC
This is now also tracked in the GCC AWT bugzilla as
Comment 11 Andrew Pinski 2004-07-29 21:16:25 UTC
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.
Comment 12 Sven de Marothy 2004-11-09 23:17:30 UTC
Created attachment 7507 [details]
Java testcase
Comment 13 Sven de Marothy 2004-11-09 23:17:53 UTC
Created attachment 7508 [details]
C test case
Comment 14 Sven de Marothy 2004-11-09 23:18:25 UTC
Created attachment 7509 [details]
C test case with FPU set to double-precision rounding
Comment 15 Sven de Marothy 2004-11-09 23:22:54 UTC
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.
Comment 16 Andrew Pinski 2004-11-09 23:32:25 UTC
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).
Comment 17 from-classpath 2004-11-10 22:07:36 UTC
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.
Comment 18 from-classpath 2004-11-11 05:45:49 UTC
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.
Comment 19 from-classpath 2004-11-11 18:23:34 UTC
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.
Comment 20 from-classpath 2004-11-11 19:40:04 UTC
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...
Comment 21 Thomas Fitzsimmons 2005-08-24 19:43:31 UTC
*** Bug 22724 has been marked as a duplicate of this bug. ***
Comment 22 Thomas Fitzsimmons 2005-08-24 23:49:22 UTC
We need to reimplement QuadCurve2D.solveQuadratic, possibly using the algorithm
described here: http://www.library.cornell.edu/nr/bookcpdf/c5-6.pdf
Comment 23 Andrew Haley 2006-08-02 10:56:07 UTC
QuadCurve2D.solveQuadratic already uses the algorithm in Numerical Recipes.

*** This bug has been marked as a duplicate of 323 ***