diff options
| -rw-r--r-- | libdimension/polynomial.c | 79 |
1 files changed, 73 insertions, 6 deletions
diff --git a/libdimension/polynomial.c b/libdimension/polynomial.c index 0ca31ab..9d2ac6a 100644 --- a/libdimension/polynomial.c +++ b/libdimension/polynomial.c @@ -320,7 +320,11 @@ static inline size_t dmnsn_solve_linear(const double poly[2], double x[1]) { x[0] = -poly[0]; - return (x[0] >= dmnsn_epsilon) ? 1 : 0; + if (x[0] >= dmnsn_epsilon) { + return 1; + } else { + return 0; + } } /** Solve a normalized quadratic polynomial algebraically. */ @@ -332,12 +336,71 @@ dmnsn_solve_quadratic(const double poly[3], double x[2]) double s = sqrt(disc); x[0] = (-poly[1] + s)/2.0; x[1] = (-poly[1] - s)/2.0; - return (x[0] >= dmnsn_epsilon) ? ((x[1] >= dmnsn_epsilon) ? 2 : 1) : 0; + if (x[1] >= dmnsn_epsilon) { + return 2; + } else if (x[0] >= dmnsn_epsilon) { + return 1; + } else { + return 0; + } } else { return 0; } } +/** Solve a normalized cubic polynomial algebraically. */ +static inline size_t +dmnsn_solve_cubic(double poly[4], double x[3]) +{ + /* Reduce to a monic trinomial (t^3 + p*t + q, t = x + b/3) */ + double b2 = poly[2]*poly[2]; + double p = poly[1] - b2/3.0; + double q = poly[0] - poly[2]*(9.0*poly[1] - 2.0*b2)/27.0; + + double det = 4.0*p*p*p + 27.0*q*q; + double bdiv3 = poly[2]/3; + + if (det <= -dmnsn_epsilon) { + /* Three real roots -- this implies p < 0 */ + double msqrtp3 = -sqrt(-p/3.0); + double theta = acos(3*q/(2*p*msqrtp3))/3.0; + + /* Store the roots in order from largest to smallest */ + x[2] = 2.0*msqrtp3*cos(theta) - bdiv3; + x[0] = -2.0*msqrtp3*cos(4.0*atan(1.0)/3.0 - theta) - bdiv3; + x[1] = -(x[0] + x[2] + poly[2]); + + if (x[2] >= dmnsn_epsilon) { + return 3; + } else if (x[1] >= dmnsn_epsilon) { + return 2; + } else if (x[0] >= dmnsn_epsilon) { + return 1; + } else { + return 0; + } + } else { + /* One real root */ + if (det < 0.0) + det = 0.0; + + double cbrtdet = cbrt(sqrt(det/108.0) + fabs(q)/2.0); + double abst = cbrtdet - p/(3.0*cbrtdet); + + if (q >= 0) { + x[0] = -abst - bdiv3; + } else { + x[0] = abst - bdiv3; + } + + if (x[0] >= dmnsn_epsilon) { + return 1; + } else { + return 0; + } + } +} + /* Uspensky's algorithm */ size_t dmnsn_solve_polynomial(const double poly[], size_t degree, double x[]) @@ -358,10 +421,11 @@ dmnsn_solve_polynomial(const double poly[], size_t degree, double x[]) /* Eliminate simple zero roots */ dmnsn_eliminate_zero_roots(p, °ree); - if (degree >= 3) { - /* Find isolating intervals for (degree - 2) roots of p[] */ - double ranges[degree - 2][2]; - size_t n = dmnsn_uspensky_bounds(p, degree, ranges, degree - 2); + static const size_t max_algebraic = 3; + if (degree > max_algebraic) { + /* Find isolating intervals for (degree - max_algebraic) roots of p[] */ + double ranges[degree - max_algebraic][2]; + size_t n = dmnsn_uspensky_bounds(p, degree, ranges, degree - max_algebraic); /* Calculate a finite upper bound for the roots of p[] */ double absmax = dmnsn_root_bound(p, degree); @@ -389,6 +453,9 @@ dmnsn_solve_polynomial(const double poly[], size_t degree, double x[]) case 2: i += dmnsn_solve_quadratic(p, x + i); break; + case 3: + i += dmnsn_solve_cubic(p, x + i); + break; } return i; |