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 "SkFloatBits.h" 8 #include "SkPathOpsTypes.h" 9 10 11 12 // from http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ 13 // FIXME: move to SkFloatBits.h 14 static bool equal_ulps(float a, float b, int epsilon) { 15 SkFloatIntUnion floatIntA, floatIntB; 16 floatIntA.fFloat = a; 17 floatIntB.fFloat = b; 18 // Different signs means they do not match. 19 if ((floatIntA.fSignBitInt < 0) != (floatIntB.fSignBitInt < 0)) { 20 // Check for equality to make sure +0 == -0 21 return a == b; 22 } 23 // Find the difference in ULPs. 24 int ulpsDiff = abs(floatIntA.fSignBitInt - floatIntB.fSignBitInt); 25 return ulpsDiff <= epsilon; 26 } 27 28 static bool less_ulps(float a, float b, int epsilon) { 29 SkFloatIntUnion floatIntA, floatIntB; 30 floatIntA.fFloat = a; 31 floatIntB.fFloat = b; 32 // Check different signs with float epsilon since we only care if they're both close to 0. 33 if ((floatIntA.fSignBitInt < 0) != (floatIntB.fSignBitInt < 0)) { 34 return a <= b + FLT_EPSILON * epsilon; 35 } 36 // Find the difference in ULPs. 37 return floatIntA.fSignBitInt <= floatIntB.fSignBitInt + epsilon; 38 } 39 40 bool AlmostEqualUlps(float a, float b) { 41 const int UlpsEpsilon = 16; 42 return equal_ulps(a, b, UlpsEpsilon); 43 } 44 45 bool RoughlyEqualUlps(float a, float b) { 46 const int UlpsEpsilon = 256; 47 return equal_ulps(a, b, UlpsEpsilon); 48 } 49 50 bool AlmostBetweenUlps(float a, float b, float c) { 51 const int UlpsEpsilon = 1; 52 return a <= c ? less_ulps(a, b, UlpsEpsilon) && less_ulps(b, c, UlpsEpsilon) 53 : less_ulps(b, a, UlpsEpsilon) && less_ulps(c, b, UlpsEpsilon); 54 } 55 56 int UlpsDistance(float a, float b) { 57 SkFloatIntUnion floatIntA, floatIntB; 58 floatIntA.fFloat = a; 59 floatIntB.fFloat = b; 60 // Different signs means they do not match. 61 if ((floatIntA.fSignBitInt < 0) != (floatIntB.fSignBitInt < 0)) { 62 // Check for equality to make sure +0 == -0 63 return a == b ? 0 : SK_MaxS32; 64 } 65 // Find the difference in ULPs. 66 return abs(floatIntA.fSignBitInt - floatIntB.fSignBitInt); 67 } 68 69 // cube root approximation using bit hack for 64-bit float 70 // adapted from Kahan's cbrt 71 static double cbrt_5d(double d) { 72 const unsigned int B1 = 715094163; 73 double t = 0.0; 74 unsigned int* pt = (unsigned int*) &t; 75 unsigned int* px = (unsigned int*) &d; 76 pt[1] = px[1] / 3 + B1; 77 return t; 78 } 79 80 // iterative cube root approximation using Halley's method (double) 81 static double cbrta_halleyd(const double a, const double R) { 82 const double a3 = a * a * a; 83 const double b = a * (a3 + R + R) / (a3 + a3 + R); 84 return b; 85 } 86 87 // cube root approximation using 3 iterations of Halley's method (double) 88 static double halley_cbrt3d(double d) { 89 double a = cbrt_5d(d); 90 a = cbrta_halleyd(a, d); 91 a = cbrta_halleyd(a, d); 92 return cbrta_halleyd(a, d); 93 } 94 95 double SkDCubeRoot(double x) { 96 if (approximately_zero_cubed(x)) { 97 return 0; 98 } 99 double result = halley_cbrt3d(fabs(x)); 100 if (x < 0) { 101 result = -result; 102 } 103 return result; 104 } 105
