davideberly
commented
1 year ago The inaccurate result when using Rational for angle GTE_C_QUARTER_PI+1e-08 was due to the biquadratic solver, a block of code mixing floating-point arithmetic and rational arithmetic. The floating-point block had a line of code that led to subtractive cancellation of significant bits, the result effectively noise. I rewrote the biquadratic solver and used an algebraic identity to reformulate the offending line of code. Now this particular case produces accurate results when using Rational. I posted the modified file RootsPolynomial.h.
I am working on the root-finder improvements, both in accuracy and speed. Once finished, I have a plan to fix the inaccuracies when using floating-point arithmetic for ellipse-ellipse intersection.
The ellipse-ellipse intersection finding code (IntrEllipse2Ellipse2.h) produced inaccurate intersection points for 2 ellipses centered at the origin and each having extents a = 2 and b = 1. For an angle t, the first ellipse has major axis (cos(t),sin(t)) and minor axis (-sin(t),cos(t)). The second ellipse has major axis (-sin(t),cos(t)) and minor axis (-cos(t),-sin(t)). For t = 1e-8, the double-only code is inaccurate due to computation of the w-value on line 701; the numerator and denominator are nearly zero and approximately the same magnitude, leading to inaccuracy of the (x,y) intersection point. When using Rational type BSRational, the intersection points are accurate.
When t = GTE_C_QUARTER_PI+1e-08, both the double-only and the Rational output is inaccurate. The Rational inaccuracy is due to rounding errors in RootsPolynomial::SolveBiquadratic. This function correctly classifies the real-valued roots, but to compute the roots there are std::sqrt calls where the inputs are Rational, converted to 'double', std::sqrt is computed, and the result is converted back to Rational. The quadratic, cubic and quartic root finders need to be improved (perhaps by an iterative algorithm for root polishing).