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
19static constexpr double BL_MATH_PI = 3.14159265358979323846; //!< pi.
20static constexpr double BL_MATH_1p5_PI = 4.71238898038468985769; //!< pi * 1.5.
21static constexpr double BL_MATH_2_PI = 6.28318530717958647692; //!< pi * 2.
22static constexpr double BL_MATH_PI_DIV_2 = 1.57079632679489661923; //!< pi / 2.
23static constexpr double BL_MATH_PI_DIV_3 = 1.04719755119659774615; //!< pi / 3.
24static constexpr double BL_MATH_PI_DIV_4 = 0.78539816339744830962; //!< pi / 4.
25
26static constexpr double BL_SQRT_0p5 = 0.70710678118654746172; //!< sqrt(0.5).
27static constexpr double BL_SQRT_2 = 1.41421356237309504880; //!< sqrt(2).
28static constexpr double BL_SQRT_3 = 1.73205080756887729353; //!< sqrt(3).
29
30static constexpr double BL_MATH_AFTER_0 = 0.49e-323; //!< First value after 0.0.
31static constexpr double BL_MATH_BEFORE_1 = 0.999999999999999889; //!< First value before 1.0.
32
33static 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.
54static constexpr double BL_MATH_KAPPA = 0.55228474983;
55
56// ============================================================================
57// [Helper Functions]
58// ============================================================================
59
60template<typename T>
61static constexpr T blSum(const T& first) { return first; }
62
63template<typename T, typename... Args>
64static constexpr T blSum(const T& first, Args&&... args) { return first + blSum(std::forward<Args>(args)...); }
65
66// ============================================================================
67// [Classification & Limits]
68// ============================================================================
69
70template<typename T> constexpr T blEpsilon() noexcept = delete;
71template<> constexpr float blEpsilon<float>() noexcept { return 1e-8f; }
72template<> constexpr double blEpsilon<double>() noexcept { return 1e-14; }
73
74static BL_INLINE bool blIsNaN(float x) noexcept { return std::isnan(x); }
75static BL_INLINE bool blIsNaN(double x) noexcept { return std::isnan(x); }
76
77static BL_INLINE bool blIsInf(float x) noexcept { return std::isinf(x); }
78static BL_INLINE bool blIsInf(double x) noexcept { return std::isinf(x); }
79
80static BL_INLINE bool blIsFinite(float x) noexcept { return std::isfinite(x); }
81static BL_INLINE bool blIsFinite(double x) noexcept { return std::isfinite(x); }
82
83template<typename T, typename... Args>
84static BL_INLINE bool blIsNaN(const T& first, Args&&... args) noexcept {
85 return blIsNaN(first) | blIsNaN(std::forward<Args>(args)...);
86}
87
88template<typename T, typename... Args>
89static BL_INLINE bool blIsInf(const T& first, Args&&... args) noexcept {
90 return blIsInf(first) | blIsInf(std::forward<Args>(args)...);
91}
92
93template<typename T, typename... Args>
94static 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
102static BL_INLINE float blCopySign(float x, float y) noexcept { return std::copysign(x, y); }
103static 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
112namespace {
113
114template<int ControlFlags>
115BL_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
120template<int ControlFlags>
121BL_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
128static BL_INLINE float blNearby(float x) noexcept { return SIMD::vcvtf128f32(bl_roundf_sse4_1<_MM_FROUND_CUR_DIRECTION>(x)); }
129static BL_INLINE double blNearby(double x) noexcept { return SIMD::vcvtd128d64(bl_roundd_sse4_1<_MM_FROUND_CUR_DIRECTION>(x)); }
130static BL_INLINE float blTrunc (float x) noexcept { return SIMD::vcvtf128f32(bl_roundf_sse4_1<_MM_FROUND_TO_ZERO >(x)); }
131static BL_INLINE double blTrunc (double x) noexcept { return SIMD::vcvtd128d64(bl_roundd_sse4_1<_MM_FROUND_TO_ZERO >(x)); }
132static BL_INLINE float blFloor (float x) noexcept { return SIMD::vcvtf128f32(bl_roundf_sse4_1<_MM_FROUND_TO_NEG_INF >(x)); }
133static BL_INLINE double blFloor (double x) noexcept { return SIMD::vcvtd128d64(bl_roundd_sse4_1<_MM_FROUND_TO_NEG_INF >(x)); }
134static BL_INLINE float blCeil (float x) noexcept { return SIMD::vcvtf128f32(bl_roundf_sse4_1<_MM_FROUND_TO_POS_INF >(x)); }
135static 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
168static 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
180static 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
198static 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
211static 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
224static 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
236static 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
255static 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
268static 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
283BL_PRAGMA_FAST_MATH_PUSH
284
285static BL_INLINE float blNearby(float x) noexcept { return rintf(x); }
286static BL_INLINE double blNearby(double x) noexcept { return rint(x); }
287static BL_INLINE float blTrunc(float x) noexcept { return truncf(x); }
288static BL_INLINE double blTrunc(double x) noexcept { return trunc(x); }
289static BL_INLINE float blFloor(float x) noexcept { return floorf(x); }
290static BL_INLINE double blFloor(double x) noexcept { return floor(x); }
291static BL_INLINE float blCeil(float x) noexcept { return ceilf(x); }
292static BL_INLINE double blCeil(double x) noexcept { return ceil(x); }
293
294BL_PRAGMA_FAST_MATH_POP
295
296#endif
297
298static BL_INLINE float blRound(float x) noexcept { float y = blFloor(x); return y + (x - y >= 0.5f ? 1.0f : 0.0f); }
299static 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
305static 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
324static 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
343static BL_INLINE int blTruncToInt(float x) noexcept { return int(x); }
344static BL_INLINE int blTruncToInt(double x) noexcept { return int(x); }
345
346#if defined(BL_TARGET_OPT_SSE4_1)
347static BL_INLINE int blFloorToInt(float x) noexcept { return _mm_cvttss_si32(bl_roundf_sse4_1<_MM_FROUND_TO_NEG_INF>(x)); }
348static BL_INLINE int blFloorToInt(double x) noexcept { return _mm_cvttsd_si32(bl_roundd_sse4_1<_MM_FROUND_TO_NEG_INF>(x)); }
349
350static BL_INLINE int blCeilToInt(float x) noexcept { return _mm_cvttss_si32(bl_roundf_sse4_1<_MM_FROUND_TO_POS_INF>(x)); }
351static BL_INLINE int blCeilToInt(double x) noexcept { return _mm_cvttsd_si32(bl_roundd_sse4_1<_MM_FROUND_TO_POS_INF>(x)); }
352#else
353static BL_INLINE int blFloorToInt(float x) noexcept { int y = blNearbyToInt(x); return y - (float(y) > x); }
354static BL_INLINE int blFloorToInt(double x) noexcept { int y = blNearbyToInt(x); return y - (double(y) > x); }
355
356static BL_INLINE int blCeilToInt(float x) noexcept { int y = blNearbyToInt(x); return y + (float(y) < x); }
357static BL_INLINE int blCeilToInt(double x) noexcept { int y = blNearbyToInt(x); return y + (double(y) < x); }
358#endif
359
360static BL_INLINE int blRoundToInt(float x) noexcept { int y = blNearbyToInt(x); return y + (float(y) - x == -0.5f); }
361static BL_INLINE int blRoundToInt(double x) noexcept { int y = blNearbyToInt(x); return y + (double(y) - x == -0.5); }
362
363static 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
382static 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
401static BL_INLINE int64_t blTruncToInt64(float x) noexcept { return int64_t(x); }
402static BL_INLINE int64_t blTruncToInt64(double x) noexcept { return int64_t(x); }
403
404static BL_INLINE int64_t blFloorToInt64(float x) noexcept { int64_t y = blTruncToInt64(x); return y - int64_t(float(y) > x); }
405static BL_INLINE int64_t blFloorToInt64(double x) noexcept { int64_t y = blTruncToInt64(x); return y - int64_t(double(y) > x); }
406
407static BL_INLINE int64_t blCeilToInt64(float x) noexcept { int64_t y = blTruncToInt64(x); return y - int64_t(float(y) < x); }
408static BL_INLINE int64_t blCeilToInt64(double x) noexcept { int64_t y = blTruncToInt64(x); return y - int64_t(double(y) < x); }
409
410static BL_INLINE int64_t blRoundToInt64(float x) noexcept { int64_t y = blNearbyToInt64(x); return y + int64_t(float(y) - x == -0.5f); }
411static 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`.
422template<typename T>
423static 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)`.
428template<typename T>
429static 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
442static BL_INLINE float blPow(float x, float y) noexcept { return powf(x, y); }
443static BL_INLINE double blPow(double x, double y) noexcept { return pow(x, y); }
444
445template<typename T> constexpr T blSquare(const T& x) noexcept { return x * x; }
446template<typename T> constexpr T blPow3(const T& x) noexcept { return x * x * x; }
447
448BL_PRAGMA_FAST_MATH_PUSH
449
450static BL_INLINE float blSqrt(float x) noexcept { return sqrtf(x); }
451static BL_INLINE double blSqrt(double x) noexcept { return sqrt(x); }
452
453BL_PRAGMA_FAST_MATH_POP
454
455static BL_INLINE float blCbrt(float x) noexcept { return cbrtf(x); }
456static BL_INLINE double blCbrt(double x) noexcept { return cbrt(x); }
457
458static BL_INLINE float blHypot(float x, float y) noexcept { return hypotf(x, y); }
459static BL_INLINE double blHypot(double x, double y) noexcept { return hypot(x, y); }
460
461// ============================================================================
462// [Trigonometry]
463// ============================================================================
464
465static BL_INLINE float blSin(float x) noexcept { return sinf(x); }
466static BL_INLINE double blSin(double x) noexcept { return sin(x); }
467
468static BL_INLINE float blCos(float x) noexcept { return cosf(x); }
469static BL_INLINE double blCos(double x) noexcept { return cos(x); }
470
471static BL_INLINE float blTan(float x) noexcept { return tanf(x); }
472static BL_INLINE double blTan(double x) noexcept { return tan(x); }
473
474static BL_INLINE float blAsin(float x) noexcept { return asinf(x); }
475static BL_INLINE double blAsin(double x) noexcept { return asin(x); }
476
477static BL_INLINE float blAcos(float x) noexcept { return acosf(x); }
478static BL_INLINE double blAcos(double x) noexcept { return acos(x); }
479
480static BL_INLINE float blAtan(float x) noexcept { return atanf(x); }
481static BL_INLINE double blAtan(double x) noexcept { return atan(x); }
482
483static BL_INLINE float blAtan2(float y, float x) noexcept { return atan2f(y, x); }
484static 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.
496template<typename V, typename T = double>
497static 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`.
502template<typename T>
503static 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.
510template<typename V, typename T = double>
511static 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`.
516template<typename T>
517static 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.
558static 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
579static 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.
584static 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`.
597BL_HIDDEN size_t blCubicRoots(double* dst, const double* poly, double tMin, double tMax) noexcept;
598
599//! \overload
600static 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`.
608BL_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).
615template<typename T>
616static BL_INLINE bool blIsBetween0And1(const T& x) noexcept { return x >= T(0) && x <= T(1); }
617
618// ============================================================================
619// [BLMath - Near]
620// ============================================================================
621
622template<typename T>
623static constexpr bool isNear(T x, T y, T eps = blEpsilon<T>()) noexcept { return blAbs(x - y) <= eps; }
624
625template<typename T>
626static constexpr bool isNearZero(T x, T eps = blEpsilon<T>()) noexcept { return blAbs(x) <= eps; }
627
628template<typename T>
629static constexpr bool isNearZeroPositive(T x, T eps = blEpsilon<T>()) noexcept { return x >= T(0) && x <= eps; }
630
631template<typename T>
632static constexpr bool isNearOne(T x, T eps = blEpsilon<T>()) noexcept { return isNear(x, T(1), eps); }
633
634#endif // BLEND2D_BLMATH_P_H
635