caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 1 | #include "CubicUtilities.h" |
| 2 | #include "Intersection_Tests.h" |
caryclark@google.com | 73ca624 | 2013-01-17 21:02:47 +0000 | [diff] [blame] | 3 | #include "QuadraticUtilities.h" |
| 4 | #include "QuarticRoot.h" |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 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 | |
caryclark@google.com | 73ca624 | 2013-01-17 21:02:47 +0000 | [diff] [blame] | 17 | |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 18 | static void quadraticTest(bool limit) { |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 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]; |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 24 | double B = rootB[bIndex]; |
| 25 | double C = rootC[cIndex]; |
| 26 | if (limit) { |
| 27 | B = (B - 6) / 12; |
| 28 | C = (C - 6) / 12; |
| 29 | } |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 30 | const double b = A * (B + C); |
| 31 | const double c = A * B * C; |
| 32 | double roots[2]; |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 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 | } |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 42 | SkASSERT(rootCount == expected); |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 43 | if (!rootCount) { |
| 44 | continue; |
| 45 | } |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 46 | SkASSERT(approximately_equal(roots[0], -B) |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 47 | || approximately_equal(roots[0], -C)); |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 48 | if (expected > 1) { |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 49 | SkASSERT(!approximately_equal(roots[0], roots[1])); |
| 50 | SkASSERT(approximately_equal(roots[1], -B) |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 51 | || approximately_equal(roots[1], -C)); |
| 52 | } |
| 53 | } |
| 54 | } |
| 55 | } |
| 56 | } |
| 57 | |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 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 | } |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 82 | SkASSERT(rootCount == expected); |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 83 | if (!rootCount) { |
| 84 | return; |
| 85 | } |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 86 | SkASSERT(approximately_equal(roots[0], -B) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 87 | || approximately_equal(roots[0], -C) |
| 88 | || approximately_equal(roots[0], -D)); |
| 89 | if (expected <= 1) { |
| 90 | return; |
| 91 | } |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 92 | SkASSERT(!approximately_equal(roots[0], roots[1])); |
| 93 | SkASSERT(approximately_equal(roots[1], -B) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 94 | || approximately_equal(roots[1], -C) |
| 95 | || approximately_equal(roots[1], -D)); |
| 96 | if (expected <= 2) { |
| 97 | return; |
| 98 | } |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 99 | SkASSERT(!approximately_equal(roots[0], roots[2]) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 100 | && !approximately_equal(roots[1], roots[2])); |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 101 | SkASSERT(approximately_equal(roots[2], -B) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 102 | || approximately_equal(roots[2], -C) |
| 103 | || approximately_equal(roots[2], -D)); |
| 104 | } |
| 105 | |
| 106 | static void cubicTest(bool limit) { |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 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) { |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 112 | testOneCubic(limit, aIndex, bIndex, cIndex, dIndex); |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 113 | } |
| 114 | } |
| 115 | } |
| 116 | } |
| 117 | } |
| 118 | |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 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(A, b, c, d, e, roots); |
| 135 | } |
| 136 | const int expected = 1 + (B != C) + (B != D && C != D) + (B != E && C != E && D != E); |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 137 | SkASSERT(rootCount == expected); |
| 138 | SkASSERT(AlmostEqualUlps(roots[0], -B) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 139 | || AlmostEqualUlps(roots[0], -C) |
| 140 | || AlmostEqualUlps(roots[0], -D) |
| 141 | || AlmostEqualUlps(roots[0], -E)); |
| 142 | if (expected <= 1) { |
| 143 | return; |
| 144 | } |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 145 | SkASSERT(!AlmostEqualUlps(roots[0], roots[1])); |
| 146 | SkASSERT(AlmostEqualUlps(roots[1], -B) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 147 | || AlmostEqualUlps(roots[1], -C) |
| 148 | || AlmostEqualUlps(roots[1], -D) |
| 149 | || AlmostEqualUlps(roots[1], -E)); |
| 150 | if (expected <= 2) { |
| 151 | return; |
| 152 | } |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 153 | SkASSERT(!AlmostEqualUlps(roots[0], roots[2]) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 154 | && !AlmostEqualUlps(roots[1], roots[2])); |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 155 | SkASSERT(AlmostEqualUlps(roots[2], -B) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 156 | || AlmostEqualUlps(roots[2], -C) |
| 157 | || AlmostEqualUlps(roots[2], -D) |
| 158 | || AlmostEqualUlps(roots[2], -E)); |
| 159 | if (expected <= 3) { |
| 160 | return; |
| 161 | } |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 162 | SkASSERT(!AlmostEqualUlps(roots[0], roots[3]) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 163 | && !AlmostEqualUlps(roots[1], roots[3]) |
| 164 | && !AlmostEqualUlps(roots[2], roots[3])); |
caryclark@google.com | aa35831 | 2013-01-29 20:28:49 +0000 | [diff] [blame^] | 165 | SkASSERT(AlmostEqualUlps(roots[3], -B) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 166 | || AlmostEqualUlps(roots[3], -C) |
| 167 | || AlmostEqualUlps(roots[3], -D) |
| 168 | || AlmostEqualUlps(roots[3], -E)); |
| 169 | } |
| 170 | |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 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) { |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 179 | testOneQuartic(aIndex, bIndex, cIndex, dIndex, eIndex); |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 180 | } |
| 181 | } |
| 182 | } |
| 183 | } |
| 184 | } |
| 185 | } |
| 186 | |
| 187 | void QuarticRoot_Test() { |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 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); |
caryclark@google.com | 235f56a | 2012-09-14 14:19:30 +0000 | [diff] [blame] | 194 | quarticTest(); |
| 195 | } |