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) || 75 approximately_zero_when_compared_to(t4 + t3 + t2 + t1 + t0, // 1 is one root 76 SkTMax(fabs(t4), SkTMax(fabs(t3), SkTMax(fabs(t2), SkTMax(fabs(t1), fabs(t0))))))); 77 // note that -C == A + B + D + E 78 int num = SkDCubic::RootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); 79 for (int i = 0; i < num; ++i) { 80 if (approximately_equal(roots[i], 1)) { 81 return num; 82 } 83 } 84 roots[num++] = 1; 85 return num; 86 } 87 return -1; 88 } 89 90 int SkQuarticRootsReal(int firstCubicRoot, const double A, const double B, const double C, 91 const double D, const double E, double s[4]) { 92 double u, v; 93 /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */ 94 const double invA = 1 / A; 95 const double a = B * invA; 96 const double b = C * invA; 97 const double c = D * invA; 98 const double d = E * invA; 99 /* substitute x = y - a/4 to eliminate cubic term: 100 x^4 + px^2 + qx + r = 0 */ 101 const double a2 = a * a; 102 const double p = -3 * a2 / 8 + b; 103 const double q = a2 * a / 8 - a * b / 2 + c; 104 const double r = -3 * a2 * a2 / 256 + a2 * b / 16 - a * c / 4 + d; 105 int num; 106 double largest = SkTMax(fabs(p), fabs(q)); 107 if (approximately_zero_when_compared_to(r, largest)) { 108 /* no absolute term: y(y^3 + py + q) = 0 */ 109 num = SkDCubic::RootsReal(1, 0, p, q, s); 110 s[num++] = 0; 111 } else { 112 /* solve the resolvent cubic ... */ 113 double cubicRoots[3]; 114 int roots = SkDCubic::RootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots); 115 int index; 116 /* ... and take one real solution ... */ 117 double z; 118 num = 0; 119 int num2 = 0; 120 for (index = firstCubicRoot; index < roots; ++index) { 121 z = cubicRoots[index]; 122 /* ... to build two quadric equations */ 123 u = z * z - r; 124 v = 2 * z - p; 125 if (approximately_zero_squared(u)) { 126 u = 0; 127 } else if (u > 0) { 128 u = sqrt(u); 129 } else { 130 continue; 131 } 132 if (approximately_zero_squared(v)) { 133 v = 0; 134 } else if (v > 0) { 135 v = sqrt(v); 136 } else { 137 continue; 138 } 139 num = SkDQuad::RootsReal(1, q < 0 ? -v : v, z - u, s); 140 num2 = SkDQuad::RootsReal(1, q < 0 ? v : -v, z + u, s + num); 141 if (!((num | num2) & 1)) { 142 break; // prefer solutions without single quad roots 143 } 144 } 145 num += num2; 146 if (!num) { 147 return 0; // no valid cubic root 148 } 149 } 150 /* resubstitute */ 151 const double sub = a / 4; 152 for (int i = 0; i < num; ++i) { 153 s[i] -= sub; 154 } 155 // eliminate duplicates 156 for (int i = 0; i < num - 1; ++i) { 157 for (int j = i + 1; j < num; ) { 158 if (AlmostDequalUlps(s[i], s[j])) { 159 if (j < --num) { 160 s[j] = s[num]; 161 } 162 } else { 163 ++j; 164 } 165 } 166 } 167 return num; 168 } 169
