remove SkFP.h and replace SkFP with SkScalar stop respecting SK_SOFTWARE_FLOAT, assume its always f…

src/core/SkGeometry.cpp

 
 /*
  * Copyright 2006 The Android Open Source Project
  *
  * Use of this source code is governed by a BSD-style license that can be
  * found in the LICENSE file.
  */
 
 
 #include "SkGeometry.h"
@@ skipping 781 matching lines @@
     for (int i = count - 1; i > 0; --i)
         for (int j = i; j > 0; --j)
             if (array[j] < array[j-1])
             {
                 T   tmp(array[j]);
                 array[j] = array[j-1];
                 array[j-1] = tmp;
             }
 }
 
-#include "SkFP.h"
-
 // newton refinement
 #if 0
 static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root)
 {
     //  x1 = x0 - f(t) / f'(t)
 
     SkFP    T = SkScalarToFloat(root);
     SkFP    N, D;
 
     // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2]
@@ skipping 76 matching lines @@
             SkASSERT(dst[j-1] < dst[j]);
         }
     }
 }
 #endif
 
 #if defined _WIN32 && _MSC_VER >= 1300  && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop
 #pragma warning ( disable : 4702 )
 #endif
 
+static SkScalar SkScalarCubeRoot(SkScalar x) {
+    return sk_float_pow(x, 0.3333333f);
+}

[caryclark, 2013/12/16 13:39:31]

Here's what I came up with for float in my experiments (I only use double in
pathops)

// cube root approximation using bit hack for 32-bit float
// adapted from Kahan's cbrt
static float cbrt_5f(float f)
{
    unsigned int* p = (unsigned int *) &f;
    *p = *p / 3 + 709921077;
    return f;
}

// iterative cube root approximation using Halley's method
static float cbrta_halleyf(const float a, const float R)
{
    const float a3 = a * a * a;
    const float b = a * (a3 + R + R) / (a3 + a3 + R);
    return b;
}

// cube root approximation using 1 iteration of Halley's method
static float halley_cbrt1f(float d)
{
    float a = cbrt_5f(d);
    return cbrta_halleyf(a, d);
}

In my tests, this was sufficiently precise. If my tests are bogus, the precision
can be improved with another iteration:

// cube root approximation using 2 iterations of Halley's method
static float halley_cbrt2f(float d)
{
    float a = cbrt_5f(d);
    a = cbrta_halleyf(a, d);
    return cbrta_halleyf(a, d);
}

+
 /*  Solve coeff(t) == 0, returning the number of roots that
     lie withing 0 < t < 1.
     coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3]
 
     Eliminates repeated roots (so that all tValues are distinct, and are always
     in increasing order.
 */
-static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3])
+static int solve_cubic_polynomial(const SkScalar coeff[4], SkScalar tValues[3])
 {
-#ifndef SK_SCALAR_IS_FLOAT
-    return 0;   // this is not yet implemented for software float
-#endif
-
     if (SkScalarNearlyZero(coeff[0]))   // we're just a quadratic
     {
         return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues);
     }
 
-    SkFP    a, b, c, Q, R;
+    SkScalar a, b, c, Q, R;
 
     {
         SkASSERT(coeff[0] != 0);
 
-        SkFP inva = SkFPInvert(coeff[0]);
-        a = SkFPMul(coeff[1], inva);
-        b = SkFPMul(coeff[2], inva);
-        c = SkFPMul(coeff[3], inva);
+        SkScalar inva = SkScalarInvert(coeff[0]);
+        a = coeff[1] * inva;
+        b = coeff[2] * inva;
+        c = coeff[3] * inva;
     }
-    Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9);
-//  R = (2*a*a*a - 9*a*b + 27*c) / 54;
-    R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2);
-    R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9));
-    R = SkFPAdd(R, SkFPMulInt(c, 27));
-    R = SkFPDivInt(R, 54);
+    Q = (a*a - b*3) / 9;
+    R = (2*a*a*a - 9*a*b + 27*c) / 54;
 
-    SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q);
-    SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3);
-    SkFP adiv3 = SkFPDivInt(a, 3);
+    SkScalar Q3 = Q * Q * Q;
+    SkScalar R2MinusQ3 = R * R - Q3;
+    SkScalar adiv3 = a / 3;
 
     SkScalar*   roots = tValues;
     SkScalar    r;
 
