DXR is a code search and navigation tool aimed at making sense of large projects. It supports full-text and regex searches as well as structural queries.

Header

Mercurial (d8847129d134)

VCS Links

Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
// 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
 *                  <math.h>, 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) ||
                approximately_zero_when_compared_to(t4 + t3 + t2 + t1 + t0,  // 1 is one root
                SkTMax(fabs(t4), SkTMax(fabs(t3), SkTMax(fabs(t2), SkTMax(fabs(t1), fabs(t0)))))));
        // 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;
    double largest = SkTMax(fabs(p), fabs(q));
    if (approximately_zero_when_compared_to(r, largest)) {
    /* 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;
}