fast-math support

Author: serge-sans-paille

include/xsimd/arch/common/xsimd_common_complex.hpp

@@ -78,11 +78,15 @@ namespace xsimd
             using batch_type = complex_batch_type_t<batch<T, A>>;
             using real_batch = typename batch_type::real_batch;
             using real_value_type = typename real_batch::value_type;
+#ifdef __FAST_MATH__
+            return { self };
+#else
             auto cond = xsimd::isinf(real(self)) || xsimd::isinf(imag(self));
             return select(cond,
                           batch_type(constants::infinity<real_batch>(),
                                      copysign(real_batch(real_value_type(0)), imag(self))),
                           batch_type(self));
+#endif
         }
 
         template <class A, class T>

include/xsimd/arch/common/xsimd_common_logical.hpp

@@ -58,10 +58,10 @@ namespace xsimd
                 else
                 {
                     // https://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel
-                    m = m - ((m >> 1) & (uint64_t) ~(uint64_t)0 / 3); // temp
-                    m = (m & (uint64_t) ~(uint64_t)0 / 15 * 3) + ((m >> 2) & (uint64_t) ~(uint64_t)0 / 15 * 3); // temp
-                    m = (m + (m >> 4)) & (uint64_t) ~(uint64_t)0 / 255 * 15; // temp
-                    return (m * ((uint64_t) ~(uint64_t)0 / 255)) >> (sizeof(uint64_t) - 1) * CHAR_BIT; // count
+                    m = m - ((m >> 1) & (uint64_t)~(uint64_t)0 / 3); // temp
+                    m = (m & (uint64_t)~(uint64_t)0 / 15 * 3) + ((m >> 2) & (uint64_t)~(uint64_t)0 / 15 * 3); // temp
+                    m = (m + (m >> 4)) & (uint64_t)~(uint64_t)0 / 255 * 15; // temp
+                    return (m * ((uint64_t)~(uint64_t)0 / 255)) >> (sizeof(uint64_t) - 1) * CHAR_BIT; // count
                 }
 #endif
             }
@@ -124,12 +124,22 @@ namespace xsimd
         template <class A>
         XSIMD_INLINE batch_bool<float, A> isinf(batch<float, A> const& self, requires_arch<common>) noexcept
         {
+#ifdef __FAST_MATH__
+            (void)self;
+            return { false };
+#else
             return abs(self) == std::numeric_limits<float>::infinity();
+#endif
         }
         template <class A>
         XSIMD_INLINE batch_bool<double, A> isinf(batch<double, A> const& self, requires_arch<common>) noexcept
         {
+#ifdef __FAST_MATH__
+            (void)self;
+            return { false };
+#else
             return abs(self) == std::numeric_limits<double>::infinity();
+#endif
         }
 
         // isfinite

include/xsimd/arch/common/xsimd_common_math.hpp

@@ -896,7 +896,9 @@ namespace xsimd
                 batch_type k = reducer_t::reduce(self, x);
                 x = reducer_t::approx(x);
                 x = select(self <= reducer_t::minlog(), batch_type(0.), ldexp(x, to_int(k)));
+#ifndef __FAST_MATH__
                 x = select(self >= reducer_t::maxlog(), constants::infinity<batch_type>(), x);
+#endif
                 return x;
             }
 
@@ -910,7 +912,9 @@ namespace xsimd
                 batch_type c = reducer_t::approx(x);
                 c = reducer_t::finalize(x, c, hi, lo);
                 c = select(self <= reducer_t::minlog(), batch_type(0.), ldexp(c, to_int(k)));
+#ifndef __FAST_MATH__
                 c = select(self >= reducer_t::maxlog(), constants::infinity<batch_type>(), c);
+#endif
                 return c;
             }
         }
@@ -1014,11 +1018,11 @@ namespace xsimd
         XSIMD_INLINE batch<T, A> expm1(batch<T, A> const& self, requires_arch<common>) noexcept
         {
             using batch_type = batch<T, A>;
-            return select(self < constants::logeps<batch_type>(),
-                          batch_type(-1.),
-                          select(self > constants::maxlog<batch_type>(),
-                                 constants::infinity<batch_type>(),
-                                 detail::expm1(self)));
+            auto x = detail::expm1(self);
+#ifndef __FAST_MATH__
+            x = select(self > constants::maxlog<batch_type>(), constants::infinity<batch_type>(), x);
+#endif
+            return select(self < constants::logeps<batch_type>(), batch_type(-1.), x);
         }
 
         template <class A, class T>
