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