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 
arguments_denormalized(float a,float b,int epsilon)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
equal_ulps(float a,float b,int epsilon,int depsilon)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 
d_equal_ulps(float a,float b,int epsilon)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 
not_equal_ulps(float a,float b,int epsilon)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 
d_not_equal_ulps(float a,float b,int epsilon)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 
less_ulps(float a,float b,int epsilon)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 
less_or_equal_ulps(float a,float b,int epsilon)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
AlmostBequalUlps(float a,float b)90 bool AlmostBequalUlps(float a, float b) {
91     const int UlpsEpsilon = 2;
92     return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
93 }
94 
AlmostPequalUlps(float a,float b)95 bool AlmostPequalUlps(float a, float b) {
96     const int UlpsEpsilon = 8;
97     return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
98 }
99 
AlmostDequalUlps(float a,float b)100 bool AlmostDequalUlps(float a, float b) {
101     const int UlpsEpsilon = 16;
102     return d_equal_ulps(a, b, UlpsEpsilon);
103 }
104 
AlmostDequalUlps(double a,double b)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 
AlmostEqualUlps(float a,float b)112 bool AlmostEqualUlps(float a, float b) {
113     const int UlpsEpsilon = 16;
114     return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
115 }
116 
NotAlmostEqualUlps(float a,float b)117 bool NotAlmostEqualUlps(float a, float b) {
118     const int UlpsEpsilon = 16;
119     return not_equal_ulps(a, b, UlpsEpsilon);
120 }
121 
NotAlmostDequalUlps(float a,float b)122 bool NotAlmostDequalUlps(float a, float b) {
123     const int UlpsEpsilon = 16;
124     return d_not_equal_ulps(a, b, UlpsEpsilon);
125 }
126 
RoughlyEqualUlps(float a,float b)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 
AlmostBetweenUlps(float a,float b,float c)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 
AlmostLessUlps(float a,float b)139 bool AlmostLessUlps(float a, float b) {
140     const int UlpsEpsilon = 16;
141     return less_ulps(a, b, UlpsEpsilon);
142 }
143 
AlmostLessOrEqualUlps(float a,float b)144 bool AlmostLessOrEqualUlps(float a, float b) {
145     const int UlpsEpsilon = 16;
146     return less_or_equal_ulps(a, b, UlpsEpsilon);
147 }
148 
UlpsDistance(float a,float b)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
cbrt_5d(double d)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)
cbrta_halleyd(const double a,const double R)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)
halley_cbrt3d(double d)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 
SkDCubeRoot(double x)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