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 |
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 | |
64 | constexpr float sk_float_degrees_to_radians(float degrees) { |
65 | return degrees * (SK_FloatPI / 180); |
66 | } |
67 | |
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 |
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 | |
84 | static inline bool sk_float_isfinite(float x) { |
85 | return SkFloatBits_IsFinite(SkFloat2Bits(x)); |
86 | } |
87 | |
88 | static inline bool sk_floats_are_finite(float a, float b) { |
89 | return sk_float_isfinite(a) && sk_float_isfinite(b); |
90 | } |
91 | |
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 | |
101 | static inline bool sk_float_isinf(float x) { |
102 | return SkFloatBits_IsInf(SkFloat2Bits(x)); |
103 | } |
104 | |
105 | static inline bool sk_float_isnan(float x) { |
106 | return !(x == x); |
107 | } |
108 | |
109 | #define sk_double_isnan(a) sk_float_isnan(a) |
110 | |
111 | #define SK_MaxS32FitsInFloat 2147483520 |
112 | #define SK_MinS32FitsInFloat -SK_MaxS32FitsInFloat |
113 | |
114 | #define SK_MaxS64FitsInFloat (SK_MaxS64 >> (63-24) << (63-24)) // 0x7fffff8000000000 |
115 | #define SK_MinS64FitsInFloat -SK_MaxS64FitsInFloat |
116 | |
117 | /** |
118 | * Return the closest int for the given float. Returns SK_MaxS32FitsInFloat for NaN. |
119 | */ |
120 | static inline int sk_float_saturate2int(float x) { |
121 | x = x < SK_MaxS32FitsInFloat ? x : SK_MaxS32FitsInFloat; |
122 | x = x > SK_MinS32FitsInFloat ? x : SK_MinS32FitsInFloat; |
123 | return (int)x; |
124 | } |
125 | |
126 | /** |
127 | * Return the closest int for the given double. Returns SK_MaxS32 for NaN. |
128 | */ |
129 | static inline int sk_double_saturate2int(double x) { |
130 | x = x < SK_MaxS32 ? x : SK_MaxS32; |
131 | x = x > SK_MinS32 ? x : SK_MinS32; |
132 | return (int)x; |
133 | } |
134 | |
135 | /** |
136 | * Return the closest int64_t for the given float. Returns SK_MaxS64FitsInFloat for NaN. |
137 | */ |
138 | static inline int64_t sk_float_saturate2int64(float x) { |
139 | x = x < SK_MaxS64FitsInFloat ? x : SK_MaxS64FitsInFloat; |
140 | x = x > SK_MinS64FitsInFloat ? x : SK_MinS64FitsInFloat; |
141 | return (int64_t)x; |
142 | } |
143 | |
144 | #define sk_float_floor2int(x) sk_float_saturate2int(sk_float_floor(x)) |
145 | #define sk_float_round2int(x) sk_float_saturate2int(sk_float_floor((x) + 0.5f)) |
146 | #define sk_float_ceil2int(x) sk_float_saturate2int(sk_float_ceil(x)) |
147 | |
148 | #define sk_float_floor2int_no_saturate(x) (int)sk_float_floor(x) |
149 | #define sk_float_round2int_no_saturate(x) (int)sk_float_floor((x) + 0.5f) |
150 | #define sk_float_ceil2int_no_saturate(x) (int)sk_float_ceil(x) |
151 | |
152 | #define sk_double_floor(x) floor(x) |
153 | #define sk_double_round(x) floor((x) + 0.5) |
154 | #define sk_double_ceil(x) ceil(x) |
155 | #define sk_double_floor2int(x) (int)floor(x) |
156 | #define sk_double_round2int(x) (int)floor((x) + 0.5) |
157 | #define sk_double_ceil2int(x) (int)ceil(x) |
158 | |
159 | // Cast double to float, ignoring any warning about too-large finite values being cast to float. |
160 | // Clang thinks this is undefined, but it's actually implementation defined to return either |
161 | // the largest float or infinity (one of the two bracketing representable floats). Good enough! |
162 | [[clang::no_sanitize("float-cast-overflow" )]] |
163 | static inline float sk_double_to_float(double x) { |
164 | return static_cast<float>(x); |
165 | } |
166 | |
167 | #define SK_FloatNaN std::numeric_limits<float>::quiet_NaN() |
168 | #define SK_FloatInfinity (+std::numeric_limits<float>::infinity()) |
169 | #define SK_FloatNegativeInfinity (-std::numeric_limits<float>::infinity()) |
170 | |
171 | #define SK_DoubleNaN std::numeric_limits<double>::quiet_NaN() |
172 | |
173 | // Returns false if any of the floats are outside of [0...1] |
174 | // Returns true if count is 0 |
175 | bool sk_floats_are_unit(const float array[], size_t count); |
176 | |
177 | static inline float sk_float_rsqrt_portable(float x) { |
178 | // Get initial estimate. |
179 | int i; |
180 | memcpy(&i, &x, 4); |
181 | i = 0x5F1FFFF9 - (i>>1); |
182 | float estimate; |
183 | memcpy(&estimate, &i, 4); |
184 | |
185 | // One step of Newton's method to refine. |
186 | const float estimate_sq = estimate*estimate; |
187 | estimate *= 0.703952253f*(2.38924456f-x*estimate_sq); |
188 | return estimate; |
189 | } |
190 | |
191 | // Fast, approximate inverse square root. |
192 | // Compare to name-brand "1.0f / sk_float_sqrt(x)". Should be around 10x faster on SSE, 2x on NEON. |
193 | static inline float sk_float_rsqrt(float x) { |
194 | // We want all this inlined, so we'll inline SIMD and just take the hit when we don't know we've got |
195 | // it at compile time. This is going to be too fast to productively hide behind a function pointer. |
196 | // |
197 | // We do one step of Newton's method to refine the estimates in the NEON and portable paths. No |
198 | // refinement is faster, but very innacurate. Two steps is more accurate, but slower than 1/sqrt. |
199 | // |
200 | // Optimized constants in the portable path courtesy of http://rrrola.wz.cz/inv_sqrt.html |
201 | #if SK_CPU_SSE_LEVEL >= SK_CPU_SSE_LEVEL_SSE1 |
202 | return _mm_cvtss_f32(_mm_rsqrt_ss(_mm_set_ss(x))); |
203 | #elif defined(SK_ARM_HAS_NEON) |
204 | // Get initial estimate. |
205 | const float32x2_t xx = vdup_n_f32(x); // Clever readers will note we're doing everything 2x. |
206 | float32x2_t estimate = vrsqrte_f32(xx); |
207 | |
208 | // One step of Newton's method to refine. |
209 | const float32x2_t estimate_sq = vmul_f32(estimate, estimate); |
210 | estimate = vmul_f32(estimate, vrsqrts_f32(xx, estimate_sq)); |
211 | return vget_lane_f32(estimate, 0); // 1 will work fine too; the answer's in both places. |
212 | #else |
213 | return sk_float_rsqrt_portable(x); |
214 | #endif |
215 | } |
216 | |
217 | // This is the number of significant digits we can print in a string such that when we read that |
218 | // string back we get the floating point number we expect. The minimum value C requires is 6, but |
219 | // most compilers support 9 |
220 | #ifdef FLT_DECIMAL_DIG |
221 | #define SK_FLT_DECIMAL_DIG FLT_DECIMAL_DIG |
222 | #else |
223 | #define SK_FLT_DECIMAL_DIG 9 |
224 | #endif |
225 | |
226 | // IEEE defines how float divide behaves for non-finite values and zero-denoms, but C does not |
227 | // so we have a helper that suppresses the possible undefined-behavior warnings. |
228 | |
229 | [[clang::no_sanitize("float-divide-by-zero" )]] |
230 | static inline float sk_ieee_float_divide(float numer, float denom) { |
231 | return numer / denom; |
232 | } |
233 | |
234 | [[clang::no_sanitize("float-divide-by-zero" )]] |
235 | static inline double sk_ieee_double_divide(double numer, double denom) { |
236 | return numer / denom; |
237 | } |
238 | |
239 | // While we clean up divide by zero, we'll replace places that do divide by zero with this TODO. |
240 | static inline float sk_ieee_float_divide_TODO_IS_DIVIDE_BY_ZERO_SAFE_HERE(float n, float d) { |
241 | return sk_ieee_float_divide(n,d); |
242 | } |
243 | static inline float sk_ieee_double_divide_TODO_IS_DIVIDE_BY_ZERO_SAFE_HERE(double n, double d) { |
244 | return sk_ieee_double_divide(n,d); |
245 | } |
246 | |
247 | static inline float sk_fmaf(float f, float m, float a) { |
248 | #if defined(FP_FAST_FMA) |
249 | return std::fmaf(f,m,a); |
250 | #else |
251 | return f*m+a; |
252 | #endif |
253 | } |
254 | |
255 | #endif |
256 | |