OLD | NEW |
| (Empty) |
1 #include "CubicUtilities.h" | |
2 #include "Intersection_Tests.h" | |
3 #include "QuadraticUtilities.h" | |
4 #include "QuarticRoot.h" | |
5 | |
6 double mulA[] = {-3, -1, 1, 3}; | |
7 size_t mulACount = sizeof(mulA) / sizeof(mulA[0]); | |
8 double rootB[] = {-9, -6, -3, -1, 0, 1, 3, 6, 9}; | |
9 size_t rootBCount = sizeof(rootB) / sizeof(rootB[0]); | |
10 double rootC[] = {-8, -6, -2, -1, 0, 1, 2, 6, 8}; | |
11 size_t rootCCount = sizeof(rootC) / sizeof(rootC[0]); | |
12 double rootD[] = {-7, -4, -1, 0, 1, 2, 5}; | |
13 size_t rootDCount = sizeof(rootD) / sizeof(rootD[0]); | |
14 double rootE[] = {-5, -1, 0, 1, 7}; | |
15 size_t rootECount = sizeof(rootE) / sizeof(rootE[0]); | |
16 | |
17 | |
18 static 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, root
s) | |
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 | |
58 static 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 | |
106 static 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 | |
119 static void testOneQuartic(size_t aIndex, size_t bIndex, size_t cIndex, size_t d
Index, | |
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 | |
171 static 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 | |
187 void 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 } | |
OLD | NEW |