1#include "CubicUtilities.h"
2#include "Intersection_Tests.h"
3#include "QuadraticUtilities.h"
4#include "QuarticRoot.h"
5
6double mulA[] = {-3, -1, 1, 3};
7size_t mulACount = sizeof(mulA) / sizeof(mulA[0]);
8double rootB[] = {-9, -6, -3, -1, 0, 1, 3, 6, 9};
9size_t rootBCount = sizeof(rootB) / sizeof(rootB[0]);
10double rootC[] = {-8, -6, -2, -1, 0, 1, 2, 6, 8};
11size_t rootCCount = sizeof(rootC) / sizeof(rootC[0]);
12double rootD[] = {-7, -4, -1, 0, 1, 2, 5};
13size_t rootDCount = sizeof(rootD) / sizeof(rootD[0]);
14double rootE[] = {-5, -1, 0, 1, 7};
15size_t rootECount = sizeof(rootE) / sizeof(rootE[0]);
16
17
18static void quadraticTest(bool limit) {
19    // (x - a)(x - b) == x^2 - (a + b)x + ab
20    for (size_t aIndex = 0; aIndex < mulACount; ++aIndex) {
21        for (size_t bIndex = 0; bIndex < rootBCount; ++bIndex) {
22            for (size_t cIndex = 0; cIndex < rootCCount; ++cIndex) {
23                const double A = mulA[aIndex];
24                double B = rootB[bIndex];
25                double C = rootC[cIndex];
26                if (limit) {
27                    B = (B - 6) / 12;
28                    C = (C - 6) / 12;
29                }
30                const double b = A * (B + C);
31                const double c = A * B * C;
32                double roots[2];
33                const int rootCount = limit ? quadraticRootsValidT(A, b, c, roots)
34                    : quadraticRootsReal(A, b, c, roots);
35                int expected;
36                if (limit) {
37                    expected = B <= 0 && B >= -1;
38                    expected += B != C && C <= 0 && C >= -1;
39                } else {
40                    expected = 1 + (B != C);
41                }
42                SkASSERT(rootCount == expected);
43                if (!rootCount) {
44                    continue;
45                }
46                SkASSERT(approximately_equal(roots[0], -B)
47                        || approximately_equal(roots[0], -C));
48                if (expected > 1) {
49                    SkASSERT(!approximately_equal(roots[0], roots[1]));
50                    SkASSERT(approximately_equal(roots[1], -B)
51                            || approximately_equal(roots[1], -C));
52                }
53            }
54        }
55    }
56}
57
58static void testOneCubic(bool limit, size_t aIndex, size_t bIndex, size_t cIndex, size_t dIndex) {
59    const double A = mulA[aIndex];
60    double B = rootB[bIndex];
61    double C = rootC[cIndex];
62    double D = rootD[dIndex];
63    if (limit) {
64        B = (B - 6) / 12;
65        C = (C - 6) / 12;
66        D = (C - 2) / 6;
67    }
68    const double b = A * (B + C + D);
69    const double c = A * (B * C + C * D + B * D);
70    const double d = A * B * C * D;
71    double roots[3];
72    const int rootCount = limit ? cubicRootsValidT(A, b, c, d, roots)
73            : cubicRootsReal(A, b, c, d, roots);
74    int expected;
75    if (limit) {
76        expected = B <= 0 && B >= -1;
77        expected += B != C && C <= 0 && C >= -1;
78        expected += B != D && C != D && D <= 0 && D >= -1;
79    } else {
80        expected = 1 + (B != C) + (B != D && C != D);
81    }
82    SkASSERT(rootCount == expected);
83    if (!rootCount) {
84        return;
85    }
86    SkASSERT(approximately_equal(roots[0], -B)
87            || approximately_equal(roots[0], -C)
88            || approximately_equal(roots[0], -D));
89    if (expected <= 1) {
90        return;
91    }
92    SkASSERT(!approximately_equal(roots[0], roots[1]));
93    SkASSERT(approximately_equal(roots[1], -B)
94            || approximately_equal(roots[1], -C)
95            || approximately_equal(roots[1], -D));
96    if (expected <= 2) {
97        return;
98    }
99    SkASSERT(!approximately_equal(roots[0], roots[2])
100            && !approximately_equal(roots[1], roots[2]));
101    SkASSERT(approximately_equal(roots[2], -B)
102            || approximately_equal(roots[2], -C)
103            || approximately_equal(roots[2], -D));
104}
105
106static void cubicTest(bool limit) {
107    // (x - a)(x - b)(x - c) == x^3 - (a + b + c)x^2 + (ab + bc + ac)x - abc
108    for (size_t aIndex = 0; aIndex < mulACount; ++aIndex) {
109        for (size_t bIndex = 0; bIndex < rootBCount; ++bIndex) {
110            for (size_t cIndex = 0; cIndex < rootCCount; ++cIndex) {
111                for (size_t dIndex = 0; dIndex < rootDCount; ++dIndex) {
112                    testOneCubic(limit, aIndex, bIndex, cIndex, dIndex);
113                }
114            }
115        }
116    }
117}
118
119static void testOneQuartic(size_t aIndex, size_t bIndex, size_t cIndex, size_t dIndex,
120        size_t eIndex) {
121    const double A = mulA[aIndex];
122    const double B = rootB[bIndex];
123    const double C = rootC[cIndex];
124    const double D = rootD[dIndex];
125    const double E = rootE[eIndex];
126    const double b = A * (B + C + D + E);
127    const double c = A * (B * C + C * D + B * D + B * E + C * E + D * E);
128    const double d = A * (B * C * D + B * C * E + B * D * E + C * D * E);
129    const double e = A * B * C * D * E;
130    double roots[4];
131    bool oneHint = approximately_zero(A + b + c + d + e);
132    int rootCount = reducedQuarticRoots(A, b, c, d, e, oneHint, roots);
133    if (rootCount < 0) {
134        rootCount = quarticRootsReal(0, A, b, c, d, e, roots);
135    }
136    const int expected = 1 + (B != C) + (B != D && C != D) + (B != E && C != E && D != E);
137    SkASSERT(rootCount == expected);
138    SkASSERT(AlmostEqualUlps(roots[0], -B)
139            || AlmostEqualUlps(roots[0], -C)
140            || AlmostEqualUlps(roots[0], -D)
141            || AlmostEqualUlps(roots[0], -E));
142    if (expected <= 1) {
143        return;
144    }
145    SkASSERT(!AlmostEqualUlps(roots[0], roots[1]));
146    SkASSERT(AlmostEqualUlps(roots[1], -B)
147            || AlmostEqualUlps(roots[1], -C)
148            || AlmostEqualUlps(roots[1], -D)
149            || AlmostEqualUlps(roots[1], -E));
150    if (expected <= 2) {
151        return;
152    }
153    SkASSERT(!AlmostEqualUlps(roots[0], roots[2])
154            && !AlmostEqualUlps(roots[1], roots[2]));
155    SkASSERT(AlmostEqualUlps(roots[2], -B)
156            || AlmostEqualUlps(roots[2], -C)
157            || AlmostEqualUlps(roots[2], -D)
158            || AlmostEqualUlps(roots[2], -E));
159    if (expected <= 3) {
160        return;
161    }
162    SkASSERT(!AlmostEqualUlps(roots[0], roots[3])
163            && !AlmostEqualUlps(roots[1], roots[3])
164            && !AlmostEqualUlps(roots[2], roots[3]));
165    SkASSERT(AlmostEqualUlps(roots[3], -B)
166            || AlmostEqualUlps(roots[3], -C)
167            || AlmostEqualUlps(roots[3], -D)
168            || AlmostEqualUlps(roots[3], -E));
169}
170
171static void quarticTest() {
172    // (x - a)(x - b)(x - c)(x - d) == x^4 - (a + b + c + d)x^3
173    //   + (ab + bc + cd + ac + bd + cd)x^2 - (abc + bcd + abd + acd) * x + abcd
174    for (size_t aIndex = 0; aIndex < mulACount; ++aIndex) {
175        for (size_t bIndex = 0; bIndex < rootBCount; ++bIndex) {
176            for (size_t cIndex = 0; cIndex < rootCCount; ++cIndex) {
177                for (size_t dIndex = 0; dIndex < rootDCount; ++dIndex) {
178                    for (size_t eIndex = 0; eIndex < rootECount; ++eIndex) {
179                        testOneQuartic(aIndex, bIndex, cIndex, dIndex, eIndex);
180                    }
181                }
182            }
183        }
184    }
185}
186
187void QuarticRoot_Test() {
188    testOneCubic(false, 0, 5, 5, 4);
189    testOneQuartic(0, 0, 2, 4, 3);
190    quadraticTest(true);
191    quadraticTest(false);
192    cubicTest(true);
193    cubicTest(false);
194    quarticTest();
195}
196