-    if (SkFPLT(R2MinusQ3, 0))   // we have 3 real roots
+    if (R2MinusQ3 < 0)   // we have 3 real roots
     {
-#ifdef SK_SCALAR_IS_FLOAT
         float theta = sk_float_acos(R / sk_float_sqrt(Q3));
         float neg2RootQ = -2 * sk_float_sqrt(Q);
 
         r = neg2RootQ * sk_float_cos(theta/3) - adiv3;
         if (is_unit_interval(r))
             *roots++ = r;
 
         r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3;
         if (is_unit_interval(r))
             *roots++ = r;
 
         r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3;
         if (is_unit_interval(r))
             *roots++ = r;
 
         SkDEBUGCODE(test_collaps_duplicates();)
 
         // now sort the roots
         int count = (int)(roots - tValues);
         SkASSERT((unsigned)count <= 3);
         bubble_sort(tValues, count);
         count = collaps_duplicates(tValues, count);
         roots = tValues + count;    // so we compute the proper count below
-#endif
     }
     else                // we have 1 real root
     {
-        SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3));
-        A = SkFPCubeRoot(A);
-        if (SkFPGT(R, 0))
-            A = SkFPNeg(A);
+        SkScalar A = SkScalarAbs(R) + SkScalarSqrt(R2MinusQ3);
+        A = SkScalarCubeRoot(A);
+        if (R > 0)
+            A = -A;
 
         if (A != 0)
-            A = SkFPAdd(A, SkFPDiv(Q, A));
-        r = SkFPToScalar(SkFPSub(A, adiv3));
+            A += Q / A;
+        r = A - adiv3;
         if (is_unit_interval(r))
             *roots++ = r;
     }
 
     return (int)(roots - tValues);
 }
 
 /*  Looking for F' dot F'' == 0
 
     A = b - a
     B = c - 2b + a
     C = d - 3c + 3b - a
 
     F' = 3Ct^2 + 6Bt + 3A
     F'' = 6Ct + 6B
 
     F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
 */
-static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4])
+static void formulate_F1DotF2(const SkScalar src[], SkScalar coeff[4])
 {
     SkScalar    a = src[2] - src[0];
     SkScalar    b = src[4] - 2 * src[2] + src[0];
     SkScalar    c = src[6] + 3 * (src[2] - src[4]) - src[0];
 
-    SkFP    A = SkScalarToFP(a);
-    SkFP    B = SkScalarToFP(b);
-    SkFP    C = SkScalarToFP(c);
-
-    coeff[0] = SkFPMul(C, C);
-    coeff[1] = SkFPMulInt(SkFPMul(B, C), 3);
-    coeff[2] = SkFPMulInt(SkFPMul(B, B), 2);
-    coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A));
-    coeff[3] = SkFPMul(A, B);
+    coeff[0] = c * c;
+    coeff[1] = 3 * b * c;
+    coeff[2] = 2 * b * b + c * a;
+    coeff[3] = a * b;
 }
 
 // EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1
 //#define kMinTValueForChopping (SK_Scalar1 / 256)
 #define kMinTValueForChopping   0
 
 /*  Looking for F' dot F'' == 0
 
     A = b - a
     B = c - 2b + a
     C = d - 3c + 3b - a
 
     F' = 3Ct^2 + 6Bt + 3A
     F'' = 6Ct + 6B
 
     F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
 */
 int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3])
 {
-    SkFP    coeffX[4], coeffY[4];
-    int     i;
+    SkScalar coeffX[4], coeffY[4];
+    int      i;
 
     formulate_F1DotF2(&src[0].fX, coeffX);
     formulate_F1DotF2(&src[0].fY, coeffY);
 
     for (i = 0; i < 4; i++)
-        coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]);
+        coeffX[i] += coeffY[i];
 
     SkScalar    t[3];
     int         count = solve_cubic_polynomial(coeffX, t);
     int         maxCount = 0;
 
     // now remove extrema where the curvature is zero (mins)
     // !!!! need a test for this !!!!
     for (i = 0; i < count; i++)
     {
         // if (not_min_curvature())
@@ skipping 591 matching lines @@
     }
     if (this->findYExtrema(&t)) {
         this->evalAt(t, &pts[count++]);
     }
     bounds->set(pts, count);
 }
 
 void SkConic::computeFastBounds(SkRect* bounds) const {
     bounds->set(fPts, 3);
 }