| // from http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c |
| /* |
| * Roots3And4.c |
| * |
| * Utility functions to find cubic and quartic roots, |
| * coefficients are passed like this: |
| * |
| * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0 |
| * |
| * The functions return the number of non-complex roots and |
| * put the values into the s array. |
| * |
| * Author: Jochen Schwarze (schwarze@isa.de) |
| * |
| * Jan 26, 1990 Version for Graphics Gems |
| * Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic |
| * (reported by Mark Podlipec), |
| * Old-style function definitions, |
| * IsZero() as a macro |
| * Nov 23, 1990 Some systems do not declare acos() and cbrt() in |
| * <math.h>, though the functions exist in the library. |
| * If large coefficients are used, EQN_EPS should be |
| * reduced considerably (e.g. to 1E-30), results will be |
| * correct but multiple roots might be reported more |
| * than once. |
| */ |
| |
| #include "SkPathOpsCubic.h" |
| #include "SkPathOpsQuad.h" |
| #include "SkQuarticRoot.h" |
| |
| int SkReducedQuarticRoots(const double t4, const double t3, const double t2, const double t1, |
| const double t0, const bool oneHint, double roots[4]) { |
| #ifdef SK_DEBUG |
| // create a string mathematica understands |
| // GDB set print repe 15 # if repeated digits is a bother |
| // set print elements 400 # if line doesn't fit |
| char str[1024]; |
| sk_bzero(str, sizeof(str)); |
| SK_SNPRINTF(str, sizeof(str), |
| "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", |
| t4, t3, t2, t1, t0); |
| SkPathOpsDebug::MathematicaIze(str, sizeof(str)); |
| #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA |
| SkDebugf("%s\n", str); |
| #endif |
| #endif |
| if (approximately_zero_when_compared_to(t4, t0) // 0 is one root |
| && approximately_zero_when_compared_to(t4, t1) |
| && approximately_zero_when_compared_to(t4, t2)) { |
| if (approximately_zero_when_compared_to(t3, t0) |
| && approximately_zero_when_compared_to(t3, t1) |
| && approximately_zero_when_compared_to(t3, t2)) { |
| return SkDQuad::RootsReal(t2, t1, t0, roots); |
| } |
| if (approximately_zero_when_compared_to(t4, t3)) { |
| return SkDCubic::RootsReal(t3, t2, t1, t0, roots); |
| } |
| } |
| if ((approximately_zero_when_compared_to(t0, t1) || approximately_zero(t1)) // 0 is one root |
| // && approximately_zero_when_compared_to(t0, t2) |
| && approximately_zero_when_compared_to(t0, t3) |
| && approximately_zero_when_compared_to(t0, t4)) { |
| int num = SkDCubic::RootsReal(t4, t3, t2, t1, roots); |
| for (int i = 0; i < num; ++i) { |
| if (approximately_zero(roots[i])) { |
| return num; |
| } |
| } |
| roots[num++] = 0; |
| return num; |
| } |
| if (oneHint) { |
| SkASSERT(approximately_zero_double(t4 + t3 + t2 + t1 + t0) || |
| approximately_zero_when_compared_to(t4 + t3 + t2 + t1 + t0, // 1 is one root |
| SkTMax(fabs(t4), SkTMax(fabs(t3), SkTMax(fabs(t2), SkTMax(fabs(t1), fabs(t0))))))); |
| // note that -C == A + B + D + E |
| int num = SkDCubic::RootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); |
| for (int i = 0; i < num; ++i) { |
| if (approximately_equal(roots[i], 1)) { |
| return num; |
| } |
| } |
| roots[num++] = 1; |
| return num; |
| } |
| return -1; |
| } |
| |
| int SkQuarticRootsReal(int firstCubicRoot, const double A, const double B, const double C, |
| const double D, const double E, double s[4]) { |
| double u, v; |
| /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */ |
| const double invA = 1 / A; |
| const double a = B * invA; |
| const double b = C * invA; |
| const double c = D * invA; |
| const double d = E * invA; |
| /* substitute x = y - a/4 to eliminate cubic term: |
| x^4 + px^2 + qx + r = 0 */ |
| const double a2 = a * a; |
| const double p = -3 * a2 / 8 + b; |
| const double q = a2 * a / 8 - a * b / 2 + c; |
| const double r = -3 * a2 * a2 / 256 + a2 * b / 16 - a * c / 4 + d; |
| int num; |
| double largest = SkTMax(fabs(p), fabs(q)); |
| if (approximately_zero_when_compared_to(r, largest)) { |
| /* no absolute term: y(y^3 + py + q) = 0 */ |
| num = SkDCubic::RootsReal(1, 0, p, q, s); |
| s[num++] = 0; |
| } else { |
| /* solve the resolvent cubic ... */ |
| double cubicRoots[3]; |
| int roots = SkDCubic::RootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots); |
| int index; |
| /* ... and take one real solution ... */ |
| double z; |
| num = 0; |
| int num2 = 0; |
| for (index = firstCubicRoot; index < roots; ++index) { |
| z = cubicRoots[index]; |
| /* ... to build two quadric equations */ |
| u = z * z - r; |
| v = 2 * z - p; |
| if (approximately_zero_squared(u)) { |
| u = 0; |
| } else if (u > 0) { |
| u = sqrt(u); |
| } else { |
| continue; |
| } |
| if (approximately_zero_squared(v)) { |
| v = 0; |
| } else if (v > 0) { |
| v = sqrt(v); |
| } else { |
| continue; |
| } |
| num = SkDQuad::RootsReal(1, q < 0 ? -v : v, z - u, s); |
| num2 = SkDQuad::RootsReal(1, q < 0 ? v : -v, z + u, s + num); |
| if (!((num | num2) & 1)) { |
| break; // prefer solutions without single quad roots |
| } |
| } |
| num += num2; |
| if (!num) { |
| return 0; // no valid cubic root |
| } |
| } |
| /* resubstitute */ |
| const double sub = a / 4; |
| for (int i = 0; i < num; ++i) { |
| s[i] -= sub; |
| } |
| // eliminate duplicates |
| for (int i = 0; i < num - 1; ++i) { |
| for (int j = i + 1; j < num; ) { |
| if (AlmostDequalUlps(s[i], s[j])) { |
| if (j < --num) { |
| s[j] = s[num]; |
| } |
| } else { |
| ++j; |
| } |
| } |
| } |
| return num; |
| } |