@@ -1245,12 +1249,20 @@ namespace xsimd
                     batch_type r1 = other(q);
                     if (any(ltza))
                     {
+#ifdef __FAST_MATH__
+                        r = negative(q, r1);
+#else
                         r = select(inf_result, constants::infinity<batch_type>(), negative(q, r1));
+#endif
                         if (all(ltza))
                             return r;
                     }
                     batch_type r2 = select(ltza, r, r1);
+#ifdef __FAST_MATH__
+                    return r2;
+#else
                     return select(a == constants::minusinfinity<batch_type>(), constants::nan<batch_type>(), select(inf_result, constants::infinity<batch_type>(), r2));
+#endif
                 }
 
             private:
@@ -1371,7 +1383,11 @@ namespace xsimd
                     }
                     batch_type r1 = other(a);
                     batch_type r2 = select(test, r, r1);
+#ifdef __FAST_MATH__
+                    return r2;
+#else
                     return select(a == constants::minusinfinity<batch_type>(), constants::nan<batch_type>(), select(inf_result, constants::infinity<batch_type>(), r2));
+#endif
                 }
 
             private:
@@ -1479,12 +1495,12 @@ namespace xsimd
             batch_type hfsq = batch_type(0.5) * f * f;
             batch_type dk = to_float(k);
             batch_type r = fma(dk, constants::log_2hi<batch_type>(), fma(s, (hfsq + R), dk * constants::log_2lo<batch_type>()) - hfsq + f);
-#ifndef XSIMD_NO_INFINITIES
-            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
+#ifdef __FAST_MATH__
+            return r;
 #else
-            batch_type zz = select(isnez, r, constants::minusinfinity<batch_type>());
-#endif
+            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
             return select(!(self >= batch_type(0.)), constants::nan<batch_type>(), zz);
+#endif
         }
 
         template <class A>
@@ -1522,12 +1538,12 @@ namespace xsimd
             batch_type t2 = z * detail::horner<batch_type, 0x3fe5555555555593ll, 0x3fd2492494229359ll, 0x3fc7466496cb03dell, 0x3fc2f112df3e5244ll>(w);
             batch_type R = t2 + t1;
             batch_type r = fma(dk, constants::log_2hi<batch_type>(), fma(s, (hfsq + R), dk * constants::log_2lo<batch_type>()) - hfsq + f);
-#ifndef XSIMD_NO_INFINITIES
-            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
+#ifdef __FAST_MATH__
+            return r;
 #else
-            batch_type zz = select(isnez, r, constants::minusinfinity<batch_type>());
-#endif
+            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
             return select(!(self >= batch_type(0.)), constants::nan<batch_type>(), zz);
+#endif
         }
 
         template <class A, class T>
@@ -1569,12 +1585,12 @@ namespace xsimd
             batch_type hfsq = batch_type(0.5) * f * f;
             batch_type dk = to_float(k);
             batch_type r = fma(fms(s, hfsq + R, hfsq) + f, constants::invlog_2<batch_type>(), dk);
-#ifndef XSIMD_NO_INFINITIES
-            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
+#ifdef __FAST_MATH__
+            return r;
 #else
-            batch_type zz = select(isnez, r, constants::minusinfinity<batch_type>());
-#endif
+            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
             return select(!(self >= batch_type(0.)), constants::nan<batch_type>(), zz);
+#endif
         }
 
         template <class A>
@@ -1617,12 +1633,12 @@ namespace xsimd
             val_lo += (dk - w1) + val_hi;
             val_hi = w1;
             batch_type r = val_lo + val_hi;
-#ifndef XSIMD_NO_INFINITIES
-            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
+#ifdef __FAST_MATH__
+            return r;
 #else
-            batch_type zz = select(isnez, r, constants::minusinfinity<batch_type>());
-#endif
+            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
             return select(!(self >= batch_type(0.)), constants::nan<batch_type>(), zz);
+#endif
         }
 
         namespace detail
@@ -1757,12 +1773,12 @@ namespace xsimd
             val_lo += (y - w1) + val_hi;
             val_hi = w1;
             batch_type r = val_lo + val_hi;
