1 // from http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c 2 /* 3 * Roots3And4.c 4 * 5 * Utility functions to find cubic and quartic roots, 6 * coefficients are passed like this: 7 * 8 * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0 9 * 10 * The functions return the number of non-complex roots and 11 * put the values into the s array. 12 * 13 * Author: Jochen Schwarze (schwarze (at) isa.de) 14 * 15 * Jan 26, 1990 Version for Graphics Gems 16 * Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic 17 * (reported by Mark Podlipec), 18 * Old-style function definitions, 19 * IsZero() as a macro 20 * Nov 23, 1990 Some systems do not declare acos() and cbrt() in 21 * <math.h>, though the functions exist in the library. 22 * If large coefficients are used, EQN_EPS should be 23 * reduced considerably (e.g. to 1E-30), results will be 24 * correct but multiple roots might be reported more 25 * than once. 26 */ 27 28 #include "SkPathOpsCubic.h" 29 #include "SkPathOpsQuad.h" 30 #include "SkQuarticRoot.h" 31 32 int SkReducedQuarticRoots(const double t4, const double t3, const double t2, const double t1, 33 const double t0, const bool oneHint, double roots[4]) { 34 #ifdef SK_DEBUG 35 // create a string mathematica understands 36 // GDB set print repe 15 # if repeated digits is a bother 37 // set print elements 400 # if line doesn't fit 38 char str[1024]; 39 sk_bzero(str, sizeof(str)); 40 SK_SNPRINTF(str, sizeof(str), 41 "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", 42 t4, t3, t2, t1, t0); 43 SkPathOpsDebug::MathematicaIze(str, sizeof(str)); 44 #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA 45 SkDebugf("%s\n", str); 46 #endif 47 #endif 48 if (approximately_zero_when_compared_to(t4, t0) // 0 is one root 49 && approximately_zero_when_compared_to(t4, t1) 50 && approximately_zero_when_compared_to(t4, t2)) { 51 if (approximately_zero_when_compared_to(t3, t0) 52 && approximately_zero_when_compared_to(t3, t1) 53 && approximately_zero_when_compared_to(t3, t2)) { 54 return SkDQuad::RootsReal(t2, t1, t0, roots); 55 } 56 if (approximately_zero_when_compared_to(t4, t3)) { 57 return SkDCubic::RootsReal(t3, t2, t1, t0, roots); 58 } 59 } 60 if ((approximately_zero_when_compared_to(t0, t1) || approximately_zero(t1)) // 0 is one root 61 // && approximately_zero_when_compared_to(t0, t2) 62 && approximately_zero_when_compared_to(t0, t3) 63 && approximately_zero_when_compared_to(t0, t4)) { 64 int num = SkDCubic::RootsReal(t4, t3, t2, t1, roots); 65 for (int i = 0; i < num; ++i) { 66 if (approximately_zero(roots[i])) { 67 return num; 68 } 69 } 70 roots[num++] = 0; 71 return num; 72 } 73 if (oneHint) { 74 SkASSERT(approximately_zero_double(t4 + t3 + t2 + t1 + t0)); // 1 is one root 75 // note that -C == A + B + D + E 76 int num = SkDCubic::RootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); 77 for (int i = 0; i < num; ++i) { 78 if (approximately_equal(roots[i], 1)) { 79 return num; 80 } 81 } 82 roots[num++] = 1; 83 return num; 84 } 85 return -1; 86 } 87 88 int SkQuarticRootsReal(int firstCubicRoot, const double A, const double B, const double C, 89 const double D, const double E, double s[4]) { 90 double u, v; 91 /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */ 92 const double invA = 1 / A; 93 const double a = B * invA; 94 const double b = C * invA; 95 const double c = D * invA; 96 const double d = E * invA; 97 /* substitute x = y - a/4 to eliminate cubic term: 98 x^4 + px^2 + qx + r = 0 */ 99 const double a2 = a * a; 100 const double p = -3 * a2 / 8 + b; 101 const double q = a2 * a / 8 - a * b / 2 + c; 102 const double r = -3 * a2 * a2 / 256 + a2 * b / 16 - a * c / 4 + d; 103 int num; 104 if (approximately_zero(r)) { 105 /* no absolute term: y(y^3 + py + q) = 0 */ 106 num = SkDCubic::RootsReal(1, 0, p, q, s); 107 s[num++] = 0; 108 } else { 109 /* solve the resolvent cubic ... */ 110 double cubicRoots[3]; 111 int roots = SkDCubic::RootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots); 112 int index; 113 /* ... and take one real solution ... */ 114 double z; 115 num = 0; 116 int num2 = 0; 117 for (index = firstCubicRoot; index < roots; ++index) { 118 z = cubicRoots[index]; 119 /* ... to build two quadric equations */ 120 u = z * z - r; 121 v = 2 * z - p; 122 if (approximately_zero_squared(u)) { 123 u = 0; 124 } else if (u > 0) { 125 u = sqrt(u); 126 } else { 127 continue; 128 } 129 if (approximately_zero_squared(v)) { 130 v = 0; 131 } else if (v > 0) { 132 v = sqrt(v); 133 } else { 134 continue; 135 } 136 num = SkDQuad::RootsReal(1, q < 0 ? -v : v, z - u, s); 137 num2 = SkDQuad::RootsReal(1, q < 0 ? v : -v, z + u, s + num); 138 if (!((num | num2) & 1)) { 139 break; // prefer solutions without single quad roots 140 } 141 } 142 num += num2; 143 if (!num) { 144 return 0; // no valid cubic root 145 } 146 } 147 /* resubstitute */ 148 const double sub = a / 4; 149 for (int i = 0; i < num; ++i) { 150 s[i] -= sub; 151 } 152 // eliminate duplicates 153 for (int i = 0; i < num - 1; ++i) { 154 for (int j = i + 1; j < num; ) { 155 if (AlmostDequalUlps(s[i], s[j])) { 156 if (j < --num) { 157 s[j] = s[num]; 158 } 159 } else { 160 ++j; 161 } 162 } 163 } 164 return num; 165 } 166
