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, 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 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 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 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 } 196
