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 static bool arguments_denormalized(float a, float b, int epsilon) { 11 float denormalizedCheck = FLT_EPSILON * epsilon / 2; 12 return fabsf(a) <= denormalizedCheck && fabsf(b) <= denormalizedCheck; 13 } 14 15 // from http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ 16 // FIXME: move to SkFloatBits.h 17 static bool equal_ulps(float a, float b, int epsilon, int depsilon) { 18 if (!SkScalarIsFinite(a) || !SkScalarIsFinite(b)) { 19 return false; 20 } 21 if (arguments_denormalized(a, b, depsilon)) { 22 return true; 23 } 24 int aBits = SkFloatAs2sCompliment(a); 25 int bBits = SkFloatAs2sCompliment(b); 26 // Find the difference in ULPs. 27 return aBits < bBits + epsilon && bBits < aBits + epsilon; 28 } 29 30 static bool d_equal_ulps(float a, float b, int epsilon) { 31 if (!SkScalarIsFinite(a) || !SkScalarIsFinite(b)) { 32 return false; 33 } 34 int aBits = SkFloatAs2sCompliment(a); 35 int bBits = SkFloatAs2sCompliment(b); 36 // Find the difference in ULPs. 37 return aBits < bBits + epsilon && bBits < aBits + epsilon; 38 } 39 40 static bool not_equal_ulps(float a, float b, int epsilon) { 41 if (!SkScalarIsFinite(a) || !SkScalarIsFinite(b)) { 42 return false; 43 } 44 if (arguments_denormalized(a, b, epsilon)) { 45 return false; 46 } 47 int aBits = SkFloatAs2sCompliment(a); 48 int bBits = SkFloatAs2sCompliment(b); 49 // Find the difference in ULPs. 50 return aBits >= bBits + epsilon || bBits >= aBits + epsilon; 51 } 52 53 static bool d_not_equal_ulps(float a, float b, int epsilon) { 54 if (!SkScalarIsFinite(a) || !SkScalarIsFinite(b)) { 55 return false; 56 } 57 int aBits = SkFloatAs2sCompliment(a); 58 int bBits = SkFloatAs2sCompliment(b); 59 // Find the difference in ULPs. 60 return aBits >= bBits + epsilon || bBits >= aBits + epsilon; 61 } 62 63 static bool less_ulps(float a, float b, int epsilon) { 64 if (!SkScalarIsFinite(a) || !SkScalarIsFinite(b)) { 65 return false; 66 } 67 if (arguments_denormalized(a, b, epsilon)) { 68 return a <= b - FLT_EPSILON * epsilon; 69 } 70 int aBits = SkFloatAs2sCompliment(a); 71 int bBits = SkFloatAs2sCompliment(b); 72 // Find the difference in ULPs. 73 return aBits <= bBits - epsilon; 74 } 75 76 static bool less_or_equal_ulps(float a, float b, int epsilon) { 77 if (!SkScalarIsFinite(a) || !SkScalarIsFinite(b)) { 78 return false; 79 } 80 if (arguments_denormalized(a, b, epsilon)) { 81 return a < b + FLT_EPSILON * epsilon; 82 } 83 int aBits = SkFloatAs2sCompliment(a); 84 int bBits = SkFloatAs2sCompliment(b); 85 // Find the difference in ULPs. 86 return aBits < bBits + epsilon; 87 } 88 89 // equality using the same error term as between 90 bool AlmostBequalUlps(float a, float b) { 91 const int UlpsEpsilon = 2; 92 return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon); 93 } 94 95 bool AlmostPequalUlps(float a, float b) { 96 const int UlpsEpsilon = 8; 97 return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon); 98 } 99 100 bool AlmostDequalUlps(float a, float b) { 101 const int UlpsEpsilon = 16; 102 return d_equal_ulps(a, b, UlpsEpsilon); 103 } 104 105 bool AlmostDequalUlps(double a, double b) { 106 if (SkScalarIsFinite(a) || SkScalarIsFinite(b)) { 107 return AlmostDequalUlps(SkDoubleToScalar(a), SkDoubleToScalar(b)); 108 } 109 return fabs(a - b) / SkTMax(fabs(a), fabs(b)) < FLT_EPSILON * 16; 110 } 111 112 bool AlmostEqualUlps(float a, float b) { 113 const int UlpsEpsilon = 16; 114 return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon); 115 } 116 117 bool NotAlmostEqualUlps(float a, float b) { 118 const int UlpsEpsilon = 16; 119 return not_equal_ulps(a, b, UlpsEpsilon); 120 } 121 122 bool NotAlmostDequalUlps(float a, float b) { 123 const int UlpsEpsilon = 16; 124 return d_not_equal_ulps(a, b, UlpsEpsilon); 125 } 126 127 bool RoughlyEqualUlps(float a, float b) { 128 const int UlpsEpsilon = 256; 129 const int DUlpsEpsilon = 1024; 130 return equal_ulps(a, b, UlpsEpsilon, DUlpsEpsilon); 131 } 132 133 bool AlmostBetweenUlps(float a, float b, float c) { 134 const int UlpsEpsilon = 2; 135 return a <= c ? less_or_equal_ulps(a, b, UlpsEpsilon) && less_or_equal_ulps(b, c, UlpsEpsilon) 136 : less_or_equal_ulps(b, a, UlpsEpsilon) && less_or_equal_ulps(c, b, UlpsEpsilon); 137 } 138 139 bool AlmostLessUlps(float a, float b) { 140 const int UlpsEpsilon = 16; 141 return less_ulps(a, b, UlpsEpsilon); 142 } 143 144 bool AlmostLessOrEqualUlps(float a, float b) { 145 const int UlpsEpsilon = 16; 146 return less_or_equal_ulps(a, b, UlpsEpsilon); 147 } 148 149 int UlpsDistance(float a, float b) { 150 if (!SkScalarIsFinite(a) || !SkScalarIsFinite(b)) { 151 return SK_MaxS32; 152 } 153 SkFloatIntUnion floatIntA, floatIntB; 154 floatIntA.fFloat = a; 155 floatIntB.fFloat = b; 156 // Different signs means they do not match. 157 if ((floatIntA.fSignBitInt < 0) != (floatIntB.fSignBitInt < 0)) { 158 // Check for equality to make sure +0 == -0 159 return a == b ? 0 : SK_MaxS32; 160 } 161 // Find the difference in ULPs. 162 return abs(floatIntA.fSignBitInt - floatIntB.fSignBitInt); 163 } 164 165 // cube root approximation using bit hack for 64-bit float 166 // adapted from Kahan's cbrt 167 static double cbrt_5d(double d) { 168 const unsigned int B1 = 715094163; 169 double t = 0.0; 170 unsigned int* pt = (unsigned int*) &t; 171 unsigned int* px = (unsigned int*) &d; 172 pt[1] = px[1] / 3 + B1; 173 return t; 174 } 175 176 // iterative cube root approximation using Halley's method (double) 177 static double cbrta_halleyd(const double a, const double R) { 178 const double a3 = a * a * a; 179 const double b = a * (a3 + R + R) / (a3 + a3 + R); 180 return b; 181 } 182 183 // cube root approximation using 3 iterations of Halley's method (double) 184 static double halley_cbrt3d(double d) { 185 double a = cbrt_5d(d); 186 a = cbrta_halleyd(a, d); 187 a = cbrta_halleyd(a, d); 188 return cbrta_halleyd(a, d); 189 } 190 191 double SkDCubeRoot(double x) { 192 if (approximately_zero_cubed(x)) { 193 return 0; 194 } 195 double result = halley_cbrt3d(fabs(x)); 196 if (x < 0) { 197 result = -result; 198 } 199 return result; 200 } 201
