OLD | NEW |
| (Empty) |
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@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 t
hat -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 d
ouble 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, cubicRo
ots); | |
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[ind
ex]); | |
169 nums[index] += quadraticRootsReal(1, q < 0 ? v : -v, z + u, tries[in
dex] + 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 } | |
OLD | NEW |