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 AlmostEqualUlps(float a, float b) {
106    const int UlpsEpsilon = 16;
107    return equal_ulps(a, b, UlpsEpsilon, UlpsEpsilon);
108}
109
110bool NotAlmostEqualUlps(float a, float b) {
111    const int UlpsEpsilon = 16;
112    return not_equal_ulps(a, b, UlpsEpsilon);
113}
114
115bool NotAlmostDequalUlps(float a, float b) {
116    const int UlpsEpsilon = 16;
117    return d_not_equal_ulps(a, b, UlpsEpsilon);
118}
119
120bool RoughlyEqualUlps(float a, float b) {
121    const int UlpsEpsilon = 256;
122    const int DUlpsEpsilon = 1024;
123    return equal_ulps(a, b, UlpsEpsilon, DUlpsEpsilon);
124}
125
126bool AlmostBetweenUlps(float a, float b, float c) {
127    const int UlpsEpsilon = 2;
128    return a <= c ? less_or_equal_ulps(a, b, UlpsEpsilon) && less_or_equal_ulps(b, c, UlpsEpsilon)
129        : less_or_equal_ulps(b, a, UlpsEpsilon) && less_or_equal_ulps(c, b, UlpsEpsilon);
130}
131
132bool AlmostLessUlps(float a, float b) {
133    const int UlpsEpsilon = 16;
134    return less_ulps(a, b, UlpsEpsilon);
135}
136
137bool AlmostLessOrEqualUlps(float a, float b) {
138    const int UlpsEpsilon = 16;
139    return less_or_equal_ulps(a, b, UlpsEpsilon);
140}
141
142int UlpsDistance(float a, float b) {
143    if (!SkScalarIsFinite(a) || !SkScalarIsFinite(b)) {
144        return SK_MaxS32;
145    }
146    SkFloatIntUnion floatIntA, floatIntB;
147    floatIntA.fFloat = a;
148    floatIntB.fFloat = b;
149    // Different signs means they do not match.
150    if ((floatIntA.fSignBitInt < 0) != (floatIntB.fSignBitInt < 0)) {
151        // Check for equality to make sure +0 == -0
152        return a == b ? 0 : SK_MaxS32;
153    }
154    // Find the difference in ULPs.
155    return abs(floatIntA.fSignBitInt - floatIntB.fSignBitInt);
156}
157
158// cube root approximation using bit hack for 64-bit float
159// adapted from Kahan's cbrt
160static double cbrt_5d(double d) {
161    const unsigned int B1 = 715094163;
162    double t = 0.0;
163    unsigned int* pt = (unsigned int*) &t;
164    unsigned int* px = (unsigned int*) &d;
165    pt[1] = px[1] / 3 + B1;
166    return t;
167}
168
169// iterative cube root approximation using Halley's method (double)
170static double cbrta_halleyd(const double a, const double R) {
171    const double a3 = a * a * a;
172    const double b = a * (a3 + R + R) / (a3 + a3 + R);
173    return b;
174}
175
176// cube root approximation using 3 iterations of Halley's method (double)
177static double halley_cbrt3d(double d) {
178    double a = cbrt_5d(d);
179    a = cbrta_halleyd(a, d);
180    a = cbrta_halleyd(a, d);
181    return cbrta_halleyd(a, d);
182}
183
184double SkDCubeRoot(double x) {
185    if (approximately_zero_cubed(x)) {
186        return 0;
187    }
188    double result = halley_cbrt3d(fabs(x));
189    if (x < 0) {
190        result = -result;
191    }
192    return result;
193}
194