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
32int 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
90int 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