math_p.h (23343B)
1 // This file is part of Blend2D project <https://blend2d.com> 2 // 3 // See blend2d.h or LICENSE.md for license and copyright information 4 // SPDX-License-Identifier: Zlib 5 6 #ifndef BLEND2D_SUPPORT_MATH_P_H_INCLUDED 7 #define BLEND2D_SUPPORT_MATH_P_H_INCLUDED 8 9 #include "../api-internal_p.h" 10 #include "../geometry.h" 11 #include "../simd/simd_p.h" 12 #include "../tables/tables_p.h" 13 #include "../support/mathconst_p.h" 14 15 //! \cond INTERNAL 16 //! \addtogroup blend2d_internal 17 //! \{ 18 19 namespace bl { 20 namespace Math { 21 namespace { 22 23 //! \name Floating Point Constants 24 //! \{ 25 26 //! Returns infinity of `T` type. 27 //! 28 //! \note `T` should be floating point. 29 template<typename T> 30 [[nodiscard]] 31 static BL_INLINE_CONSTEXPR T inf() noexcept { return std::numeric_limits<T>::infinity(); } 32 33 //! Returns a quiet NaN of `T` type. 34 //! 35 //! \note `T` should be floating point. 36 template<typename T> 37 [[nodiscard]] 38 static BL_INLINE_CONSTEXPR T nan() noexcept { return std::numeric_limits<T>::quiet_NaN(); } 39 40 template<typename T> constexpr T epsilon() noexcept = delete; 41 template<> BL_INLINE_CONSTEXPR float epsilon<float>() noexcept { return 1e-8f; } 42 template<> BL_INLINE_CONSTEXPR double epsilon<double>() noexcept { return 1e-14; } 43 44 //! \} 45 46 //! \name Floating Point Testing 47 //! \{ 48 49 static BL_INLINE_NODEBUG bool is_nan(float x) noexcept { return std::isnan(x); } 50 static BL_INLINE_NODEBUG bool is_nan(double x) noexcept { return std::isnan(x); } 51 52 template<typename T, typename... Args> 53 static BL_INLINE_NODEBUG bool is_nan(const T& first, Args&&... args) noexcept { 54 return BLInternal::bool_or(is_nan(first), (is_nan(forward<Args>(args)))...); 55 } 56 57 static BL_INLINE_NODEBUG bool is_inf(float x) noexcept { return std::isinf(x); } 58 static BL_INLINE_NODEBUG bool is_inf(double x) noexcept { return std::isinf(x); } 59 60 template<typename T, typename... Args> 61 static BL_INLINE_NODEBUG bool is_inf(const T& first, Args&&... args) noexcept { 62 return BLInternal::bool_or(is_inf(first), (is_inf(forward<Args>(args)))...); 63 } 64 65 static BL_INLINE_NODEBUG bool is_finite(const float& x) noexcept { return std::isfinite(x); } 66 static BL_INLINE_NODEBUG bool is_finite(const double& x) noexcept { return std::isfinite(x); } 67 68 static BL_INLINE_NODEBUG bool is_finite(const BLPoint& p) noexcept; 69 static BL_INLINE_NODEBUG bool is_finite(const BLBox& b) noexcept; 70 static BL_INLINE_NODEBUG bool is_finite(const BLRect& r) noexcept; 71 72 template<typename T, typename... Args> 73 static BL_INLINE_NODEBUG bool is_finite(T first, Args&&... args) noexcept { 74 return BLInternal::bool_and(is_finite(first), (is_finite(BLInternal::forward<Args>(args)))...); 75 } 76 77 static BL_INLINE_NODEBUG bool is_finite(const BLPoint& p) noexcept { return is_finite(p.x, p.y); } 78 static BL_INLINE_NODEBUG bool is_finite(const BLBox& b) noexcept { return is_finite(b.x0, b.y0, b.x1, b.y1); } 79 static BL_INLINE_NODEBUG bool is_finite(const BLRect& r) noexcept { return is_finite(r.x, r.y, r.w, r.h); } 80 81 static BL_INLINE_NODEBUG bool is_nan(const BLPoint& p) noexcept { return is_nan(p.x, p.y); } 82 83 template<typename T> 84 static BL_INLINE_CONSTEXPR bool is_near(T x, T y, T eps = epsilon<T>()) noexcept { return bl_abs(x - y) <= eps; } 85 86 template<typename T> 87 static BL_INLINE_CONSTEXPR bool is_near_zero(T x, T eps = epsilon<T>()) noexcept { return bl_abs(x) <= eps; } 88 89 template<typename T> 90 static BL_INLINE_CONSTEXPR bool is_near_zero_positive(T x, T eps = epsilon<T>()) noexcept { return x >= T(0) && x <= eps; } 91 92 template<typename T> 93 static BL_INLINE_CONSTEXPR bool is_near_one(T x, T eps = epsilon<T>()) noexcept { return is_near(x, T(1), eps); } 94 95 //! Check if `x` is within [0, 1] range (inclusive). 96 template<typename T> 97 static BL_INLINE bool is_between_0_and_1(const T& x) noexcept { return bool_and(x >= T(0), x <= T(1)); } 98 99 //! \} 100 101 //! \name Sum of Arguments 102 //! \{ 103 104 template<typename T> 105 static BL_INLINE_CONSTEXPR T sum(const T& first) { return first; } 106 107 template<typename T, typename... Args> 108 static BL_INLINE_CONSTEXPR T sum(const T& first, Args&&... args) { return first + sum(BLInternal::forward<Args>(args)...); } 109 110 //! \} 111 112 //! \name Miscellaneous Functions 113 //! \{ 114 115 static BL_INLINE_NODEBUG float copy_sign(float x, float y) noexcept { return std::copysign(x, y); } 116 static BL_INLINE_NODEBUG double copy_sign(double x, double y) noexcept { return std::copysign(x, y); } 117 static BL_INLINE_NODEBUG BLPoint copy_sign(const BLPoint& a, const BLPoint& b) noexcept { return BLPoint(copy_sign(a.x, b.x), copy_sign(a.y, b.y)); } 118 119 static BL_INLINE_NODEBUG float cut_off(float x, uint32_t bits) noexcept { 120 uint32_t msk = (uint32_t(1) << bits) - 1u; 121 return bl_bit_cast<float>(bl_bit_cast<uint32_t>(x) & ~msk); 122 } 123 124 static BL_INLINE_NODEBUG double cut_off(double x, uint32_t bits) noexcept { 125 uint64_t msk = (uint64_t(1) << bits) - 1u; 126 return bl_bit_cast<double>(bl_bit_cast<uint64_t>(x) & ~msk); 127 } 128 129 //! \} 130 131 //! \name FMA or Mul+Add (depending on target / compile flags) 132 //! \{ 133 134 static BL_INLINE_NODEBUG float madd(float x, float y, float a) noexcept { return x * y + a; } 135 static BL_INLINE_NODEBUG double madd(double x, double y, double a) noexcept { return x * y + a; } 136 137 //! \} 138 139 //! \name Rounding 140 //! \{ 141 142 #if defined(BL_TARGET_OPT_SSE4_1) 143 144 namespace { 145 146 template<int ControlFlags> 147 BL_INLINE __m128 bl_roundf_sse4_1(float x) noexcept { 148 __m128 y = SIMD::cast_from_f32(x).v; 149 return _mm_round_ss(y, y, ControlFlags | _MM_FROUND_NO_EXC); 150 } 151 152 template<int ControlFlags> 153 BL_INLINE __m128d bl_roundd_sse4_1(double x) noexcept { 154 __m128d y = SIMD::cast_from_f64(x).v; 155 return _mm_round_sd(y, y, ControlFlags | _MM_FROUND_NO_EXC); 156 } 157 158 } // {anonymous} 159 160 static BL_INLINE float nearby(float x) noexcept { return _mm_cvtss_f32(bl_roundf_sse4_1<_MM_FROUND_CUR_DIRECTION>(x)); } 161 static BL_INLINE double nearby(double x) noexcept { return _mm_cvtsd_f64(bl_roundd_sse4_1<_MM_FROUND_CUR_DIRECTION>(x)); } 162 static BL_INLINE float trunc(float x) noexcept { return _mm_cvtss_f32(bl_roundf_sse4_1<_MM_FROUND_TO_ZERO>(x)); } 163 static BL_INLINE double trunc(double x) noexcept { return _mm_cvtsd_f64(bl_roundd_sse4_1<_MM_FROUND_TO_ZERO>(x)); } 164 static BL_INLINE float floor(float x) noexcept { return _mm_cvtss_f32(bl_roundf_sse4_1<_MM_FROUND_TO_NEG_INF>(x)); } 165 static BL_INLINE double floor(double x) noexcept { return _mm_cvtsd_f64(bl_roundd_sse4_1<_MM_FROUND_TO_NEG_INF>(x)); } 166 static BL_INLINE float ceil(float x) noexcept { return _mm_cvtss_f32(bl_roundf_sse4_1<_MM_FROUND_TO_POS_INF>(x)); } 167 static BL_INLINE double ceil(double x) noexcept { return _mm_cvtsd_f64(bl_roundd_sse4_1<_MM_FROUND_TO_POS_INF>(x)); } 168 169 #elif defined(BL_TARGET_OPT_SSE2) 170 171 // Rounding is very expensive on pre-SSE4.1 X86 hardware as it requires to alter rounding bits of FPU/SSE states. 172 // The only method which is cheap is `rint()`, which uses the current FPU/SSE rounding mode to round a floating 173 // point. The code below can be used to implement `roundeven()`, which can be then used to implement any other 174 // rounding operation. Blend2D implementation then assumes that `nearby` is round to even as it's the default 175 // mode CPU is setup to and in general software like math libraries expect that this mode is used. 176 // 177 // Single Precision 178 // ---------------- 179 // 180 // ``` 181 // float roundeven(float x) { 182 // float magic = x >= 0 ? pow(2, 22) : pow(2, 22) + pow(2, 21); 183 // return x >= magic ? x : x + magic - magic; 184 // } 185 // ``` 186 // 187 // Double Precision 188 // ---------------- 189 // 190 // ``` 191 // double roundeven(double x) { 192 // double magic = x >= 0 ? pow(2, 52) : pow(2, 52) + pow(2, 51); 193 // return x >= magic ? x : x + magic - magic; 194 // } 195 // ``` 196 197 static BL_INLINE float nearby(float x) noexcept { 198 using namespace SIMD; 199 200 Vec4xF32 src = cast_from_f32(x); 201 Vec4xF32 cvt = Vec4xF32{_mm_cvt_si2ss(src.v, _mm_cvt_ss2si(src.v))}; 202 Vec4xF32 src_abs = src & common_table.p_7FFFFFFF7FFFFFFF.as<Vec4xF32>(); 203 Vec4xF32 mask = cmp_lt_f32x1(src_abs, common_table.f32_round_magic.as<Vec4xF32>()); 204 Vec4xF32 result = blendv_bits(src, cvt, mask); 205 206 return cast_to_f32(result); 207 } 208 209 static BL_INLINE float trunc(float x) noexcept { 210 using namespace SIMD; 211 212 Vec4xF32 src = cast_from_f32(x); 213 Vec4xF32 cvt = Vec4xF32{_mm_cvt_si2ss(src.v, _mm_cvtt_ss2si(src.v))}; 214 Vec4xF32 src_abs = src & common_table.p_7FFFFFFF7FFFFFFF.as<Vec4xF32>(); 215 Vec4xF32 mask = cmp_lt_f32x1(src_abs, common_table.f32_round_magic.as<Vec4xF32>()); 216 Vec4xF32 result = blendv_bits(src, cvt, mask); 217 218 return cast_to_f32(result); 219 } 220 221 static BL_INLINE float floor(float x) noexcept { 222 using namespace SIMD; 223 224 Vec4xF32 src = cast_from_f32(x); 225 Vec4xF32 magic = slli_u32<32 - 9>(srli_u32<31>(src)) | common_table.f32_round_magic.as<Vec4xF32>(); 226 227 Vec4xF32 mask = cmp_ge_f32x1(src, magic); 228 Vec4xF32 rounded = sub_f32x1(add_f32x1(src, magic), magic); 229 Vec4xF32 maybeone = cmp_lt_f32x1(src, rounded) & common_table.f32_1.as<Vec4xF32>(); 230 231 return cast_to_f32(blendv_bits(sub_f32x1(rounded, maybeone), src, mask)); 232 } 233 234 static BL_INLINE float ceil(float x) noexcept { 235 using namespace SIMD; 236 237 Vec4xF32 src = cast_from_f32(x); 238 Vec4xF32 magic = slli_u32<32 - 9>(srli_u32<31>(src)) | common_table.f32_round_magic.as<Vec4xF32>(); 239 240 Vec4xF32 mask = cmp_ge_f32x1(src, magic); 241 Vec4xF32 rounded = sub_f32x1(add_f32x1(src, magic), magic); 242 Vec4xF32 maybeone = cmp_gt_f32x1(src, rounded) & common_table.f32_1.as<Vec4xF32>(); 243 244 return cast_to_f32(blendv_bits(add_f32x1(rounded, maybeone), src, mask)); 245 } 246 247 static BL_INLINE double nearby(double x) noexcept { 248 using namespace SIMD; 249 250 Vec2xF64 src = cast_from_f64(x); 251 Vec2xF64 mask = cmp_lt_f64x1(src, common_table.f64_round_magic.as<Vec2xF64>()); 252 Vec2xF64 magic = mask & (common_table.f64_round_magic.as<Vec2xF64>() | slli_u64<64 - 13>(srli_u64<63>(src))); 253 Vec2xF64 rounded = sub_f64x1(add_f64x1(src, magic), magic); 254 255 return cast_to_f64(rounded); 256 } 257 258 static BL_INLINE double trunc(double x) noexcept { 259 using namespace SIMD; 260 261 Vec2xF64 src = cast_from_f64(x); 262 Vec2xF64 msk_abs = common_table.p_7FFFFFFFFFFFFFFF.as<Vec2xF64>(); 263 Vec2xF64 src_abs = src & msk_abs; 264 265 Vec2xF64 sign = andnot(msk_abs, src); 266 Vec2xF64 magic = common_table.f64_round_magic.as<Vec2xF64>(); 267 268 Vec2xF64 mask = cmp_ge_f64x1(src_abs, magic) | sign; 269 Vec2xF64 rounded = sub_f64x1(add_f64x1(src_abs, magic), magic); 270 Vec2xF64 maybeone = cmp_lt_f64x1(src_abs, rounded) & common_table.f64_1.as<Vec2xF64>(); 271 272 return cast_to_f64(blendv_bits(sub_f64x1(rounded, maybeone), src, mask)); 273 } 274 275 static BL_INLINE double floor(double x) noexcept { 276 using namespace SIMD; 277 278 Vec2xF64 src = cast_from_f64(x); 279 Vec2xF64 magic = slli_u64<64 - 13>(srli_u64<63>(src)) | common_table.f64_round_magic.as<Vec2xF64>(); 280 281 Vec2xF64 mask = cmp_ge_f64x1(src, magic); 282 Vec2xF64 rounded = sub_f64x1(add_f64x1(src, magic), magic); 283 Vec2xF64 maybeone = cmp_lt_f64x1(src, rounded) & common_table.f64_1.as<Vec2xF64>(); 284 285 return cast_to_f64(blendv_bits(sub_f64x1(rounded, maybeone), src, mask)); 286 } 287 288 static BL_INLINE double ceil(double x) noexcept { 289 using namespace SIMD; 290 291 Vec2xF64 src = cast_from_f64(x); 292 Vec2xF64 magic = slli_u64<64 - 13>(srli_u64<63>(src)) | common_table.f64_round_magic.as<Vec2xF64>(); 293 294 Vec2xF64 mask = cmp_ge_f64x1(src, magic); 295 Vec2xF64 rounded = sub_f64x1(add_f64x1(src, magic), magic); 296 Vec2xF64 maybeone = cmp_gt_f64x1(src, rounded) & common_table.f64_1.as<Vec2xF64>(); 297 298 return cast_to_f64(blendv_bits(add_f64x1(rounded, maybeone), src, mask)); 299 } 300 301 #else 302 303 BL_PRAGMA_FAST_MATH_PUSH 304 305 static BL_INLINE_NODEBUG float nearby(float x) noexcept { return ::rintf(x); } 306 static BL_INLINE_NODEBUG double nearby(double x) noexcept { return ::rint(x); } 307 static BL_INLINE_NODEBUG float trunc(float x) noexcept { return ::truncf(x); } 308 static BL_INLINE_NODEBUG double trunc(double x) noexcept { return ::trunc(x); } 309 static BL_INLINE_NODEBUG float floor(float x) noexcept { return ::floorf(x); } 310 static BL_INLINE_NODEBUG double floor(double x) noexcept { return ::floor(x); } 311 static BL_INLINE_NODEBUG float ceil(float x) noexcept { return ::ceilf(x); } 312 static BL_INLINE_NODEBUG double ceil(double x) noexcept { return ::ceil(x); } 313 314 BL_PRAGMA_FAST_MATH_POP 315 316 #endif 317 318 static BL_INLINE float round(float x) noexcept { float y = floor(x); return y + (x - y >= 0.5f ? 1.0f : 0.0f); } 319 static BL_INLINE double round(double x) noexcept { double y = floor(x); return y + (x - y >= 0.5 ? 1.0 : 0.0); } 320 321 //! \} 322 323 //! \name Rounding to Integer 324 //! \{ 325 326 static BL_INLINE int nearby_to_int(float x) noexcept { 327 #if defined(BL_TARGET_OPT_SSE) 328 return SIMD::cvt_f32_to_scalar_i32(SIMD::cast_from_f32(x)); 329 #else 330 return int(lrintf(x)); 331 #endif 332 } 333 334 static BL_INLINE int nearby_to_int(double x) noexcept { 335 #if defined(BL_TARGET_OPT_SSE2) 336 return SIMD::cvt_f64_to_scalar_i32(SIMD::cast_from_f64(x)); 337 #else 338 return int(lrint(x)); 339 #endif 340 } 341 342 static BL_INLINE int trunc_to_int(float x) noexcept { return int(x); } 343 static BL_INLINE int trunc_to_int(double x) noexcept { return int(x); } 344 345 static BL_INLINE BLBoxI trunc_to_int(const BLBox& box) noexcept { 346 return BLBoxI(trunc_to_int(box.x0), 347 trunc_to_int(box.y0), 348 trunc_to_int(box.x1), 349 trunc_to_int(box.y1)); 350 } 351 352 #if defined(BL_TARGET_OPT_SSE4_1) 353 static BL_INLINE int floor_to_int(float x) noexcept { return _mm_cvttss_si32(bl_roundf_sse4_1<_MM_FROUND_TO_NEG_INF>(x)); } 354 static BL_INLINE int floor_to_int(double x) noexcept { return _mm_cvttsd_si32(bl_roundd_sse4_1<_MM_FROUND_TO_NEG_INF>(x)); } 355 356 static BL_INLINE int ceil_to_int(float x) noexcept { return _mm_cvttss_si32(bl_roundf_sse4_1<_MM_FROUND_TO_POS_INF>(x)); } 357 static BL_INLINE int ceil_to_int(double x) noexcept { return _mm_cvttsd_si32(bl_roundd_sse4_1<_MM_FROUND_TO_POS_INF>(x)); } 358 #else 359 static BL_INLINE int floor_to_int(float x) noexcept { int y = nearby_to_int(x); return y - (float(y) > x); } 360 static BL_INLINE int floor_to_int(double x) noexcept { int y = nearby_to_int(x); return y - (double(y) > x); } 361 362 static BL_INLINE int ceil_to_int(float x) noexcept { int y = nearby_to_int(x); return y + (float(y) < x); } 363 static BL_INLINE int ceil_to_int(double x) noexcept { int y = nearby_to_int(x); return y + (double(y) < x); } 364 #endif 365 366 static BL_INLINE int round_to_int(float x) noexcept { int y = nearby_to_int(x); return y + (float(y) - x == -0.5f); } 367 static BL_INLINE int round_to_int(double x) noexcept { int y = nearby_to_int(x); return y + (double(y) - x == -0.5); } 368 369 static BL_INLINE int64_t nearby_to_int64(float x) noexcept { 370 #if BL_TARGET_ARCH_X86 == 64 371 return SIMD::cvt_f32_to_scalar_i64(SIMD::cast_from_f32(x)); 372 #elif BL_TARGET_ARCH_X86 == 32 && defined(__GNUC__) 373 int64_t y; 374 __asm__ __volatile__("flds %1\n" "fistpq %0\n" : "=m" (y) : "m" (x)); 375 return y; 376 #elif BL_TARGET_ARCH_X86 == 32 && defined(_MSC_VER) 377 int64_t y; 378 __asm { 379 fld dword ptr [x] 380 fistp qword ptr [y] 381 } 382 return y; 383 #else 384 return int64_t(llrintf(x)); 385 #endif 386 } 387 388 static BL_INLINE int64_t nearby_to_int64(double x) noexcept { 389 #if BL_TARGET_ARCH_X86 == 64 390 return SIMD::cvt_f64_to_scalar_i64(SIMD::cast_from_f64(x)); 391 #elif BL_TARGET_ARCH_X86 == 32 && defined(__GNUC__) 392 int64_t y; 393 __asm__ __volatile__("fldl %1\n" "fistpq %0\n" : "=m" (y) : "m" (x)); 394 return y; 395 #elif BL_TARGET_ARCH_X86 == 32 && defined(_MSC_VER) 396 int64_t y; 397 __asm { 398 fld qword ptr [x] 399 fistp qword ptr [y] 400 } 401 return y; 402 #else 403 return int64_t(llrint(x)); 404 #endif 405 } 406 407 static BL_INLINE_NODEBUG int64_t trunc_to_int64(float x) noexcept { return int64_t(x); } 408 static BL_INLINE_NODEBUG int64_t trunc_to_int64(double x) noexcept { return int64_t(x); } 409 410 static BL_INLINE_NODEBUG int64_t floor_to_int64(float x) noexcept { int64_t y = trunc_to_int64(x); return y - int64_t(float(y) > x); } 411 static BL_INLINE_NODEBUG int64_t floor_to_int64(double x) noexcept { int64_t y = trunc_to_int64(x); return y - int64_t(double(y) > x); } 412 413 static BL_INLINE_NODEBUG int64_t ceil_to_int64(float x) noexcept { int64_t y = trunc_to_int64(x); return y - int64_t(float(y) < x); } 414 static BL_INLINE_NODEBUG int64_t ceil_to_int64(double x) noexcept { int64_t y = trunc_to_int64(x); return y - int64_t(double(y) < x); } 415 416 static BL_INLINE_NODEBUG int64_t round_to_int64(float x) noexcept { int64_t y = nearby_to_int64(x); return y + int64_t(float(y) - x == -0.5f); } 417 static BL_INLINE_NODEBUG int64_t round_to_int64(double x) noexcept { int64_t y = nearby_to_int64(x); return y + int64_t(double(y) - x == -0.5); } 418 419 //! \} 420 421 //! \name Fraction & Repeat 422 //! \{ 423 424 //! Returns a fractional part of `x`. 425 //! 426 //! \note Fractional part returned is always equal or greater than zero. The 427 //! implementation is compatible to many shader implementations defined as 428 //! `frac(x) == x - floor(x)`, which would return `0.25` for `-1.75`. 429 template<typename T> 430 static BL_INLINE_NODEBUG T frac(T x) noexcept { return x - floor(x); } 431 432 //! Repeats the given value `x` in `y`, returning a value that is always equal 433 //! to or greater than zero and lesser than `y`. The return of `repeat(x, 1.0)` 434 //! should be identical to the return of `frac(x)`. 435 template<typename T> 436 static BL_INLINE T repeat(T x, T y) noexcept { 437 T a = x; 438 if (a >= y || a <= -y) 439 a = std::fmod(a, y); 440 if (a < T(0)) 441 a += y; 442 return a; 443 } 444 445 //! \} 446 447 //! \name Power Functions 448 //! \{ 449 450 template<typename T> BL_INLINE_CONSTEXPR T square(const T& x) noexcept { return x * x; } 451 template<typename T> BL_INLINE_CONSTEXPR T cube(const T& x) noexcept { return x * x * x; } 452 453 static BL_INLINE_NODEBUG float pow(float x, float y) noexcept { return ::powf(x, y); } 454 static BL_INLINE_NODEBUG double pow(double x, double y) noexcept { return ::pow(x, y); } 455 456 BL_PRAGMA_FAST_MATH_PUSH 457 458 static BL_INLINE_NODEBUG float sqrt(float x) noexcept { return ::sqrtf(x); } 459 static BL_INLINE_NODEBUG double sqrt(double x) noexcept { return ::sqrt(x); } 460 461 BL_PRAGMA_FAST_MATH_POP 462 463 static BL_INLINE_NODEBUG BLPoint sqrt(const BLPoint& p) noexcept { return BLPoint(sqrt(p.x), sqrt(p.y)); } 464 465 static BL_INLINE_NODEBUG float cbrt(float x) noexcept { return ::cbrtf(x); } 466 static BL_INLINE_NODEBUG double cbrt(double x) noexcept { return ::cbrt(x); } 467 468 static BL_INLINE_NODEBUG float hypot(float x, float y) noexcept { return ::hypotf(x, y); } 469 static BL_INLINE_NODEBUG double hypot(double x, double y) noexcept { return ::hypot(x, y); } 470 471 //! \} 472 473 //! \name Trigonometric Functions 474 //! \{ 475 476 static BL_INLINE_NODEBUG float sin(float x) noexcept { return ::sinf(x); } 477 static BL_INLINE_NODEBUG double sin(double x) noexcept { return ::sin(x); } 478 479 static BL_INLINE_NODEBUG float cos(float x) noexcept { return ::cosf(x); } 480 static BL_INLINE_NODEBUG double cos(double x) noexcept { return ::cos(x); } 481 482 static BL_INLINE_NODEBUG float tan(float x) noexcept { return ::tanf(x); } 483 static BL_INLINE_NODEBUG double tan(double x) noexcept { return ::tan(x); } 484 485 static BL_INLINE_NODEBUG float asin(float x) noexcept { return ::asinf(x); } 486 static BL_INLINE_NODEBUG double asin(double x) noexcept { return ::asin(x); } 487 488 static BL_INLINE_NODEBUG float acos(float x) noexcept { return ::acosf(x); } 489 static BL_INLINE_NODEBUG double acos(double x) noexcept { return ::acos(x); } 490 491 static BL_INLINE_NODEBUG float atan(float x) noexcept { return ::atanf(x); } 492 static BL_INLINE_NODEBUG double atan(double x) noexcept { return ::atan(x); } 493 494 static BL_INLINE_NODEBUG float atan2(float y, float x) noexcept { return ::atan2f(y, x); } 495 static BL_INLINE_NODEBUG double atan2(double y, double x) noexcept { return ::atan2(y, x); } 496 497 //! \} 498 499 //! \name Linear Interpolation 500 //! \{ 501 502 //! Linear interpolation of `a` and `b` at `t`. 503 //! 504 //! Returns `(a - t * a) + t * b`. 505 //! 506 //! \note This function should work with most geometric types Blend2D offers 507 //! that use double precision, however, it's not compatible with integral types. 508 template<typename V, typename T = double> 509 static BL_INLINE_NODEBUG V lerp(const V& a, const V& b, const T& t) noexcept { 510 return (a - t * a) + t * b; 511 } 512 513 // Linear interpolation of `a` and `b` at `t=0.5`. 514 template<typename T> 515 static BL_INLINE_NODEBUG T lerp(const T& a, const T& b) noexcept { 516 return 0.5 * a + 0.5 * b; 517 } 518 519 //! Alternative LERP implementation that is faster, but won't handle pathological 520 //! inputs. It should only be used in places in which it's known that such inputs 521 //! cannot happen. 522 template<typename V, typename T = double> 523 static BL_INLINE_NODEBUG V fast_lerp(const V& a, const V& b, const T& t) noexcept { 524 return a + t * (b - a); 525 } 526 527 //! Alternative LERP implementation at `t=0.5`. 528 template<typename T> 529 static BL_INLINE_NODEBUG T fast_lerp(const T& a, const T& b) noexcept { 530 return 0.5 * (a + b); 531 } 532 533 //! \} 534 535 //! \name Quadratic Roots 536 //! \{ 537 538 //! Solve a quadratic polynomial `Ax^2 + Bx + C = 0` and store the result in `dst`. 539 //! 540 //! Returns the number of roots found within [t_min, t_max] - `0` to `2`. 541 //! 542 //! Resources: 543 //! - http://stackoverflow.com/questions/4503849/quadratic-equation-in-ada/4504415#4504415 544 //! - http://people.csail.mit.edu/bkph/articles/Quadratics.pdf 545 //! 546 //! The standard equation: 547 //! 548 //! ``` 549 //! x0 = (-b + sqrt(delta)) / 2a 550 //! x1 = (-b - sqrt(delta)) / 2a 551 //! ``` 552 //! 553 //! When 4*a*c < b*b, computing x0 involves subtracting close numbers, and makes 554 //! you lose accuracy, so use the following instead: 555 //! 556 //! ``` 557 //! x0 = 2c / (-b - sqrt(delta)) 558 //! x1 = 2c / (-b + sqrt(delta)) 559 //! ``` 560 //! 561 //! Which yields a better x0, but whose x1 has the same problem as x0 had above. 562 //! The correct way to compute the roots is therefore: 563 //! 564 //! ``` 565 //! q = -0.5 * (b + sign(b) * sqrt(delta)) 566 //! x0 = q / a 567 //! x1 = c / q 568 //! ``` 569 //! 570 //! \note This is a branchless version designed to be easily inlinable. 571 static BL_INLINE size_t quad_roots(double dst[2], double a, double b, double c, double t_min, double t_max) noexcept { 572 double d = bl_max(b * b - 4.0 * a * c, 0.0); 573 double s = sqrt(d); 574 double q = -0.5 * (b + copy_sign(s, b)); 575 576 double t0 = q / a; 577 double t1 = c / q; 578 579 double x0 = bl_min(t0, t1); 580 double x1 = bl_max(t1, t0); 581 582 dst[0] = x0; 583 size_t n = size_t((x0 >= t_min) & (x0 <= t_max)); 584 585 dst[n] = x1; 586 n += size_t((x1 > x0) & (x1 >= t_min) & (x1 <= t_max)); 587 588 return n; 589 } 590 591 //! \overload 592 static BL_INLINE size_t quad_roots(double dst[2], const double poly[3], double t_min, double t_max) noexcept { 593 return quad_roots(dst, poly[0], poly[1], poly[2], t_min, t_max); 594 } 595 596 //! Like `quad_roots()`, but always returns two roots and doesn't sort them. 597 static BL_INLINE size_t simplified_quad_roots(double dst[2], double a, double b, double c) noexcept { 598 double d = bl_max(b * b - 4.0 * a * c, 0.0); 599 double s = sqrt(d); 600 double q = -0.5 * (b + copy_sign(s, b)); 601 602 dst[0] = q / a; 603 dst[1] = c / q; 604 return 2; 605 } 606 607 static BL_INLINE size_t simplified_quad_roots(BLPoint dst[2], const BLPoint& a, const BLPoint& b, const BLPoint& c) noexcept { 608 BLPoint d = bl_max(b * b - 4.0 * a * c, 0.0); 609 BLPoint s = sqrt(d); 610 BLPoint q = -0.5 * (b + copy_sign(s, b)); 611 612 dst[0] = q / a; 613 dst[1] = c / q; 614 615 return 2; 616 } 617 618 //! \} 619 620 } // {anonymous} 621 622 //! \name Cubic Roots 623 //! \{ 624 625 //! Solve a cubic polynomial and store the result in `dst`. 626 //! 627 //! Returns the number of roots found within [t_min, t_max] - `0` to `3`. 628 BL_HIDDEN size_t cubic_roots(double* dst, const double* poly, double t_min, double t_max) noexcept; 629 630 static BL_INLINE size_t cubic_roots(double dst[3], double a, double b, double c, double d, double t_min, double t_max) noexcept { 631 double poly[4] = { a, b, c, d }; 632 return cubic_roots(dst, poly, t_min, t_max); 633 } 634 635 //! \} 636 637 } // {Math} 638 } // {bl} 639 640 //! \} 641 //! \endcond 642 643 #endif // BLEND2D_SUPPORT_MATH_P_H_INCLUDED