gfx/skia/trunk/src/pathops/SkQuarticRoot.cpp@b8a032363ba2 (annotated)
gfx/skia/trunk/src/pathops/SkQuarticRoot.cpp
Thu, 22 Jan 2015 13:21:57 +0100
- author
- Michael Schloh von Bennewitz <michael@schloh.com>
- date
- Thu, 22 Jan 2015 13:21:57 +0100
- branch
- TOR_BUG_9701
- changeset 15
- b8a032363ba2
- permissions
- -rw-r--r--
Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6
| michael@0 | 1 | // from http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c |
| michael@0 | 2 | /* |
| michael@0 | 3 | * Roots3And4.c |
| michael@0 | 4 | * |
| michael@0 | 5 | * Utility functions to find cubic and quartic roots, |
| michael@0 | 6 | * coefficients are passed like this: |
| michael@0 | 7 | * |
| michael@0 | 8 | * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0 |
| michael@0 | 9 | * |
| michael@0 | 10 | * The functions return the number of non-complex roots and |
| michael@0 | 11 | * put the values into the s array. |
| michael@0 | 12 | * |
| michael@0 | 13 | * Author: Jochen Schwarze (schwarze@isa.de) |
| michael@0 | 14 | * |
| michael@0 | 15 | * Jan 26, 1990 Version for Graphics Gems |
| michael@0 | 16 | * Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic |
| michael@0 | 17 | * (reported by Mark Podlipec), |
| michael@0 | 18 | * Old-style function definitions, |
| michael@0 | 19 | * IsZero() as a macro |
| michael@0 | 20 | * Nov 23, 1990 Some systems do not declare acos() and cbrt() in |
| michael@0 | 21 | * <math.h>, though the functions exist in the library. |
| michael@0 | 22 | * If large coefficients are used, EQN_EPS should be |
| michael@0 | 23 | * reduced considerably (e.g. to 1E-30), results will be |
| michael@0 | 24 | * correct but multiple roots might be reported more |
| michael@0 | 25 | * than once. |
| michael@0 | 26 | */ |
| michael@0 | 27 | |
| michael@0 | 28 | #include "SkPathOpsCubic.h" |
| michael@0 | 29 | #include "SkPathOpsQuad.h" |
| michael@0 | 30 | #include "SkQuarticRoot.h" |
| michael@0 | 31 | |
| michael@0 | 32 | int SkReducedQuarticRoots(const double t4, const double t3, const double t2, const double t1, |
| michael@0 | 33 | const double t0, const bool oneHint, double roots[4]) { |
| michael@0 | 34 | #ifdef SK_DEBUG |
| michael@0 | 35 | // create a string mathematica understands |
| michael@0 | 36 | // GDB set print repe 15 # if repeated digits is a bother |
| michael@0 | 37 | // set print elements 400 # if line doesn't fit |
| michael@0 | 38 | char str[1024]; |
| michael@0 | 39 | sk_bzero(str, sizeof(str)); |
| michael@0 | 40 | SK_SNPRINTF(str, sizeof(str), |
| michael@0 | 41 | "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", |
| michael@0 | 42 | t4, t3, t2, t1, t0); |
| michael@0 | 43 | SkPathOpsDebug::MathematicaIze(str, sizeof(str)); |
| michael@0 | 44 | #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA |
| michael@0 | 45 | SkDebugf("%s\n", str); |
| michael@0 | 46 | #endif |
| michael@0 | 47 | #endif |
| michael@0 | 48 | if (approximately_zero_when_compared_to(t4, t0) // 0 is one root |
| michael@0 | 49 | && approximately_zero_when_compared_to(t4, t1) |
| michael@0 | 50 | && approximately_zero_when_compared_to(t4, t2)) { |
| michael@0 | 51 | if (approximately_zero_when_compared_to(t3, t0) |
| michael@0 | 52 | && approximately_zero_when_compared_to(t3, t1) |
| michael@0 | 53 | && approximately_zero_when_compared_to(t3, t2)) { |
| michael@0 | 54 | return SkDQuad::RootsReal(t2, t1, t0, roots); |
| michael@0 | 55 | } |
| michael@0 | 56 | if (approximately_zero_when_compared_to(t4, t3)) { |
| michael@0 | 57 | return SkDCubic::RootsReal(t3, t2, t1, t0, roots); |
| michael@0 | 58 | } |
| michael@0 | 59 | } |
| michael@0 | 60 | if ((approximately_zero_when_compared_to(t0, t1) || approximately_zero(t1)) // 0 is one root |
| michael@0 | 61 | // && approximately_zero_when_compared_to(t0, t2) |
| michael@0 | 62 | && approximately_zero_when_compared_to(t0, t3) |
| michael@0 | 63 | && approximately_zero_when_compared_to(t0, t4)) { |
| michael@0 | 64 | int num = SkDCubic::RootsReal(t4, t3, t2, t1, roots); |
| michael@0 | 65 | for (int i = 0; i < num; ++i) { |
| michael@0 | 66 | if (approximately_zero(roots[i])) { |
| michael@0 | 67 | return num; |
| michael@0 | 68 | } |
| michael@0 | 69 | } |
| michael@0 | 70 | roots[num++] = 0; |
| michael@0 | 71 | return num; |
| michael@0 | 72 | } |
| michael@0 | 73 | if (oneHint) { |
| michael@0 | 74 | SkASSERT(approximately_zero_double(t4 + t3 + t2 + t1 + t0)); // 1 is one root |
| michael@0 | 75 | // note that -C == A + B + D + E |
| michael@0 | 76 | int num = SkDCubic::RootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); |
| michael@0 | 77 | for (int i = 0; i < num; ++i) { |
| michael@0 | 78 | if (approximately_equal(roots[i], 1)) { |
| michael@0 | 79 | return num; |
| michael@0 | 80 | } |
| michael@0 | 81 | } |
| michael@0 | 82 | roots[num++] = 1; |
| michael@0 | 83 | return num; |
| michael@0 | 84 | } |
| michael@0 | 85 | return -1; |
| michael@0 | 86 | } |
| michael@0 | 87 | |
| michael@0 | 88 | int SkQuarticRootsReal(int firstCubicRoot, const double A, const double B, const double C, |
| michael@0 | 89 | const double D, const double E, double s[4]) { |
| michael@0 | 90 | double u, v; |
| michael@0 | 91 | /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */ |
| michael@0 | 92 | const double invA = 1 / A; |
| michael@0 | 93 | const double a = B * invA; |
| michael@0 | 94 | const double b = C * invA; |
| michael@0 | 95 | const double c = D * invA; |
| michael@0 | 96 | const double d = E * invA; |
| michael@0 | 97 | /* substitute x = y - a/4 to eliminate cubic term: |
| michael@0 | 98 | x^4 + px^2 + qx + r = 0 */ |
| michael@0 | 99 | const double a2 = a * a; |
| michael@0 | 100 | const double p = -3 * a2 / 8 + b; |
| michael@0 | 101 | const double q = a2 * a / 8 - a * b / 2 + c; |
| michael@0 | 102 | const double r = -3 * a2 * a2 / 256 + a2 * b / 16 - a * c / 4 + d; |
| michael@0 | 103 | int num; |
| michael@0 | 104 | if (approximately_zero(r)) { |
| michael@0 | 105 | /* no absolute term: y(y^3 + py + q) = 0 */ |
| michael@0 | 106 | num = SkDCubic::RootsReal(1, 0, p, q, s); |
| michael@0 | 107 | s[num++] = 0; |
| michael@0 | 108 | } else { |
| michael@0 | 109 | /* solve the resolvent cubic ... */ |
| michael@0 | 110 | double cubicRoots[3]; |
| michael@0 | 111 | int roots = SkDCubic::RootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots); |
| michael@0 | 112 | int index; |
| michael@0 | 113 | /* ... and take one real solution ... */ |
| michael@0 | 114 | double z; |
| michael@0 | 115 | num = 0; |
| michael@0 | 116 | int num2 = 0; |
| michael@0 | 117 | for (index = firstCubicRoot; index < roots; ++index) { |
| michael@0 | 118 | z = cubicRoots[index]; |
| michael@0 | 119 | /* ... to build two quadric equations */ |
| michael@0 | 120 | u = z * z - r; |
| michael@0 | 121 | v = 2 * z - p; |
| michael@0 | 122 | if (approximately_zero_squared(u)) { |
| michael@0 | 123 | u = 0; |
| michael@0 | 124 | } else if (u > 0) { |
| michael@0 | 125 | u = sqrt(u); |
| michael@0 | 126 | } else { |
| michael@0 | 127 | continue; |
| michael@0 | 128 | } |
| michael@0 | 129 | if (approximately_zero_squared(v)) { |
| michael@0 | 130 | v = 0; |
| michael@0 | 131 | } else if (v > 0) { |
| michael@0 | 132 | v = sqrt(v); |
| michael@0 | 133 | } else { |
| michael@0 | 134 | continue; |
| michael@0 | 135 | } |
| michael@0 | 136 | num = SkDQuad::RootsReal(1, q < 0 ? -v : v, z - u, s); |
| michael@0 | 137 | num2 = SkDQuad::RootsReal(1, q < 0 ? v : -v, z + u, s + num); |
| michael@0 | 138 | if (!((num | num2) & 1)) { |
| michael@0 | 139 | break; // prefer solutions without single quad roots |
| michael@0 | 140 | } |
| michael@0 | 141 | } |
| michael@0 | 142 | num += num2; |
| michael@0 | 143 | if (!num) { |
| michael@0 | 144 | return 0; // no valid cubic root |
| michael@0 | 145 | } |
| michael@0 | 146 | } |
| michael@0 | 147 | /* resubstitute */ |
| michael@0 | 148 | const double sub = a / 4; |
| michael@0 | 149 | for (int i = 0; i < num; ++i) { |
| michael@0 | 150 | s[i] -= sub; |
| michael@0 | 151 | } |
| michael@0 | 152 | // eliminate duplicates |
| michael@0 | 153 | for (int i = 0; i < num - 1; ++i) { |
| michael@0 | 154 | for (int j = i + 1; j < num; ) { |
| michael@0 | 155 | if (AlmostDequalUlps(s[i], s[j])) { |
| michael@0 | 156 | if (j < --num) { |
| michael@0 | 157 | s[j] = s[num]; |
| michael@0 | 158 | } |
| michael@0 | 159 | } else { |
| michael@0 | 160 | ++j; |
| michael@0 | 161 | } |
| michael@0 | 162 | } |
| michael@0 | 163 | } |
| michael@0 | 164 | return num; |
| michael@0 | 165 | } |