-#ifndef XSIMD_NO_INFINITIES
-            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
+#ifdef __FAST_MATH__
+            return r;
 #else
-            batch_type zz = select(isnez, r, constants::minusinfinity<batch_type>());
-#endif
+            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
             return select(!(self >= batch_type(0.)), constants::nan<batch_type>(), zz);
+#endif
         }
 
         template <class A, class T>
@@ -1805,12 +1821,12 @@ namespace xsimd
             /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
             batch_type c = select(batch_bool_cast<float>(k >= i_type(2)), batch_type(1.) - (uf - self), self - (uf - batch_type(1.))) / uf;
             batch_type r = fma(dk, constants::log_2hi<batch_type>(), fma(s, (hfsq + R), dk * constants::log_2lo<batch_type>() + c) - hfsq + f);
-#ifndef XSIMD_NO_INFINITIES
-            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
+#ifdef __FAST_MATH__
+            return r;
 #else
-            batch_type zz = select(isnez, r, constants::minusinfinity<batch_type>());
-#endif
+            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
             return select(!(uf >= batch_type(0.)), constants::nan<batch_type>(), zz);
+#endif
         }
 
         template <class A>
@@ -1838,12 +1854,12 @@ namespace xsimd
             batch_type R = t2 + t1;
             batch_type dk = to_float(k);
             batch_type r = fma(dk, constants::log_2hi<batch_type>(), fma(s, hfsq + R, dk * constants::log_2lo<batch_type>() + c) - hfsq + f);
-#ifndef XSIMD_NO_INFINITIES
-            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
+#ifdef __FAST_MATH__
+            return r;
 #else
-            batch_type zz = select(isnez, r, constants::minusinfinity<batch_type>());
-#endif
+            batch_type zz = select(isnez, select(self == constants::infinity<batch_type>(), constants::infinity<batch_type>(), r), constants::minusinfinity<batch_type>());
             return select(!(uf >= batch_type(0.)), constants::nan<batch_type>(), zz);
+#endif
         }
 
         template <class A, class T>
@@ -1980,13 +1996,21 @@ namespace xsimd
                 static XSIMD_INLINE batch_type next(const batch_type& b) noexcept
                 {
                     batch_type n = ::xsimd::bitwise_cast<T>(::xsimd::bitwise_cast<int_type>(b) + int_type(1));
+#ifdef __FAST_MATH__
+                    return n;
+#else
                     return select(b == constants::infinity<batch_type>(), b, n);
+#endif
                 }
 
                 static XSIMD_INLINE batch_type prev(const batch_type& b) noexcept
                 {
                     batch_type p = ::xsimd::bitwise_cast<T>(::xsimd::bitwise_cast<int_type>(b) - int_type(1));
+#ifdef __FAST_MATH__
+                    return p;
+#else
                     return select(b == constants::minusinfinity<batch_type>(), b, p);
+#endif
                 }
             };
         }
@@ -2355,10 +2379,12 @@ namespace xsimd
                 y *= v;
                 y = select(test, y, y * v);
                 y *= constants::sqrt_2pi<batch_type>() * w;
-#ifndef XSIMD_NO_INFINITIES
+#ifdef __FAST_MATH__
+                return y;
+#else
                 y = select(isinf(x), x, y);
-#endif
                 return select(x > stirlinglargelim, constants::infinity<batch_type>(), y);
+#endif
             }
 
             /* origin: boost/simd/arch/common/detail/common/gamma_kernel.hpp */
@@ -2501,7 +2527,11 @@ namespace xsimd
             }
             batch_type r1 = detail::tgamma_other(self, test);
             batch_type r2 = select(test, r, r1);
+#ifdef __FAST_MATH__
+            return r2;
+#else
             return select(self == batch_type(0.), copysign(constants::infinity<batch_type>(), self), select(nan_result, constants::nan<batch_type>(), r2));
+#endif
         }
 
     }

include/xsimd/arch/common/xsimd_common_trigo.hpp

@@ -310,10 +310,13 @@ namespace xsimd
             num = x2 + num * num;
             real_batch den = y - one;
             den = x2 + den * den;
