Eigen  5.0.1
 
Loading...
Searching...
No Matches
MathFunctions.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2006-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
5// Copyright (c) 2021, NVIDIA CORPORATION. All rights reserved.
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10// SPDX-License-Identifier: MPL-2.0
11
12#ifndef EIGEN_MATHFUNCTIONS_H
13#define EIGEN_MATHFUNCTIONS_H
14
15// TODO: consider moving these constants to NumTraits.
16// Source: WolframAlpha
17#define EIGEN_PI 3.141592653589793238462643383279502884197169399375105820974944592307816406L
18#define EIGEN_LOG2E 1.442695040888963407359924681001892137426645954152985934135449406931109219L
19#define EIGEN_LN2 0.693147180559945309417232121458176568075500134360255254120680009493393621L
20
21// IWYU pragma: private
22#include "./InternalHeaderCheck.h"
23
24namespace Eigen {
25
26namespace internal {
27
48
49template <typename T, typename dummy = void>
50struct global_math_functions_filtering_base {
51 using type = T;
52};
53
54template <typename T>
55struct global_math_functions_filtering_base<T,
56 void_t<typename T::Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl>> {
57 using type = typename T::Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl;
58};
59
60#define EIGEN_MATHFUNC_IMPL(func, scalar) \
61 Eigen::internal::func##_impl<typename Eigen::internal::global_math_functions_filtering_base<scalar>::type>
62
63/****************************************************************************
64 * Implementation of real *
65 ****************************************************************************/
66
67template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
68struct real_default_impl {
69 using RealScalar = typename NumTraits<Scalar>::Real;
70 EIGEN_DEVICE_FUNC static constexpr RealScalar run(const Scalar& x) { return x; }
71};
72
73template <typename Scalar>
74struct real_default_impl<Scalar, true> {
75 using RealScalar = typename NumTraits<Scalar>::Real;
76 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
77 using std::real;
78 return real(x);
79 }
80};
81
82template <typename Scalar>
83struct real_impl : real_default_impl<Scalar> {};
84
85#if defined(EIGEN_GPU_COMPILE_PHASE)
86template <typename T>
87struct real_impl<std::complex<T>> {
88 using RealScalar = T;
89 EIGEN_DEVICE_FUNC static inline T run(const std::complex<T>& x) { return x.real(); }
90};
91#endif
92
93/****************************************************************************
94 * Implementation of imag *
95 ****************************************************************************/
96
97template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
98struct imag_default_impl {
99 using RealScalar = typename NumTraits<Scalar>::Real;
100 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar&) { return RealScalar(0); }
101};
102
103template <typename Scalar>
104struct imag_default_impl<Scalar, true> {
105 using RealScalar = typename NumTraits<Scalar>::Real;
106 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
107 using std::imag;
108 return imag(x);
109 }
110};
111
112template <typename Scalar>
113struct imag_impl : imag_default_impl<Scalar> {};
114
115#if defined(EIGEN_GPU_COMPILE_PHASE)
116template <typename T>
117struct imag_impl<std::complex<T>> {
118 using RealScalar = T;
119 EIGEN_DEVICE_FUNC static inline T run(const std::complex<T>& x) { return x.imag(); }
120};
121#endif
122
123/****************************************************************************
124 * Implementation of real_ref *
125 ****************************************************************************/
126
127template <typename Scalar>
128struct real_ref_impl {
129 using RealScalar = typename NumTraits<Scalar>::Real;
130 EIGEN_DEVICE_FUNC static inline RealScalar& run(Scalar& x) { return reinterpret_cast<RealScalar*>(&x)[0]; }
131 EIGEN_DEVICE_FUNC static inline const RealScalar& run(const Scalar& x) {
132 return reinterpret_cast<const RealScalar*>(&x)[0];
133 }
134};
135
136/****************************************************************************
137 * Implementation of imag_ref *
138 ****************************************************************************/
139
140template <typename Scalar, bool IsComplex>
141struct imag_ref_default_impl {
142 using RealScalar = typename NumTraits<Scalar>::Real;
143 EIGEN_DEVICE_FUNC static inline RealScalar& run(Scalar& x) { return reinterpret_cast<RealScalar*>(&x)[1]; }
144 EIGEN_DEVICE_FUNC static inline const RealScalar& run(const Scalar& x) {
145 return reinterpret_cast<const RealScalar*>(&x)[1];
146 }
147};
148
149template <typename Scalar>
150struct imag_ref_default_impl<Scalar, false> {
151 using RealScalar = typename NumTraits<Scalar>::Real;
152 EIGEN_DEVICE_FUNC constexpr static inline RealScalar run(Scalar&) { return RealScalar(0); }
153 EIGEN_DEVICE_FUNC constexpr static inline RealScalar run(const Scalar&) { return RealScalar(0); }
154};
155
156template <typename Scalar>
157struct imag_ref_impl : imag_ref_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
158
159} // namespace internal
160
161namespace numext {
162
163template <typename Scalar>
164EIGEN_DEVICE_FUNC inline decltype(auto) real(const Scalar& x) {
165 return EIGEN_MATHFUNC_IMPL(real, Scalar)::run(x);
166}
167
168template <typename Scalar>
169EIGEN_DEVICE_FUNC inline auto real_ref(const Scalar& x)
170 -> internal::add_const_on_value_type_t<decltype(internal::real_ref_impl<Scalar>::run(x))> {
171 return internal::real_ref_impl<Scalar>::run(x);
172}
173
174template <typename Scalar>
175EIGEN_DEVICE_FUNC inline auto real_ref(Scalar& x)
176 -> decltype(EIGEN_MATHFUNC_IMPL(real_ref, std::remove_const_t<Scalar>)::run(x)) {
177 return EIGEN_MATHFUNC_IMPL(real_ref, std::remove_const_t<Scalar>)::run(x);
178}
179
180template <typename Scalar>
181EIGEN_DEVICE_FUNC inline decltype(auto) imag(const Scalar& x) {
182 return EIGEN_MATHFUNC_IMPL(imag, Scalar)::run(x);
183}
184
185template <typename Scalar>
186EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar select(const Scalar& mask, const Scalar& a, const Scalar& b) {
187 return numext::is_exactly_zero(mask) ? b : a;
188}
189
190} // namespace numext
191
192namespace internal {
193
194/****************************************************************************
195 * Implementation of conj *
196 ****************************************************************************/
197
198template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
199struct conj_default_impl {
200 EIGEN_DEVICE_FUNC static constexpr Scalar run(const Scalar& x) { return x; }
201};
202
203template <typename Scalar>
204struct conj_default_impl<Scalar, true> {
205 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& x) {
206 using std::conj;
207 return conj(x);
208 }
209};
210
211template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
212struct conj_impl : conj_default_impl<Scalar, IsComplex> {};
213
214/****************************************************************************
215 * Implementation of mul *
216 ****************************************************************************/
217
218// Unsigned operands narrower than int promote to int, where a product whose wrapped value is representable
219// can still overflow: (unsigned short)0xffff squared is 0xfffe0001 > INT_MAX. Multiplying in unsigned int,
220// wider than every such operand, wraps instead. Mixed operand types, signed types, and bool keep operator*.
221template <typename Lhs, typename Rhs,
222 bool NarrowUnsigned = std::is_same<Lhs, Rhs>::value && std::is_unsigned<Lhs>::value &&
223 !std::is_same<Lhs, bool>::value && (sizeof(Lhs) < sizeof(int))>
224struct mul_impl {
225 EIGEN_DEVICE_FUNC static constexpr EIGEN_ALWAYS_INLINE auto run(const Lhs& a, const Rhs& b) { return a * b; }
226};
227
228template <typename Scalar>
229struct mul_impl<Scalar, Scalar, true> {
230 EIGEN_DEVICE_FUNC static constexpr EIGEN_ALWAYS_INLINE Scalar run(const Scalar& a, const Scalar& b) {
231 return static_cast<Scalar>(static_cast<unsigned int>(a) * static_cast<unsigned int>(b));
232 }
233};
234
235// MSVC's STL evaluates T * complex<T> as complex<T>(a) * b, whose 0 * inf terms turn 2 * (inf, 2) into
236// (inf, NaN). Scale component-wise, as complex<T> * T, libstdc++, libc++, and the GPU overloads do.
237template <typename RealScalar>
238struct mul_impl<RealScalar, std::complex<RealScalar>, false> {
239 EIGEN_DEVICE_FUNC static constexpr EIGEN_ALWAYS_INLINE std::complex<RealScalar> run(
240 const RealScalar& a, const std::complex<RealScalar>& b) {
241 return std::complex<RealScalar>(a * b.real(), a * b.imag());
242 }
243};
244
248template <typename Lhs, typename Rhs>
249EIGEN_DEVICE_FUNC constexpr EIGEN_ALWAYS_INLINE auto mul(const Lhs& a, const Rhs& b) {
250 return mul_impl<Lhs, Rhs>::run(a, b);
251}
252
253/****************************************************************************
254 * Implementation of abs2 *
255 ****************************************************************************/
256
257template <typename Scalar, bool IsComplex>
258struct abs2_impl_default {
259 using RealScalar = typename NumTraits<Scalar>::Real;
260 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) { return internal::mul(x, x); }
261};
262
263template <typename Scalar>
264struct abs2_impl_default<Scalar, true> // IsComplex
265{
266 using RealScalar = typename NumTraits<Scalar>::Real;
267 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
268 return numext::real(x) * numext::real(x) + numext::imag(x) * numext::imag(x);
269 }
270};
271
272template <typename Scalar>
273struct abs2_impl {
274 using RealScalar = typename NumTraits<Scalar>::Real;
275 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
276 return abs2_impl_default<Scalar, NumTraits<Scalar>::IsComplex>::run(x);
277 }
278};
279
280/****************************************************************************
281 * Implementation of sqrt/rsqrt *
282 ****************************************************************************/
283
284template <typename Scalar>
285struct sqrt_impl {
286 EIGEN_DEVICE_FUNC static EIGEN_ALWAYS_INLINE Scalar run(const Scalar& x) {
287 EIGEN_USING_STD(sqrt);
288 return sqrt(x);
289 }
290};
291
292// Complex sqrt defined in MathFunctionsImpl.h.
293template <typename ComplexT>
294EIGEN_DEVICE_FUNC constexpr ComplexT complex_sqrt(const ComplexT& a_x);
295
296// Custom implementation is faster than `std::sqrt`, works on
297// GPU, and correctly handles special cases (unlike MSVC).
298template <typename T>
299struct sqrt_impl<std::complex<T>> {
300 EIGEN_DEVICE_FUNC static EIGEN_ALWAYS_INLINE std::complex<T> run(const std::complex<T>& x) { return complex_sqrt(x); }
301};
302
303// Default implementation relies on numext::sqrt, at bottom of file.
304template <typename T>
305struct rsqrt_impl;
306
307// Complex rsqrt defined in MathFunctionsImpl.h.
308template <typename ComplexT>
309EIGEN_DEVICE_FUNC constexpr ComplexT complex_rsqrt(const ComplexT& a_x);
310
311template <typename T>
312struct rsqrt_impl<std::complex<T>> {
313 EIGEN_DEVICE_FUNC static EIGEN_ALWAYS_INLINE std::complex<T> run(const std::complex<T>& x) {
314 return complex_rsqrt(x);
315 }
316};
317
318/****************************************************************************
319 * Implementation of norm1 *
320 ****************************************************************************/
321
322template <typename Scalar, bool IsComplex>
323struct norm1_default_impl;
324
325template <typename Scalar>
326struct norm1_default_impl<Scalar, true> {
327 using RealScalar = typename NumTraits<Scalar>::Real;
328 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
329 EIGEN_USING_STD(abs);
330 return abs(numext::real(x)) + abs(numext::imag(x));
331 }
332};
333
334template <typename Scalar>
335struct norm1_default_impl<Scalar, false> {
336 using RealScalar = typename NumTraits<Scalar>::Real;
337 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
338 EIGEN_USING_STD(abs);
339 return abs(x);
340 }
341};
342
343template <typename Scalar>
344struct norm1_impl : norm1_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
345
346/****************************************************************************
347 * Implementation of hypot *
348 ****************************************************************************/
349
350template <typename Scalar>
351struct hypot_impl;
352
353/****************************************************************************
354 * Implementation of cast *
355 ****************************************************************************/
356
357template <typename OldType, typename NewType, typename EnableIf = void>
358struct cast_impl {
359 EIGEN_DEVICE_FUNC static inline NewType run(const OldType& x) { return static_cast<NewType>(x); }
360};
361
362template <typename OldType>
363struct cast_impl<OldType, bool> {
364 EIGEN_DEVICE_FUNC static inline bool run(const OldType& x) { return x != OldType(0); }
365};
366
367// Casting from S -> Complex<T> leads to an implicit conversion from S to T,
368// generating warnings on clang. Here we explicitly cast the real component.
369template <typename OldType, typename NewType>
370struct cast_impl<OldType, NewType, std::enable_if_t<!NumTraits<OldType>::IsComplex && NumTraits<NewType>::IsComplex>> {
371 EIGEN_DEVICE_FUNC static inline NewType run(const OldType& x) {
372 using NewReal = typename NumTraits<NewType>::Real;
373 return static_cast<NewType>(static_cast<NewReal>(x));
374 }
375};
376
377// Returns NewType directly to avoid unintended intermediate conversions.
378
379template <typename OldType, typename NewType>
380EIGEN_DEVICE_FUNC inline NewType cast(const OldType& x) {
381 return cast_impl<OldType, NewType>::run(x);
382}
383
384/****************************************************************************
385 * Implementation of arg *
386 ****************************************************************************/
387
388// Visual Studio 2017 has a bug where arg(float) returns 0 for negative inputs.
389// This seems to be fixed in VS 2019.
390#if (!EIGEN_COMP_MSVC || EIGEN_COMP_MSVC >= 1920)
391// std::arg is only defined for types of std::complex, or integer types or float/double/long double
392template <typename Scalar, bool HasStdImpl = NumTraits<Scalar>::IsComplex || std::is_integral<Scalar>::value ||
393 std::is_same<Scalar, float>::value ||
394 std::is_same<Scalar, double>::value ||
395 std::is_same<Scalar, long double>::value>
396struct arg_default_impl;
397
398template <typename Scalar>
399struct arg_default_impl<Scalar, true> {
400 using RealScalar = typename NumTraits<Scalar>::Real;
401 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
402 // There is no official ::arg on device in CUDA/HIP, so we always need to use std::arg.
403 using std::arg;
404 return static_cast<RealScalar>(arg(x));
405 }
406};
407
408// Must be non-complex floating-point type (e.g. half/bfloat16).
409template <typename Scalar>
410struct arg_default_impl<Scalar, false> {
411 using RealScalar = typename NumTraits<Scalar>::Real;
412 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
413 return (x < Scalar(0)) ? RealScalar(EIGEN_PI) : RealScalar(0);
414 }
415};
416#else
417template <typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
418struct arg_default_impl {
419 using RealScalar = typename NumTraits<Scalar>::Real;
420 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
421 return (x < RealScalar(0)) ? RealScalar(EIGEN_PI) : RealScalar(0);
422 }
423};
424
425template <typename Scalar>
426struct arg_default_impl<Scalar, true> {
427 using RealScalar = typename NumTraits<Scalar>::Real;
428 EIGEN_DEVICE_FUNC static inline RealScalar run(const Scalar& x) {
429 EIGEN_USING_STD(arg);
430 return arg(x);
431 }
432};
433#endif
434template <typename Scalar>
435struct arg_impl : arg_default_impl<Scalar> {};
436
437/****************************************************************************
438 * Implementation of expm1 *
439 ****************************************************************************/
440
441// This implementation is based on GSL Math's expm1.
442namespace std_fallback {
443// fallback expm1 implementation in case there is no expm1(Scalar) function in namespace of Scalar,
444// or that there is no suitable std::expm1 function available. Implementation
445// attributed to Kahan. See: http://www.plunk.org/~hatch/rightway.php.
446template <typename Scalar>
447EIGEN_DEVICE_FUNC inline Scalar expm1(const Scalar& x) {
448 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
449 using RealScalar = typename NumTraits<Scalar>::Real;
450
451 EIGEN_USING_STD(exp);
452 Scalar u = exp(x);
453 if (numext::equal_strict(u, Scalar(1))) {
454 return x;
455 }
456 Scalar um1 = u - RealScalar(1);
457 if (numext::equal_strict(um1, Scalar(-1))) {
458 return RealScalar(-1);
459 }
460
461 EIGEN_USING_STD(log);
462 Scalar logu = log(u);
463 return numext::equal_strict(u, logu) ? u : (u - RealScalar(1)) * x / logu;
464}
465} // namespace std_fallback
466
467template <typename Scalar>
468struct expm1_impl {
469 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& x) {
470 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
471 EIGEN_USING_STD(expm1);
472 return expm1(x);
473 }
474};
475
476/****************************************************************************
477 * Implementation of log *
478 ****************************************************************************/
479
480// Complex log defined in MathFunctionsImpl.h.
481template <typename ComplexT>
482EIGEN_DEVICE_FUNC constexpr ComplexT complex_log(const ComplexT& z);
483
484template <typename Scalar>
485struct log_impl {
486 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& x) {
487 EIGEN_USING_STD(log);
488 return static_cast<Scalar>(log(x));
489 }
490};
491
492template <typename Scalar>
493struct log_impl<std::complex<Scalar>> {
494 EIGEN_DEVICE_FUNC static inline std::complex<Scalar> run(const std::complex<Scalar>& z) { return complex_log(z); }
495};
496
497/****************************************************************************
498 * Implementation of log1p *
499 ****************************************************************************/
500
501namespace std_fallback {
502// fallback log1p implementation in case there is no log1p(Scalar) function in namespace of Scalar,
503// or that there is no suitable std::log1p function available
504template <typename Scalar>
505EIGEN_DEVICE_FUNC inline Scalar log1p(const Scalar& x) {
506 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
507 using RealScalar = typename NumTraits<Scalar>::Real;
508 EIGEN_USING_STD(log);
509 Scalar x1p = RealScalar(1) + x;
510 Scalar log_1p = log_impl<Scalar>::run(x1p);
511 const bool is_small = numext::equal_strict(x1p, Scalar(1));
512 // The correction below is 0 * inf or inf / inf once log(1 + x) is not finite; that log is the answer then.
513 EIGEN_USING_STD(isfinite);
514 if (!((isfinite)(numext::real(log_1p)) && (isfinite)(numext::imag(log_1p)))) return log_1p;
515 return is_small ? x : x * (log_1p / (x1p - RealScalar(1)));
516}
517} // namespace std_fallback
518
519template <typename Scalar>
520struct log1p_impl {
521 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
522
523 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& x) {
524 EIGEN_USING_STD(log1p);
525 return log1p(x);
526 }
527};
528
529// Specialization for complex types that are not supported by std::log1p.
530template <typename RealScalar>
531struct log1p_impl<std::complex<RealScalar>> {
532 EIGEN_STATIC_ASSERT_NON_INTEGER(RealScalar)
533
534 EIGEN_DEVICE_FUNC static inline std::complex<RealScalar> run(const std::complex<RealScalar>& x) {
535 return std_fallback::log1p(x);
536 }
537};
538
539// numext::exp2 of a complex argument; defined after the numext functions it is built from.
540template <typename RealScalar>
541struct complex_exp2_impl;
542
543/****************************************************************************
544 * Implementation of pow *
545 ****************************************************************************/
546
547template <typename ScalarX, typename ScalarY,
548 bool IsInteger = NumTraits<ScalarX>::IsInteger && NumTraits<ScalarY>::IsInteger>
549struct pow_impl {
550 using result_type =
551 typename ScalarBinaryOpTraits<ScalarX, ScalarY, internal::scalar_pow_op<ScalarX, ScalarY>>::ReturnType;
552 static EIGEN_DEVICE_FUNC inline result_type run(const ScalarX& x, const ScalarY& y) {
553 EIGEN_USING_STD(pow);
554 return pow(x, y);
555 }
556};
557
558template <typename ScalarX, typename ScalarY>
559struct pow_impl<ScalarX, ScalarY, true> {
560 using result_type = ScalarX;
561 static EIGEN_DEVICE_FUNC inline ScalarX run(ScalarX x, ScalarY y) {
562 ScalarX res(1);
563 eigen_assert(!NumTraits<ScalarY>::IsSigned || y >= 0);
564 if (y & 1) res = internal::mul(res, x);
565 y >>= 1;
566 while (y) {
567 x = internal::mul(x, x);
568 if (y & 1) res = internal::mul(res, x);
569 y >>= 1;
570 }
571 return res;
572 }
573};
574
575constexpr int floor_log2(unsigned int value) {
576 int result = 0;
577 while (value > 1) {
578 value >>= 1;
579 ++result;
580 }
581 return result;
582}
583
584template <typename BitsType, typename EnableIf = void>
585struct count_bits_impl {
586 static_assert(std::is_integral<BitsType>::value && std::is_unsigned<BitsType>::value,
587 "BitsType must be an unsigned integer");
588 static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) {
589 int n = CHAR_BIT * sizeof(BitsType);
590 int shift = n / 2;
591 while (bits > 0 && shift > 0) {
592 BitsType y = bits >> shift;
593 if (y > 0) {
594 n -= shift;
595 bits = y;
596 }
597 shift /= 2;
598 }
599 if (shift == 0) {
600 --n;
601 }
602 return n;
603 }
604
605 static EIGEN_DEVICE_FUNC inline int ctz(BitsType bits) {
606 int n = CHAR_BIT * sizeof(BitsType);
607 int shift = n / 2;
608 while (bits > 0 && shift > 0) {
609 BitsType y = bits << shift;
610 if (y > 0) {
611 n -= shift;
612 bits = y;
613 }
614 shift /= 2;
615 }
616 if (shift == 0) {
617 --n;
618 }
619 return n;
620 }
621
622 static EIGEN_DEVICE_FUNC inline int popcount(BitsType bits) {
623 int n = 0;
624 while (bits) {
625 bits &= bits - 1;
626 ++n;
627 }
628 return n;
629 }
630};
631
632// Count leading zeros.
633template <typename BitsType>
634EIGEN_DEVICE_FUNC inline int clz(BitsType bits) {
635 return count_bits_impl<BitsType>::clz(bits);
636}
637
638// Count trailing zeros.
639template <typename BitsType>
640EIGEN_DEVICE_FUNC inline int ctz(BitsType bits) {
641 return count_bits_impl<BitsType>::ctz(bits);
642}
643
644// Count set bits (population count).
645template <typename BitsType>
646EIGEN_DEVICE_FUNC inline int popcount(BitsType bits) {
647 return count_bits_impl<BitsType>::popcount(bits);
648}
649
650#if EIGEN_COMP_GNUC || EIGEN_COMP_CLANG
651
652template <typename BitsType>
653struct count_bits_impl<
654 BitsType, std::enable_if_t<std::is_integral<BitsType>::value && sizeof(BitsType) <= sizeof(unsigned int)>> {
655 static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT);
656 static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) {
657 static constexpr int kLeadingBitsOffset = (sizeof(unsigned int) - sizeof(BitsType)) * CHAR_BIT;
658 return bits == 0 ? kNumBits : __builtin_clz(static_cast<unsigned int>(bits)) - kLeadingBitsOffset;
659 }
660
661 static EIGEN_DEVICE_FUNC inline int ctz(BitsType bits) {
662 return bits == 0 ? kNumBits : __builtin_ctz(static_cast<unsigned int>(bits));
663 }
664
665 static EIGEN_DEVICE_FUNC inline int popcount(BitsType bits) {
666 return __builtin_popcount(static_cast<unsigned int>(bits));
667 }
668};
669
670template <typename BitsType>
671struct count_bits_impl<BitsType,
672 std::enable_if_t<std::is_integral<BitsType>::value && sizeof(unsigned int) < sizeof(BitsType) &&
673 sizeof(BitsType) <= sizeof(unsigned long)>> {
674 static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT);
675 static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) {
676 static constexpr int kLeadingBitsOffset = (sizeof(unsigned long) - sizeof(BitsType)) * CHAR_BIT;
677 return bits == 0 ? kNumBits : __builtin_clzl(static_cast<unsigned long>(bits)) - kLeadingBitsOffset;
678 }
679
680 static EIGEN_DEVICE_FUNC inline int ctz(BitsType bits) {
681 return bits == 0 ? kNumBits : __builtin_ctzl(static_cast<unsigned long>(bits));
682 }
683
684 static EIGEN_DEVICE_FUNC inline int popcount(BitsType bits) {
685 return __builtin_popcountl(static_cast<unsigned long>(bits));
686 }
687};
688
689template <typename BitsType>
690struct count_bits_impl<BitsType,
691 std::enable_if_t<std::is_integral<BitsType>::value && sizeof(unsigned long) < sizeof(BitsType) &&
692 sizeof(BitsType) <= sizeof(unsigned long long)>> {
693 static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT);
694 static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) {
695 static constexpr int kLeadingBitsOffset = (sizeof(unsigned long long) - sizeof(BitsType)) * CHAR_BIT;
696 return bits == 0 ? kNumBits : __builtin_clzll(static_cast<unsigned long long>(bits)) - kLeadingBitsOffset;
697 }
698
699 static EIGEN_DEVICE_FUNC inline int ctz(BitsType bits) {
700 return bits == 0 ? kNumBits : __builtin_ctzll(static_cast<unsigned long long>(bits));
701 }
702
703 static EIGEN_DEVICE_FUNC inline int popcount(BitsType bits) {
704 return __builtin_popcountll(static_cast<unsigned long long>(bits));
705 }
706};
707
708#elif EIGEN_COMP_MSVC
709
710// `__popcnt`/`__popcnt64` require `POPCNT` hardware support, which MSVC cannot guarantee
711// at its default baseline (unlike `_BitScanReverse`/`_BitScanForward`, which lower to
712// baseline `bsr`/`bsf`). Fall back to a portable count when building without SSE4.2
713// enabled.
714template <typename BitsType>
715EIGEN_DEVICE_FUNC inline int popcount_fallback(BitsType bits) {
716 int n = 0;
717 while (bits) {
718 bits &= bits - 1;
719 ++n;
720 }
721 return n;
722}
723
724template <typename BitsType>
725struct count_bits_impl<
726 BitsType, std::enable_if_t<std::is_integral<BitsType>::value && sizeof(BitsType) <= sizeof(unsigned long)>> {
727 static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT);
728 static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) {
729 unsigned long out;
730 _BitScanReverse(&out, static_cast<unsigned long>(bits));
731 return bits == 0 ? kNumBits : (kNumBits - 1) - static_cast<int>(out);
732 }
733
734 static EIGEN_DEVICE_FUNC inline int ctz(BitsType bits) {
735 unsigned long out;
736 _BitScanForward(&out, static_cast<unsigned long>(bits));
737 return bits == 0 ? kNumBits : static_cast<int>(out);
738 }
739
740 static EIGEN_DEVICE_FUNC inline int popcount(BitsType bits) {
741#if defined(EIGEN_VECTORIZE_SSE4_2)
742 return static_cast<int>(__popcnt(static_cast<unsigned int>(bits)));
743#else
744 return popcount_fallback(bits);
745#endif
746 }
747};
748
749#ifdef _WIN64
750
751template <typename BitsType>
752struct count_bits_impl<BitsType,
753 std::enable_if_t<std::is_integral<BitsType>::value && sizeof(unsigned long) < sizeof(BitsType) &&
754 sizeof(BitsType) <= sizeof(__int64)>> {
755 static constexpr int kNumBits = static_cast<int>(sizeof(BitsType) * CHAR_BIT);
756 static EIGEN_DEVICE_FUNC inline int clz(BitsType bits) {
757 unsigned long out;
758 _BitScanReverse64(&out, static_cast<unsigned __int64>(bits));
759 return bits == 0 ? kNumBits : (kNumBits - 1) - static_cast<int>(out);
760 }
761
762 static EIGEN_DEVICE_FUNC inline int ctz(BitsType bits) {
763 unsigned long out;
764 _BitScanForward64(&out, static_cast<unsigned __int64>(bits));
765 return bits == 0 ? kNumBits : static_cast<int>(out);
766 }
767
768 static EIGEN_DEVICE_FUNC inline int popcount(BitsType bits) {
769#if defined(EIGEN_VECTORIZE_SSE4_2)
770 return static_cast<int>(__popcnt64(static_cast<unsigned __int64>(bits)));
771#else
772 return popcount_fallback(bits);
773#endif
774 }
775};
776
777#endif // _WIN64
778
779#endif // EIGEN_COMP_GNUC || EIGEN_COMP_CLANG
780
781template <typename BitsType>
782struct log_2_impl {
783 static constexpr int kTotalBits = sizeof(BitsType) * CHAR_BIT;
784 static EIGEN_DEVICE_FUNC inline int run_ceil(const BitsType& x) {
785 const int n = kTotalBits - clz(x);
786 bool power_of_two = (x & (x - 1)) == 0;
787 return x == 0 ? 0 : power_of_two ? (n - 1) : n;
788 }
789 static EIGEN_DEVICE_FUNC inline int run_floor(const BitsType& x) {
790 const int n = kTotalBits - clz(x);
791 return x == 0 ? 0 : n - 1;
792 }
793};
794
795template <typename BitsType>
796int log2_ceil(const BitsType& x) {
797 return log_2_impl<BitsType>::run_ceil(x);
798}
799
800template <typename BitsType>
801int log2_floor(const BitsType& x) {
802 return log_2_impl<BitsType>::run_floor(x);
803}
804
805// Implementation of is* functions
806
807template <typename T>
808EIGEN_DEVICE_FUNC std::enable_if_t<!(std::numeric_limits<T>::has_infinity || std::numeric_limits<T>::has_quiet_NaN ||
809 std::numeric_limits<T>::has_signaling_NaN),
810 bool>
811isfinite_impl(const T&) {
812 return true;
813}
814
815template <typename T>
816EIGEN_DEVICE_FUNC std::enable_if_t<(std::numeric_limits<T>::has_infinity || std::numeric_limits<T>::has_quiet_NaN ||
817 std::numeric_limits<T>::has_signaling_NaN) &&
818 (!NumTraits<T>::IsComplex),
819 bool>
820isfinite_impl(const T& x) {
821 EIGEN_USING_STD(isfinite);
822 return isfinite EIGEN_NOT_A_MACRO(x);
823}
824
825template <typename T>
826EIGEN_DEVICE_FUNC std::enable_if_t<!std::numeric_limits<T>::has_infinity, bool> isinf_impl(const T&) {
827 return false;
828}
829
830template <typename T>
831EIGEN_DEVICE_FUNC std::enable_if_t<(std::numeric_limits<T>::has_infinity && !NumTraits<T>::IsComplex), bool> isinf_impl(
832 const T& x) {
833 EIGEN_USING_STD(isinf);
834 return isinf EIGEN_NOT_A_MACRO(x);
835}
836
837template <typename T>
838EIGEN_DEVICE_FUNC
839 std::enable_if_t<!(std::numeric_limits<T>::has_quiet_NaN || std::numeric_limits<T>::has_signaling_NaN), bool>
840 isnan_impl(const T&) {
841 return false;
842}
843
844template <typename T>
845EIGEN_DEVICE_FUNC std::enable_if_t<
846 (std::numeric_limits<T>::has_quiet_NaN || std::numeric_limits<T>::has_signaling_NaN) && (!NumTraits<T>::IsComplex),
847 bool>
848isnan_impl(const T& x) {
849 EIGEN_USING_STD(isnan);
850 return isnan EIGEN_NOT_A_MACRO(x);
851}
852
853// The following overload are defined at the end of this file
854template <typename T>
855EIGEN_DEVICE_FUNC bool isfinite_impl(const std::complex<T>& x);
856template <typename T>
857EIGEN_DEVICE_FUNC bool isnan_impl(const std::complex<T>& x);
858template <typename T>
859EIGEN_DEVICE_FUNC bool isinf_impl(const std::complex<T>& x);
860template <typename T>
861EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS T ptanh_float(const T& a_x);
862
863/****************************************************************************
864 * Implementation of sign *
865 ****************************************************************************/
866template <typename Scalar, bool IsComplex = (NumTraits<Scalar>::IsComplex != 0),
867 bool IsInteger = (NumTraits<Scalar>::IsInteger != 0)>
868struct sign_impl {
869 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& a) { return Scalar((a > Scalar(0)) - (a < Scalar(0))); }
870};
871
872template <typename Scalar>
873struct sign_impl<Scalar, false, false> {
874 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& a) {
875 return (isnan_impl<Scalar>)(a) ? a : Scalar((a > Scalar(0)) - (a < Scalar(0)));
876 }
877};
878
879template <typename Scalar, bool IsInteger>
880struct sign_impl<Scalar, true, IsInteger> {
881 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& a) {
882 using real_type = typename NumTraits<Scalar>::Real;
883 EIGEN_USING_STD(abs);
884 real_type aa = abs(a);
885 if (aa == real_type(0)) return Scalar(0);
886 // Divide rather than multiply by 1/aa: the reciprocal overflows for subnormal aa and is itself
887 // subnormal, hence inexact, for aa near the top of the range. psign_impl for complex packets
888 // divides for the same reason.
889 return Scalar(numext::real(a) / aa, numext::imag(a) / aa);
890 }
891};
892
893// The sign function for bool is the identity.
894template <>
895struct sign_impl<bool, false, true> {
896 EIGEN_DEVICE_FUNC static inline bool run(const bool& a) { return a; }
897};
898
899template <typename Scalar, bool IsComplex = (NumTraits<Scalar>::IsComplex != 0),
900 bool IsInteger = (NumTraits<Scalar>::IsInteger != 0)>
901struct copysign_impl {
902 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& a, const Scalar& b) {
903 EIGEN_USING_STD(copysign);
904 return Scalar(copysign(a, b));
905 }
906};
907
908template <typename Scalar, bool IsInteger>
909struct copysign_impl<Scalar, true, IsInteger> {
910 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& a, const Scalar& b) {
911 EIGEN_USING_STD(copysign);
912 return Scalar(copysign(numext::real(a), numext::real(b)), copysign(numext::imag(a), numext::imag(b)));
913 }
914};
915
916template <typename Scalar>
917struct copysign_impl<Scalar, false, true> {
918 EIGEN_DEVICE_FUNC static inline Scalar run(const Scalar& a, const Scalar& b) {
919 EIGEN_IF_CONSTEXPR (!NumTraits<Scalar>::IsSigned) return a;
920 const Scalar abs_a = a < Scalar(0) ? -a : a;
921 return b < Scalar(0) ? -abs_a : abs_a;
922 }
923};
924
925template <>
926struct copysign_impl<bool, false, true> {
927 EIGEN_DEVICE_FUNC static inline bool run(const bool& a, const bool&) { return a; }
928};
929
930// suppress "unary minus operator applied to unsigned type, result still unsigned" warnings on MSVC
931// note: `0 - a` is distinct from `-a` when Scalar is a floating point type and `a` is zero
932
933template <typename Scalar, bool IsInteger = NumTraits<Scalar>::IsInteger>
934struct negate_impl {
935 static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Scalar run(const Scalar& a) { return -a; }
936};
937
938template <typename Scalar>
939struct negate_impl<Scalar, true> {
940 EIGEN_STATIC_ASSERT((!std::is_same<Scalar, bool>::value), NEGATE IS NOT DEFINED FOR BOOLEAN TYPES)
941 static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Scalar run(const Scalar& a) { return Scalar(0) - a; }
942};
943
944template <typename Scalar, bool IsInteger = NumTraits<typename unpacket_traits<Scalar>::type>::IsInteger>
945struct nearest_integer_impl {
946 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_floor(const Scalar& x) {
947 EIGEN_USING_STD(floor) return floor(x);
948 }
949 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_ceil(const Scalar& x) {
950 EIGEN_USING_STD(ceil) return ceil(x);
951 }
952 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_rint(const Scalar& x) {
953 EIGEN_USING_STD(rint) return rint(x);
954 }
955 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_round(const Scalar& x) {
956 EIGEN_USING_STD(round) return round(x);
957 }
958 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_trunc(const Scalar& x) {
959 EIGEN_USING_STD(trunc) return trunc(x);
960 }
961};
962template <typename Scalar>
963struct nearest_integer_impl<Scalar, true> {
964 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_floor(const Scalar& x) { return x; }
965 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_ceil(const Scalar& x) { return x; }
966 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_rint(const Scalar& x) { return x; }
967 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_round(const Scalar& x) { return x; }
968 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run_trunc(const Scalar& x) { return x; }
969};
970
971// Extra namespace to prevent leaking std::fma into Eigen::internal.
972namespace has_fma_detail {
973
974template <typename T, typename EnableIf = void>
975struct has_fma_impl : public std::false_type {};
976
977using std::fma;
978
979template <typename T>
980struct has_fma_impl<
981 T, std::enable_if_t<std::is_same<T, decltype(fma(std::declval<T>(), std::declval<T>(), std::declval<T>()))>::value>>
982 : public std::true_type {};
983
984} // namespace has_fma_detail
985
986template <typename T>
987struct has_fma : public has_fma_detail::has_fma_impl<T> {};
988
989// Default implementation.
990template <typename T, typename Enable = void>
991struct fma_impl {
992 static_assert(has_fma<T>::value, "No function fma(...) for type. Please provide an implementation.");
993};
994
995// STD or ADL version if it exists.
996template <typename T>
997struct fma_impl<T, std::enable_if_t<has_fma<T>::value>> {
998 static EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T run(const T& a, const T& b, const T& c) {
999 using std::fma;
1000 return fma(a, b, c);
1001 }
1002};
1003
1004#if defined(EIGEN_GPUCC)
1005template <>
1006struct has_fma<float> : public std::true_type {};
1007
1008template <>
1009struct fma_impl<float, void> {
1010 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float run(const float& a, const float& b, const float& c) {
1011 return ::fmaf(a, b, c);
1012 }
1013};
1014
1015template <>
1016struct has_fma<double> : public std::true_type {};
1017
1018template <>
1019struct fma_impl<double, void> {
1020 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double run(const double& a, const double& b, const double& c) {
1021 return ::fma(a, b, c);
1022 }
1023};
1024#endif
1025
1026// Basic multiply-add.
1027template <typename Scalar, typename EnableIf = void>
1028struct madd_impl {
1029 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(const Scalar& x, const Scalar& y, const Scalar& z) {
1030 return x * y + z;
1031 }
1032};
1033
1034#if EIGEN_SCALAR_MADD_USE_FMA
1035template <typename Scalar>
1036struct madd_impl<Scalar, std::enable_if_t<has_fma<Scalar>::value>> {
1037 static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(const Scalar& x, const Scalar& y, const Scalar& z) {
1038 return fma_impl<Scalar>::run(x, y, z);
1039 }
1040};
1041#endif
1042
1043// Representation properties shared by operations that inspect an IEEE-style binary floating-point value.
1044template <typename Scalar>
1045struct binary_floating_point_traits {
1046 using Bits = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type;
1047 static constexpr int kFractionBits = std::numeric_limits<Scalar>::digits - 1;
1048 static constexpr int kExponentBits = int(sizeof(Scalar) * CHAR_BIT) - std::numeric_limits<Scalar>::digits;
1049 static constexpr int kExponentBias = std::numeric_limits<Scalar>::max_exponent - 1;
1050 static constexpr Bits kSignBit = Bits(1) << (sizeof(Scalar) * CHAR_BIT - 1);
1051 static constexpr Bits kExponentUnit = Bits(1) << kFractionBits;
1052 static constexpr Bits kFractionMask = kExponentUnit - 1;
1053 static constexpr Bits kExponentMask = ((Bits(1) << kExponentBits) - 1) << kFractionBits;
1054 static constexpr Bits kMaxFiniteExponentBits = kExponentMask - kExponentUnit;
1055
1056 EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Bits bits(const Scalar& value) { return numext::bit_cast<Bits>(value); }
1057
1058 EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Bits magnitude(const Scalar& value) { return bits(value) & ~kSignBit; }
1059};
1060
1061template <typename Scalar>
1062EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_exactly_zero_no_flush_impl(const Scalar& value, false_type) {
1063 return numext::is_exactly_zero(value);
1064}
1065
1066template <typename Scalar>
1067EIGEN_DEVICE_FUNC EIGEN_DONT_INLINE bool is_zero_magnitude_bits(
1068 typename binary_floating_point_traits<Scalar>::Bits magnitude) {
1069 // Keep the integer comparison out of line: Clang can otherwise fold it into a
1070 // floating-point comparison that treats subnormals as zero under DAZ/FZ.
1071 return magnitude == 0;
1072}
1073
1074// Floating-point comparisons can treat subnormals as zero under DAZ/FTZ.
1075template <typename Scalar>
1076struct sign_impl_binary {
1077 EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Scalar run(const Scalar& a) {
1078 using Binary = binary_floating_point_traits<Scalar>;
1079 using Bits = typename Binary::Bits;
1080 const Bits bits = Binary::bits(a);
1081 const Bits magnitude = bits & ~Binary::kSignBit;
1082 // magnitude | -magnitude has its high bit set iff magnitude is nonzero.
1083 const Bits nonzero = (magnitude | (Bits(0) - magnitude)) >> (sizeof(Bits) * CHAR_BIT - 1);
1084 const Bits nonzero_mask = Bits(0) - nonzero;
1085 const Bits nan_mask = Bits(0) - Bits(magnitude > Binary::kExponentMask);
1086 const Bits signed_one = (bits & Binary::kSignBit) | Binary::bits(Scalar(1));
1087 return numext::bit_cast<Scalar>((bits & nan_mask) | (signed_one & nonzero_mask & ~nan_mask));
1088 }
1089};
1090
1091template <>
1092struct sign_impl<float, false, false> : sign_impl_binary<float> {};
1093
1094template <>
1095struct sign_impl<double, false, false> : sign_impl_binary<double> {};
1096
1097template <typename Scalar>
1098EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_exactly_zero_no_flush_impl(const Scalar& value, true_type) {
1099 return is_zero_magnitude_bits<Scalar>(binary_floating_point_traits<Scalar>::magnitude(value));
1100}
1101
1102// |x - y|, as the packet path computes it: abs keeps the +0 of (+0, -0) and the NaN of (inf, inf).
1103template <typename T, bool IsInteger = NumTraits<T>::IsInteger>
1104struct absdiff_impl {
1105 EIGEN_DEVICE_FUNC static EIGEN_ALWAYS_INLINE T run(const T& x, const T& y) {
1106 EIGEN_USING_STD(abs);
1107 return abs(x - y);
1108 }
1109};
1110
1111// Unsigned x - y wraps, so integers order the operands.
1112template <typename T>
1113struct absdiff_impl<T, true> {
1114 EIGEN_DEVICE_FUNC static EIGEN_ALWAYS_INLINE T run(const T& x, const T& y) { return x > y ? x - y : y - x; }
1115};
1116
1117} // end namespace internal
1118
1119/****************************************************************************
1120 * Generic math functions *
1121 ****************************************************************************/
1122
1123namespace numext {
1124
1125template <typename X>
1126EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool is_exactly_zero_no_flush(const X& x) {
1127 constexpr bool HasBinaryEncoding = (std::is_same<X, float>::value || std::is_same<X, double>::value) &&
1128 std::numeric_limits<X>::is_iec559 && std::numeric_limits<X>::radix == 2 &&
1129 (sizeof(X) == sizeof(uint32_t) || sizeof(X) == sizeof(uint64_t));
1130 return internal::is_exactly_zero_no_flush_impl(x, internal::bool_constant<HasBinaryEncoding>());
1131}
1132
1133#if !defined(EIGEN_GPUCC)
1134template <typename T>
1135EIGEN_DEVICE_FUNC constexpr EIGEN_ALWAYS_INLINE T mini(const T& x, const T& y) {
1136 EIGEN_USING_STD(min)
1137 return min EIGEN_NOT_A_MACRO(x, y);
1138}
1139
1140template <typename T>
1141EIGEN_DEVICE_FUNC constexpr EIGEN_ALWAYS_INLINE T maxi(const T& x, const T& y) {
1142 EIGEN_USING_STD(max)
1143 return max EIGEN_NOT_A_MACRO(x, y);
1144}
1145#else
1146template <typename T>
1147EIGEN_DEVICE_FUNC constexpr EIGEN_ALWAYS_INLINE T mini(const T& x, const T& y) {
1148 return y < x ? y : x;
1149}
1150#if !defined(EIGEN_CONSTEXPR_ARE_DEVICE_FUNC)
1151// Without relaxed constexpr, numeric GPU scalars keep fmin/fmax's number-preferring
1152// NaN behavior. With relaxed constexpr, they use the constexpr ternary overloads,
1153// matching std::min/std::max behavior and supporting custom less-comparable scalars.
1154template <>
1155EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float mini(const float& x, const float& y) {
1156 return fminf(x, y);
1157}
1158template <>
1159EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double mini(const double& x, const double& y) {
1160 return fmin(x, y);
1161}
1162
1163#ifndef EIGEN_GPU_COMPILE_PHASE
1164template <>
1165EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE long double mini(const long double& x, const long double& y) {
1166#if defined(EIGEN_HIPCC)
1167 // no "fminl" on HIP yet
1168 return (x < y) ? x : y;
1169#else
1170 return fminl(x, y);
1171#endif
1172}
1173#endif
1174#endif
1175
1176template <typename T>
1177EIGEN_DEVICE_FUNC constexpr EIGEN_ALWAYS_INLINE T maxi(const T& x, const T& y) {
1178 return x < y ? y : x;
1179}
1180#if !defined(EIGEN_CONSTEXPR_ARE_DEVICE_FUNC)
1181template <>
1182EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float maxi(const float& x, const float& y) {
1183 return fmaxf(x, y);
1184}
1185template <>
1186EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double maxi(const double& x, const double& y) {
1187 return fmax(x, y);
1188}
1189#ifndef EIGEN_GPU_COMPILE_PHASE
1190template <>
1191EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE long double maxi(const long double& x, const long double& y) {
1192#if defined(EIGEN_HIPCC)
1193 // no "fmaxl" on HIP yet
1194 return (x > y) ? x : y;
1195#else
1196 return fmaxl(x, y);
1197#endif
1198}
1199#endif
1200#endif
1201#endif
1202
1203#if defined(SYCL_DEVICE_ONLY)
1204
1205#define SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_BINARY(NAME, FUNC) \
1206 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_char) \
1207 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_short) \
1208 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_int) \
1209 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_long)
1210#define SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_UNARY(NAME, FUNC) \
1211 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_char) \
1212 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_short) \
1213 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_int) \
1214 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_long)
1215#define SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_BINARY(NAME, FUNC) \
1216 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_uchar) \
1217 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_ushort) \
1218 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_uint) \
1219 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_ulong)
1220#define SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_UNARY(NAME, FUNC) \
1221 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_uchar) \
1222 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_ushort) \
1223 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_uint) \
1224 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_ulong)
1225#define SYCL_SPECIALIZE_INTEGER_TYPES_BINARY(NAME, FUNC) \
1226 SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_BINARY(NAME, FUNC) \
1227 SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_BINARY(NAME, FUNC)
1228#define SYCL_SPECIALIZE_INTEGER_TYPES_UNARY(NAME, FUNC) \
1229 SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_UNARY(NAME, FUNC) \
1230 SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_UNARY(NAME, FUNC)
1231#define SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(NAME, FUNC) \
1232 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_float) \
1233 SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, cl::sycl::cl_double)
1234#define SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(NAME, FUNC) \
1235 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_float) \
1236 SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, cl::sycl::cl_double)
1237#define SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE(NAME, FUNC, RET_TYPE) \
1238 SYCL_SPECIALIZE_GEN_UNARY_FUNC(NAME, FUNC, RET_TYPE, cl::sycl::cl_float) \
1239 SYCL_SPECIALIZE_GEN_UNARY_FUNC(NAME, FUNC, RET_TYPE, cl::sycl::cl_double)
1240
1241#define SYCL_SPECIALIZE_GEN_UNARY_FUNC(NAME, FUNC, RET_TYPE, ARG_TYPE) \
1242 template <> \
1243 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE RET_TYPE NAME(const ARG_TYPE& x) { \
1244 return cl::sycl::FUNC(x); \
1245 }
1246
1247#define SYCL_SPECIALIZE_UNARY_FUNC(NAME, FUNC, TYPE) SYCL_SPECIALIZE_GEN_UNARY_FUNC(NAME, FUNC, TYPE, TYPE)
1248
1249#define SYCL_SPECIALIZE_GEN1_BINARY_FUNC(NAME, FUNC, RET_TYPE, ARG_TYPE1, ARG_TYPE2) \
1250 template <> \
1251 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE RET_TYPE NAME(const ARG_TYPE1& x, const ARG_TYPE2& y) { \
1252 return cl::sycl::FUNC(x, y); \
1253 }
1254
1255#define SYCL_SPECIALIZE_GEN2_BINARY_FUNC(NAME, FUNC, RET_TYPE, ARG_TYPE) \
1256 SYCL_SPECIALIZE_GEN1_BINARY_FUNC(NAME, FUNC, RET_TYPE, ARG_TYPE, ARG_TYPE)
1257
1258#define SYCL_SPECIALIZE_BINARY_FUNC(NAME, FUNC, TYPE) SYCL_SPECIALIZE_GEN2_BINARY_FUNC(NAME, FUNC, TYPE, TYPE)
1259
1260SYCL_SPECIALIZE_INTEGER_TYPES_BINARY(mini, min)
1261SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(mini, fmin)
1262SYCL_SPECIALIZE_INTEGER_TYPES_BINARY(maxi, max)
1263SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(maxi, fmax)
1264
1265#endif
1266
1267template <typename Scalar>
1268EIGEN_DEVICE_FUNC inline typename NumTraits<Scalar>::Real arg(const Scalar& x) {
1269 return EIGEN_MATHFUNC_IMPL(arg, Scalar)::run(x);
1270}
1271
1272template <typename Scalar>
1273EIGEN_DEVICE_FUNC inline auto imag_ref(const Scalar& x)
1274 -> internal::add_const_on_value_type_t<decltype(internal::imag_ref_impl<Scalar>::run(x))> {
1275 return internal::imag_ref_impl<Scalar>::run(x);
1276}
1277
1278template <typename Scalar>
1279EIGEN_DEVICE_FUNC inline auto imag_ref(Scalar& x)
1280 -> decltype(EIGEN_MATHFUNC_IMPL(imag_ref, std::remove_const_t<Scalar>)::run(x)) {
1281 return EIGEN_MATHFUNC_IMPL(imag_ref, std::remove_const_t<Scalar>)::run(x);
1282}
1283
1284template <typename Scalar>
1285EIGEN_DEVICE_FUNC inline Scalar conj(const Scalar& x) {
1286 return EIGEN_MATHFUNC_IMPL(conj, Scalar)::run(x);
1287}
1288
1289template <typename Scalar>
1290EIGEN_DEVICE_FUNC inline Scalar sign(const Scalar& x) {
1291 return EIGEN_MATHFUNC_IMPL(sign, Scalar)::run(x);
1292}
1293
1294template <typename Scalar>
1295EIGEN_DEVICE_FUNC inline Scalar copysign(const Scalar& x, const Scalar& y) {
1296 return EIGEN_MATHFUNC_IMPL(copysign, Scalar)::run(x, y);
1297}
1298
1299template <typename Scalar>
1300EIGEN_DEVICE_FUNC inline Scalar negate(const Scalar& x) {
1301 return EIGEN_MATHFUNC_IMPL(negate, Scalar)::run(x);
1302}
1303
1304template <typename Scalar>
1305EIGEN_DEVICE_FUNC inline decltype(auto) abs2(const Scalar& x) {
1306 return EIGEN_MATHFUNC_IMPL(abs2, Scalar)::run(x);
1307}
1308
1309EIGEN_DEVICE_FUNC inline bool abs2(bool x) { return x; }
1310
1311template <typename T>
1312EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T absdiff(const T& x, const T& y) {
1313 return internal::absdiff_impl<T>::run(x, y);
1314}
1315template <>
1316EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float absdiff(const float& x, const float& y) {
1317 return fabsf(x - y);
1318}
1319template <>
1320EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double absdiff(const double& x, const double& y) {
1321 return fabs(x - y);
1322}
1323
1324// HIP and CUDA do not support long double.
1325#ifndef EIGEN_GPU_COMPILE_PHASE
1326template <>
1327EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE long double absdiff(const long double& x, const long double& y) {
1328 return fabsl(x - y);
1329}
1330#endif
1331
1332template <typename Scalar>
1333EIGEN_DEVICE_FUNC inline typename NumTraits<Scalar>::Real norm1(const Scalar& x) {
1334 return EIGEN_MATHFUNC_IMPL(norm1, Scalar)::run(x);
1335}
1336
1337template <typename Scalar>
1338EIGEN_DEVICE_FUNC inline typename NumTraits<Scalar>::Real hypot(const Scalar& x, const Scalar& y) {
1339 return EIGEN_MATHFUNC_IMPL(hypot, Scalar)::run(x, y);
1340}
1341
1342#if defined(SYCL_DEVICE_ONLY)
1343SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(hypot, hypot)
1344#endif
1345
1346template <typename Scalar>
1347EIGEN_DEVICE_FUNC inline Scalar log1p(const Scalar& x) {
1348 return EIGEN_MATHFUNC_IMPL(log1p, Scalar)::run(x);
1349}
1350
1351#if defined(SYCL_DEVICE_ONLY)
1352SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(log1p, log1p)
1353#endif
1354
1355#if defined(EIGEN_GPUCC)
1356template <>
1357EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float log1p(const float& x) {
1358 return ::log1pf(x);
1359}
1360
1361template <>
1362EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double log1p(const double& x) {
1363 return ::log1p(x);
1364}
1365#endif
1366
1367template <typename ScalarX, typename ScalarY>
1368EIGEN_DEVICE_FUNC inline typename internal::pow_impl<ScalarX, ScalarY>::result_type pow(const ScalarX& x,
1369 const ScalarY& y) {
1370 return internal::pow_impl<ScalarX, ScalarY>::run(x, y);
1371}
1372
1373#if defined(SYCL_DEVICE_ONLY)
1374SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(pow, pow)
1375#endif
1376
1377template <typename T>
1378EIGEN_DEVICE_FUNC bool(isnan)(const T& x) {
1379 return internal::isnan_impl(x);
1380}
1381template <typename T>
1382EIGEN_DEVICE_FUNC bool(isinf)(const T& x) {
1383 return internal::isinf_impl(x);
1384}
1385template <typename T>
1386EIGEN_DEVICE_FUNC bool(isfinite)(const T& x) {
1387 return internal::isfinite_impl(x);
1388}
1389
1390#if defined(SYCL_DEVICE_ONLY)
1391SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE(isnan, isnan, bool)
1392SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE(isinf, isinf, bool)
1393SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE(isfinite, isfinite, bool)
1394#endif
1395
1396template <typename Scalar>
1397EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar rint(const Scalar& x) {
1398 return internal::nearest_integer_impl<Scalar>::run_rint(x);
1399}
1400
1401template <typename Scalar>
1402EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar round(const Scalar& x) {
1403 return internal::nearest_integer_impl<Scalar>::run_round(x);
1404}
1405
1406template <typename Scalar>
1407EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar(floor)(const Scalar& x) {
1408 return internal::nearest_integer_impl<Scalar>::run_floor(x);
1409}
1410
1411template <typename Scalar>
1412EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar(ceil)(const Scalar& x) {
1413 return internal::nearest_integer_impl<Scalar>::run_ceil(x);
1414}
1415
1416template <typename Scalar>
1417EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar(trunc)(const Scalar& x) {
1418 return internal::nearest_integer_impl<Scalar>::run_trunc(x);
1419}
1420
1421#if defined(SYCL_DEVICE_ONLY)
1422SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(round, round)
1423SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(floor, floor)
1424SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(ceil, ceil)
1425SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(trunc, trunc)
1426#endif
1427
1428#if defined(EIGEN_GPUCC)
1429template <>
1430EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float floor(const float& x) {
1431 return ::floorf(x);
1432}
1433template <>
1434EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double floor(const double& x) {
1435 return ::floor(x);
1436}
1437template <>
1438EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float ceil(const float& x) {
1439 return ::ceilf(x);
1440}
1441template <>
1442EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double ceil(const double& x) {
1443 return ::ceil(x);
1444}
1445template <>
1446EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float trunc(const float& x) {
1447 return ::truncf(x);
1448}
1449template <>
1450EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double trunc(const double& x) {
1451 return ::trunc(x);
1452}
1453#endif
1454
1455// Integer division with rounding up.
1456// T is assumed to be an integer type with a>=0, and b>0
1457template <typename T>
1458EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE constexpr T div_ceil(T a, T b) {
1459 using UnsignedT = std::make_unsigned_t<T>;
1460 EIGEN_STATIC_ASSERT((NumTraits<T>::IsInteger), THIS FUNCTION IS FOR INTEGER TYPES)
1461 // Note: explicitly declaring a and b as non-negative values allows the compiler to use better optimizations
1462 const UnsignedT ua = UnsignedT(a);
1463 const UnsignedT ub = UnsignedT(b);
1464 // Note: This form is used because it cannot overflow.
1465 return ua == 0 ? 0 : (ua - 1) / ub + 1;
1466}
1467
1468// Integer round down to nearest power of b
1469// T is assumed to be an integer type with a>=0, and b>0
1470template <typename T, typename U>
1471EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE constexpr T round_down(T a, U b) {
1472 using UnsignedT = std::make_unsigned_t<T>;
1473 using UnsignedU = std::make_unsigned_t<U>;
1474 EIGEN_STATIC_ASSERT((NumTraits<T>::IsInteger), THIS FUNCTION IS FOR INTEGER TYPES)
1475 EIGEN_STATIC_ASSERT((NumTraits<U>::IsInteger), THIS FUNCTION IS FOR INTEGER TYPES)
1476 // Note: explicitly declaring a and b as non-negative values allows the compiler to use better optimizations
1477 const UnsignedT ua = UnsignedT(a);
1478 const UnsignedU ub = UnsignedU(b);
1479 return ub * (ua / ub);
1480}
1481
1482template <typename T>
1483EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T log2(T x) {
1484 EIGEN_USING_STD(log2);
1485 return log2(x);
1486}
1487
1490constexpr int log2(int x) {
1491 unsigned int v(x);
1492 constexpr int table[32] = {0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
1493 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31};
1494 v |= v >> 1;
1495 v |= v >> 2;
1496 v |= v >> 4;
1497 v |= v >> 8;
1498 v |= v >> 16;
1499 return table[(v * 0x07C4ACDDU) >> 27];
1500}
1501
1511template <typename Scalar>
1512EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE Scalar sqrt(const Scalar& x) {
1513 return EIGEN_MATHFUNC_IMPL(sqrt, Scalar)::run(x);
1514}
1515
1516// Boolean specialization, avoids implicit float to bool conversion (-Wimplicit-conversion-floating-point-to-bool).
1517template <>
1518EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_DEVICE_FUNC bool sqrt<bool>(const bool& x) {
1519 return x;
1520}
1521
1522#if defined(SYCL_DEVICE_ONLY)
1523SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(sqrt, sqrt)
1524#endif
1525
1527template <typename T>
1528EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::enable_if_t<!NumTraits<T>::IsComplex, T> cbrt(const T& x) {
1529 EIGEN_USING_STD(cbrt);
1530 return static_cast<T>(cbrt(x));
1531}
1532
1533template <typename T>
1534EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::enable_if_t<NumTraits<T>::IsComplex, T> cbrt(const T& x) {
1535 EIGEN_USING_STD(pow);
1536 return pow(x, typename NumTraits<T>::Real(1.0 / 3.0));
1537}
1538
1540template <typename T>
1541EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T rsqrt(const T& x) {
1542 return internal::rsqrt_impl<T>::run(x);
1543}
1544
1545template <typename T>
1546EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T log(const T& x) {
1547 return internal::log_impl<T>::run(x);
1548}
1549
1550#if defined(SYCL_DEVICE_ONLY)
1551SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(log, log)
1552#endif
1553
1554#if defined(EIGEN_GPUCC)
1555template <>
1556EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float log(const float& x) {
1557 return ::logf(x);
1558}
1559
1560template <>
1561EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double log(const double& x) {
1562 return ::log(x);
1563}
1564#endif
1565
1566template <typename T>
1567EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
1568 std::enable_if_t<NumTraits<T>::IsSigned && !NumTraits<T>::IsComplex, typename NumTraits<T>::Real>
1569 abs(const T& x) {
1570 EIGEN_USING_STD(abs);
1571 return abs(x);
1572}
1573
1574template <typename T>
1575EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::enable_if_t<NumTraits<T>::IsComplex, typename NumTraits<T>::Real> abs(
1576 const T& x) {
1577 return numext::hypot(numext::real(x), numext::imag(x));
1578}
1579
1580template <typename T>
1581EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
1582 std::enable_if_t<!NumTraits<T>::IsSigned && !NumTraits<T>::IsComplex, typename NumTraits<T>::Real>
1583 abs(const T& x) {
1584 return x;
1585}
1586
1587#if defined(SYCL_DEVICE_ONLY)
1588SYCL_SPECIALIZE_INTEGER_TYPES_UNARY(abs, abs)
1589SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(abs, fabs)
1590#endif
1591
1592#if defined(EIGEN_GPUCC)
1593template <>
1594EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float abs(const float& x) {
1595 return ::fabsf(x);
1596}
1597
1598template <>
1599EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double abs(const double& x) {
1600 return ::fabs(x);
1601}
1602
1603template <>
1604EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float abs(const std::complex<float>& x) {
1605 return ::hypotf(x.real(), x.imag());
1606}
1607
1608template <>
1609EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double abs(const std::complex<double>& x) {
1610 return ::hypot(x.real(), x.imag());
1611}
1612#endif
1613
1614template <typename Scalar, bool IsInteger = NumTraits<Scalar>::IsInteger, bool IsSigned = NumTraits<Scalar>::IsSigned>
1615struct signbit_impl;
1616template <typename Scalar>
1617struct signbit_impl<Scalar, false, true> {
1618 static constexpr size_t Size = sizeof(Scalar);
1619 static constexpr size_t Shift = (CHAR_BIT * Size) - 1;
1620 using intSize_t = typename get_integer_by_size<Size>::signed_type;
1621 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static Scalar run(const Scalar& x) {
1622 intSize_t a = bit_cast<intSize_t, Scalar>(x);
1623 a = a >> Shift;
1624 Scalar result = bit_cast<Scalar, intSize_t>(a);
1625 return result;
1626 }
1627};
1628template <typename Scalar>
1629struct signbit_impl<Scalar, true, true> {
1630 static constexpr size_t Size = sizeof(Scalar);
1631 static constexpr size_t Shift = (CHAR_BIT * Size) - 1;
1632 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static constexpr Scalar run(const Scalar& x) { return x >> Shift; }
1633};
1634template <typename Scalar>
1635struct signbit_impl<Scalar, true, false> {
1636 EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static constexpr Scalar run(const Scalar&) { return Scalar(0); }
1637};
1638template <typename Scalar>
1639EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE static constexpr Scalar signbit(const Scalar& x) {
1640 return signbit_impl<Scalar>::run(x);
1641}
1642
1643template <typename T>
1644EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T exp(const T& x) {
1645 EIGEN_USING_STD(exp);
1646 return exp(x);
1647}
1648
1649// MSVC through at least 19.51 mishandles some finite-real, non-finite-imaginary inputs to std::exp(complex).
1650#if EIGEN_COMP_MSVC
1651template <typename RealScalar>
1652EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::complex<RealScalar> exp(const std::complex<RealScalar>& x) {
1653 // Evaluate eagerly to preserve the floating-point exceptions raised by std::exp before correcting its result.
1654 const std::complex<RealScalar> result = std::exp(x);
1655 // If z is (x,±∞) (for any finite x), the result is (NaN,NaN) and FE_INVALID is raised.
1656 // If z is (x,NaN) (for any finite x), the result is (NaN,NaN) and FE_INVALID may be raised.
1657 if ((isfinite)(real_ref(x)) && !(isfinite)(imag_ref(x))) {
1658 return std::complex<RealScalar>(NumTraits<RealScalar>::quiet_NaN(), NumTraits<RealScalar>::quiet_NaN());
1659 }
1660#if EIGEN_COMP_MSVC < 1931
1661 // Before 19.31, MSVC also mishandles a positive-infinite real part with a non-finite imaginary part.
1662 // If z is (+∞,±∞), the result is (±∞,NaN) and FE_INVALID is raised (the sign of the real part is unspecified)
1663 // If z is (+∞,NaN), the result is (±∞,NaN) (the sign of the real part is unspecified)
1664 if ((real_ref(x) == NumTraits<RealScalar>::infinity() && !(isfinite)(imag_ref(x)))) {
1665 return std::complex<RealScalar>(NumTraits<RealScalar>::infinity(), NumTraits<RealScalar>::quiet_NaN());
1666 }
1667#endif
1668 return result;
1669}
1670#endif
1671
1672#if defined(SYCL_DEVICE_ONLY)
1673SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(exp, exp)
1674#endif
1675
1676#if defined(EIGEN_GPUCC)
1677template <>
1678EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float exp(const float& x) {
1679 return ::expf(x);
1680}
1681
1682template <>
1683EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double exp(const double& x) {
1684 return ::exp(x);
1685}
1686
1687template <>
1688EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::complex<float> exp(const std::complex<float>& x) {
1689 float com = ::expf(x.real());
1690 float res_real = com * ::cosf(x.imag());
1691 float res_imag = com * ::sinf(x.imag());
1692 return std::complex<float>(res_real, res_imag);
1693}
1694
1695template <>
1696EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::complex<double> exp(const std::complex<double>& x) {
1697 double com = ::exp(x.real());
1698 double res_real = com * ::cos(x.imag());
1699 double res_imag = com * ::sin(x.imag());
1700 return std::complex<double>(res_real, res_imag);
1701}
1702#endif
1703
1704template <typename T>
1705EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T exp2(const T& x) {
1706 EIGEN_USING_STD(exp2);
1707 return exp2(x);
1708}
1709
1710// There is no std::exp2 for complex arguments.
1711template <typename RealScalar>
1712EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE std::complex<RealScalar> exp2(const std::complex<RealScalar>& x) {
1713 return internal::complex_exp2_impl<RealScalar>::run(x);
1714}
1715
1716#if defined(SYCL_DEVICE_ONLY)
1717SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(exp2, exp2)
1718#endif
1719
1720#if defined(EIGEN_GPUCC)
1721template <>
1722EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float exp2(const float& x) {
1723 return ::exp2f(x);
1724}
1725
1726template <>
1727EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double exp2(const double& x) {
1728 return ::exp2(x);
1729}
1730#endif
1731
1732// Exact scaling by 2^exponent, including denormals and unrepresentable scale factors.
1733// Results outside the finite range saturate to zero or infinity.
1734template <typename T>
1735EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T ldexp(const T& x, int exponent) {
1736 EIGEN_USING_STD(ldexp);
1737 return static_cast<T>(ldexp(x, exponent));
1738}
1739
1740#if defined(SYCL_DEVICE_ONLY)
1741template <>
1742EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE cl::sycl::cl_float ldexp(const cl::sycl::cl_float& x, int exponent) {
1743 return cl::sycl::ldexp(x, exponent);
1744}
1745template <>
1746EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE cl::sycl::cl_double ldexp(const cl::sycl::cl_double& x, int exponent) {
1747 return cl::sycl::ldexp(x, exponent);
1748}
1749#endif
1750
1751#if defined(EIGEN_GPUCC)
1752template <>
1753EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float ldexp(const float& x, int exponent) {
1754 return ::ldexpf(x, exponent);
1755}
1756
1757template <>
1758EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double ldexp(const double& x, int exponent) {
1759 return ::ldexp(x, exponent);
1760}
1761#endif
1762
1763template <typename Scalar>
1764EIGEN_DEVICE_FUNC inline Scalar expm1(const Scalar& x) {
1765 return EIGEN_MATHFUNC_IMPL(expm1, Scalar)::run(x);
1766}
1767
1768#if defined(SYCL_DEVICE_ONLY)
1769SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(expm1, expm1)
1770#endif
1771
1772#if defined(EIGEN_GPUCC)
1773template <>
1774EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float expm1(const float& x) {
1775 return ::expm1f(x);
1776}
1777
1778template <>
1779EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double expm1(const double& x) {
1780 return ::expm1(x);
1781}
1782#endif
1783
1784template <typename T>
1785EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T cos(const T& x) {
1786 EIGEN_USING_STD(cos);
1787 return cos(x);
1788}
1789
1790#if defined(SYCL_DEVICE_ONLY)
1791SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(cos, cos)
1792#endif
1793
1794#if defined(EIGEN_GPUCC)
1795template <>
1796EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float cos(const float& x) {
1797 return ::cosf(x);
1798}
1799
1800template <>
1801EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double cos(const double& x) {
1802 return ::cos(x);
1803}
1804#endif
1805
1806template <typename T>
1807EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T sin(const T& x) {
1808 EIGEN_USING_STD(sin);
1809 return sin(x);
1810}
1811
1812#if defined(SYCL_DEVICE_ONLY)
1813SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(sin, sin)
1814#endif
1815
1816#if defined(EIGEN_GPUCC)
1817template <>
1818EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float sin(const float& x) {
1819 return ::sinf(x);
1820}
1821
1822template <>
1823EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double sin(const double& x) {
1824 return ::sin(x);
1825}
1826#endif
1827
1828template <typename T>
1829EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T tan(const T& x) {
1830 EIGEN_USING_STD(tan);
1831 return tan(x);
1832}
1833
1834#if defined(SYCL_DEVICE_ONLY)
1835SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(tan, tan)
1836#endif
1837
1838#if defined(EIGEN_GPUCC)
1839template <>
1840EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float tan(const float& x) {
1841 return ::tanf(x);
1842}
1843
1844template <>
1845EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double tan(const double& x) {
1846 return ::tan(x);
1847}
1848#endif
1849
1850template <typename T>
1851EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T acos(const T& x) {
1852 EIGEN_USING_STD(acos);
1853 return acos(x);
1854}
1855
1856template <typename T>
1857EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T acosh(const T& x) {
1858 EIGEN_USING_STD(acosh);
1859 return static_cast<T>(acosh(x));
1860}
1861
1862#if defined(SYCL_DEVICE_ONLY)
1863SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(acos, acos)
1864SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(acosh, acosh)
1865#endif
1866
1867#if defined(EIGEN_GPUCC)
1868template <>
1869EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float acos(const float& x) {
1870 return ::acosf(x);
1871}
1872
1873template <>
1874EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double acos(const double& x) {
1875 return ::acos(x);
1876}
1877#endif
1878
1879template <typename T>
1880EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T asin(const T& x) {
1881 EIGEN_USING_STD(asin);
1882 return asin(x);
1883}
1884
1885template <typename T>
1886EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T asinh(const T& x) {
1887 EIGEN_USING_STD(asinh);
1888 return static_cast<T>(asinh(x));
1889}
1890
1891#if defined(SYCL_DEVICE_ONLY)
1892SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(asin, asin)
1893SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(asinh, asinh)
1894#endif
1895
1896#if defined(EIGEN_GPUCC)
1897template <>
1898EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float asin(const float& x) {
1899 return ::asinf(x);
1900}
1901
1902template <>
1903EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double asin(const double& x) {
1904 return ::asin(x);
1905}
1906#endif
1907
1908template <typename T>
1909EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T atan(const T& x) {
1910 EIGEN_USING_STD(atan);
1911 return static_cast<T>(atan(x));
1912}
1913
1914template <typename T, std::enable_if_t<!NumTraits<T>::IsComplex, int> = 0>
1915EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T atan2(const T& y, const T& x) {
1916 EIGEN_USING_STD(atan2);
1917 return static_cast<T>(atan2(y, x));
1918}
1919
1920template <typename T>
1921EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T atanh(const T& x) {
1922 EIGEN_USING_STD(atanh);
1923 return static_cast<T>(atanh(x));
1924}
1925
1926#if defined(SYCL_DEVICE_ONLY)
1927SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(atan, atan)
1928SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(atanh, atanh)
1929#endif
1930
1931#if defined(EIGEN_GPUCC)
1932template <>
1933EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float atan(const float& x) {
1934 return ::atanf(x);
1935}
1936
1937template <>
1938EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double atan(const double& x) {
1939 return ::atan(x);
1940}
1941#endif
1942
1943template <typename T>
1944EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T cosh(const T& x) {
1945 EIGEN_USING_STD(cosh);
1946 return static_cast<T>(cosh(x));
1947}
1948
1949#if defined(SYCL_DEVICE_ONLY)
1950SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(cosh, cosh)
1951#endif
1952
1953#if defined(EIGEN_GPUCC)
1954template <>
1955EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float cosh(const float& x) {
1956 return ::coshf(x);
1957}
1958
1959template <>
1960EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double cosh(const double& x) {
1961 return ::cosh(x);
1962}
1963#endif
1964
1965template <typename T>
1966EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T sinh(const T& x) {
1967 EIGEN_USING_STD(sinh);
1968 return static_cast<T>(sinh(x));
1969}
1970
1971#if defined(SYCL_DEVICE_ONLY)
1972SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(sinh, sinh)
1973#endif
1974
1975#if defined(EIGEN_GPUCC)
1976template <>
1977EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float sinh(const float& x) {
1978 return ::sinhf(x);
1979}
1980
1981template <>
1982EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double sinh(const double& x) {
1983 return ::sinh(x);
1984}
1985#endif
1986
1987template <typename T>
1988EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T tanh(const T& x) {
1989 EIGEN_USING_STD(tanh);
1990 return tanh(x);
1991}
1992
1993#if (!defined(EIGEN_GPUCC)) && EIGEN_FAST_MATH && !defined(SYCL_DEVICE_ONLY)
1994EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float tanh(float x) { return internal::ptanh_float(x); }
1995#endif
1996
1997#if defined(SYCL_DEVICE_ONLY)
1998SYCL_SPECIALIZE_FLOATING_TYPES_UNARY(tanh, tanh)
1999#endif
2000
2001#if defined(EIGEN_GPUCC)
2002template <>
2003EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float tanh(const float& x) {
2004 return ::tanhf(x);
2005}
2006
2007template <>
2008EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double tanh(const double& x) {
2009 return ::tanh(x);
2010}
2011#endif
2012
2013template <typename T>
2014EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T fmod(const T& a, const T& b) {
2015 EIGEN_USING_STD(fmod);
2016 return fmod(a, b);
2017}
2018
2019#if defined(SYCL_DEVICE_ONLY)
2020SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(fmod, fmod)
2021#endif
2022
2023#if defined(EIGEN_GPUCC)
2024template <>
2025EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float fmod(const float& a, const float& b) {
2026 return ::fmodf(a, b);
2027}
2028
2029template <>
2030EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double fmod(const double& a, const double& b) {
2031 return ::fmod(a, b);
2032}
2033#endif
2034
2035template <typename T>
2036EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T nextafter(const T& from, const T& to) {
2037 EIGEN_USING_STD(nextafter);
2038 return nextafter(from, to);
2039}
2040
2046template <typename Scalar,
2047 std::enable_if_t<!std::is_same<Scalar, long double>::value && std::numeric_limits<Scalar>::is_iec559 &&
2048 std::numeric_limits<Scalar>::radix == 2 &&
2049 (sizeof(Scalar) == sizeof(uint16_t) || sizeof(Scalar) == sizeof(uint32_t) ||
2050 sizeof(Scalar) == sizeof(uint64_t)),
2051 int> = 0>
2052EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar ceil_power_of_two(const Scalar& value) {
2053 using Binary = internal::binary_floating_point_traits<Scalar>;
2054 using Bits = typename Binary::Bits;
2055 const Bits valueBits = Binary::bits(value);
2056 const Bits resultBits = Bits((valueBits + Binary::kFractionMask) & Binary::kExponentMask);
2057 return bit_cast<Scalar>(resultBits);
2058}
2059
2060// long double has no portable object representation: depending on the ABI it can be binary64, padded x87 extended
2061// precision, IEEE binary128, or IBM double-double. Construct the result arithmetically instead of inspecting its bits.
2062#if !defined(EIGEN_GPU_COMPILE_PHASE)
2063EIGEN_STRONG_INLINE long double ceil_power_of_two(const long double& value) {
2064 EIGEN_USING_STD(frexp);
2065 EIGEN_USING_STD(ldexp);
2066 int exponent = 0;
2067 const long double fraction = frexp(value, &exponent);
2068 if (numext::is_exactly_zero(fraction)) return 0.0L;
2069 if (numext::equal_strict(fraction, 0.5L)) --exponent;
2070 if (exponent < std::numeric_limits<long double>::min_exponent - 1)
2071 exponent = std::numeric_limits<long double>::min_exponent - 1;
2072 return ldexp(1.0L, exponent);
2073}
2074#endif
2075
2076#if defined(SYCL_DEVICE_ONLY)
2077SYCL_SPECIALIZE_FLOATING_TYPES_BINARY(nextafter, nextafter)
2078#endif
2079
2080#if defined(EIGEN_GPUCC)
2081template <>
2082EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE float nextafter(const float& from, const float& to) {
2083 return ::nextafterf(from, to);
2084}
2085
2086template <>
2087EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE double nextafter(const double& from, const double& to) {
2088 return ::nextafter(from, to);
2089}
2090#endif
2091
2092#if defined(SYCL_DEVICE_ONLY)
2093#undef SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_BINARY
2094#undef SYCL_SPECIALIZE_SIGNED_INTEGER_TYPES_UNARY
2095#undef SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_BINARY
2096#undef SYCL_SPECIALIZE_UNSIGNED_INTEGER_TYPES_UNARY
2097#undef SYCL_SPECIALIZE_INTEGER_TYPES_BINARY
2098#undef SYCL_SPECIALIZE_INTEGER_TYPES_UNARY
2099#undef SYCL_SPECIALIZE_FLOATING_TYPES_BINARY
2100#undef SYCL_SPECIALIZE_FLOATING_TYPES_UNARY
2101#undef SYCL_SPECIALIZE_FLOATING_TYPES_UNARY_FUNC_RET_TYPE
2102#undef SYCL_SPECIALIZE_GEN_UNARY_FUNC
2103#undef SYCL_SPECIALIZE_UNARY_FUNC
2104#undef SYCL_SPECIALIZE_GEN1_BINARY_FUNC
2105#undef SYCL_SPECIALIZE_GEN2_BINARY_FUNC
2106#undef SYCL_SPECIALIZE_BINARY_FUNC
2107#endif
2108
2109// A logical shift is a shift of the bit pattern, so it goes through the unsigned type of the same
2110// width whatever Scalar's signedness. A Scalar narrower than int promotes to int for the shift
2111// itself, which the explicit truncation undoes: for the left shift the promoted value can carry set
2112// bits above Scalar's width, and dropping them is the operation rather than an accident of the cast.
2113template <typename Scalar, typename Enable = std::enable_if_t<std::is_integral<Scalar>::value>>
2114EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar logical_shift_left(const Scalar& a, int n) {
2115 using UnsignedScalar = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type;
2116 return bit_cast<Scalar, UnsignedScalar>(static_cast<UnsignedScalar>(bit_cast<UnsignedScalar, Scalar>(a) << n));
2117}
2118
2119template <typename Scalar, typename Enable = std::enable_if_t<std::is_integral<Scalar>::value>>
2120EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar logical_shift_right(const Scalar& a, int n) {
2121 using UnsignedScalar = typename numext::get_integer_by_size<sizeof(Scalar)>::unsigned_type;
2122 return bit_cast<Scalar, UnsignedScalar>(static_cast<UnsignedScalar>(bit_cast<UnsignedScalar, Scalar>(a) >> n));
2123}
2124
2125// An arithmetic shift propagates the sign bit, so it coincides with the logical shift when Scalar is
2126// unsigned and has none. Shifting through the signed type unconditionally would sign-extend an
2127// ordinary value bit, which is what every backend's parithmetic_shift_right on unsigned packets
2128// avoids, leaving the scalar and vectorized paths of one expression disagreeing. Scalar's own
2129// operator>> already selects on its signedness, so no reinterpretation is needed; the cast only
2130// undoes the integral promotion a Scalar narrower than int is subject to.
2131template <typename Scalar, typename Enable = std::enable_if_t<std::is_integral<Scalar>::value>>
2132EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar arithmetic_shift_right(const Scalar& a, int n) {
2133 return static_cast<Scalar>(a >> n);
2134}
2135
2136template <typename Scalar>
2137EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar fma(const Scalar& x, const Scalar& y, const Scalar& z) {
2138 return internal::fma_impl<Scalar>::run(x, y, z);
2139}
2140
2141// Multiply-add.
2142template <typename Scalar>
2143EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar madd(const Scalar& x, const Scalar& y, const Scalar& z) {
2144 return internal::madd_impl<Scalar>::run(x, y, z);
2145}
2146
2147} // end namespace numext
2148
2149namespace internal {
2150
2151template <typename T>
2152EIGEN_DEVICE_FUNC bool isfinite_impl(const std::complex<T>& x) {
2153 EIGEN_IF_CONSTEXPR ((std::is_floating_point<T>::value)) {
2154 // Eager classification lets coefficient-wise loops remain branch-free.
2155 return static_cast<unsigned int>((numext::isfinite)(numext::real(x))) & (numext::isfinite)(numext::imag(x));
2156 } else {
2157 return (numext::isfinite)(numext::real(x)) && (numext::isfinite)(numext::imag(x));
2158 }
2159}
2160
2161template <typename T>
2162EIGEN_DEVICE_FUNC bool isnan_impl(const std::complex<T>& x) {
2163 EIGEN_IF_CONSTEXPR ((std::is_floating_point<T>::value)) {
2164 return static_cast<unsigned int>((numext::isnan)(numext::real(x))) | (numext::isnan)(numext::imag(x));
2165 } else {
2166 return (numext::isnan)(numext::real(x)) || (numext::isnan)(numext::imag(x));
2167 }
2168}
2169
2170template <typename T>
2171EIGEN_DEVICE_FUNC bool isinf_impl(const std::complex<T>& x) {
2172 EIGEN_IF_CONSTEXPR ((std::is_floating_point<T>::value)) {
2173 const bool has_inf = static_cast<unsigned int>((numext::isinf)(numext::real(x))) | (numext::isinf)(numext::imag(x));
2174 const bool has_nan = (numext::isnan)(x);
2175 return has_inf & !has_nan;
2176 } else {
2177 return ((numext::isinf)(numext::real(x)) || (numext::isinf)(numext::imag(x))) && (!(numext::isnan)(x));
2178 }
2179}
2180
2181/****************************************************************************
2182 * Implementation of fuzzy comparisons *
2183 ****************************************************************************/
2184
2185template <typename Scalar, bool IsComplex, bool IsInteger>
2186struct scalar_fuzzy_default_impl {};
2187
2188template <typename Scalar>
2189struct scalar_fuzzy_default_impl<Scalar, false, false> {
2190 using RealScalar = typename NumTraits<Scalar>::Real;
2191 template <typename OtherScalar>
2192 EIGEN_DEVICE_FUNC static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y,
2193 const RealScalar& prec) {
2194 return numext::abs(x) <= numext::abs(y) * prec;
2195 }
2196 EIGEN_DEVICE_FUNC static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec) {
2197 return numext::abs(x - y) <= numext::mini(numext::abs(x), numext::abs(y)) * prec;
2198 }
2199 EIGEN_DEVICE_FUNC static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar& prec) {
2200 return x <= y || isApprox(x, y, prec);
2201 }
2202};
2203
2204template <typename Scalar>
2205struct scalar_fuzzy_default_impl<Scalar, false, true> {
2206 using RealScalar = typename NumTraits<Scalar>::Real;
2207 template <typename OtherScalar>
2208 EIGEN_DEVICE_FUNC static inline bool isMuchSmallerThan(const Scalar& x, const Scalar&, const RealScalar&) {
2209 return x == Scalar(0);
2210 }
2211 EIGEN_DEVICE_FUNC static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar&) { return x == y; }
2212 EIGEN_DEVICE_FUNC static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar&) {
2213 return x <= y;
2214 }
2215};
2216
2217template <typename Scalar>
2218struct scalar_fuzzy_default_impl<Scalar, true, false> {
2219 using RealScalar = typename NumTraits<Scalar>::Real;
2220 template <typename OtherScalar>
2221 EIGEN_DEVICE_FUNC static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y,
2222 const RealScalar& prec) {
2223 return numext::abs2(x) <= numext::abs2(y) * prec * prec;
2224 }
2225 EIGEN_DEVICE_FUNC static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec) {
2226 return numext::abs2(x - y) <= numext::mini(numext::abs2(x), numext::abs2(y)) * prec * prec;
2227 }
2228};
2229
2230template <typename Scalar>
2231struct scalar_fuzzy_impl
2232 : scalar_fuzzy_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
2233
2234template <typename Scalar, typename OtherScalar>
2235EIGEN_DEVICE_FUNC inline bool isMuchSmallerThan(
2236 const Scalar& x, const OtherScalar& y,
2237 const typename NumTraits<Scalar>::Real& precision = NumTraits<Scalar>::dummy_precision()) {
2238 return scalar_fuzzy_impl<Scalar>::template isMuchSmallerThan<OtherScalar>(x, y, precision);
2239}
2240
2241template <typename Scalar>
2242EIGEN_DEVICE_FUNC inline bool isApprox(
2243 const Scalar& x, const Scalar& y,
2244 const typename NumTraits<Scalar>::Real& precision = NumTraits<Scalar>::dummy_precision()) {
2245 return scalar_fuzzy_impl<Scalar>::isApprox(x, y, precision);
2246}
2247
2248template <typename Scalar>
2249EIGEN_DEVICE_FUNC inline bool isApproxOrLessThan(
2250 const Scalar& x, const Scalar& y,
2251 const typename NumTraits<Scalar>::Real& precision = NumTraits<Scalar>::dummy_precision()) {
2252 return scalar_fuzzy_impl<Scalar>::isApproxOrLessThan(x, y, precision);
2253}
2254
2255/******************************************
2256*** The special case of the bool type ***
2257******************************************/
2258
2259template <>
2260struct scalar_fuzzy_impl<bool> {
2261 using RealScalar = bool;
2262
2263 template <typename OtherScalar>
2264 EIGEN_DEVICE_FUNC static inline bool isMuchSmallerThan(const bool& x, const bool&, const bool&) {
2265 return !x;
2266 }
2267
2268 EIGEN_DEVICE_FUNC static inline bool isApprox(bool x, bool y, bool) { return x == y; }
2269
2270 EIGEN_DEVICE_FUNC static inline bool isApproxOrLessThan(const bool& x, const bool& y, const bool&) {
2271 return (!x) || y;
2272 }
2273};
2274
2275} // end namespace internal
2276
2277// Default implementations that rely on other numext implementations
2278namespace internal {
2279
2280// Specialization for complex types that are not supported by std::expm1.
2281template <typename RealScalar>
2282struct expm1_impl<std::complex<RealScalar>> {
2283 EIGEN_STATIC_ASSERT_NON_INTEGER(RealScalar)
2284
2285 EIGEN_DEVICE_FUNC static inline std::complex<RealScalar> run(const std::complex<RealScalar>& x) {
2286 RealScalar xr = x.real();
2287 RealScalar xi = x.imag();
2288 // expm1(z) = exp(z) - 1
2289 // = exp(x + i * y) - 1
2290 // = exp(x) * (cos(y) + i * sin(y)) - 1
2291 // = exp(x) * cos(y) - 1 + i * exp(x) * sin(y)
2292 // Imag(expm1(z)) = exp(x) * sin(y)
2293 // Real(expm1(z)) = exp(x) * cos(y) - 1
2294 // = exp(x) * cos(y) - 1.
2295 // = expm1(x) + exp(x) * (cos(y) - 1)
2296 // = expm1(x) + exp(x) * (2 * sin(y / 2) ** 2)
2297 RealScalar erm1 = numext::expm1<RealScalar>(xr);
2298 RealScalar er = erm1 + RealScalar(1.);
2299 // C99 Annex G for exp, less one: (+inf, inf or NaN) -> (inf, NaN) and (-inf, inf or NaN) -> (-1, 0).
2300 // xi * 0 is that NaN, raising FE_INVALID for xi = inf as Annex G asks.
2301 if ((numext::isinf)(xr) && !(numext::isfinite)(xi)) {
2302 return xr > RealScalar(0) ? std::complex<RealScalar>(xr, xi * RealScalar(0))
2303 : std::complex<RealScalar>(RealScalar(-1), RealScalar(0));
2304 }
2305 const RealScalar s = numext::sin(xi);
2306 const RealScalar c = numext::cos(xi);
2307 // 1 - cos(y) = 2 sin(y / 2)^2, taken as sin(y)^2 / (1 + cos(y)) where the difference would cancel.
2308 const RealScalar one_minus_c = c > RealScalar(0) ? s * s / (RealScalar(1) + c) : RealScalar(1) - c;
2309 // exp(x) overflows where exp(x) cos(y) and exp(x) sin(y) need not; multiply exp(x / 2) in twice then, or
2310 // exp(x / 3) three times where that overflows too, and exp(x) sin(y) is finite only for denormal y.
2311 if ((numext::isinf)(er) && (numext::isfinite)(xr)) {
2312 const RealScalar h = numext::exp(xr / RealScalar(2));
2313 if (!(numext::isinf)(h)) {
2314 return std::complex<RealScalar>((h * c) * h - RealScalar(1), numext::is_exactly_zero(xi) ? xi : (h * s) * h);
2315 }
2316 const RealScalar h3 = numext::exp(xr / RealScalar(3));
2317 return std::complex<RealScalar>(((h3 * c) * h3) * h3, numext::is_exactly_zero(xi) ? xi : ((h3 * s) * h3) * h3);
2318 }
2319 // For x < 0 both terms of expm1(x) - exp(x) (1 - cos(y)) are <= 0, and the result is -1 exactly once exp(x) is
2320 // negligible. For x >= 0, expm1(x) cos(y) - (1 - cos(y)) has the smaller terms, and no inf - inf at x = inf.
2321 const RealScalar real_part = xr < RealScalar(0) ? erm1 - er * one_minus_c : erm1 * c - one_minus_c;
2322 // exp(x) * sin(0) is NaN for infinite exp(x); keep the exact zero instead.
2323 return std::complex<RealScalar>(real_part, numext::is_exactly_zero(xi) ? xi : er * s);
2324 }
2325};
2326
2327// 2^(a + ib) = 2^a (cos(t) + i sin(t)) with t = b ln(2). The modulus is the real exp2, so 2^n is exact. Rounding t
2328// would cost |t| eps in the phase, so b ln(2) = t_hi + t_lo is kept in double-word precision, ln(2) = ln2_hi + ln2_lo:
2329// cos(t) = cos(t_hi) + (cos(t_hi) (cos(t_lo) - 1) - sin(t_hi) sin(t_lo)),
2330// sin(t) = sin(t_hi) + (sin(t_hi) (cos(t_lo) - 1) + cos(t_hi) sin(t_lo)).
2331// Special values are those of exp(z) in C99 Annex G, since ln(2) > 0.
2332template <typename RealScalar>
2333struct complex_exp2_impl {
2334 EIGEN_STATIC_ASSERT_NON_INTEGER(RealScalar)
2335
2336 EIGEN_DEVICE_FUNC static inline std::complex<RealScalar> run(const std::complex<RealScalar>& x) {
2337 using Complex = std::complex<RealScalar>;
2338 const RealScalar a = x.real();
2339 const RealScalar b = x.imag();
2340 // (a, +-0) -> (2^a, +-0), also for a = +-inf and NaN.
2341 if (numext::is_exactly_zero(b)) return Complex(numext::exp2(a), b);
2342 // (+inf, inf or NaN) -> (inf, NaN) and (-inf, inf or NaN) -> (0, 0), where 2^a cos(t) would be NaN.
2343 // b * 0 is that NaN, raising FE_INVALID for b = inf as Annex G asks.
2344 if ((numext::isinf)(a) && !(numext::isfinite)(b)) {
2345 return a > RealScalar(0) ? Complex(a, b * RealScalar(0)) : Complex(RealScalar(0), RealScalar(0));
2346 }
2347 // ln2_lo = (double(ln(2)) - ln2_hi) + (ln(2) - double(ln(2))), the first difference exact in Wide.
2348 using Wide =
2349 std::conditional_t<(NumTraits<RealScalar>::digits() > NumTraits<double>::digits()), RealScalar, double>;
2350 const RealScalar ln2_hi = RealScalar(EIGEN_LN2);
2351 const RealScalar ln2_lo = RealScalar((Wide(static_cast<double>(EIGEN_LN2)) - Wide(ln2_hi)) +
2352 Wide(2.3190468138462996154948554638754786504e-17L));
2353 const RealScalar t_hi = b * ln2_hi;
2354 const RealScalar t_lo = numext::fma(b, ln2_hi, -t_hi) + b * ln2_lo;
2355 // |t_lo| <= ulp(t_hi) / 2, so the Taylor terms below are exact to rounding unless |t| exceeds about 2 / sqrt(eps).
2356 // The phase is within an ulp of b ln(2) while |t| < 1 / eps.
2357 RealScalar sin_lo = t_lo;
2358 RealScalar cos_lo_m1 = -t_lo * t_lo / RealScalar(2);
2359 if (numext::abs(t_lo) > numext::sqrt(NumTraits<RealScalar>::epsilon())) {
2360 sin_lo = numext::sin(t_lo);
2361 const RealScalar sin_half = numext::sin(t_lo / RealScalar(2));
2362 cos_lo_m1 = RealScalar(-2) * sin_half * sin_half;
2363 }
2364 const RealScalar cos_hi = numext::cos(t_hi);
2365 const RealScalar sin_hi = numext::sin(t_hi);
2366 const RealScalar c = cos_hi + (cos_hi * cos_lo_m1 - sin_hi * sin_lo);
2367 const RealScalar s = sin_hi + (sin_hi * cos_lo_m1 + cos_hi * sin_lo);
2368 // 2^a overflows or is denormal where 2^a cos(t) and 2^a sin(t) need not be; apply 2^e after the products then.
2369 // Past 2 e_max, a part is finite only for denormal t.
2370 const int e_max = NumTraits<RealScalar>::max_exponent() - 1;
2371 const int e = a > RealScalar(e_max) ? (a > RealScalar(2 * e_max) ? 2 * e_max : e_max)
2372 : (a < RealScalar(NumTraits<RealScalar>::min_exponent()) ? -e_max : 0);
2373 const RealScalar m = numext::exp2(a - RealScalar(e));
2374 if (numext::abs(t_hi) < (numext::numeric_limits<RealScalar>::min)()) {
2375 // A denormal t has lost bits that 2^a can bring back. There cos(t) = 1 and sin(t) = t, taken from b 2^digits.
2376 const int p = NumTraits<RealScalar>::digits();
2377 const RealScalar b_p = numext::ldexp(b, p);
2378 return Complex(numext::ldexp(m, e), numext::ldexp(m * numext::fma(b_p, ln2_hi, b_p * ln2_lo), e - p));
2379 }
2380 if (e == 0) return Complex(m * c, m * s);
2381 return Complex(numext::ldexp(m * c, e), numext::ldexp(m * s, e));
2382 }
2383};
2384
2385template <typename T>
2386struct rsqrt_impl {
2387// C4804: unsafe use of type 'bool' in operation. Unavoidable when instantiated with T=bool.
2388#if EIGEN_COMP_MSVC
2389#pragma warning(push)
2390#pragma warning(disable : 4804)
2391#endif
2392 EIGEN_DEVICE_FUNC static EIGEN_ALWAYS_INLINE T run(const T& x) { return T(1) / numext::sqrt(x); }
2393#if EIGEN_COMP_MSVC
2394#pragma warning(pop)
2395#endif
2396};
2397
2398#if defined(EIGEN_GPU_COMPILE_PHASE)
2399template <typename T>
2400struct conj_impl<std::complex<T>, true> {
2401 EIGEN_DEVICE_FUNC static inline std::complex<T> run(const std::complex<T>& x) {
2402 return std::complex<T>(numext::real(x), -numext::imag(x));
2403 }
2404};
2405#endif
2406
2407// Complex multiply and division operators.
2408// Note that these do not handle the case if inf+NaNi, which is considered an infinity.
2409// This is for consistency with our standard pmul, pdiv implementations.
2410template <typename T>
2411EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> complex_multiply(const std::complex<T>& a,
2412 const std::complex<T>& b) {
2413 const T a_real = numext::real(a);
2414 const T a_imag = numext::imag(a);
2415 const T b_real = numext::real(b);
2416 const T b_imag = numext::imag(b);
2417 return std::complex<T>(a_real * b_real - a_imag * b_imag, a_imag * b_real + a_real * b_imag);
2418}
2419
2420template <typename T>
2421EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> complex_divide_fast(const std::complex<T>& a,
2422 const std::complex<T>& b) {
2423 const T a_real = numext::real(a);
2424 const T a_imag = numext::imag(a);
2425 const T b_real = numext::real(b);
2426 const T b_imag = numext::imag(b);
2427 const T norm = (b_real * b_real + b_imag * b_imag);
2428 return std::complex<T>((a_real * b_real + a_imag * b_imag) / norm, (a_imag * b_real - a_real * b_imag) / norm);
2429}
2430
2431template <typename T>
2432EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> complex_divide_smith(const std::complex<T>& a,
2433 const std::complex<T>& b) {
2434 const T a_real = numext::real(a);
2435 const T a_imag = numext::imag(a);
2436 const T b_real = numext::real(b);
2437 const T b_imag = numext::imag(b);
2438 // Smith's complex division (https://arxiv.org/pdf/1210.4539.pdf),
2439 // guards against over/under-flow.
2440 const bool scale_imag = numext::abs(b_imag) <= numext::abs(b_real);
2441 const T rscale = scale_imag ? T(1) : b_real / b_imag;
2442 const T iscale = scale_imag ? b_imag / b_real : T(1);
2443 const T denominator = b_real * rscale + b_imag * iscale;
2444 return std::complex<T>((a_real * rscale + a_imag * iscale) / denominator,
2445 (a_imag * rscale - a_real * iscale) / denominator);
2446}
2447
2448template <typename T>
2449EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> complex_divide(const std::complex<T>& a,
2450 const std::complex<T>& b) {
2451#if EIGEN_FAST_MATH
2452 return complex_divide_fast(a, b);
2453#else
2454 return complex_divide_smith(a, b);
2455#endif
2456}
2457
2458} // end namespace internal
2459
2460namespace numext {
2461
2470template <typename LhsScalar, typename RhsScalar>
2471EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE auto divide(const LhsScalar& a, const RhsScalar& b) -> decltype(a / b) {
2472 return a / b;
2473}
2474
2475template <typename T>
2476EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> divide(const std::complex<T>& a, const std::complex<T>& b) {
2477 return internal::complex_divide_smith(a, b);
2478}
2479
2480template <typename T>
2481EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> divide(const T& a, const std::complex<T>& b) {
2482 return internal::complex_divide_smith(std::complex<T>(a, T(0)), b);
2483}
2484
2485template <typename T>
2486EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> divide(const std::complex<T>& a, const T& b) {
2487 return a / b;
2488}
2489
2490} // end namespace numext
2491
2492} // end namespace Eigen
2493
2494#endif // EIGEN_MATHFUNCTIONS_H