1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/gfx/skia/trunk/src/pathops/SkQuarticRoot.cpp Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,165 @@
1.4 +// from http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c
1.5 +/*
1.6 + * Roots3And4.c
1.7 + *
1.8 + * Utility functions to find cubic and quartic roots,
1.9 + * coefficients are passed like this:
1.10 + *
1.11 + * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
1.12 + *
1.13 + * The functions return the number of non-complex roots and
1.14 + * put the values into the s array.
1.15 + *
1.16 + * Author: Jochen Schwarze (schwarze@isa.de)
1.17 + *
1.18 + * Jan 26, 1990 Version for Graphics Gems
1.19 + * Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic
1.20 + * (reported by Mark Podlipec),
1.21 + * Old-style function definitions,
1.22 + * IsZero() as a macro
1.23 + * Nov 23, 1990 Some systems do not declare acos() and cbrt() in
1.24 + * <math.h>, though the functions exist in the library.
1.25 + * If large coefficients are used, EQN_EPS should be
1.26 + * reduced considerably (e.g. to 1E-30), results will be
1.27 + * correct but multiple roots might be reported more
1.28 + * than once.
1.29 + */
1.30 +
1.31 +#include "SkPathOpsCubic.h"
1.32 +#include "SkPathOpsQuad.h"
1.33 +#include "SkQuarticRoot.h"
1.34 +
1.35 +int SkReducedQuarticRoots(const double t4, const double t3, const double t2, const double t1,
1.36 + const double t0, const bool oneHint, double roots[4]) {
1.37 +#ifdef SK_DEBUG
1.38 + // create a string mathematica understands
1.39 + // GDB set print repe 15 # if repeated digits is a bother
1.40 + // set print elements 400 # if line doesn't fit
1.41 + char str[1024];
1.42 + sk_bzero(str, sizeof(str));
1.43 + SK_SNPRINTF(str, sizeof(str),
1.44 + "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]",
1.45 + t4, t3, t2, t1, t0);
1.46 + SkPathOpsDebug::MathematicaIze(str, sizeof(str));
1.47 +#if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA
1.48 + SkDebugf("%s\n", str);
1.49 +#endif
1.50 +#endif
1.51 + if (approximately_zero_when_compared_to(t4, t0) // 0 is one root
1.52 + && approximately_zero_when_compared_to(t4, t1)
1.53 + && approximately_zero_when_compared_to(t4, t2)) {
1.54 + if (approximately_zero_when_compared_to(t3, t0)
1.55 + && approximately_zero_when_compared_to(t3, t1)
1.56 + && approximately_zero_when_compared_to(t3, t2)) {
1.57 + return SkDQuad::RootsReal(t2, t1, t0, roots);
1.58 + }
1.59 + if (approximately_zero_when_compared_to(t4, t3)) {
1.60 + return SkDCubic::RootsReal(t3, t2, t1, t0, roots);
1.61 + }
1.62 + }
1.63 + if ((approximately_zero_when_compared_to(t0, t1) || approximately_zero(t1)) // 0 is one root
1.64 + // && approximately_zero_when_compared_to(t0, t2)
1.65 + && approximately_zero_when_compared_to(t0, t3)
1.66 + && approximately_zero_when_compared_to(t0, t4)) {
1.67 + int num = SkDCubic::RootsReal(t4, t3, t2, t1, roots);
1.68 + for (int i = 0; i < num; ++i) {
1.69 + if (approximately_zero(roots[i])) {
1.70 + return num;
1.71 + }
1.72 + }
1.73 + roots[num++] = 0;
1.74 + return num;
1.75 + }
1.76 + if (oneHint) {
1.77 + SkASSERT(approximately_zero_double(t4 + t3 + t2 + t1 + t0)); // 1 is one root
1.78 + // note that -C == A + B + D + E
1.79 + int num = SkDCubic::RootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots);
1.80 + for (int i = 0; i < num; ++i) {
1.81 + if (approximately_equal(roots[i], 1)) {
1.82 + return num;
1.83 + }
1.84 + }
1.85 + roots[num++] = 1;
1.86 + return num;
1.87 + }
1.88 + return -1;
1.89 +}
1.90 +
1.91 +int SkQuarticRootsReal(int firstCubicRoot, const double A, const double B, const double C,
1.92 + const double D, const double E, double s[4]) {
1.93 + double u, v;
1.94 + /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
1.95 + const double invA = 1 / A;
1.96 + const double a = B * invA;
1.97 + const double b = C * invA;
1.98 + const double c = D * invA;
1.99 + const double d = E * invA;
1.100 + /* substitute x = y - a/4 to eliminate cubic term:
1.101 + x^4 + px^2 + qx + r = 0 */
1.102 + const double a2 = a * a;
1.103 + const double p = -3 * a2 / 8 + b;
1.104 + const double q = a2 * a / 8 - a * b / 2 + c;
1.105 + const double r = -3 * a2 * a2 / 256 + a2 * b / 16 - a * c / 4 + d;
1.106 + int num;
1.107 + if (approximately_zero(r)) {
1.108 + /* no absolute term: y(y^3 + py + q) = 0 */
1.109 + num = SkDCubic::RootsReal(1, 0, p, q, s);
1.110 + s[num++] = 0;
1.111 + } else {
1.112 + /* solve the resolvent cubic ... */
1.113 + double cubicRoots[3];
1.114 + int roots = SkDCubic::RootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots);
1.115 + int index;
1.116 + /* ... and take one real solution ... */
1.117 + double z;
1.118 + num = 0;
1.119 + int num2 = 0;
1.120 + for (index = firstCubicRoot; index < roots; ++index) {
1.121 + z = cubicRoots[index];
1.122 + /* ... to build two quadric equations */
1.123 + u = z * z - r;
1.124 + v = 2 * z - p;
1.125 + if (approximately_zero_squared(u)) {
1.126 + u = 0;
1.127 + } else if (u > 0) {
1.128 + u = sqrt(u);
1.129 + } else {
1.130 + continue;
1.131 + }
1.132 + if (approximately_zero_squared(v)) {
1.133 + v = 0;
1.134 + } else if (v > 0) {
1.135 + v = sqrt(v);
1.136 + } else {
1.137 + continue;
1.138 + }
1.139 + num = SkDQuad::RootsReal(1, q < 0 ? -v : v, z - u, s);
1.140 + num2 = SkDQuad::RootsReal(1, q < 0 ? v : -v, z + u, s + num);
1.141 + if (!((num | num2) & 1)) {
1.142 + break; // prefer solutions without single quad roots
1.143 + }
1.144 + }
1.145 + num += num2;
1.146 + if (!num) {
1.147 + return 0; // no valid cubic root
1.148 + }
1.149 + }
1.150 + /* resubstitute */
1.151 + const double sub = a / 4;
1.152 + for (int i = 0; i < num; ++i) {
1.153 + s[i] -= sub;
1.154 + }
1.155 + // eliminate duplicates
1.156 + for (int i = 0; i < num - 1; ++i) {
1.157 + for (int j = i + 1; j < num; ) {
1.158 + if (AlmostDequalUlps(s[i], s[j])) {
1.159 + if (j < --num) {
1.160 + s[j] = s[num];
1.161 + }
1.162 + } else {
1.163 + ++j;
1.164 + }
1.165 + }
1.166 + }
1.167 + return num;
1.168 +}
