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
11static 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
18static 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
28static 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
41static 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
48static 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
58static 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
71static 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
78static 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
88static 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
99bool AlmostBequalUlps(float a, float b) {
100    const int UlpsEpsilon = 2;
101    return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
102}
103
104bool AlmostPequalUlps(float a, float b) {
105    const int UlpsEpsilon = 8;
106    return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
107}
108
109bool AlmostDequalUlps(float a, float b) {
110    const int UlpsEpsilon = 16;
111    return d_equal_ulps(a, b, UlpsEpsilon);
112}
113
114bool AlmostDequalUlps(double a, double b) {
115    return AlmostDequalUlps(SkDoubleToScalar(a), SkDoubleToScalar(b));
116}
117
118bool AlmostEqualUlps(float a, float b) {
119    const int UlpsEpsilon = 16;
120    return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
121}
122
123bool AlmostEqualUlps_Pin(float a, float b) {
124    const int UlpsEpsilon = 16;
125    return equal_ulps_pin(a, b, UlpsEpsilon, UlpsEpsilon);
126}
127
128bool NotAlmostEqualUlps(float a, float b) {
129    const int UlpsEpsilon = 16;
130    return not_equal_ulps(a, b, UlpsEpsilon);
131}
132
133bool NotAlmostEqualUlps_Pin(float a, float b) {
134    const int UlpsEpsilon = 16;
135    return not_equal_ulps_pin(a, b, UlpsEpsilon);
136}
137
138bool NotAlmostDequalUlps(float a, float b) {
139    const int UlpsEpsilon = 16;
140    return d_not_equal_ulps(a, b, UlpsEpsilon);
141}
142
143bool 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
149bool 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
155bool AlmostLessUlps(float a, float b) {
156    const int UlpsEpsilon = 16;
157    return less_ulps(a, b, UlpsEpsilon);
158}
159
160bool AlmostLessOrEqualUlps(float a, float b) {
161    const int UlpsEpsilon = 16;
162    return less_or_equal_ulps(a, b, UlpsEpsilon);
163}
164
165int 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
180static 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)
190static 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)
197static 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
204double 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
215SkOpGlobalState::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