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
10static 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
17static 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
30static 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
40static 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
53static 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
63static 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
76static 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
90bool AlmostBequalUlps(float a, float b) {
91    const int UlpsEpsilon = 2;
92    return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
93}
94
95bool AlmostPequalUlps(float a, float b) {
96    const int UlpsEpsilon = 8;
97    return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
98}
99
100bool AlmostDequalUlps(float a, float b) {
101    const int UlpsEpsilon = 16;
102    return d_equal_ulps(a, b, UlpsEpsilon);
103}
104
105bool 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
112bool AlmostEqualUlps(float a, float b) {
113    const int UlpsEpsilon = 16;
114    return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
115}
116
117bool NotAlmostEqualUlps(float a, float b) {
118    const int UlpsEpsilon = 16;
119    return not_equal_ulps(a, b, UlpsEpsilon);
120}
121
122bool NotAlmostDequalUlps(float a, float b) {
123    const int UlpsEpsilon = 16;
124    return d_not_equal_ulps(a, b, UlpsEpsilon);
125}
126
127bool 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
133bool 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
139bool AlmostLessUlps(float a, float b) {
140    const int UlpsEpsilon = 16;
141    return less_ulps(a, b, UlpsEpsilon);
142}
143
144bool AlmostLessOrEqualUlps(float a, float b) {
145    const int UlpsEpsilon = 16;
146    return less_or_equal_ulps(a, b, UlpsEpsilon);
147}
148
149int 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
167static 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)
177static 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)
184static 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
191double 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