diff -r 000000000000 -r 6474c204b198 gfx/skia/trunk/src/pathops/SkQuarticRoot.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/gfx/skia/trunk/src/pathops/SkQuarticRoot.cpp Wed Dec 31 06:09:35 2014 +0100 @@ -0,0 +1,165 @@ +// from http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c +/* + * Roots3And4.c + * + * Utility functions to find cubic and quartic roots, + * coefficients are passed like this: + * + * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0 + * + * The functions return the number of non-complex roots and + * put the values into the s array. + * + * Author: Jochen Schwarze (schwarze@isa.de) + * + * Jan 26, 1990 Version for Graphics Gems + * Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic + * (reported by Mark Podlipec), + * Old-style function definitions, + * IsZero() as a macro + * Nov 23, 1990 Some systems do not declare acos() and cbrt() in + * , though the functions exist in the library. + * If large coefficients are used, EQN_EPS should be + * reduced considerably (e.g. to 1E-30), results will be + * correct but multiple roots might be reported more + * than once. + */ + +#include "SkPathOpsCubic.h" +#include "SkPathOpsQuad.h" +#include "SkQuarticRoot.h" + +int SkReducedQuarticRoots(const double t4, const double t3, const double t2, const double t1, + const double t0, const bool oneHint, double roots[4]) { +#ifdef SK_DEBUG + // create a string mathematica understands + // GDB set print repe 15 # if repeated digits is a bother + // set print elements 400 # if line doesn't fit + char str[1024]; + sk_bzero(str, sizeof(str)); + SK_SNPRINTF(str, sizeof(str), + "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", + t4, t3, t2, t1, t0); + SkPathOpsDebug::MathematicaIze(str, sizeof(str)); +#if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA + SkDebugf("%s\n", str); +#endif +#endif + if (approximately_zero_when_compared_to(t4, t0) // 0 is one root + && approximately_zero_when_compared_to(t4, t1) + && approximately_zero_when_compared_to(t4, t2)) { + if (approximately_zero_when_compared_to(t3, t0) + && approximately_zero_when_compared_to(t3, t1) + && approximately_zero_when_compared_to(t3, t2)) { + return SkDQuad::RootsReal(t2, t1, t0, roots); + } + if (approximately_zero_when_compared_to(t4, t3)) { + return SkDCubic::RootsReal(t3, t2, t1, t0, roots); + } + } + if ((approximately_zero_when_compared_to(t0, t1) || approximately_zero(t1)) // 0 is one root + // && approximately_zero_when_compared_to(t0, t2) + && approximately_zero_when_compared_to(t0, t3) + && approximately_zero_when_compared_to(t0, t4)) { + int num = SkDCubic::RootsReal(t4, t3, t2, t1, roots); + for (int i = 0; i < num; ++i) { + if (approximately_zero(roots[i])) { + return num; + } + } + roots[num++] = 0; + return num; + } + if (oneHint) { + SkASSERT(approximately_zero_double(t4 + t3 + t2 + t1 + t0)); // 1 is one root + // note that -C == A + B + D + E + int num = SkDCubic::RootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); + for (int i = 0; i < num; ++i) { + if (approximately_equal(roots[i], 1)) { + return num; + } + } + roots[num++] = 1; + return num; + } + return -1; +} + +int SkQuarticRootsReal(int firstCubicRoot, const double A, const double B, const double C, + const double D, const double E, double s[4]) { + double u, v; + /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */ + const double invA = 1 / A; + const double a = B * invA; + const double b = C * invA; + const double c = D * invA; + const double d = E * invA; + /* substitute x = y - a/4 to eliminate cubic term: + x^4 + px^2 + qx + r = 0 */ + const double a2 = a * a; + const double p = -3 * a2 / 8 + b; + const double q = a2 * a / 8 - a * b / 2 + c; + const double r = -3 * a2 * a2 / 256 + a2 * b / 16 - a * c / 4 + d; + int num; + if (approximately_zero(r)) { + /* no absolute term: y(y^3 + py + q) = 0 */ + num = SkDCubic::RootsReal(1, 0, p, q, s); + s[num++] = 0; + } else { + /* solve the resolvent cubic ... */ + double cubicRoots[3]; + int roots = SkDCubic::RootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots); + int index; + /* ... and take one real solution ... */ + double z; + num = 0; + int num2 = 0; + for (index = firstCubicRoot; index < roots; ++index) { + z = cubicRoots[index]; + /* ... to build two quadric equations */ + u = z * z - r; + v = 2 * z - p; + if (approximately_zero_squared(u)) { + u = 0; + } else if (u > 0) { + u = sqrt(u); + } else { + continue; + } + if (approximately_zero_squared(v)) { + v = 0; + } else if (v > 0) { + v = sqrt(v); + } else { + continue; + } + num = SkDQuad::RootsReal(1, q < 0 ? -v : v, z - u, s); + num2 = SkDQuad::RootsReal(1, q < 0 ? v : -v, z + u, s + num); + if (!((num | num2) & 1)) { + break; // prefer solutions without single quad roots + } + } + num += num2; + if (!num) { + return 0; // no valid cubic root + } + } + /* resubstitute */ + const double sub = a / 4; + for (int i = 0; i < num; ++i) { + s[i] -= sub; + } + // eliminate duplicates + for (int i = 0; i < num - 1; ++i) { + for (int j = i + 1; j < num; ) { + if (AlmostDequalUlps(s[i], s[j])) { + if (j < --num) { + s[j] = s[num]; + } + } else { + ++j; + } + } + } + return num; +}