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 "SkPathOpsCubic.h"
29 #include "SkPathOpsQuad.h"
30 #include "SkQuarticRoot.h"
31
SkReducedQuarticRoots(const double t4,const double t3,const double t2,const double t1,const double t0,const bool oneHint,double roots[4])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
SkQuarticRootsReal(int firstCubicRoot,const double A,const double B,const double C,const double D,const double E,double s[4])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