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