1 /*
2 * Copyright 2006 The Android Open Source Project
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
8 #ifndef SkFloatingPoint_DEFINED
9 #define SkFloatingPoint_DEFINED
10
11 #include "include/core/SkTypes.h"
12 #include "include/private/SkFloatBits.h"
13 #include "include/private/SkSafe_math.h"
14 #include <float.h>
15 #include <math.h>
16 #include <cmath>
17 #include <cstring>
18 #include <limits>
19
20
21 #if SK_CPU_SSE_LEVEL >= SK_CPU_SSE_LEVEL_SSE1
22 #include <xmmintrin.h>
23 #elif defined(SK_ARM_HAS_NEON)
24 #include <arm_neon.h>
25 #endif
26
27 // For _POSIX_VERSION
28 #if defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
29 #include <unistd.h>
30 #endif
31
32 constexpr float SK_FloatSqrt2 = 1.41421356f;
33 constexpr float SK_FloatPI = 3.14159265f;
34 constexpr double SK_DoublePI = 3.14159265358979323846264338327950288;
35
36 // C++98 cmath std::pow seems to be the earliest portable way to get float pow.
37 // However, on Linux including cmath undefines isfinite.
38 // http://gcc.gnu.org/bugzilla/show_bug.cgi?id=14608
sk_float_pow(float base,float exp)39 static inline float sk_float_pow(float base, float exp) {
40 return powf(base, exp);
41 }
42
43 #define sk_float_sqrt(x) sqrtf(x)
44 #define sk_float_sin(x) sinf(x)
45 #define sk_float_cos(x) cosf(x)
46 #define sk_float_tan(x) tanf(x)
47 #define sk_float_floor(x) floorf(x)
48 #define sk_float_ceil(x) ceilf(x)
49 #define sk_float_trunc(x) truncf(x)
50 #ifdef SK_BUILD_FOR_MAC
51 # define sk_float_acos(x) static_cast<float>(acos(x))
52 # define sk_float_asin(x) static_cast<float>(asin(x))
53 #else
54 # define sk_float_acos(x) acosf(x)
55 # define sk_float_asin(x) asinf(x)
56 #endif
57 #define sk_float_atan2(y,x) atan2f(y,x)
58 #define sk_float_abs(x) fabsf(x)
59 #define sk_float_copysign(x, y) copysignf(x, y)
60 #define sk_float_mod(x,y) fmodf(x,y)
61 #define sk_float_exp(x) expf(x)
62 #define sk_float_log(x) logf(x)
63
sk_float_degrees_to_radians(float degrees)64 constexpr float sk_float_degrees_to_radians(float degrees) {
65 return degrees * (SK_FloatPI / 180);
66 }
67
sk_float_radians_to_degrees(float radians)68 constexpr float sk_float_radians_to_degrees(float radians) {
69 return radians * (180 / SK_FloatPI);
70 }
71
72 #define sk_float_round(x) sk_float_floor((x) + 0.5f)
73
74 // can't find log2f on android, but maybe that just a tool bug?
75 #ifdef SK_BUILD_FOR_ANDROID
sk_float_log2(float x)76 static inline float sk_float_log2(float x) {
77 const double inv_ln_2 = 1.44269504088896;
78 return (float)(log(x) * inv_ln_2);
79 }
80 #else
81 #define sk_float_log2(x) log2f(x)
82 #endif
83
sk_float_isfinite(float x)84 static inline bool sk_float_isfinite(float x) {
85 return SkFloatBits_IsFinite(SkFloat2Bits(x));
86 }
87
sk_floats_are_finite(float a,float b)88 static inline bool sk_floats_are_finite(float a, float b) {
89 return sk_float_isfinite(a) && sk_float_isfinite(b);
90 }
91
sk_floats_are_finite(const float array[],int count)92 static inline bool sk_floats_are_finite(const float array[], int count) {
93 float prod = 0;
94 for (int i = 0; i < count; ++i) {
95 prod *= array[i];
96 }
97 // At this point, prod will either be NaN or 0
98 return prod == 0; // if prod is NaN, this check will return false
99 }
100
sk_float_isinf(float x)101 static inline bool sk_float_isinf(float x) {
102 return SkFloatBits_IsInf(SkFloat2Bits(x));
103 }
104
105 #ifdef SK_BUILD_FOR_WIN
106 #define sk_float_isnan(x) _isnan(x)
107 #elif defined(__clang__) || defined(__GNUC__)
108 #define sk_float_isnan(x) __builtin_isnan(x)
109 #else
110 #define sk_float_isnan(x) isnan(x)
111 #endif
112
113 #define sk_double_isnan(a) sk_float_isnan(a)
114
115 #define SK_MaxS32FitsInFloat 2147483520
116 #define SK_MinS32FitsInFloat -SK_MaxS32FitsInFloat
117
118 #define SK_MaxS64FitsInFloat (SK_MaxS64 >> (63-24) << (63-24)) // 0x7fffff8000000000
119 #define SK_MinS64FitsInFloat -SK_MaxS64FitsInFloat
120
121 /**
122 * Return the closest int for the given float. Returns SK_MaxS32FitsInFloat for NaN.
123 */
sk_float_saturate2int(float x)124 static inline int sk_float_saturate2int(float x) {
125 x = SkTMin<float>(x, SK_MaxS32FitsInFloat);
126 x = SkTMax<float>(x, SK_MinS32FitsInFloat);
127 return (int)x;
128 }
129
130 /**
131 * Return the closest int for the given double. Returns SK_MaxS32 for NaN.
132 */
sk_double_saturate2int(double x)133 static inline int sk_double_saturate2int(double x) {
134 x = SkTMin<double>(x, SK_MaxS32);
135 x = SkTMax<double>(x, SK_MinS32);
136 return (int)x;
137 }
138
139 /**
140 * Return the closest int64_t for the given float. Returns SK_MaxS64FitsInFloat for NaN.
141 */
sk_float_saturate2int64(float x)142 static inline int64_t sk_float_saturate2int64(float x) {
143 x = SkTMin<float>(x, SK_MaxS64FitsInFloat);
144 x = SkTMax<float>(x, SK_MinS64FitsInFloat);
145 return (int64_t)x;
146 }
147
148 #define sk_float_floor2int(x) sk_float_saturate2int(sk_float_floor(x))
149 #define sk_float_round2int(x) sk_float_saturate2int(sk_float_floor((x) + 0.5f))
150 #define sk_float_ceil2int(x) sk_float_saturate2int(sk_float_ceil(x))
151
152 #define sk_float_floor2int_no_saturate(x) (int)sk_float_floor(x)
153 #define sk_float_round2int_no_saturate(x) (int)sk_float_floor((x) + 0.5f)
154 #define sk_float_ceil2int_no_saturate(x) (int)sk_float_ceil(x)
155
156 #define sk_double_floor(x) floor(x)
157 #define sk_double_round(x) floor((x) + 0.5)
158 #define sk_double_ceil(x) ceil(x)
159 #define sk_double_floor2int(x) (int)floor(x)
160 #define sk_double_round2int(x) (int)floor((x) + 0.5)
161 #define sk_double_ceil2int(x) (int)ceil(x)
162
163 // Cast double to float, ignoring any warning about too-large finite values being cast to float.
164 // Clang thinks this is undefined, but it's actually implementation defined to return either
165 // the largest float or infinity (one of the two bracketing representable floats). Good enough!
166 #if defined(__clang__) && (__clang_major__ * 1000 + __clang_minor__) >= 3007
167 __attribute__((no_sanitize("float-cast-overflow")))
168 #endif
sk_double_to_float(double x)169 static inline float sk_double_to_float(double x) {
170 return static_cast<float>(x);
171 }
172
173 #define SK_FloatNaN std::numeric_limits<float>::quiet_NaN()
174 #define SK_FloatInfinity (+std::numeric_limits<float>::infinity())
175 #define SK_FloatNegativeInfinity (-std::numeric_limits<float>::infinity())
176
177 #define SK_DoubleNaN std::numeric_limits<double>::quiet_NaN()
178
179 // Returns false if any of the floats are outside of [0...1]
180 // Returns true if count is 0
181 bool sk_floats_are_unit(const float array[], size_t count);
182
sk_float_rsqrt_portable(float x)183 static inline float sk_float_rsqrt_portable(float x) {
184 // Get initial estimate.
185 int i;
186 memcpy(&i, &x, 4);
187 i = 0x5F1FFFF9 - (i>>1);
188 float estimate;
189 memcpy(&estimate, &i, 4);
190
191 // One step of Newton's method to refine.
192 const float estimate_sq = estimate*estimate;
193 estimate *= 0.703952253f*(2.38924456f-x*estimate_sq);
194 return estimate;
195 }
196
197 // Fast, approximate inverse square root.
198 // Compare to name-brand "1.0f / sk_float_sqrt(x)". Should be around 10x faster on SSE, 2x on NEON.
sk_float_rsqrt(float x)199 static inline float sk_float_rsqrt(float x) {
200 // We want all this inlined, so we'll inline SIMD and just take the hit when we don't know we've got
201 // it at compile time. This is going to be too fast to productively hide behind a function pointer.
202 //
203 // We do one step of Newton's method to refine the estimates in the NEON and portable paths. No
204 // refinement is faster, but very innacurate. Two steps is more accurate, but slower than 1/sqrt.
205 //
206 // Optimized constants in the portable path courtesy of http://rrrola.wz.cz/inv_sqrt.html
207 #if SK_CPU_SSE_LEVEL >= SK_CPU_SSE_LEVEL_SSE1
208 return _mm_cvtss_f32(_mm_rsqrt_ss(_mm_set_ss(x)));
209 #elif defined(SK_ARM_HAS_NEON)
210 // Get initial estimate.
211 const float32x2_t xx = vdup_n_f32(x); // Clever readers will note we're doing everything 2x.
212 float32x2_t estimate = vrsqrte_f32(xx);
213
214 // One step of Newton's method to refine.
215 const float32x2_t estimate_sq = vmul_f32(estimate, estimate);
216 estimate = vmul_f32(estimate, vrsqrts_f32(xx, estimate_sq));
217 return vget_lane_f32(estimate, 0); // 1 will work fine too; the answer's in both places.
218 #else
219 return sk_float_rsqrt_portable(x);
220 #endif
221 }
222
223 // This is the number of significant digits we can print in a string such that when we read that
224 // string back we get the floating point number we expect. The minimum value C requires is 6, but
225 // most compilers support 9
226 #ifdef FLT_DECIMAL_DIG
227 #define SK_FLT_DECIMAL_DIG FLT_DECIMAL_DIG
228 #else
229 #define SK_FLT_DECIMAL_DIG 9
230 #endif
231
232 // IEEE defines how float divide behaves for non-finite values and zero-denoms, but C does not
233 // so we have a helper that suppresses the possible undefined-behavior warnings.
234
235 #ifdef __clang__
236 __attribute__((no_sanitize("float-divide-by-zero")))
237 #endif
sk_ieee_float_divide(float numer,float denom)238 static inline float sk_ieee_float_divide(float numer, float denom) {
239 return numer / denom;
240 }
241
242 #ifdef __clang__
243 __attribute__((no_sanitize("float-divide-by-zero")))
244 #endif
sk_ieee_double_divide(double numer,double denom)245 static inline double sk_ieee_double_divide(double numer, double denom) {
246 return numer / denom;
247 }
248
249 // While we clean up divide by zero, we'll replace places that do divide by zero with this TODO.
sk_ieee_float_divide_TODO_IS_DIVIDE_BY_ZERO_SAFE_HERE(float n,float d)250 static inline float sk_ieee_float_divide_TODO_IS_DIVIDE_BY_ZERO_SAFE_HERE(float n, float d) {
251 return sk_ieee_float_divide(n,d);
252 }
sk_ieee_double_divide_TODO_IS_DIVIDE_BY_ZERO_SAFE_HERE(double n,double d)253 static inline float sk_ieee_double_divide_TODO_IS_DIVIDE_BY_ZERO_SAFE_HERE(double n, double d) {
254 return sk_ieee_double_divide(n,d);
255 }
256
sk_fmaf(float f,float m,float a)257 static inline float sk_fmaf(float f, float m, float a) {
258 #if defined(FP_FAST_FMA)
259 return std::fmaf(f,m,a);
260 #else
261 return f*m+a;
262 #endif
263 }
264
265 #endif
266