1/***************************************************************************
2 * Copyright (c) Johan Mabille, Sylvain Corlay, Wolf Vollprecht and *
3 * Martin Renou *
4 * Copyright (c) QuantStack *
5 * Copyright (c) Serge Guelton *
6 * *
7 * Distributed under the terms of the BSD 3-Clause License. *
8 * *
9 * The full license is in the file LICENSE, distributed with this software. *
10 ****************************************************************************/
11
12#ifndef XSIMD_SCALAR_HPP
13#define XSIMD_SCALAR_HPP
14
15#include <cassert>
16#include <cmath>
17#include <complex>
18#include <cstdint>
19#include <cstring>
20#include <limits>
21#include <type_traits>
22
23#ifdef XSIMD_ENABLE_XTL_COMPLEX
24#include "xtl/xcomplex.hpp"
25#endif
26
27namespace xsimd
28{
29 template <class T, class A>
30 class batch;
31 template <class T, class A>
32 class batch_bool;
33
34 using std::abs;
35
36 using std::acos;
37 using std::acosh;
38 using std::arg;
39 using std::asin;
40 using std::asinh;
41 using std::atan;
42 using std::atan2;
43 using std::atanh;
44 using std::cbrt;
45 using std::ceil;
46 using std::conj;
47 using std::copysign;
48 using std::cos;
49 using std::cosh;
50 using std::erf;
51 using std::erfc;
52 using std::exp;
53 using std::exp2;
54 using std::expm1;
55 using std::fabs;
56 using std::fdim;
57 using std::floor;
58 using std::fmax;
59 using std::fmin;
60 using std::fmod;
61 using std::hypot;
62 using std::ldexp;
63 using std::lgamma;
64 using std::log;
65 using std::log10;
66 using std::log1p;
67 using std::log2;
68 using std::modf;
69 using std::nearbyint;
70 using std::nextafter;
71 using std::norm;
72 using std::polar;
73 using std::proj;
74 using std::remainder;
75 using std::rint;
76 using std::round;
77 using std::sin;
78 using std::sinh;
79 using std::sqrt;
80 using std::tan;
81 using std::tanh;
82 using std::tgamma;
83 using std::trunc;
84
85#ifndef _WIN32
86 using std::isfinite;
87 using std::isinf;
88 using std::isnan;
89#else
90
91 // Windows defines catch all templates
92 template <class T>
93 inline typename std::enable_if<std::is_floating_point<T>::value, bool>::type
94 isfinite(T var) noexcept
95 {
96 return std::isfinite(var);
97 }
98
99 template <class T>
100 inline typename std::enable_if<std::is_integral<T>::value, bool>::type
101 isfinite(T var) noexcept
102 {
103 return isfinite(double(var));
104 }
105
106 template <class T>
107 inline typename std::enable_if<std::is_floating_point<T>::value, bool>::type
108 isinf(T var) noexcept
109 {
110 return std::isinf(var);
111 }
112
113 template <class T>
114 inline typename std::enable_if<std::is_integral<T>::value, bool>::type
115 isinf(T var) noexcept
116 {
117 return isinf(double(var));
118 }
119
120 template <class T>
121 inline typename std::enable_if<std::is_floating_point<T>::value, bool>::type
122 isnan(T var) noexcept
123 {
124 return std::isnan(var);
125 }
126
127 template <class T>
128 inline typename std::enable_if<std::is_integral<T>::value, bool>::type
129 isnan(T var) noexcept
130 {
131 return isnan(double(var));
132 }
133#endif
134
135 template <class T, class Tp>
136 inline auto add(T const& x, Tp const& y) noexcept -> decltype(x + y)
137 {
138 return x + y;
139 }
140
141 template <class T>
142 inline typename std::enable_if<std::is_integral<T>::value, T>::type
143 bitwise_and(T x, T y) noexcept
144 {
145 return x & y;
146 }
147
148 inline float bitwise_and(float x, float y) noexcept
149 {
150 uint32_t ix, iy;
151 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(float));
152 std::memcpy(dest: (void*)&iy, src: (void*)&y, n: sizeof(float));
153 uint32_t ir = bitwise_and(x: ix, y: iy);
154 float r;
155 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(float));
156 return r;
157 }
158
159 inline double bitwise_and(double x, double y) noexcept
160 {
161 uint64_t ix, iy;
162 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(double));
163 std::memcpy(dest: (void*)&iy, src: (void*)&y, n: sizeof(double));
164 uint64_t ir = bitwise_and(x: ix, y: iy);
165 double r;
166 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(double));
167 return r;
168 }
169
170 template <class T>
171 inline typename std::enable_if<std::is_integral<T>::value, T>::type
172 bitwise_andnot(T x, T y) noexcept
173 {
174 return x & ~y;
175 }
176
177 inline float bitwise_andnot(float x, float y) noexcept
178 {
179 uint32_t ix, iy;
180 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(float));
181 std::memcpy(dest: (void*)&iy, src: (void*)&y, n: sizeof(float));
182 uint32_t ir = bitwise_andnot(x: ix, y: iy);
183 float r;
184 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(float));
185 return r;
186 }
187
188 inline double bitwise_andnot(double x, double y) noexcept
189 {
190 uint64_t ix, iy;
191 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(double));
192 std::memcpy(dest: (void*)&iy, src: (void*)&y, n: sizeof(double));
193 uint64_t ir = bitwise_andnot(x: ix, y: iy);
194 double r;
195 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(double));
196 return r;
197 }
198
199 template <class T>
200 inline typename std::enable_if<std::is_integral<T>::value, T>::type
201 bitwise_not(T x) noexcept
202 {
203 return ~x;
204 }
205
206 inline float bitwise_not(float x) noexcept
207 {
208 uint32_t ix;
209 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(float));
210 uint32_t ir = bitwise_not(x: ix);
211 float r;
212 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(float));
213 return r;
214 }
215
216 inline double bitwise_not(double x) noexcept
217 {
218 uint64_t ix;
219 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(double));
220 uint64_t ir = bitwise_not(x: ix);
221 double r;
222 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(double));
223 return r;
224 }
225
226 template <class T>
227 inline typename std::enable_if<std::is_integral<T>::value, T>::type
228 bitwise_or(T x, T y) noexcept
229 {
230 return x | y;
231 }
232
233 inline float bitwise_or(float x, float y) noexcept
234 {
235 uint32_t ix, iy;
236 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(float));
237 std::memcpy(dest: (void*)&iy, src: (void*)&y, n: sizeof(float));
238 uint32_t ir = bitwise_or(x: ix, y: iy);
239 float r;
240 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(float));
241 return r;
242 }
243
244 inline double bitwise_or(double x, double y) noexcept
245 {
246 uint64_t ix, iy;
247 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(double));
248 std::memcpy(dest: (void*)&iy, src: (void*)&y, n: sizeof(double));
249 uint64_t ir = bitwise_or(x: ix, y: iy);
250 double r;
251 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(double));
252 return r;
253 }
254
255 template <class T>
256 inline typename std::enable_if<std::is_integral<T>::value, T>::type
257 bitwise_xor(T x, T y) noexcept
258 {
259 return x ^ y;
260 }
261
262 inline float bitwise_xor(float x, float y) noexcept
263 {
264 uint32_t ix, iy;
265 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(float));
266 std::memcpy(dest: (void*)&iy, src: (void*)&y, n: sizeof(float));
267 uint32_t ir = bitwise_xor(x: ix, y: iy);
268 float r;
269 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(float));
270 return r;
271 }
272
273 inline double bitwise_xor(double x, double y) noexcept
274 {
275 uint64_t ix, iy;
276 std::memcpy(dest: (void*)&ix, src: (void*)&x, n: sizeof(double));
277 std::memcpy(dest: (void*)&iy, src: (void*)&y, n: sizeof(double));
278 uint64_t ir = bitwise_xor(x: ix, y: iy);
279 double r;
280 std::memcpy(dest: (void*)&r, src: (void*)&ir, n: sizeof(double));
281 return r;
282 }
283
284 template <class T, class Tp>
285 inline auto div(T const& x, Tp const& y) noexcept -> decltype(x / y)
286 {
287 return x / y;
288 }
289
290 template <class T, class Tp>
291 inline auto mod(T const& x, Tp const& y) noexcept -> decltype(x % y)
292 {
293 return x % y;
294 }
295
296 template <class T, class Tp>
297 inline auto mul(T const& x, Tp const& y) noexcept -> decltype(x * y)
298 {
299 return x * y;
300 }
301
302 template <class T>
303 inline auto neg(T const& x) noexcept -> decltype(-x)
304 {
305 return -x;
306 }
307
308 template <class T>
309 inline auto pos(T const& x) noexcept -> decltype(+x)
310 {
311 return +x;
312 }
313
314 inline float reciprocal(float const& x) noexcept
315 {
316 return 1.f / x;
317 }
318
319 inline double reciprocal(double const& x) noexcept
320 {
321 return 1. / x;
322 }
323
324#ifdef XSIMD_ENABLE_NUMPY_COMPLEX
325 template <class T>
326 inline bool isnan(std::complex<T> var) noexcept
327 {
328 return std::isnan(std::real(var)) || std::isnan(std::imag(var));
329 }
330
331 template <class T>
332 inline bool isinf(std::complex<T> var) noexcept
333 {
334 return std::isinf(std::real(var)) || std::isinf(std::imag(var));
335 }
336#endif
337
338#ifdef XSIMD_ENABLE_XTL_COMPLEX
339 using xtl::abs;
340 using xtl::acos;
341 using xtl::acosh;
342 using xtl::asin;
343 using xtl::asinh;
344 using xtl::atan;
345 using xtl::atanh;
346 using xtl::cos;
347 using xtl::cosh;
348 using xtl::exp;
349 using xtl::log;
350 using xtl::log10;
351 using xtl::norm;
352 using xtl::pow;
353 using xtl::proj;
354 using xtl::sin;
355 using xtl::sinh;
356 using xtl::sqrt;
357 using xtl::tan;
358 using xtl::tanh;
359#endif
360
361 template <typename T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
362 inline T clip(const T& val, const T& low, const T& hi) noexcept
363 {
364 assert(low <= hi && "ordered clipping bounds");
365 return low > val ? low : (hi < val ? hi : val);
366 }
367
368 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
369 inline bool is_flint(const T& x) noexcept
370 {
371 return std::isnan(x - x) ? false : (x - std::trunc(x)) == T(0);
372 }
373
374 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
375 inline bool is_even(const T& x) noexcept
376 {
377 return is_flint(x * T(0.5));
378 }
379
380 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
381 inline bool is_odd(const T& x) noexcept
382 {
383 return is_even(x - 1.);
384 }
385
386 inline int32_t nearbyint_as_int(float var) noexcept
387 {
388 return static_cast<int32_t>(std::nearbyint(x: var));
389 }
390
391 inline int64_t nearbyint_as_int(double var) noexcept
392 {
393 return static_cast<int64_t>(std::nearbyint(x: var));
394 }
395
396 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
397 inline bool eq(const T& x0, const T& x1) noexcept
398 {
399 return x0 == x1;
400 }
401
402 template <class T>
403 inline bool eq(const std::complex<T>& x0, const std::complex<T>& x1) noexcept
404 {
405 return x0 == x1;
406 }
407
408 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
409 inline bool ge(const T& x0, const T& x1) noexcept
410 {
411 return x0 >= x1;
412 }
413
414 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
415 inline bool gt(const T& x0, const T& x1) noexcept
416 {
417 return x0 > x1;
418 }
419
420 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
421 inline bool le(const T& x0, const T& x1) noexcept
422 {
423 return x0 <= x1;
424 }
425
426 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
427 inline bool lt(const T& x0, const T& x1) noexcept
428 {
429 return x0 < x1;
430 }
431
432 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
433 inline bool neq(const T& x0, const T& x1) noexcept
434 {
435 return x0 != x1;
436 }
437
438 template <class T>
439 inline bool neq(const std::complex<T>& x0, const std::complex<T>& x1) noexcept
440 {
441 return !(x0 == x1);
442 }
443
444#if defined(_GNU_SOURCE) && !defined(__APPLE__) && !defined(__MINGW32__) && !defined(__ANDROID__)
445 inline float exp10(const float& x) noexcept
446 {
447 return ::exp10f(x: x);
448 }
449 inline double exp10(const double& x) noexcept
450 {
451 return ::exp10(x: x);
452 }
453#endif
454
455 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
456 inline T exp10(const T& x) noexcept
457 {
458 // FIXME: very inefficient
459 return std::pow(T(10), x);
460 }
461
462 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
463 inline auto rsqrt(const T& x) noexcept -> decltype(std::sqrt(x))
464 {
465 using float_type = decltype(std::sqrt(x));
466 return static_cast<float_type>(1) / std::sqrt(x);
467 }
468
469 namespace detail
470 {
471 template <class C>
472 inline C expm1_complex_scalar_impl(const C& val) noexcept
473 {
474 using T = typename C::value_type;
475 T isin = std::sin(val.imag());
476 T rem1 = std::expm1(val.real());
477 T re = rem1 + T(1.);
478 T si = std::sin(val.imag() * T(0.5));
479 return std::complex<T>(rem1 - T(2.) * re * si * si, re * isin);
480 }
481 }
482
483 template <class T>
484 inline std::complex<T> expm1(const std::complex<T>& val) noexcept
485 {
486 return detail::expm1_complex_scalar_impl(val);
487 }
488
489#ifdef XSIMD_ENABLE_XTL_COMPLEX
490 template <class T, bool i3ec>
491 inline xtl::xcomplex<T, T, i3ec> expm1(const xtl::xcomplex<T, T, i3ec>& val) noexcept
492 {
493 return detail::expm1_complex_scalar_impl(val);
494 }
495#endif
496
497 namespace detail
498 {
499 template <class C>
500 inline C log1p_complex_scalar_impl(const C& val) noexcept
501 {
502 using T = typename C::value_type;
503 C u = C(1.) + val;
504 return u == C(1.) ? val : (u.real() <= T(0.) ? log(u) : log(u) * val / (u - C(1.)));
505 }
506 }
507
508 template <class T>
509 inline std::complex<T> log1p(const std::complex<T>& val) noexcept
510 {
511 return detail::log1p_complex_scalar_impl(val);
512 }
513
514 template <class T>
515 inline std::complex<T> log2(const std::complex<T>& val) noexcept
516 {
517 return log(val) / std::log(T(2));
518 }
519
520 template <typename T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
521 inline T sadd(const T& lhs, const T& rhs) noexcept
522 {
523 if (std::numeric_limits<T>::is_signed)
524 {
525 if ((lhs > 0) && (rhs > std::numeric_limits<T>::max() - lhs))
526 {
527 return std::numeric_limits<T>::max();
528 }
529 else if ((lhs < 0) && (rhs < std::numeric_limits<T>::lowest() - lhs))
530 {
531 return std::numeric_limits<T>::lowest();
532 }
533 else
534 {
535 return lhs + rhs;
536 }
537 }
538 else
539 {
540 if (rhs > std::numeric_limits<T>::max() - lhs)
541 {
542 return std::numeric_limits<T>::max();
543 }
544 else
545 {
546 return lhs + rhs;
547 }
548 }
549 }
550
551 template <typename T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
552 inline T ssub(const T& lhs, const T& rhs) noexcept
553 {
554 if (std::numeric_limits<T>::is_signed)
555 {
556 return sadd(lhs, (T)-rhs);
557 }
558 else
559 {
560 if (lhs < rhs)
561 {
562 return std::numeric_limits<T>::lowest();
563 }
564 else
565 {
566 return lhs - rhs;
567 }
568 }
569 }
570
571 namespace detail
572 {
573 template <class T>
574 struct value_type_or_type_helper
575 {
576 using type = T;
577 };
578 template <class T, class A>
579 struct value_type_or_type_helper<batch<T, A>>
580 {
581 using type = T;
582 };
583
584 template <class T>
585 using value_type_or_type = typename value_type_or_type_helper<T>::type;
586
587 template <class T0, class T1>
588 inline typename std::enable_if<std::is_integral<T1>::value, T0>::type
589 ipow(const T0& x, const T1& n) noexcept
590 {
591 static_assert(std::is_integral<T1>::value, "second argument must be an integer");
592 T0 a = x;
593 T1 b = n;
594 bool const recip = b < 0;
595 T0 r(static_cast<value_type_or_type<T0>>(1));
596 while (1)
597 {
598 if (b & 1)
599 {
600 r *= a;
601 }
602 b /= 2;
603 if (b == 0)
604 {
605 break;
606 }
607 a *= a;
608 }
609 return recip ? static_cast<T0>(1) / r : r;
610 }
611 }
612
613 template <class T0, class T1>
614 inline typename std::enable_if<std::is_integral<T1>::value, T0>::type
615 pow(const T0& x, const T1& n) noexcept
616 {
617 return detail::ipow(x, n);
618 }
619
620 template <class T0, class T1>
621 inline auto
622 pow(const T0& t0, const T1& t1) noexcept
623 -> typename std::enable_if<std::is_scalar<T0>::value && std::is_floating_point<T1>::value, decltype(std::pow(t0, t1))>::type
624 {
625 return std::pow(t0, t1);
626 }
627
628 template <class T0, class T1>
629 inline typename std::enable_if<std::is_integral<T1>::value, std::complex<T0>>::type
630 pow(const std::complex<T0>& t0, const T1& t1) noexcept
631 {
632 return detail::ipow(t0, t1);
633 }
634
635 template <class T0, class T1>
636 inline typename std::enable_if<!std::is_integral<T1>::value, std::complex<T0>>::type
637 pow(const std::complex<T0>& t0, const T1& t1) noexcept
638 {
639 return std::pow(t0, t1);
640 }
641
642 template <class T0, class T1>
643 inline auto
644 pow(const T0& t0, const std::complex<T1>& t1) noexcept
645 -> typename std::enable_if<std::is_scalar<T0>::value, decltype(std::pow(t0, t1))>::type
646 {
647 return std::pow(t0, t1);
648 }
649
650 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
651 inline bool bitofsign(T const& x) noexcept
652 {
653 return x < T(0);
654 }
655
656 template <class T>
657 inline auto signbit(T const& v) noexcept -> decltype(bitofsign(v))
658 {
659 return bitofsign(v);
660 }
661
662 inline double sign(bool const& v) noexcept
663 {
664 return v;
665 }
666
667 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
668 inline T sign(const T& v) noexcept
669 {
670 return v < T(0) ? T(-1.) : v == T(0) ? T(0.)
671 : T(1.);
672 }
673
674 namespace detail
675 {
676 template <class C>
677 inline C sign_complex_scalar_impl(const C& v) noexcept
678 {
679 using value_type = typename C::value_type;
680 if (v.real())
681 {
682 return C(sign(v.real()), value_type(0));
683 }
684 else
685 {
686 return C(sign(v.imag()), value_type(0));
687 }
688 }
689 }
690
691 template <class T>
692 inline std::complex<T> sign(const std::complex<T>& v) noexcept
693 {
694 return detail::sign_complex_scalar_impl(v);
695 }
696
697#ifdef XSIMD_ENABLE_XTL_COMPLEX
698 template <class T, bool i3ec>
699 inline xtl::xcomplex<T, T, i3ec> sign(const xtl::xcomplex<T, T, i3ec>& v) noexcept
700 {
701 return detail::sign_complex_scalar_impl(v);
702 }
703#endif
704
705 inline double signnz(bool const&) noexcept
706 {
707 return 1;
708 }
709
710 template <class T, class = typename std::enable_if<std::is_scalar<T>::value>::type>
711 inline T signnz(const T& v) noexcept
712 {
713 return v < T(0) ? T(-1.) : T(1.);
714 }
715
716 template <class T, class Tp>
717 inline auto sub(T const& x, Tp const& y) noexcept -> decltype(x - y)
718 {
719 return x - y;
720 }
721
722#ifdef XSIMD_ENABLE_XTL_COMPLEX
723 template <class T, bool i3ec>
724 inline xtl::xcomplex<T, T, i3ec> log2(const xtl::xcomplex<T, T, i3ec>& val) noexcept
725 {
726 return log(val) / log(T(2));
727 }
728#endif
729
730#ifdef XSIMD_ENABLE_XTL_COMPLEX
731 template <class T, bool i3ec>
732 inline xtl::xcomplex<T, T, i3ec> log1p(const xtl::xcomplex<T, T, i3ec>& val) noexcept
733 {
734 return detail::log1p_complex_scalar_impl(val);
735 }
736#endif
737
738 template <class T0, class T1>
739 inline auto min(T0 const& self, T1 const& other) noexcept
740 -> typename std::enable_if<std::is_scalar<T0>::value && std::is_scalar<T1>::value,
741 typename std::decay<decltype(self > other ? other : self)>::type>::type
742 {
743 return self > other ? other : self;
744 }
745
746 // numpy defines minimum operator on complex using lexical comparison
747 template <class T0, class T1>
748 inline std::complex<typename std::common_type<T0, T1>::type>
749 min(std::complex<T0> const& self, std::complex<T1> const& other) noexcept
750 {
751 return (self.real() < other.real()) ? (self) : (self.real() == other.real() ? (self.imag() < other.imag() ? self : other) : other);
752 }
753
754 template <class T0, class T1>
755 inline auto max(T0 const& self, T1 const& other) noexcept
756 -> typename std::enable_if<std::is_scalar<T0>::value && std::is_scalar<T1>::value,
757 typename std::decay<decltype(self > other ? other : self)>::type>::type
758 {
759 return self < other ? other : self;
760 }
761
762 // numpy defines maximum operator on complex using lexical comparison
763 template <class T0, class T1>
764 inline std::complex<typename std::common_type<T0, T1>::type>
765 max(std::complex<T0> const& self, std::complex<T1> const& other) noexcept
766 {
767 return (self.real() > other.real()) ? (self) : (self.real() == other.real() ? (self.imag() > other.imag() ? self : other) : other);
768 }
769
770 template <class T>
771 inline typename std::enable_if<std::is_integral<T>::value, T>::type fma(const T& a, const T& b, const T& c) noexcept
772 {
773 return a * b + c;
774 }
775
776 template <class T>
777 inline typename std::enable_if<std::is_floating_point<T>::value, T>::type fma(const T& a, const T& b, const T& c) noexcept
778 {
779 return std::fma(a, b, c);
780 }
781
782 template <class T>
783 inline typename std::enable_if<std::is_scalar<T>::value, T>::type fms(const T& a, const T& b, const T& c) noexcept
784 {
785 return a * b - c;
786 }
787
788 namespace detail
789 {
790 template <class C>
791 inline C fma_complex_scalar_impl(const C& a, const C& b, const C& c) noexcept
792 {
793 return { fms(a.real(), b.real(), fms(a.imag(), b.imag(), c.real())),
794 fma(a.real(), b.imag(), fma(a.imag(), b.real(), c.imag())) };
795 }
796 }
797
798 template <class T>
799 inline std::complex<T> fma(const std::complex<T>& a, const std::complex<T>& b, const std::complex<T>& c) noexcept
800 {
801 return detail::fma_complex_scalar_impl(a, b, c);
802 }
803
804#ifdef XSIMD_ENABLE_XTL_COMPLEX
805 template <class T, bool i3ec>
806 inline xtl::xcomplex<T, T, i3ec> fma(const xtl::xcomplex<T, T, i3ec>& a, const xtl::xcomplex<T, T, i3ec>& b, const xtl::xcomplex<T, T, i3ec>& c) noexcept
807 {
808 return detail::fma_complex_scalar_impl(a, b, c);
809 }
810#endif
811
812 namespace detail
813 {
814 template <class C>
815 inline C fms_complex_scalar_impl(const C& a, const C& b, const C& c) noexcept
816 {
817 return { fms(a.real(), b.real(), fma(a.imag(), b.imag(), c.real())),
818 fma(a.real(), b.imag(), fms(a.imag(), b.real(), c.imag())) };
819 }
820 }
821
822 template <class T>
823 inline std::complex<T> fms(const std::complex<T>& a, const std::complex<T>& b, const std::complex<T>& c) noexcept
824 {
825 return detail::fms_complex_scalar_impl(a, b, c);
826 }
827
828#ifdef XSIMD_ENABLE_XTL_COMPLEX
829 template <class T, bool i3ec>
830 inline xtl::xcomplex<T, T, i3ec> fms(const xtl::xcomplex<T, T, i3ec>& a, const xtl::xcomplex<T, T, i3ec>& b, const xtl::xcomplex<T, T, i3ec>& c) noexcept
831 {
832 return detail::fms_complex_scalar_impl(a, b, c);
833 }
834#endif
835
836 template <class T>
837 inline typename std::enable_if<std::is_integral<T>::value, T>::type fnma(const T& a, const T& b, const T& c) noexcept
838 {
839 return -(a * b) + c;
840 }
841
842 template <class T>
843 inline typename std::enable_if<std::is_floating_point<T>::value, T>::type fnma(const T& a, const T& b, const T& c) noexcept
844 {
845 return std::fma(-a, b, c);
846 }
847
848 namespace detail
849 {
850 template <class C>
851 inline C fnma_complex_scalar_impl(const C& a, const C& b, const C& c) noexcept
852 {
853 return { fms(a.imag(), b.imag(), fms(a.real(), b.real(), c.real())),
854 -fma(a.real(), b.imag(), fms(a.imag(), b.real(), c.imag())) };
855 }
856 }
857
858 template <class T>
859 inline std::complex<T> fnma(const std::complex<T>& a, const std::complex<T>& b, const std::complex<T>& c) noexcept
860 {
861 return detail::fnma_complex_scalar_impl(a, b, c);
862 }
863
864#ifdef XSIMD_ENABLE_XTL_COMPLEX
865 template <class T, bool i3ec>
866 inline xtl::xcomplex<T, T, i3ec> fnma(const xtl::xcomplex<T, T, i3ec>& a, const xtl::xcomplex<T, T, i3ec>& b, const xtl::xcomplex<T, T, i3ec>& c) noexcept
867 {
868 return detail::fnma_complex_scalar_impl(a, b, c);
869 }
870#endif
871
872 template <class T>
873 inline typename std::enable_if<std::is_integral<T>::value, T>::type fnms(const T& a, const T& b, const T& c) noexcept
874 {
875 return -(a * b) - c;
876 }
877
878 template <class T>
879 inline typename std::enable_if<std::is_floating_point<T>::value, T>::type fnms(const T& a, const T& b, const T& c) noexcept
880 {
881 return -std::fma(a, b, c);
882 }
883
884 namespace detail
885 {
886 template <class C>
887 inline C fnms_complex_scalar_impl(const C& a, const C& b, const C& c) noexcept
888 {
889 return { fms(a.imag(), b.imag(), fma(a.real(), b.real(), c.real())),
890 -fma(a.real(), b.imag(), fma(a.imag(), b.real(), c.imag())) };
891 }
892 }
893
894 template <class T>
895 inline std::complex<T> fnms(const std::complex<T>& a, const std::complex<T>& b, const std::complex<T>& c) noexcept
896 {
897 return detail::fnms_complex_scalar_impl(a, b, c);
898 }
899
900#ifdef XSIMD_ENABLE_XTL_COMPLEX
901 template <class T, bool i3ec>
902 inline xtl::xcomplex<T, T, i3ec> fnms(const xtl::xcomplex<T, T, i3ec>& a, const xtl::xcomplex<T, T, i3ec>& b, const xtl::xcomplex<T, T, i3ec>& c) noexcept
903 {
904 return detail::fnms_complex_scalar_impl(a, b, c);
905 }
906#endif
907
908 namespace detail
909 {
910#define XSIMD_HASSINCOS_TRAIT(func) \
911 template <class S> \
912 struct has##func \
913 { \
914 template <class T> \
915 static auto get(T* ptr) -> decltype(func(std::declval<T>(), std::declval<T*>(), std::declval<T*>()), std::true_type {}); \
916 static std::false_type get(...); \
917 static constexpr bool value = decltype(get((S*)nullptr))::value; \
918 }
919
920#define XSIMD_HASSINCOS(func, T) has##func<T>::value
921
922 XSIMD_HASSINCOS_TRAIT(sincos);
923 XSIMD_HASSINCOS_TRAIT(sincosf);
924 XSIMD_HASSINCOS_TRAIT(__sincos);
925 XSIMD_HASSINCOS_TRAIT(__sincosf);
926
927 struct generic_sincosf
928 {
929 template <class T>
930 typename std::enable_if<XSIMD_HASSINCOS(sincosf, T), void>::type
931 operator()(float val, T& s, T& c)
932 {
933 sincosf(val, &s, &c);
934 }
935
936 template <class T>
937 typename std::enable_if<!XSIMD_HASSINCOS(sincosf, T) && XSIMD_HASSINCOS(__sincosf, T), void>::type
938 operator()(float val, T& s, T& c)
939 {
940 __sincosf(val, &s, &c);
941 }
942
943 template <class T>
944 typename std::enable_if<!XSIMD_HASSINCOS(sincosf, T) && !XSIMD_HASSINCOS(__sincosf, T), void>::type
945 operator()(float val, T& s, T& c)
946 {
947 s = std::sin(x: val);
948 c = std::cos(x: val);
949 }
950 };
951
952 struct generic_sincos
953 {
954 template <class T>
955 typename std::enable_if<XSIMD_HASSINCOS(sincos, T), void>::type
956 operator()(double val, T& s, T& c)
957 {
958 sincos(val, &s, &c);
959 }
960
961 template <class T>
962 typename std::enable_if<!XSIMD_HASSINCOS(sincos, T) && XSIMD_HASSINCOS(__sincos, T), void>::type
963 operator()(double val, T& s, T& c)
964 {
965 __sincos(val, &s, &c);
966 }
967
968 template <class T>
969 typename std::enable_if<!XSIMD_HASSINCOS(sincos, T) && !XSIMD_HASSINCOS(__sincos, T), void>::type
970 operator()(double val, T& s, T& c)
971 {
972 s = std::sin(x: val);
973 c = std::cos(x: val);
974 }
975 };
976
977#undef XSIMD_HASSINCOS_TRAIT
978#undef XSIMD_HASSINCOS
979 }
980
981 inline std::pair<float, float> sincos(float val) noexcept
982 {
983 float s, c;
984 detail::generic_sincosf {}(val, s, c);
985 return std::make_pair(x&: s, y&: c);
986 }
987
988 inline std::pair<double, double> sincos(double val) noexcept
989 {
990 double s, c;
991 detail::generic_sincos {}(val, s, c);
992 return std::make_pair(x&: s, y&: c);
993 }
994
995 template <class T>
996 inline std::pair<std::complex<T>, std::complex<T>>
997 sincos(const std::complex<T>& val) noexcept
998 {
999 return std::make_pair(std::sin(val), std::cos(val));
1000 }
1001
1002#ifdef XSIMD_ENABLE_XTL_COMPLEX
1003 template <class T>
1004 inline std::pair<xtl::xcomplex<T>, xtl::xcomplex<T>> sincos(const xtl::xcomplex<T>& val) noexcept
1005 {
1006 return std::make_pair(sin(val), cos(val));
1007 }
1008#endif
1009
1010 template <class T, class _ = typename std::enable_if<std::is_floating_point<T>::value, void>::type>
1011 inline T frexp(T const& val, int& exp) noexcept
1012 {
1013 return std::frexp(val, &exp);
1014 }
1015
1016 template <class T>
1017 inline T select(bool cond, T const& true_br, T const& false_br) noexcept
1018 {
1019 return cond ? true_br : false_br;
1020 }
1021
1022}
1023
1024#endif
1025