OLD | NEW |
(Empty) | |
| 1 /* |
| 2 * Copyright 2012 Google Inc. |
| 3 * |
| 4 * Use of this source code is governed by a BSD-style license that can be |
| 5 * found in the LICENSE file. |
| 6 */ |
| 7 #include "SkPathOpsTypes.h" |
| 8 |
| 9 const int UlpsEpsilon = 16; |
| 10 |
| 11 // from http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-num
bers-2012-edition/ |
| 12 union SkPathOpsUlpsFloat { |
| 13 int32_t fInt; |
| 14 float fFloat; |
| 15 |
| 16 SkPathOpsUlpsFloat(float num = 0.0f) : fFloat(num) {} |
| 17 bool negative() const { return fInt < 0; } |
| 18 }; |
| 19 |
| 20 bool AlmostEqualUlps(float A, float B) { |
| 21 SkPathOpsUlpsFloat uA(A); |
| 22 SkPathOpsUlpsFloat uB(B); |
| 23 // Different signs means they do not match. |
| 24 if (uA.negative() != uB.negative()) |
| 25 { |
| 26 // Check for equality to make sure +0 == -0 |
| 27 return A == B; |
| 28 } |
| 29 // Find the difference in ULPs. |
| 30 int ulpsDiff = abs(uA.fInt - uB.fInt); |
| 31 return ulpsDiff <= UlpsEpsilon; |
| 32 } |
| 33 |
| 34 // cube root approximation using bit hack for 64-bit float |
| 35 // adapted from Kahan's cbrt |
| 36 static double cbrt_5d(double d) { |
| 37 const unsigned int B1 = 715094163; |
| 38 double t = 0.0; |
| 39 unsigned int* pt = (unsigned int*) &t; |
| 40 unsigned int* px = (unsigned int*) &d; |
| 41 pt[1] = px[1] / 3 + B1; |
| 42 return t; |
| 43 } |
| 44 |
| 45 // iterative cube root approximation using Halley's method (double) |
| 46 static double cbrta_halleyd(const double a, const double R) { |
| 47 const double a3 = a * a * a; |
| 48 const double b = a * (a3 + R + R) / (a3 + a3 + R); |
| 49 return b; |
| 50 } |
| 51 |
| 52 // cube root approximation using 3 iterations of Halley's method (double) |
| 53 static double halley_cbrt3d(double d) { |
| 54 double a = cbrt_5d(d); |
| 55 a = cbrta_halleyd(a, d); |
| 56 a = cbrta_halleyd(a, d); |
| 57 return cbrta_halleyd(a, d); |
| 58 } |
| 59 |
| 60 double SkDCubeRoot(double x) { |
| 61 if (approximately_zero_cubed(x)) { |
| 62 return 0; |
| 63 } |
| 64 double result = halley_cbrt3d(fabs(x)); |
| 65 if (x < 0) { |
| 66 result = -result; |
| 67 } |
| 68 return result; |
| 69 } |
OLD | NEW |