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 <math.h> 29 #include "CubicUtilities.h" 30 #include "QuadraticUtilities.h" 31 #include "QuarticRoot.h" 32 33 int reducedQuarticRoots(const double t4, const double t3, const double t2, const double t1, 34 const double t0, const bool oneHint, double roots[4]) { 35 #ifdef SK_DEBUG 36 // create a string mathematica understands 37 // GDB set print repe 15 # if repeated digits is a bother 38 // set print elements 400 # if line doesn't fit 39 char str[1024]; 40 bzero(str, sizeof(str)); 41 sprintf(str, "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 mathematica_ize(str, sizeof(str)); 44 #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA 45 SkDebugf("%s\n", str); 46 #endif 47 #endif 48 #if 0 && SK_DEBUG 49 bool t4Or = approximately_zero_when_compared_to(t4, t0) // 0 is one root 50 || approximately_zero_when_compared_to(t4, t1) 51 || approximately_zero_when_compared_to(t4, t2); 52 bool t4And = approximately_zero_when_compared_to(t4, t0) // 0 is one root 53 && approximately_zero_when_compared_to(t4, t1) 54 && approximately_zero_when_compared_to(t4, t2); 55 if (t4Or != t4And) { 56 SkDebugf("%s t4 or and\n", __FUNCTION__); 57 } 58 bool t3Or = approximately_zero_when_compared_to(t3, t0) 59 || approximately_zero_when_compared_to(t3, t1) 60 || approximately_zero_when_compared_to(t3, t2); 61 bool t3And = approximately_zero_when_compared_to(t3, t0) 62 && approximately_zero_when_compared_to(t3, t1) 63 && approximately_zero_when_compared_to(t3, t2); 64 if (t3Or != t3And) { 65 SkDebugf("%s t3 or and\n", __FUNCTION__); 66 } 67 bool t0Or = approximately_zero_when_compared_to(t0, t1) // 0 is one root 68 && approximately_zero_when_compared_to(t0, t2) 69 && approximately_zero_when_compared_to(t0, t3) 70 && approximately_zero_when_compared_to(t0, t4); 71 bool t0And = approximately_zero_when_compared_to(t0, t1) // 0 is one root 72 && approximately_zero_when_compared_to(t0, t2) 73 && approximately_zero_when_compared_to(t0, t3) 74 && approximately_zero_when_compared_to(t0, t4); 75 if (t0Or != t0And) { 76 SkDebugf("%s t0 or and\n", __FUNCTION__); 77 } 78 #endif 79 if (approximately_zero_when_compared_to(t4, t0) // 0 is one root 80 && approximately_zero_when_compared_to(t4, t1) 81 && approximately_zero_when_compared_to(t4, t2)) { 82 if (approximately_zero_when_compared_to(t3, t0) 83 && approximately_zero_when_compared_to(t3, t1) 84 && approximately_zero_when_compared_to(t3, t2)) { 85 return quadraticRootsReal(t2, t1, t0, roots); 86 } 87 if (approximately_zero_when_compared_to(t4, t3)) { 88 return cubicRootsReal(t3, t2, t1, t0, roots); 89 } 90 } 91 if ((approximately_zero_when_compared_to(t0, t1) || approximately_zero(t1))// 0 is one root 92 // && approximately_zero_when_compared_to(t0, t2) 93 && approximately_zero_when_compared_to(t0, t3) 94 && approximately_zero_when_compared_to(t0, t4)) { 95 int num = cubicRootsReal(t4, t3, t2, t1, roots); 96 for (int i = 0; i < num; ++i) { 97 if (approximately_zero(roots[i])) { 98 return num; 99 } 100 } 101 roots[num++] = 0; 102 return num; 103 } 104 if (oneHint) { 105 SkASSERT(approximately_zero(t4 + t3 + t2 + t1 + t0)); // 1 is one root 106 int num = cubicRootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); // note that -C==A+B+D+E 107 for (int i = 0; i < num; ++i) { 108 if (approximately_equal(roots[i], 1)) { 109 return num; 110 } 111 } 112 roots[num++] = 1; 113 return num; 114 } 115 return -1; 116 } 117 118 int quarticRootsReal(int firstCubicRoot, const double A, const double B, const double C, 119 const double D, const double E, double s[4]) { 120 double u, v; 121 /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */ 122 const double invA = 1 / A; 123 const double a = B * invA; 124 const double b = C * invA; 125 const double c = D * invA; 126 const double d = E * invA; 127 /* substitute x = y - a/4 to eliminate cubic term: 128 x^4 + px^2 + qx + r = 0 */ 129 const double a2 = a * a; 130 const double p = -3 * a2 / 8 + b; 131 const double q = a2 * a / 8 - a * b / 2 + c; 132 const double r = -3 * a2 * a2 / 256 + a2 * b / 16 - a * c / 4 + d; 133 int num; 134 if (approximately_zero(r)) { 135 /* no absolute term: y(y^3 + py + q) = 0 */ 136 num = cubicRootsReal(1, 0, p, q, s); 137 s[num++] = 0; 138 } else { 139 /* solve the resolvent cubic ... */ 140 double cubicRoots[3]; 141 int roots = cubicRootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots); 142 int index; 143 #if 0 && SK_DEBUG // enable to verify that any cubic root is as good as any other 144 double tries[3][4]; 145 int nums[3]; 146 for (index = 0; index < roots; ++index) { 147 /* ... and take one real solution ... */ 148 const double z = cubicRoots[index]; 149 /* ... to build two quadric equations */ 150 u = z * z - r; 151 v = 2 * z - p; 152 if (approximately_zero_squared(u)) { 153 u = 0; 154 } else if (u > 0) { 155 u = sqrt(u); 156 } else { 157 SkDebugf("%s u=%1.9g <0\n", __FUNCTION__, u); 158 continue; 159 } 160 if (approximately_zero_squared(v)) { 161 v = 0; 162 } else if (v > 0) { 163 v = sqrt(v); 164 } else { 165 SkDebugf("%s v=%1.9g <0\n", __FUNCTION__, v); 166 continue; 167 } 168 nums[index] = quadraticRootsReal(1, q < 0 ? -v : v, z - u, tries[index]); 169 nums[index] += quadraticRootsReal(1, q < 0 ? v : -v, z + u, tries[index] + nums[index]); 170 /* resubstitute */ 171 const double sub = a / 4; 172 for (int i = 0; i < nums[index]; ++i) { 173 tries[index][i] -= sub; 174 } 175 } 176 for (index = 0; index < roots; ++index) { 177 SkDebugf("%s", __FUNCTION__); 178 for (int idx2 = 0; idx2 < nums[index]; ++idx2) { 179 SkDebugf(" %1.9g", tries[index][idx2]); 180 } 181 SkDebugf("\n"); 182 } 183 #endif 184 /* ... and take one real solution ... */ 185 double z; 186 num = 0; 187 int num2 = 0; 188 for (index = firstCubicRoot; index < roots; ++index) { 189 z = cubicRoots[index]; 190 /* ... to build two quadric equations */ 191 u = z * z - r; 192 v = 2 * z - p; 193 if (approximately_zero_squared(u)) { 194 u = 0; 195 } else if (u > 0) { 196 u = sqrt(u); 197 } else { 198 continue; 199 } 200 if (approximately_zero_squared(v)) { 201 v = 0; 202 } else if (v > 0) { 203 v = sqrt(v); 204 } else { 205 continue; 206 } 207 num = quadraticRootsReal(1, q < 0 ? -v : v, z - u, s); 208 num2 = quadraticRootsReal(1, q < 0 ? v : -v, z + u, s + num); 209 if (!((num | num2) & 1)) { 210 break; // prefer solutions without single quad roots 211 } 212 } 213 num += num2; 214 if (!num) { 215 return 0; // no valid cubic root 216 } 217 } 218 /* resubstitute */ 219 const double sub = a / 4; 220 for (int i = 0; i < num; ++i) { 221 s[i] -= sub; 222 } 223 // eliminate duplicates 224 for (int i = 0; i < num - 1; ++i) { 225 for (int j = i + 1; j < num; ) { 226 if (AlmostEqualUlps(s[i], s[j])) { 227 if (j < --num) { 228 s[j] = s[num]; 229 } 230 } else { 231 ++j; 232 } 233 } 234 } 235 return num; 236 } 237