-            batch_type res = select((x == real_batch(0.)) && (y == real_batch(1.)),
-                                    batch_type(real_batch(0.), constants::infinity<real_batch>()),
-                                    batch_type(w, 0.25 * log(num / den)));
-            return res;
+#ifdef __FAST_MATH__
+            return batch_type(w, 0.25 * log(num / den));
+#else
+            return select((x == real_batch(0.)) && (y == real_batch(1.)),
+                          batch_type(real_batch(0.), constants::infinity<real_batch>()),
+                          batch_type(w, 0.25 * log(num / den)));
+#endif
         }
 
         // atanh
@@ -583,12 +586,14 @@ namespace xsimd
                         for (std::size_t i = 0; i < size; ++i)
                         {
                             double arg = args[i];
+#ifndef __FAST_MATH__
                             if (arg == std::numeric_limits<value_type>::infinity())
                             {
                                 tmp[i] = 0.;
                                 txr[i] = std::numeric_limits<value_type>::quiet_NaN();
                             }
                             else
+#endif
                             {
                                 double y[2];
                                 std::int32_t n = ::xsimd::detail::__ieee754_rem_pio2(arg, y);
@@ -841,11 +846,15 @@ namespace xsimd
             using batch_type = batch<std::complex<T>, A>;
             using real_batch = typename batch_type::real_batch;
             real_batch d = cos(2 * z.real()) + cosh(2 * z.imag());
-            batch_type winf(constants::infinity<real_batch>(), constants::infinity<real_batch>());
             real_batch wreal = sin(2 * z.real()) / d;
             real_batch wimag = sinh(2 * z.imag());
+#ifdef __FAST_MATH__
+            return batch_type(wreal, real_batch(1.)), batch_type(wreal, wimag / d);
+#else
+            batch_type winf(constants::infinity<real_batch>(), constants::infinity<real_batch>());
             batch_type wres = select(isinf(wimag), batch_type(wreal, real_batch(1.)), batch_type(wreal, wimag / d));
             return select(d == real_batch(0.), winf, wres);
+#endif
         }
 
         // tanh

include/xsimd/arch/xsimd_constants.hpp

@@ -61,7 +61,10 @@ namespace xsimd
 #pragma GCC push_options
 #pragma GCC optimize("signed-zeros")
 #endif
+#ifndef __FAST_MATH__
         XSIMD_DEFINE_CONSTANT(infinity, (std::numeric_limits<float>::infinity()), (std::numeric_limits<double>::infinity()))
+        XSIMD_DEFINE_CONSTANT(minusinfinity, (-infinity<float>()), (-infinity<double>()))
+#endif
         XSIMD_DEFINE_CONSTANT(invlog_2, 1.442695040888963407359924681001892137426645954152986f, 1.442695040888963407359924681001892137426645954152986)
         XSIMD_DEFINE_CONSTANT_HEX(invlog_2hi, 0x3fb8b000, 0x3ff7154765200000)
         XSIMD_DEFINE_CONSTANT_HEX(invlog_2lo, 0xb9389ad4, 0x3de705fc2eefa200)
@@ -83,7 +86,6 @@ namespace xsimd
         XSIMD_DEFINE_CONSTANT(minlog, -88.3762626647949f, -708.3964185322641)
         XSIMD_DEFINE_CONSTANT(minlog2, -127.0f, -1023.)
         XSIMD_DEFINE_CONSTANT(minlog10, -37.89999771118164f, -308.2547155599167)
-        XSIMD_DEFINE_CONSTANT(minusinfinity, (-infinity<float>()), (-infinity<double>()))
         XSIMD_DEFINE_CONSTANT_HEX(nan, 0xffffffff, 0xffffffffffffffff)
         XSIMD_DEFINE_CONSTANT_HEX(oneosqrteps, 0x453504f3, 0x4190000000000000)
         XSIMD_DEFINE_CONSTANT_HEX(oneotwoeps, 0x4a800000, 0x4320000000000000)

include/xsimd/types/xsimd_api.hpp

@@ -1114,6 +1114,7 @@ namespace xsimd
         return kernel::incr_if<A>(x, mask, A {});
     }
 
+#ifndef __FAST_MATH__
     /**
      * @ingroup batch_constant
      *
@@ -1128,6 +1129,7 @@ namespace xsimd
         detail::static_check_supported_config<T, A>();
         return B(std::numeric_limits<T>::infinity());
     }
+#endif
 
     /**
      * @ingroup batch_data_transfer