1 | // [Blend2D] |
2 | // 2D Vector Graphics Powered by a JIT Compiler. |
3 | // |
4 | // [License] |
5 | // Zlib - See LICENSE.md file in the package. |
6 | |
7 | #ifndef BLEND2D_BLMATH_P_H |
8 | #define BLEND2D_BLMATH_P_H |
9 | |
10 | #include "./blapi-internal_p.h" |
11 | #include "./blsimd_p.h" |
12 | #include "./blsupport_p.h" |
13 | #include "./bltables_p.h" |
14 | |
15 | // ============================================================================ |
16 | // [Global Constants] |
17 | // ============================================================================ |
18 | |
19 | static constexpr double BL_MATH_PI = 3.14159265358979323846; //!< pi. |
20 | static constexpr double BL_MATH_1p5_PI = 4.71238898038468985769; //!< pi * 1.5. |
21 | static constexpr double BL_MATH_2_PI = 6.28318530717958647692; //!< pi * 2. |
22 | static constexpr double BL_MATH_PI_DIV_2 = 1.57079632679489661923; //!< pi / 2. |
23 | static constexpr double BL_MATH_PI_DIV_3 = 1.04719755119659774615; //!< pi / 3. |
24 | static constexpr double BL_MATH_PI_DIV_4 = 0.78539816339744830962; //!< pi / 4. |
25 | |
26 | static constexpr double BL_SQRT_0p5 = 0.70710678118654746172; //!< sqrt(0.5). |
27 | static constexpr double BL_SQRT_2 = 1.41421356237309504880; //!< sqrt(2). |
28 | static constexpr double BL_SQRT_3 = 1.73205080756887729353; //!< sqrt(3). |
29 | |
30 | static constexpr double BL_MATH_AFTER_0 = 0.49e-323; //!< First value after 0.0. |
31 | static constexpr double BL_MATH_BEFORE_1 = 0.999999999999999889; //!< First value before 1.0. |
32 | |
33 | static constexpr double BL_MATH_ANGLE_EPSILON = 1e-8; |
34 | |
35 | //! Constant that is used to approximate elliptic arcs with cubic curves. |
36 | //! Since it's an approximation there are various approaches that can be |
37 | //! used to calculate the best value. The most used KAPPA is: |
38 | //! |
39 | //! k = (4/3) * (sqrt(2) - 1) ~= 0.55228474983 |
40 | //! |
41 | //! which has a maximum error of 0.00027253. However, according to this post |
42 | //! |
43 | //! http://spencermortensen.com/articles/bezier-circle/ |
44 | //! |
45 | //! the maximum error can be further reduced by 28% if we change the |
46 | //! approximation constraint to have the maximum radial distance from the |
47 | //! circle to the curve as small as possible. The an alternative constant |
48 | //! |
49 | //! k = 1/2 +- sqrt(12 - 20*c - 3*c^2)/(4 - 6*c) ~= 0.551915024494. |
50 | //! |
51 | //! can be used to reduce the maximum error to 0.00019608. We don't use the |
52 | //! alternative, because we need to caculate the KAPPA for arcs that are not |
53 | //! 90deg, in that case the KAPPA must be calculated for such angles. |
54 | static constexpr double BL_MATH_KAPPA = 0.55228474983; |
55 | |
56 | // ============================================================================ |
57 | // [Helper Functions] |
58 | // ============================================================================ |
59 | |
60 | template<typename T> |
61 | static constexpr T blSum(const T& first) { return first; } |
62 | |
63 | template<typename T, typename... Args> |
64 | static constexpr T blSum(const T& first, Args&&... args) { return first + blSum(std::forward<Args>(args)...); } |
65 | |
66 | // ============================================================================ |
67 | // [Classification & Limits] |
68 | // ============================================================================ |
69 | |
70 | template<typename T> constexpr T blEpsilon() noexcept = delete; |
71 | template<> constexpr float blEpsilon<float>() noexcept { return 1e-8f; } |
72 | template<> constexpr double blEpsilon<double>() noexcept { return 1e-14; } |
73 | |
74 | static BL_INLINE bool blIsNaN(float x) noexcept { return std::isnan(x); } |
75 | static BL_INLINE bool blIsNaN(double x) noexcept { return std::isnan(x); } |
76 | |
77 | static BL_INLINE bool blIsInf(float x) noexcept { return std::isinf(x); } |
78 | static BL_INLINE bool blIsInf(double x) noexcept { return std::isinf(x); } |
79 | |
80 | static BL_INLINE bool blIsFinite(float x) noexcept { return std::isfinite(x); } |
81 | static BL_INLINE bool blIsFinite(double x) noexcept { return std::isfinite(x); } |
82 | |
83 | template<typename T, typename... Args> |
84 | static BL_INLINE bool blIsNaN(const T& first, Args&&... args) noexcept { |
85 | return blIsNaN(first) | blIsNaN(std::forward<Args>(args)...); |
86 | } |
87 | |
88 | template<typename T, typename... Args> |
89 | static BL_INLINE bool blIsInf(const T& first, Args&&... args) noexcept { |
90 | return blIsInf(first) | blIsInf(std::forward<Args>(args)...); |
91 | } |
92 | |
93 | template<typename T, typename... Args> |
94 | static BL_INLINE bool blIsFinite(const T& first, Args&&... args) noexcept { |
95 | return blIsFinite(first) & blIsFinite(std::forward<Args>(args)...); |
96 | } |
97 | |
98 | // ============================================================================ |
99 | // [Miscellaneous] |
100 | // ============================================================================ |
101 | |
102 | static BL_INLINE float blCopySign(float x, float y) noexcept { return std::copysign(x, y); } |
103 | static BL_INLINE double blCopySign(double x, double y) noexcept { return std::copysign(x, y); } |
104 | |
105 | // ============================================================================ |
106 | // [Rounding] |
107 | // ============================================================================ |
108 | |
109 | |
110 | #if defined(BL_TARGET_OPT_SSE4_1) |
111 | |
112 | namespace { |
113 | |
114 | template<int ControlFlags> |
115 | BL_INLINE __m128 bl_roundf_sse4_1(float x) noexcept { |
116 | __m128 y = SIMD::vcvtf32f128(x); |
117 | return _mm_round_ss(y, y, ControlFlags | _MM_FROUND_NO_EXC); |
118 | } |
119 | |
120 | template<int ControlFlags> |
121 | BL_INLINE __m128d bl_roundd_sse4_1(double x) noexcept { |
122 | __m128d y = SIMD::vcvtd64d128(x); |
123 | return _mm_round_sd(y, y, ControlFlags | _MM_FROUND_NO_EXC); |
124 | } |
125 | |
126 | } // {anonymous} |
127 | |
128 | static BL_INLINE float blNearby(float x) noexcept { return SIMD::vcvtf128f32(bl_roundf_sse4_1<_MM_FROUND_CUR_DIRECTION>(x)); } |
129 | static BL_INLINE double blNearby(double x) noexcept { return SIMD::vcvtd128d64(bl_roundd_sse4_1<_MM_FROUND_CUR_DIRECTION>(x)); } |
130 | static BL_INLINE float blTrunc (float x) noexcept { return SIMD::vcvtf128f32(bl_roundf_sse4_1<_MM_FROUND_TO_ZERO >(x)); } |
131 | static BL_INLINE double blTrunc (double x) noexcept { return SIMD::vcvtd128d64(bl_roundd_sse4_1<_MM_FROUND_TO_ZERO >(x)); } |
132 | static BL_INLINE float blFloor (float x) noexcept { return SIMD::vcvtf128f32(bl_roundf_sse4_1<_MM_FROUND_TO_NEG_INF >(x)); } |
133 | static BL_INLINE double blFloor (double x) noexcept { return SIMD::vcvtd128d64(bl_roundd_sse4_1<_MM_FROUND_TO_NEG_INF >(x)); } |
134 | static BL_INLINE float blCeil (float x) noexcept { return SIMD::vcvtf128f32(bl_roundf_sse4_1<_MM_FROUND_TO_POS_INF >(x)); } |
135 | static BL_INLINE double blCeil (double x) noexcept { return SIMD::vcvtd128d64(bl_roundd_sse4_1<_MM_FROUND_TO_POS_INF >(x)); } |
136 | |
137 | #elif defined(BL_TARGET_OPT_SSE2) |
138 | |
139 | // Rounding is very expensive on pre-SSE4.1 X86 hardware as it requires to alter |
140 | // rounding bits of FPU/SSE states. The only method which is cheap is `rint()`, |
141 | // which uses the current FPU/SSE rounding mode to round a floating point. The |
142 | // code below can be used to implement `roundeven()`, which can be then used to |
143 | // implement any other rounding operation. Blend2D implementation then assumes |
144 | // that `blNearby` is round to even as it's the default mode CPU is setup to. |
145 | // |
146 | // Single Precision |
147 | // ---------------- |
148 | // |
149 | // ``` |
150 | // float roundeven(float x) { |
151 | // float maxn = pow(2, 23); // 8388608 |
152 | // float magic = pow(2, 23) + pow(2, 22); // 12582912 |
153 | // return x >= maxn ? x : x + magic - magic; |
154 | // } |
155 | // ``` |
156 | // |
157 | // Double Precision |
158 | // ---------------- |
159 | // |
160 | // ``` |
161 | // double roundeven(double x) { |
162 | // double maxn = pow(2, 52); // 4503599627370496 |
163 | // double magic = pow(2, 52) + pow(2, 51); // 6755399441055744 |
164 | // return x >= maxn ? x : x + magic - magic; |
165 | // } |
166 | // ``` |
167 | |
168 | static BL_INLINE float blNearby(float x) noexcept { |
169 | using namespace SIMD; |
170 | |
171 | F128 src = vcvtf32f128(x); |
172 | F128 magic = v_const_as<F128>(blCommonTable.f128_round_magic); |
173 | |
174 | F128 mask = vcmpgess(src, v_const_as<F128>(blCommonTable.f128_round_max)); |
175 | F128 rounded = vsubss(vaddss(src, magic), magic); |
176 | |
177 | return vcvtf128f32(vblendmask(rounded, src, mask)); |
178 | } |
179 | |
180 | static BL_INLINE float blTrunc(float x) noexcept { |
181 | using namespace SIMD; |
182 | |
183 | F128 src = vcvtf32f128(x); |
184 | |
185 | F128 msk_abs = v_const_as<F128>(blCommonTable.f128_abs); |
186 | F128 src_abs = vand(src, msk_abs); |
187 | |
188 | F128 sign = vandnot_a(msk_abs, src); |
189 | F128 magic = v_const_as<F128>(blCommonTable.f128_round_magic); |
190 | |
191 | F128 mask = vor(vcmpgess(src_abs, v_const_as<F128>(blCommonTable.f128_round_max)), sign); |
192 | F128 rounded = vsubss(vaddss(src_abs, magic), magic); |
193 | F128 maybeone = vand(vcmpltss(src_abs, rounded), v_const_as<F128>(blCommonTable.f128_1)); |
194 | |
195 | return vcvtf128f32(vblendmask(vsubss(rounded, maybeone), src, mask)); |
196 | } |
197 | |
198 | static BL_INLINE float blFloor(float x) noexcept { |
199 | using namespace SIMD; |
200 | |
201 | F128 src = vcvtf32f128(x); |
202 | F128 magic = v_const_as<F128>(blCommonTable.f128_round_magic); |
203 | |
204 | F128 mask = vcmpgess(src, v_const_as<F128>(blCommonTable.f128_round_max)); |
205 | F128 rounded = vsubss(vaddss(src, magic), magic); |
206 | F128 maybeone = vand(vcmpltss(src, rounded), v_const_as<F128>(blCommonTable.f128_1)); |
207 | |
208 | return vcvtf128f32(vblendmask(vsubss(rounded, maybeone), src, mask)); |
209 | } |
210 | |
211 | static BL_INLINE float blCeil(float x) noexcept { |
212 | using namespace SIMD; |
213 | |
214 | F128 src = SIMD::vcvtf32f128(x); |
215 | F128 magic = v_const_as<F128>(blCommonTable.f128_round_magic); |
216 | |
217 | F128 mask = vcmpgess(src, v_const_as<F128>(blCommonTable.f128_round_max)); |
218 | F128 rounded = vsubss(vaddss(src, magic), magic); |
219 | F128 maybeone = vand(vcmpgtss(src, rounded), v_const_as<F128>(blCommonTable.f128_1)); |
220 | |
221 | return vcvtf128f32(vblendmask(vaddss(rounded, maybeone), src, mask)); |
222 | } |
223 | |
224 | static BL_INLINE double blNearby(double x) noexcept { |
225 | using namespace SIMD; |
226 | |
227 | D128 src = vcvtd64d128(x); |
228 | D128 magic = v_const_as<D128>(blCommonTable.d128_round_magic); |
229 | |
230 | D128 mask = vcmpgesd(src, v_const_as<D128>(blCommonTable.d128_round_max)); |
231 | D128 rounded = vsubsd(vaddsd(src, magic), magic); |
232 | |
233 | return vcvtd128d64(vblendmask(rounded, src, mask)); |
234 | } |
235 | |
236 | static BL_INLINE double blTrunc(double x) noexcept { |
237 | using namespace SIMD; |
238 | |
239 | static const uint64_t kSepMask[1] = { 0x7FFFFFFFFFFFFFFFu }; |
240 | |
241 | D128 src = vcvtd64d128(x); |
242 | D128 msk_abs = _mm_load_sd(reinterpret_cast<const double*>(kSepMask)); |
243 | D128 src_abs = vand(src, msk_abs); |
244 | |
245 | D128 sign = vandnot_a(msk_abs, src); |
246 | D128 magic = v_const_as<D128>(blCommonTable.d128_round_magic); |
247 | |
248 | D128 mask = vor(vcmpgesd(src_abs, v_const_as<D128>(blCommonTable.d128_round_max)), sign); |
249 | D128 rounded = vsubsd(vaddsd(src_abs, magic), magic); |
250 | D128 maybeone = vand(vcmpltsd(src_abs, rounded), v_const_as<D128>(blCommonTable.d128_1)); |
251 | |
252 | return vcvtd128d64(vblendmask(vsubsd(rounded, maybeone), src, mask)); |
253 | } |
254 | |
255 | static BL_INLINE double blFloor(double x) noexcept { |
256 | using namespace SIMD; |
257 | |
258 | D128 src = vcvtd64d128(x); |
259 | D128 magic = v_const_as<D128>(blCommonTable.d128_round_magic); |
260 | |
261 | D128 mask = vcmpgesd(src, v_const_as<D128>(blCommonTable.d128_round_max)); |
262 | D128 rounded = vsubsd(vaddsd(src, magic), magic); |
263 | D128 maybeone = vand(vcmpltsd(src, rounded), v_const_as<D128>(blCommonTable.d128_1)); |
264 | |
265 | return vcvtd128d64(vblendmask(vsubsd(rounded, maybeone), src, mask)); |
266 | } |
267 | |
268 | static BL_INLINE double blCeil(double x) noexcept { |
269 | using namespace SIMD; |
270 | |
271 | D128 src = vcvtd64d128(x); |
272 | D128 magic = v_const_as<D128>(blCommonTable.d128_round_magic); |
273 | |
274 | D128 mask = vcmpgesd(src, v_const_as<D128>(blCommonTable.d128_round_max)); |
275 | D128 rounded = vsubsd(vaddsd(src, magic), magic); |
276 | D128 maybeone = vand(vcmpgtsd(src, rounded), v_const_as<D128>(blCommonTable.d128_1)); |
277 | |
278 | return vcvtd128d64(vblendmask(vaddsd(rounded, maybeone), src, mask)); |
279 | } |
280 | |
281 | #else |
282 | |
283 | BL_PRAGMA_FAST_MATH_PUSH |
284 | |
285 | static BL_INLINE float blNearby(float x) noexcept { return rintf(x); } |
286 | static BL_INLINE double blNearby(double x) noexcept { return rint(x); } |
287 | static BL_INLINE float blTrunc(float x) noexcept { return truncf(x); } |
288 | static BL_INLINE double blTrunc(double x) noexcept { return trunc(x); } |
289 | static BL_INLINE float blFloor(float x) noexcept { return floorf(x); } |
290 | static BL_INLINE double blFloor(double x) noexcept { return floor(x); } |
291 | static BL_INLINE float blCeil(float x) noexcept { return ceilf(x); } |
292 | static BL_INLINE double blCeil(double x) noexcept { return ceil(x); } |
293 | |
294 | BL_PRAGMA_FAST_MATH_POP |
295 | |
296 | #endif |
297 | |
298 | static BL_INLINE float blRound(float x) noexcept { float y = blFloor(x); return y + (x - y >= 0.5f ? 1.0f : 0.0f); } |
299 | static BL_INLINE double blRound(double x) noexcept { double y = blFloor(x); return y + (x - y >= 0.5 ? 1.0 : 0.0); } |
300 | |
301 | // ============================================================================ |
302 | // [Rounding to Integer] |
303 | // ============================================================================ |
304 | |
305 | static BL_INLINE int blNearbyToInt(float x) noexcept { |
306 | #if defined(BL_TARGET_OPT_SSE) |
307 | return _mm_cvtss_si32(SIMD::vcvtf32f128(x)); |
308 | #elif BL_TARGET_ARCH_X86 == 32 && defined(__GNUC__) |
309 | int y; |
310 | __asm__ __volatile__("flds %1\n" "fistpl %0\n" : "=m" (y) : "m" (x)); |
311 | return y; |
312 | #elif BL_TARGET_ARCH_X86 == 32 && defined(_MSC_VER) |
313 | int y; |
314 | __asm { |
315 | fld dword ptr [x] |
316 | fistp dword ptr [y] |
317 | } |
318 | return y; |
319 | #else |
320 | return int(lrintf(x)); |
321 | #endif |
322 | } |
323 | |
324 | static BL_INLINE int blNearbyToInt(double x) noexcept { |
325 | #if defined(BL_TARGET_OPT_SSE2) |
326 | return _mm_cvtsd_si32(SIMD::vcvtd64d128(x)); |
327 | #elif BL_TARGET_ARCH_X86 == 32 && defined(__GNUC__) |
328 | int y; |
329 | __asm__ __volatile__("fldl %1\n" "fistpl %0\n" : "=m" (y) : "m" (x)); |
330 | return y; |
331 | #elif BL_TARGET_ARCH_X86 == 32 && defined(_MSC_VER) |
332 | int y; |
333 | __asm { |
334 | fld qword ptr [x] |
335 | fistp dword ptr [y] |
336 | } |
337 | return y; |
338 | #else |
339 | return int(lrint(x)); |
340 | #endif |
341 | } |
342 | |
343 | static BL_INLINE int blTruncToInt(float x) noexcept { return int(x); } |
344 | static BL_INLINE int blTruncToInt(double x) noexcept { return int(x); } |
345 | |
346 | #if defined(BL_TARGET_OPT_SSE4_1) |
347 | static BL_INLINE int blFloorToInt(float x) noexcept { return _mm_cvttss_si32(bl_roundf_sse4_1<_MM_FROUND_TO_NEG_INF>(x)); } |
348 | static BL_INLINE int blFloorToInt(double x) noexcept { return _mm_cvttsd_si32(bl_roundd_sse4_1<_MM_FROUND_TO_NEG_INF>(x)); } |
349 | |
350 | static BL_INLINE int blCeilToInt(float x) noexcept { return _mm_cvttss_si32(bl_roundf_sse4_1<_MM_FROUND_TO_POS_INF>(x)); } |
351 | static BL_INLINE int blCeilToInt(double x) noexcept { return _mm_cvttsd_si32(bl_roundd_sse4_1<_MM_FROUND_TO_POS_INF>(x)); } |
352 | #else |
353 | static BL_INLINE int blFloorToInt(float x) noexcept { int y = blNearbyToInt(x); return y - (float(y) > x); } |
354 | static BL_INLINE int blFloorToInt(double x) noexcept { int y = blNearbyToInt(x); return y - (double(y) > x); } |
355 | |
356 | static BL_INLINE int blCeilToInt(float x) noexcept { int y = blNearbyToInt(x); return y + (float(y) < x); } |
357 | static BL_INLINE int blCeilToInt(double x) noexcept { int y = blNearbyToInt(x); return y + (double(y) < x); } |
358 | #endif |
359 | |
360 | static BL_INLINE int blRoundToInt(float x) noexcept { int y = blNearbyToInt(x); return y + (float(y) - x == -0.5f); } |
361 | static BL_INLINE int blRoundToInt(double x) noexcept { int y = blNearbyToInt(x); return y + (double(y) - x == -0.5); } |
362 | |
363 | static BL_INLINE int64_t blNearbyToInt64(float x) noexcept { |
364 | #if BL_TARGET_ARCH_X86 == 64 |
365 | return _mm_cvtss_si64(SIMD::vcvtf32f128(x)); |
366 | #elif BL_TARGET_ARCH_X86 == 32 && defined(__GNUC__) |
367 | int64_t y; |
368 | __asm__ __volatile__("flds %1\n" "fistpq %0\n" : "=m" (y) : "m" (x)); |
369 | return y; |
370 | #elif BL_TARGET_ARCH_X86 == 32 && defined(_MSC_VER) |
371 | int64_t y; |
372 | __asm { |
373 | fld dword ptr [x] |
374 | fistp qword ptr [y] |
375 | } |
376 | return y; |
377 | #else |
378 | return int64_t(llrintf(x)); |
379 | #endif |
380 | } |
381 | |
382 | static BL_INLINE int64_t blNearbyToInt64(double x) noexcept { |
383 | #if BL_TARGET_ARCH_X86 == 64 |
384 | return _mm_cvtsd_si64(_mm_set_sd(x)); |
385 | #elif BL_TARGET_ARCH_X86 == 32 && defined(__GNUC__) |
386 | int64_t y; |
387 | __asm__ __volatile__("fldl %1\n" "fistpq %0\n" : "=m" (y) : "m" (x)); |
388 | return y; |
389 | #elif BL_TARGET_ARCH_X86 == 32 && defined(_MSC_VER) |
390 | int64_t y; |
391 | __asm { |
392 | fld qword ptr [x] |
393 | fistp qword ptr [y] |
394 | } |
395 | return y; |
396 | #else |
397 | return int64_t(llrint(x)); |
398 | #endif |
399 | } |
400 | |
401 | static BL_INLINE int64_t blTruncToInt64(float x) noexcept { return int64_t(x); } |
402 | static BL_INLINE int64_t blTruncToInt64(double x) noexcept { return int64_t(x); } |
403 | |
404 | static BL_INLINE int64_t blFloorToInt64(float x) noexcept { int64_t y = blTruncToInt64(x); return y - int64_t(float(y) > x); } |
405 | static BL_INLINE int64_t blFloorToInt64(double x) noexcept { int64_t y = blTruncToInt64(x); return y - int64_t(double(y) > x); } |
406 | |
407 | static BL_INLINE int64_t blCeilToInt64(float x) noexcept { int64_t y = blTruncToInt64(x); return y - int64_t(float(y) < x); } |
408 | static BL_INLINE int64_t blCeilToInt64(double x) noexcept { int64_t y = blTruncToInt64(x); return y - int64_t(double(y) < x); } |
409 | |
410 | static BL_INLINE int64_t blRoundToInt64(float x) noexcept { int64_t y = blNearbyToInt64(x); return y + int64_t(float(y) - x == -0.5f); } |
411 | static BL_INLINE int64_t blRoundToInt64(double x) noexcept { int64_t y = blNearbyToInt64(x); return y + int64_t(double(y) - x == -0.5); } |
412 | |
413 | // ============================================================================ |
414 | // [Fraction / Repeat] |
415 | // ============================================================================ |
416 | |
417 | //! Returns a fractional part of `x`. |
418 | //! |
419 | //! \note Fractional part returned is always equal or greater than zero. The |
420 | //! implementation is compatible to many shader implementations defined as |
421 | //! `frac(x) == x - floor(x)`, which would return `0.25` for `-1.75`. |
422 | template<typename T> |
423 | static BL_INLINE T blFrac(T x) noexcept { return x - blFloor(x); } |
424 | |
425 | //! Repeats the given value `x` in `y`, returning a value that is always equal |
426 | //! to or greater than zero and lesser than `y`. The return of `repeat(x, 1.0)` |
427 | //! should be identical to the return of `frac(x)`. |
428 | template<typename T> |
429 | static BL_INLINE T blRepeat(T x, T y) noexcept { |
430 | T a = x; |
431 | if (a >= y || a <= -y) |
432 | a = std::fmod(a, y); |
433 | if (a < T(0)) |
434 | a += y; |
435 | return a; |
436 | } |
437 | |
438 | // ============================================================================ |
439 | // [Power] |
440 | // ============================================================================ |
441 | |
442 | static BL_INLINE float blPow(float x, float y) noexcept { return powf(x, y); } |
443 | static BL_INLINE double blPow(double x, double y) noexcept { return pow(x, y); } |
444 | |
445 | template<typename T> constexpr T blSquare(const T& x) noexcept { return x * x; } |
446 | template<typename T> constexpr T blPow3(const T& x) noexcept { return x * x * x; } |
447 | |
448 | BL_PRAGMA_FAST_MATH_PUSH |
449 | |
450 | static BL_INLINE float blSqrt(float x) noexcept { return sqrtf(x); } |
451 | static BL_INLINE double blSqrt(double x) noexcept { return sqrt(x); } |
452 | |
453 | BL_PRAGMA_FAST_MATH_POP |
454 | |
455 | static BL_INLINE float blCbrt(float x) noexcept { return cbrtf(x); } |
456 | static BL_INLINE double blCbrt(double x) noexcept { return cbrt(x); } |
457 | |
458 | static BL_INLINE float blHypot(float x, float y) noexcept { return hypotf(x, y); } |
459 | static BL_INLINE double blHypot(double x, double y) noexcept { return hypot(x, y); } |
460 | |
461 | // ============================================================================ |
462 | // [Trigonometry] |
463 | // ============================================================================ |
464 | |
465 | static BL_INLINE float blSin(float x) noexcept { return sinf(x); } |
466 | static BL_INLINE double blSin(double x) noexcept { return sin(x); } |
467 | |
468 | static BL_INLINE float blCos(float x) noexcept { return cosf(x); } |
469 | static BL_INLINE double blCos(double x) noexcept { return cos(x); } |
470 | |
471 | static BL_INLINE float blTan(float x) noexcept { return tanf(x); } |
472 | static BL_INLINE double blTan(double x) noexcept { return tan(x); } |
473 | |
474 | static BL_INLINE float blAsin(float x) noexcept { return asinf(x); } |
475 | static BL_INLINE double blAsin(double x) noexcept { return asin(x); } |
476 | |
477 | static BL_INLINE float blAcos(float x) noexcept { return acosf(x); } |
478 | static BL_INLINE double blAcos(double x) noexcept { return acos(x); } |
479 | |
480 | static BL_INLINE float blAtan(float x) noexcept { return atanf(x); } |
481 | static BL_INLINE double blAtan(double x) noexcept { return atan(x); } |
482 | |
483 | static BL_INLINE float blAtan2(float y, float x) noexcept { return atan2f(y, x); } |
484 | static BL_INLINE double blAtan2(double y, double x) noexcept { return atan2(y, x); } |
485 | |
486 | // ============================================================================ |
487 | // [Linear Interpolation] |
488 | // ============================================================================ |
489 | |
490 | //! Linear interpolation of `a` and `b` at `t`. |
491 | //! |
492 | //! Returns `(a - t * a) + t * b`. |
493 | //! |
494 | //! \note This function should work with most geometric types Blend2D offers |
495 | //! that use double precision, however, it's not compatible with integral types. |
496 | template<typename V, typename T = double> |
497 | static BL_INLINE V blLerp(const V& a, const V& b, const T& t) noexcept { |
498 | return (a - t * a) + t * b; |
499 | } |
500 | |
501 | // Linear interpolation of `a` and `b` at `t=0.5`. |
502 | template<typename T> |
503 | static BL_INLINE T blLerp(const T& a, const T& b) noexcept { |
504 | return 0.5 * a + 0.5 * b; |
505 | } |
506 | |
507 | //! Alternative LERP implementation that is faster, but won't handle pathological |
508 | //! inputs. It should only be used in places in which it's known that such inputs |
509 | //! cannot happen. |
510 | template<typename V, typename T = double> |
511 | static BL_INLINE V blFastLerp(const V& a, const V& b, const T& t) noexcept { |
512 | return a + t * (b - a); |
513 | } |
514 | |
515 | //! Alternative LERP implementation at `t=0.5`. |
516 | template<typename T> |
517 | static BL_INLINE T blFastLerp(const T& a, const T& b) noexcept { |
518 | return 0.5 * (a + b); |
519 | } |
520 | |
521 | // ============================================================================ |
522 | // [Roots] |
523 | // ============================================================================ |
524 | |
525 | //! Solve a quadratic polynomial `Ax^2 + Bx + C = 0` and store the result in `dst`. |
526 | //! |
527 | //! Returns the number of roots found within [tMin, tMax] - `0` to `2`. |
528 | //! |
529 | //! Resources: |
530 | //! - http://stackoverflow.com/questions/4503849/quadratic-equation-in-ada/4504415#4504415 |
531 | //! - http://people.csail.mit.edu/bkph/articles/Quadratics.pdf |
532 | //! |
533 | //! The standard equation: |
534 | //! |
535 | //! ``` |
536 | //! x0 = (-b + sqrt(delta)) / 2a |
537 | //! x1 = (-b - sqrt(delta)) / 2a |
538 | //! ``` |
539 | //! |
540 | //! When 4*a*c < b*b, computing x0 involves subtracting close numbers, and makes |
541 | //! you lose accuracy, so use the following instead: |
542 | //! |
543 | //! ``` |
544 | //! x0 = 2c / (-b - sqrt(delta)) |
545 | //! x1 = 2c / (-b + sqrt(delta)) |
546 | //! ``` |
547 | //! |
548 | //! Which yields a better x0, but whose x1 has the same problem as x0 had above. |
549 | //! The correct way to compute the roots is therefore: |
550 | //! |
551 | //! ``` |
552 | //! q = -0.5 * (b + sign(b) * sqrt(delta)) |
553 | //! x0 = q / a |
554 | //! x1 = c / q |
555 | //! ``` |
556 | //! |
557 | //! \note This is a branchless version designed to be easily inlineable. |
558 | static BL_INLINE size_t blQuadRoots(double dst[2], double a, double b, double c, double tMin, double tMax) noexcept { |
559 | double d = blMax(b * b - 4.0 * a * c, 0.0); |
560 | double s = blSqrt(d); |
561 | double q = -0.5 * (b + blCopySign(s, b)); |
562 | |
563 | double t0 = q / a; |
564 | double t1 = c / q; |
565 | |
566 | double x0 = blMin(t0, t1); |
567 | double x1 = blMax(t1, t0); |
568 | |
569 | dst[0] = x0; |
570 | size_t n = size_t((x0 >= tMin) & (x0 <= tMax)); |
571 | |
572 | dst[n] = x1; |
573 | n += size_t((x1 > x0) & (x1 >= tMin) & (x1 <= tMax)); |
574 | |
575 | return n; |
576 | } |
577 | |
578 | //! \overload |
579 | static BL_INLINE size_t blQuadRoots(double dst[2], const double poly[3], double tMin, double tMax) noexcept { |
580 | return blQuadRoots(dst, poly[0], poly[1], poly[2], tMin, tMax); |
581 | } |
582 | |
583 | //! Like `blQuadRoots()`, but always returns two roots and doesn't sort them. |
584 | static BL_INLINE size_t blSimplifiedQuadRoots(double dst[2], double a, double b, double c) noexcept { |
585 | double d = blMax(b * b - 4.0 * a * c, 0.0); |
586 | double s = blSqrt(d); |
587 | double q = -0.5 * (b + blCopySign(s, b)); |
588 | |
589 | dst[0] = q / a; |
590 | dst[1] = c / q; |
591 | return 2; |
592 | } |
593 | |
594 | //! Solve a cubic polynomial and store the result in `dst`. |
595 | //! |
596 | //! Returns the number of roots found within [tMin, tMax] - `0` to `3`. |
597 | BL_HIDDEN size_t blCubicRoots(double* dst, const double* poly, double tMin, double tMax) noexcept; |
598 | |
599 | //! \overload |
600 | static BL_INLINE size_t blCubicRoots(double dst[3], double a, double b, double c, double d, double tMin, double tMax) noexcept { |
601 | double poly[4] = { a, b, c, d }; |
602 | return ::blCubicRoots(dst, poly, tMin, tMax); |
603 | } |
604 | |
605 | //! Solve an Nth degree polynomial and store the result in `dst`. |
606 | //! |
607 | //! Returns the number of roots found [tMin, tMax] - `0` to `n - 1`. |
608 | BL_HIDDEN size_t blPolyRoots(double* dst, const double* poly, int n, double tMin, double tMax) noexcept; |
609 | |
610 | // ============================================================================ |
611 | // [blIsBetween0And1] |
612 | // ============================================================================ |
613 | |
614 | //! Check if `x` is within [0, 1] range (inclusive). |
615 | template<typename T> |
616 | static BL_INLINE bool blIsBetween0And1(const T& x) noexcept { return x >= T(0) && x <= T(1); } |
617 | |
618 | // ============================================================================ |
619 | // [BLMath - Near] |
620 | // ============================================================================ |
621 | |
622 | template<typename T> |
623 | static constexpr bool isNear(T x, T y, T eps = blEpsilon<T>()) noexcept { return blAbs(x - y) <= eps; } |
624 | |
625 | template<typename T> |
626 | static constexpr bool isNearZero(T x, T eps = blEpsilon<T>()) noexcept { return blAbs(x) <= eps; } |
627 | |
628 | template<typename T> |
629 | static constexpr bool isNearZeroPositive(T x, T eps = blEpsilon<T>()) noexcept { return x >= T(0) && x <= eps; } |
630 | |
631 | template<typename T> |
632 | static constexpr bool isNearOne(T x, T eps = blEpsilon<T>()) noexcept { return isNear(x, T(1), eps); } |
633 | |
634 | #endif // BLEND2D_BLMATH_P_H |
635 | |