blob: 46bb4f50247a2e670db7bbe2833360d9cffb18e8 [file] [log] [blame]
caryclark@google.com235f56a2012-09-14 14:19:30 +00001// 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
skia.committer@gmail.com055c7c22012-09-15 02:01:41 +000017 * (reported by Mark Podlipec),
18 * Old-style function definitions,
19 * IsZero() as a macro
caryclark@google.com235f56a2012-09-14 14:19:30 +000020 * 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"
caryclark@google.com73ca6242013-01-17 21:02:47 +000030#include "QuadraticUtilities.h"
caryclark@google.com235f56a2012-09-14 14:19:30 +000031#include "QuarticRoot.h"
32
caryclark@google.com9f602912013-01-24 21:47:16 +000033int reducedQuarticRoots(const double t4, const double t3, const double t2, const double t1,
34 const double t0, const bool oneHint, double roots[4]) {
mtklein753b8702014-08-20 07:38:46 -070035#ifdef SK_DEBUG
caryclark@google.com73ca6242013-01-17 21:02:47 +000036 // create a string mathematica understands
caryclark@google.com9f602912013-01-24 21:47:16 +000037 // GDB set print repe 15 # if repeated digits is a bother
38 // set print elements 400 # if line doesn't fit
caryclark@google.com73ca6242013-01-17 21:02:47 +000039 char str[1024];
40 bzero(str, sizeof(str));
caryclark@google.com9f602912013-01-24 21:47:16 +000041 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);
caryclark@google.com5e0500f2013-02-20 12:51:37 +000043 mathematica_ize(str, sizeof(str));
caryclark@google.com4aaaaea2013-02-28 16:12:39 +000044#if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA
caryclark@google.com5e0500f2013-02-20 12:51:37 +000045 SkDebugf("%s\n", str);
46#endif
caryclark@google.com73ca6242013-01-17 21:02:47 +000047#endif
caryclark@google.comc83c70e2013-02-22 21:50:07 +000048#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
caryclark@google.combeda3892013-02-07 13:13:41 +000079 if (approximately_zero_when_compared_to(t4, t0) // 0 is one root
80 && approximately_zero_when_compared_to(t4, t1)
caryclark@google.com5e0500f2013-02-20 12:51:37 +000081 && approximately_zero_when_compared_to(t4, t2)) {
caryclark@google.combeda3892013-02-07 13:13:41 +000082 if (approximately_zero_when_compared_to(t3, t0)
83 && approximately_zero_when_compared_to(t3, t1)
84 && approximately_zero_when_compared_to(t3, t2)) {
caryclark@google.com9f602912013-01-24 21:47:16 +000085 return quadraticRootsReal(t2, t1, t0, roots);
86 }
caryclark@google.com5e0500f2013-02-20 12:51:37 +000087 if (approximately_zero_when_compared_to(t4, t3)) {
88 return cubicRootsReal(t3, t2, t1, t0, roots);
89 }
caryclark@google.comd1688742012-09-18 20:08:37 +000090 }
caryclark@google.comc83c70e2013-02-22 21:50:07 +000091 if ((approximately_zero_when_compared_to(t0, t1) || approximately_zero(t1))// 0 is one root
92 // && approximately_zero_when_compared_to(t0, t2)
caryclark@google.comf9502d72013-02-04 14:06:49 +000093 && approximately_zero_when_compared_to(t0, t3)
94 && approximately_zero_when_compared_to(t0, t4)) {
caryclark@google.com9f602912013-01-24 21:47:16 +000095 int num = cubicRootsReal(t4, t3, t2, t1, roots);
caryclark@google.comd1688742012-09-18 20:08:37 +000096 for (int i = 0; i < num; ++i) {
caryclark@google.com9f602912013-01-24 21:47:16 +000097 if (approximately_zero(roots[i])) {
caryclark@google.comd1688742012-09-18 20:08:37 +000098 return num;
99 }
100 }
caryclark@google.com9f602912013-01-24 21:47:16 +0000101 roots[num++] = 0;
caryclark@google.comd1688742012-09-18 20:08:37 +0000102 return num;
103 }
caryclark@google.com9f602912013-01-24 21:47:16 +0000104 if (oneHint) {
caryclark@google.comaa358312013-01-29 20:28:49 +0000105 SkASSERT(approximately_zero(t4 + t3 + t2 + t1 + t0)); // 1 is one root
caryclark@google.com9f602912013-01-24 21:47:16 +0000106 int num = cubicRootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); // note that -C==A+B+D+E
caryclark@google.com6aea33f2012-10-09 14:11:58 +0000107 for (int i = 0; i < num; ++i) {
caryclark@google.com9f602912013-01-24 21:47:16 +0000108 if (approximately_equal(roots[i], 1)) {
caryclark@google.com6aea33f2012-10-09 14:11:58 +0000109 return num;
110 }
111 }
caryclark@google.com9f602912013-01-24 21:47:16 +0000112 roots[num++] = 1;
caryclark@google.com6aea33f2012-10-09 14:11:58 +0000113 return num;
114 }
caryclark@google.com9f602912013-01-24 21:47:16 +0000115 return -1;
caryclark@google.comd1688742012-09-18 20:08:37 +0000116}
caryclark@google.com235f56a2012-09-14 14:19:30 +0000117
caryclark@google.com5e0500f2013-02-20 12:51:37 +0000118int quarticRootsReal(int firstCubicRoot, const double A, const double B, const double C,
119 const double D, const double E, double s[4]) {
caryclark@google.comd1688742012-09-18 20:08:37 +0000120 double u, v;
caryclark@google.com235f56a2012-09-14 14:19:30 +0000121 /* 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:
skia.committer@gmail.com055c7c22012-09-15 02:01:41 +0000128 x^4 + px^2 + qx + r = 0 */
caryclark@google.com235f56a2012-09-14 14:19:30 +0000129 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;
caryclark@google.com73ca6242013-01-17 21:02:47 +0000133 int num;
caryclark@google.com235f56a2012-09-14 14:19:30 +0000134 if (approximately_zero(r)) {
skia.committer@gmail.com055c7c22012-09-15 02:01:41 +0000135 /* no absolute term: y(y^3 + py + q) = 0 */
caryclark@google.com9f602912013-01-24 21:47:16 +0000136 num = cubicRootsReal(1, 0, p, q, s);
caryclark@google.com235f56a2012-09-14 14:19:30 +0000137 s[num++] = 0;
138 } else {
139 /* solve the resolvent cubic ... */
caryclark@google.com9f602912013-01-24 21:47:16 +0000140 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
caryclark@google.com73ca6242013-01-17 21:02:47 +0000184 /* ... and take one real solution ... */
caryclark@google.com9f602912013-01-24 21:47:16 +0000185 double z;
186 num = 0;
187 int num2 = 0;
caryclark@google.com5e0500f2013-02-20 12:51:37 +0000188 for (index = firstCubicRoot; index < roots; ++index) {
caryclark@google.com9f602912013-01-24 21:47:16 +0000189 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 }
caryclark@google.com235f56a2012-09-14 14:19:30 +0000212 }
caryclark@google.com9f602912013-01-24 21:47:16 +0000213 num += num2;
214 if (!num) {
215 return 0; // no valid cubic root
caryclark@google.com235f56a2012-09-14 14:19:30 +0000216 }
caryclark@google.com9f602912013-01-24 21:47:16 +0000217 }
218 /* resubstitute */
219 const double sub = a / 4;
220 for (int i = 0; i < num; ++i) {
221 s[i] -= sub;
caryclark@google.com235f56a2012-09-14 14:19:30 +0000222 }
223 // eliminate duplicates
caryclark@google.com73ca6242013-01-17 21:02:47 +0000224 for (int i = 0; i < num - 1; ++i) {
caryclark@google.com235f56a2012-09-14 14:19:30 +0000225 for (int j = i + 1; j < num; ) {
caryclark@google.com9f602912013-01-24 21:47:16 +0000226 if (AlmostEqualUlps(s[i], s[j])) {
caryclark@google.com235f56a2012-09-14 14:19:30 +0000227 if (j < --num) {
228 s[j] = s[num];
229 }
230 } else {
231 ++j;
232 }
233 }
234 }
caryclark@google.com235f56a2012-09-14 14:19:30 +0000235 return num;
236}