arch/common/xsimd_common_trigo.hpp

0001 /***************************************************************************
0002  * Copyright (c) Johan Mabille, Sylvain Corlay, Wolf Vollprecht and         *
0003  * Martin Renou                                                             *
0004  * Copyright (c) QuantStack                                                 *
0005  * Copyright (c) Serge Guelton                                              *
0006  *                                                                          *
0007  * Distributed under the terms of the BSD 3-Clause License.                 *
0008  *                                                                          *
0009  * The full license is in the file LICENSE, distributed with this software. *
0010  ****************************************************************************/
0011
0012 #ifndef XSIMD_COMMON_TRIGO_HPP
0013 #define XSIMD_COMMON_TRIGO_HPP
0014
0015 #include "./xsimd_common_details.hpp"
0016
0017 #include <array>
0018
0019 namespace xsimd
0020 {
0021
0022     namespace kernel
0023     {
0024         /* origin: boost/simd/arch/common/detail/simd/trig_base.hpp */
0025         /*
0026          * ====================================================
0027          * copyright 2016 NumScale SAS
0028          *
0029          * Distributed under the Boost Software License, Version 1.0.
0030          * (See copy at http://boost.org/LICENSE_1_0.txt)
0031          * ====================================================
0032          */
0033
0034         using namespace types;
0035
0036         // acos
0037         template <class A, class T>
0038         XSIMD_INLINE batch<T, A> acos(batch<T, A> const& self, requires_arch<common>) noexcept
0039         {
0040             using batch_type = batch<T, A>;
0041             batch_type x = abs(self);
0042             auto x_larger_05 = x > batch_type(0.5);
0043             x = select(x_larger_05, sqrt(fma(batch_type(-0.5), x, batch_type(0.5))), self);
0044             x = asin(x);
0045             x = select(x_larger_05, x + x, x);
0046             x = select(self < batch_type(-0.5), constants::pi<batch_type>() - x, x);
0047             return select(x_larger_05, x, constants::pio2<batch_type>() - x);
0048         }
0049         template <class A, class T>
0050         XSIMD_INLINE batch<std::complex<T>, A> acos(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
0051         {
0052             using batch_type = batch<std::complex<T>, A>;
0053             using real_batch = typename batch_type::real_batch;
0054             batch_type tmp = asin(z);
0055             return { constants::pio2<real_batch>() - tmp.real(), -tmp.imag() };
0056         }
0057
0058         // acosh
0059         /* origin: boost/simd/arch/common/simd/function/acosh.hpp */
0060         /*
0061          * ====================================================
0062          * copyright 2016 NumScale SAS
0063          *
0064          * Distributed under the Boost Software License, Version 1.0.
0065          * (See copy at http://boost.org/LICENSE_1_0.txt)
0066          * ====================================================
0067          */
0068         template <class A, class T>
0069         XSIMD_INLINE batch<T, A> acosh(batch<T, A> const& self, requires_arch<common>) noexcept
0070         {
0071             using batch_type = batch<T, A>;
0072             batch_type x = self - batch_type(1.);
0073             auto test = x > constants::oneotwoeps<batch_type>();
0074             batch_type z = select(test, self, x + sqrt(x + x + x * x));
0075             batch_type l1pz = log1p(z);
0076             return select(test, l1pz + constants::log_2<batch_type>(), l1pz);
0077         }
0078         template <class A, class T>
0079         XSIMD_INLINE batch<std::complex<T>, A> acosh(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
0080         {
0081             using batch_type = batch<std::complex<T>, A>;
0082             batch_type w = acos(z);
0083             w = batch_type(-w.imag(), w.real());
0084             return w;
0085         }
0086
0087         // asin
0088         template <class A>
0089         XSIMD_INLINE batch<float, A> asin(batch<float, A> const& self, requires_arch<common>) noexcept
0090         {
0091             using batch_type = batch<float, A>;
0092             batch_type x = abs(self);
0093             batch_type sign = bitofsign(self);
0094             auto x_larger_05 = x > batch_type(0.5);
0095             batch_type z = select(x_larger_05, batch_type(0.5) * (batch_type(1.) - x), x * x);
0096             x = select(x_larger_05, sqrt(z), x);
0097             batch_type z1 = detail::horner<batch_type,
0098                                            0x3e2aaae4,
0099                                            0x3d9980f6,
0100                                            0x3d3a3ec7,
0101                                            0x3cc617e3,
0102                                            0x3d2cb352>(z);
0103             z1 = fma(z1, z * x, x);
0104             z = select(x_larger_05, constants::pio2<batch_type>() - (z1 + z1), z1);
0105             return z ^ sign;
0106         }
0107         template <class A>
0108         XSIMD_INLINE batch<double, A> asin(batch<double, A> const& self, requires_arch<common>) noexcept
0109         {
0110             using batch_type = batch<double, A>;
0111             batch_type x = abs(self);
0112             auto small_cond = x < constants::sqrteps<batch_type>();
0113             batch_type ct1 = batch_type(bit_cast<double>(int64_t(0x3fe4000000000000)));
0114             batch_type zz1 = batch_type(1.) - x;
0115             batch_type vp = zz1 * detail::horner<batch_type, 0x403c896240f3081dull, 0xc03991aaac01ab68ull, 0x401bdff5baf33e6aull, 0xbfe2079259f9290full, 0x3f684fc3988e9f08ull>(zz1) / detail::horner1<batch_type, 0x40756709b0b644beull, 0xc077fe08959063eeull, 0x40626219af6a7f42ull, 0xc035f2a2b6bf5d8cull>(zz1);
0116             zz1 = sqrt(zz1 + zz1);
0117             batch_type z = constants::pio4<batch_type>() - zz1;
0118             zz1 = fms(zz1, vp, constants::pio_2lo<batch_type>());
0119             z = z - zz1;
0120             zz1 = z + constants::pio4<batch_type>();
0121             batch_type zz2 = self * self;
0122             z = zz2 * detail::horner<batch_type, 0xc020656c06ceafd5ull, 0x40339007da779259ull, 0xc0304331de27907bull, 0x4015c74b178a2dd9ull, 0xbfe34341333e5c16ull, 0x3f716b9b0bd48ad3ull>(zz2) / detail::horner1<batch_type, 0xc04898220a3607acull, 0x4061705684ffbf9dull, 0xc06265bb6d3576d7ull, 0x40519fc025fe9054ull, 0xc02d7b590b5e0eabull>(zz2);
0123             zz2 = fma(x, z, x);
0124             return select(x > batch_type(1.), constants::nan<batch_type>(),
0125                           select(small_cond, x,
0126                                  select(x > ct1, zz1, zz2))
0127                               ^ bitofsign(self));
0128         }
0129         template <class A, class T>
0130         XSIMD_INLINE batch<std::complex<T>, A> asin(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
0131         {
0132             using batch_type = batch<std::complex<T>, A>;
0133             using real_batch = typename batch_type::real_batch;
0134             real_batch x = z.real();
0135             real_batch y = z.imag();
0136
0137             batch_type ct(-y, x);
0138             batch_type zz(real_batch(1.) - (x - y) * (x + y), -2 * x * y);
0139             zz = log(ct + sqrt(zz));
0140             batch_type resg(zz.imag(), -zz.real());
0141
0142             return select(y == real_batch(0.),
0143                           select(fabs(x) > real_batch(1.),
0144                                  batch_type(constants::pio2<real_batch>(), real_batch(0.)),
0145                                  batch_type(asin(x), real_batch(0.))),
0146                           resg);
0147         }
0148
0149         // asinh
0150         /* origin: boost/simd/arch/common/simd/function/asinh.hpp */
0151         /*
0152          * ====================================================
0153          * copyright 2016 NumScale SAS
0154          *
0155          * Distributed under the Boost Software License, Version 1.0.
0156          * (See copy at http://boost.org/LICENSE_1_0.txt)
0157          * ====================================================
0158          */
0159         namespace detail
0160         {
0161             template <class A, class T, class = std::enable_if_t<std::is_integral<T>::value>>
0162             XSIMD_INLINE batch<T, A>
0163             average(const batch<T, A>& x1, const batch<T, A>& x2) noexcept
0164             {
0165                 return (x1 & x2) + ((x1 ^ x2) >> 1);
0166             }
0167
0168             template <class A, class T>
0169             XSIMD_INLINE batch<T, A>
0170             averagef(const batch<T, A>& x1, const batch<T, A>& x2) noexcept
0171             {
0172                 using batch_type = batch<T, A>;
0173                 return fma(x1, batch_type(0.5), x2 * batch_type(0.5));
0174             }
0175             template <class A>
0176             XSIMD_INLINE batch<float, A> average(batch<float, A> const& x1, batch<float, A> const& x2) noexcept
0177             {
0178                 return averagef(x1, x2);
0179             }
0180             template <class A>
0181             XSIMD_INLINE batch<double, A> average(batch<double, A> const& x1, batch<double, A> const& x2) noexcept
0182             {
0183                 return averagef(x1, x2);
0184             }
0185         }
0186         template <class A>
0187         XSIMD_INLINE batch<float, A> asinh(batch<float, A> const& self, requires_arch<common>) noexcept
0188         {
0189             using batch_type = batch<float, A>;
0190             batch_type x = abs(self);
0191             auto lthalf = x < batch_type(0.5);
0192             batch_type x2 = x * x;
0193             batch_type bts = bitofsign(self);
0194             batch_type z(0.);
0195             if (any(lthalf))
0196             {
0197                 z = detail::horner<batch_type,
0198                                    0x3f800000,
0199                                    0xbe2aa9ad,
0200                                    0x3d9949b1,
0201                                    0xbd2ee581,
0202                                    0x3ca4d6e6>(x2)
0203                     * x;
0204                 if (all(lthalf))
0205                     return z ^ bts;
0206             }
0207             batch_type tmp = select(x > constants::oneosqrteps<batch_type>(), x, detail::average(x, hypot(batch_type(1.), x)));
0208 #ifndef XSIMD_NO_NANS
0209             return select(isnan(self), constants::nan<batch_type>(), select(lthalf, z, log(tmp) + constants::log_2<batch_type>()) ^ bts);
0210 #else
0211             return select(lthalf, z, log(tmp) + constants::log_2<batch_type>()) ^ bts;
0212 #endif
0213         }
0214         template <class A>
0215         XSIMD_INLINE batch<double, A> asinh(batch<double, A> const& self, requires_arch<common>) noexcept
0216         {
0217             using batch_type = batch<double, A>;
0218             batch_type x = abs(self);
0219             auto test = x > constants::oneosqrteps<batch_type>();
0220             batch_type z = select(test, x - batch_type(1.), x + x * x / (batch_type(1.) + hypot(batch_type(1.), x)));
0221 #ifndef XSIMD_NO_INFINITIES
0222             z = select(x == constants::infinity<batch_type>(), x, z);
0223 #endif
0224             batch_type l1pz = log1p(z);
0225             z = select(test, l1pz + constants::log_2<batch_type>(), l1pz);
0226             return bitofsign(self) ^ z;
0227         }
0228         template <class A, class T>
0229         XSIMD_INLINE batch<std::complex<T>, A> asinh(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
0230         {
0231             using batch_type = batch<std::complex<T>, A>;
0232             batch_type w = asin(batch_type(-z.imag(), z.real()));
0233             w = batch_type(w.imag(), -w.real());
0234             return w;
0235         }
0236
0237         // atan
0238         namespace detail
0239         {
0240             template <class A>
0241             static XSIMD_INLINE batch<float, A> kernel_atan(const batch<float, A>& x, const batch<float, A>& recx) noexcept
0242             {
0243                 using batch_type = batch<float, A>;
0244                 const auto flag1 = x < constants::tan3pio8<batch_type>();
0245                 const auto flag2 = (x >= batch_type(bit_cast<float>((uint32_t)0x3ed413cd))) && flag1;
0246                 batch_type yy = select(flag1, batch_type(0.), constants::pio2<batch_type>());
0247                 yy = select(flag2, constants::pio4<batch_type>(), yy);
0248                 batch_type xx = select(flag1, x, -recx);
0249                 xx = select(flag2, (x - batch_type(1.)) / (x + batch_type(1.)), xx);
0250                 const batch_type z = xx * xx;
0251                 batch_type z1 = detail::horner<batch_type,
0252                                                0xbeaaaa2aul,
0253                                                0x3e4c925ful,
0254                                                0xbe0e1b85ul,
0255                                                0x3da4f0d1ul>(z);
0256                 z1 = fma(xx, z1 * z, xx);
0257                 z1 = select(flag2, z1 + constants::pio_4lo<batch_type>(), z1);
0258                 z1 = select(!flag1, z1 + constants::pio_2lo<batch_type>(), z1);
0259                 return yy + z1;
0260             }
0261             template <class A>
0262             static XSIMD_INLINE batch<double, A> kernel_atan(const batch<double, A>& x, const batch<double, A>& recx) noexcept
0263             {
0264                 using batch_type = batch<double, A>;
0265                 const auto flag1 = x < constants::tan3pio8<batch_type>();
0266                 const auto flag2 = (x >= constants::tanpio8<batch_type>()) && flag1;
0267                 batch_type yy = select(flag1, batch_type(0.), constants::pio2<batch_type>());
0268                 yy = select(flag2, constants::pio4<batch_type>(), yy);
0269                 batch_type xx = select(flag1, x, -recx);
0270                 xx = select(flag2, (x - batch_type(1.)) / (x + batch_type(1.)), xx);
0271                 batch_type z = xx * xx;
0272                 z *= detail::horner<batch_type,
0273                                     0xc0503669fd28ec8eull,
0274                                     0xc05eb8bf2d05ba25ull,
0275                                     0xc052c08c36880273ull,
0276                                     0xc03028545b6b807aull,
0277                                     0xbfec007fa1f72594ull>(z)
0278                     / detail::horner1<batch_type,
0279                                       0x4068519efbbd62ecull,
0280                                       0x407e563f13b049eaull,
0281                                       0x407b0e18d2e2be3bull,
0282                                       0x4064a0dd43b8fa25ull,
0283                                       0x4038dbc45b14603cull>(z);
0284                 z = fma(xx, z, xx);
0285                 z = select(flag2, z + constants::pio_4lo<batch_type>(), z);
0286                 z = z + select(flag1, batch_type(0.), constants::pio_2lo<batch_type>());
0287                 return yy + z;
0288             }
0289         }
0290         template <class A, class T>
0291         XSIMD_INLINE batch<T, A> atan(batch<T, A> const& self, requires_arch<common>) noexcept
0292         {
0293             using batch_type = batch<T, A>;
0294             const batch_type absa = abs(self);
0295             const batch_type x = detail::kernel_atan(absa, batch_type(1.) / absa);
0296             return x ^ bitofsign(self);
0297         }
0298         template <class A, class T>
0299         XSIMD_INLINE batch<std::complex<T>, A> atan(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
0300         {
0301             using batch_type = batch<std::complex<T>, A>;
0302             using real_batch = typename batch_type::real_batch;
0303             real_batch x = z.real();
0304             real_batch y = z.imag();
0305             real_batch x2 = x * x;
0306             real_batch one(1.);
0307             real_batch a = one - x2 - (y * y);
0308             real_batch w = 0.5 * atan2(2. * x, a);
0309             real_batch num = y + one;
0310             num = x2 + num * num;
0311             real_batch den = y - one;
0312             den = x2 + den * den;
0313 #ifdef __FAST_MATH__
0314             return batch_type(w, 0.25 * log(num / den));
0315 #else
0316             return select((x == real_batch(0.)) && (y == real_batch(1.)),
0317                           batch_type(real_batch(0.), constants::infinity<real_batch>()),
0318                           batch_type(w, 0.25 * log(num / den)));
0319 #endif
0320         }
0321
0322         // atanh
0323         /* origin: boost/simd/arch/common/simd/function/acosh.hpp */
0324         /*
0325          * ====================================================
0326          * copyright 2016 NumScale SAS
0327          *
0328          * Distributed under the Boost Software License, Version 1.0.
0329          * (See copy at http://boost.org/LICENSE_1_0.txt)
0330          * ====================================================
0331          */
0332         template <class A, class T>
0333         XSIMD_INLINE batch<T, A> atanh(batch<T, A> const& self, requires_arch<common>) noexcept
0334         {
0335             using batch_type = batch<T, A>;
0336             batch_type x = abs(self);
0337             batch_type t = x + x;
0338             batch_type z = batch_type(1.) - x;
0339             auto test = x < batch_type(0.5);
0340             batch_type tmp = select(test, x, t) / z;
0341             return bitofsign(self) ^ (batch_type(0.5) * log1p(select(test, fma(t, tmp, t), tmp)));
0342         }
0343         template <class A, class T>
0344         XSIMD_INLINE batch<std::complex<T>, A> atanh(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
0345         {
0346             using batch_type = batch<std::complex<T>, A>;
0347             batch_type w = atan(batch_type(-z.imag(), z.real()));
0348             w = batch_type(w.imag(), -w.real());
0349             return w;
0350         }
0351
0352         // atan2
0353         template <class A, class T>
0354         XSIMD_INLINE batch<T, A> atan2(batch<T, A> const& self, batch<T, A> const& other, requires_arch<common>) noexcept
0355         {
0356             using batch_type = batch<T, A>;
0357 #ifdef __FAST_MATH__
0358             const batch_type q = abs(self / other);
0359             const batch_type q_p = abs(other / self);
0360 #else
0361             const batch_type q = abs(self / other);
0362             const batch_type q_p = 1. / q;
0363 #endif
0364             const batch_type z = detail::kernel_atan(q, q_p);
0365             return select(other > batch_type(0.), z, constants::pi<batch_type>() - z) * signnz(self);
0366         }
0367
0368         // cos
0369         namespace detail
0370         {
0371             template <class T, class A>
0372             XSIMD_INLINE batch<T, A> quadrant(const batch<T, A>& x) noexcept
0373             {
0374                 return x & batch<T, A>(3);
0375             }
0376
0377             template <class A>
0378             XSIMD_INLINE batch<float, A> quadrant(const batch<float, A>& x) noexcept
0379             {
0380                 return to_float(quadrant(to_int(x)));
0381             }
0382
0383             template <class A>
0384             XSIMD_INLINE batch<double, A> quadrant(const batch<double, A>& x) noexcept
0385             {
0386                 using batch_type = batch<double, A>;
0387                 batch_type a = x * batch_type(0.25);
0388                 return (a - floor(a)) * batch_type(4.);
0389             }
0390             /* origin: boost/simd/arch/common/detail/simd/f_trig_evaluation.hpp */
0391             /*
0392              * ====================================================
0393              * copyright 2016 NumScale SAS
0394              *
0395              * Distributed under the Boost Software License, Version 1.0.
0396              * (See copy at http://boost.org/LICENSE_1_0.txt)
0397              * ====================================================
0398              */
0399
0400             template <class A>
0401             XSIMD_INLINE batch<float, A> cos_eval(const batch<float, A>& z) noexcept
0402             {
0403                 using batch_type = batch<float, A>;
0404                 batch_type y = detail::horner<batch_type,
0405                                               0x3d2aaaa5,
0406                                               0xbab60619,
0407                                               0x37ccf5ce>(z);
0408                 return batch_type(1.) + fma(z, batch_type(-0.5), y * z * z);
0409             }
0410
0411             template <class A>
0412             XSIMD_INLINE batch<float, A> sin_eval(const batch<float, A>& z, const batch<float, A>& x) noexcept
0413             {
0414                 using batch_type = batch<float, A>;
0415                 batch_type y = detail::horner<batch_type,
0416                                               0xbe2aaaa2,
0417                                               0x3c08839d,
0418                                               0xb94ca1f9>(z);
0419                 return fma(y * z, x, x);
0420             }
0421
0422             template <class A>
0423             static XSIMD_INLINE batch<float, A> base_tancot_eval(const batch<float, A>& z) noexcept
0424             {
0425                 using batch_type = batch<float, A>;
0426                 batch_type zz = z * z;
0427                 batch_type y = detail::horner<batch_type,
0428                                               0x3eaaaa6f,
0429                                               0x3e0896dd,
0430                                               0x3d5ac5c9,
0431                                               0x3cc821b5,
0432                                               0x3b4c779c,
0433                                               0x3c19c53b>(zz);
0434                 return fma(y, zz * z, z);
0435             }
0436
0437             template <class A, class BB>
0438             static XSIMD_INLINE batch<float, A> tan_eval(const batch<float, A>& z, const BB& test) noexcept
0439             {
0440                 using batch_type = batch<float, A>;
0441                 batch_type y = base_tancot_eval(z);
0442                 return select(test, y, -batch_type(1.) / y);
0443             }
0444
0445             template <class A, class BB>
0446             static XSIMD_INLINE batch<float, A> cot_eval(const batch<float, A>& z, const BB& test) noexcept
0447             {
0448                 using batch_type = batch<float, A>;
0449                 batch_type y = base_tancot_eval(z);
0450                 return select(test, batch_type(1.) / y, -y);
0451             }
0452
0453             /* origin: boost/simd/arch/common/detail/simd/d_trig_evaluation.hpp */
0454             /*
0455              * ====================================================
0456              * copyright 2016 NumScale SAS
0457              *
0458              * Distributed under the Boost Software License, Version 1.0.
0459              * (See copy at http://boost.org/LICENSE_1_0.txt)
0460              * ====================================================
0461              */
0462             template <class A>
0463             static XSIMD_INLINE batch<double, A> cos_eval(const batch<double, A>& z) noexcept
0464             {
0465                 using batch_type = batch<double, A>;
0466                 batch_type y = detail::horner<batch_type,
0467                                               0x3fe0000000000000ull,
0468                                               0xbfa5555555555551ull,
0469                                               0x3f56c16c16c15d47ull,
0470                                               0xbefa01a019ddbcd9ull,
0471                                               0x3e927e4f8e06d9a5ull,
0472                                               0xbe21eea7c1e514d4ull,
0473                                               0x3da8ff831ad9b219ull>(z);
0474                 return batch_type(1.) - y * z;
0475             }
0476
0477             template <class A>
0478             static XSIMD_INLINE batch<double, A> sin_eval(const batch<double, A>& z, const batch<double, A>& x) noexcept
0479             {
0480                 using batch_type = batch<double, A>;
0481                 batch_type y = detail::horner<batch_type,
0482                                               0xbfc5555555555548ull,
0483                                               0x3f8111111110f7d0ull,
0484                                               0xbf2a01a019bfdf03ull,
0485                                               0x3ec71de3567d4896ull,
0486                                               0xbe5ae5e5a9291691ull,
0487                                               0x3de5d8fd1fcf0ec1ull>(z);
0488                 return fma(y * z, x, x);
0489             }
0490
0491             template <class A>
0492             static XSIMD_INLINE batch<double, A> base_tancot_eval(const batch<double, A>& z) noexcept
0493             {
0494                 using batch_type = batch<double, A>;
0495                 batch_type zz = z * z;
0496                 batch_type num = detail::horner<batch_type,
0497                                                 0xc1711fead3299176ull,
0498                                                 0x413199eca5fc9dddull,
0499                                                 0xc0c992d8d24f3f38ull>(zz);
0500                 batch_type den = detail::horner1<batch_type,
0501                                                  0xc189afe03cbe5a31ull,
0502                                                  0x4177d98fc2ead8efull,
0503                                                  0xc13427bc582abc96ull,
0504                                                  0x40cab8a5eeb36572ull>(zz);
0505                 return fma(z, (zz * (num / den)), z);
0506             }
0507
0508             template <class A, class BB>
0509             static XSIMD_INLINE batch<double, A> tan_eval(const batch<double, A>& z, const BB& test) noexcept
0510             {
0511                 using batch_type = batch<double, A>;
0512                 batch_type y = base_tancot_eval(z);
0513                 return select(test, y, -batch_type(1.) / y);
0514             }
0515
0516             template <class A, class BB>
0517             static XSIMD_INLINE batch<double, A> cot_eval(const batch<double, A>& z, const BB& test) noexcept
0518             {
0519                 using batch_type = batch<double, A>;
0520                 batch_type y = base_tancot_eval(z);
0521                 return select(test, batch_type(1.) / y, -y);
0522             }
0523             /* origin: boost/simd/arch/common/detail/simd/trig_reduction.hpp */
0524             /*
0525              * ====================================================
0526              * copyright 2016 NumScale SAS
0527              *
0528              * Distributed under the Boost Software License, Version 1.0.
0529              * (See copy at http://boost.org/LICENSE_1_0.txt)
0530              * ====================================================
0531              */
0532
0533             struct trigo_radian_tag
0534             {
0535             };
0536             struct trigo_pi_tag
0537             {
0538             };
0539
0540             template <class B, class Tag = trigo_radian_tag>
0541             struct trigo_reducer
0542             {
0543                 static XSIMD_INLINE B reduce(const B& x, B& xr) noexcept
0544                 {
0545                     if (all(x <= constants::pio4<B>()))
0546                     {
0547                         xr = x;
0548                         return B(0.);
0549                     }
0550                     else if (all(x <= constants::pio2<B>()))
0551                     {
0552                         auto test = x > constants::pio4<B>();
0553                         xr = x - constants::pio2_1<B>();
0554                         detail::reassociation_barrier(xr, "ordered pio2 subtraction");
0555                         xr -= constants::pio2_2<B>();
0556                         detail::reassociation_barrier(xr, "ordered pio2 subtraction");
0557                         xr -= constants::pio2_3<B>();
0558                         detail::reassociation_barrier(xr, "ordered pio2 subtraction");
0559                         xr = select(test, xr, x);
0560                         return select(test, B(1.), B(0.));
0561                     }
0562                     else if (all(x <= constants::twentypi<B>()))
0563                     {
0564                         B xi = nearbyint(x * constants::twoopi<B>());
0565                         detail::reassociation_barrier(xi, "preserve quadrant selection");
0566                         xr = fnma(xi, constants::pio2_1<B>(), x);
0567                         detail::reassociation_barrier(xr, "compensated range reduction");
0568                         xr -= xi * constants::pio2_2<B>();
0569                         detail::reassociation_barrier(xr, "compensated range reduction");
0570                         xr -= xi * constants::pio2_3<B>();
0571                         detail::reassociation_barrier(xr, "compensated range reduction");
0572                         return quadrant(xi);
0573                     }
0574                     else if (all(x <= constants::mediumpi<B>()))
0575                     {
0576                         B fn = nearbyint(x * constants::twoopi<B>());
0577                         detail::reassociation_barrier(fn, "multi-term range reduction");
0578                         B r = x - fn * constants::pio2_1<B>();
0579                         detail::reassociation_barrier(r, "multi-term range reduction");
0580                         B w = fn * constants::pio2_1t<B>();
0581                         B t = r;
0582                         w = fn * constants::pio2_2<B>();
0583                         r = t - w;
0584                         detail::reassociation_barrier(r, "multi-term range reduction");
0585                         w = fn * constants::pio2_2t<B>() - ((t - r) - w);
0586                         t = r;
0587                         w = fn * constants::pio2_3<B>();
0588                         r = t - w;
0589                         detail::reassociation_barrier(r, "multi-term range reduction");
0590                         w = fn * constants::pio2_3t<B>() - ((t - r) - w);
0591                         xr = r - w;
0592                         detail::reassociation_barrier(xr, "multi-term range reduction");
0593                         return quadrant(fn);
0594                     }
0595                     else
0596                     {
0597                         static constexpr std::size_t size = B::size;
0598                         using value_type = typename B::value_type;
0599                         alignas(B::arch_type::alignment()) std::array<value_type, size> tmp;
0600                         alignas(B::arch_type::alignment()) std::array<value_type, size> txr;
0601                         alignas(B::arch_type::alignment()) std::array<value_type, size> args;
0602                         x.store_aligned(args.data());
0603
0604                         for (std::size_t i = 0; i < size; ++i)
0605                         {
0606                             double arg = args[i];
0607 #ifndef __FAST_MATH__
0608                             if (arg == std::numeric_limits<value_type>::infinity())
0609                             {
0610                                 tmp[i] = 0.;
0611                                 txr[i] = std::numeric_limits<value_type>::quiet_NaN();
0612                             }
0613                             else
0614 #endif
0615                             {
0616                                 double y[2];
0617                                 std::int32_t n = ::xsimd::detail::__ieee754_rem_pio2(arg, y);
0618                                 tmp[i] = value_type(n & 3);
0619                                 txr[i] = value_type(y[0]);
0620                             }
0621                         }
0622                         xr = B::load_aligned(&txr[0]);
0623                         B res = B::load_aligned(&tmp[0]);
0624                         return res;
0625                     }
0626                 }
0627             };
0628
0629             template <class B>
0630             struct trigo_reducer<B, trigo_pi_tag>
0631             {
0632                 static XSIMD_INLINE B reduce(const B& x, B& xr) noexcept
0633                 {
0634                     B xi = nearbyint(x * B(2.));
0635                     B x2 = x - xi * B(0.5);
0636                     xr = x2 * constants::pi<B>();
0637                     return quadrant(xi);
0638                 }
0639             };
0640
0641         }
0642         template <class A, class T>
0643         XSIMD_INLINE batch<T, A> cos(batch<T, A> const& self, requires_arch<common>) noexcept
0644         {
0645             using batch_type = batch<T, A>;
0646             const batch_type x = abs(self);
0647             batch_type xr = constants::nan<batch_type>();
0648             const batch_type n = detail::trigo_reducer<batch_type>::reduce(x, xr);
0649             auto tmp = select(n >= batch_type(2.), batch_type(1.), batch_type(0.));
0650             auto swap_bit = fma(batch_type(-2.), tmp, n);
0651             auto sign_bit = select((swap_bit ^ tmp) != batch_type(0.), constants::signmask<batch_type>(), batch_type(0.));
0652             const batch_type z = xr * xr;
0653             const batch_type se = detail::sin_eval(z, xr);
0654             const batch_type ce = detail::cos_eval(z);
0655             const batch_type z1 = select(swap_bit != batch_type(0.), se, ce);
0656             return z1 ^ sign_bit;
0657         }
0658
0659         template <class A, class T>
0660         XSIMD_INLINE batch<std::complex<T>, A> cos(batch<std::complex<T>, A> const& z, requires_arch<common>) noexcept
0661         {
0662             return { cos(z.real()) * cosh(z.imag()), -sin(z.real()) * sinh(z.imag()) };
0663         }
0664
0665         // cosh
0666
0667         /* origin: boost/simd/arch/common/simd/function/cosh.hpp */
0668         /*
0669          * ====================================================
0670          * copyright 2016 NumScale SAS
0671          *
0672          * Distributed under the Boost Software License, Version 1.0.
0673          * (See copy at http://boost.org/LICENSE_1_0.txt)
0674          * ====================================================
0675          */
0676
0677         template <class A, class T>
0678         XSIMD_INLINE batch<T, A> cosh(batch<T, A> const& self, requires_arch<common>) noexcept
0679         {
0680             using batch_type = batch<T, A>;
0681             batch_type x = abs(self);
0682             auto test1 = x > (constants::maxlog<batch_type>() - constants::log_2<batch_type>());
0683             batch_type fac = select(test1, batch_type(0.5), batch_type(1.));
0684             batch_type tmp = exp(x * fac);
0685             batch_type tmp1 = batch_type(0.5) * tmp;
0686             return select(test1, tmp1 * tmp, detail::average(tmp, batch_type(1.) / tmp));
0687         }
0688         template <class A, class T>
0689         XSIMD_INLINE batch<std::complex<T>, A> cosh(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
0690         {
0691             auto x = z.real();
0692             auto y = z.imag();
0693             return { cosh(x) * cos(y), sinh(x) * sin(y) };
0694         }
0695
0696         // sin
0697         namespace detail
0698         {
0699             template <class A, class T, class Tag = trigo_radian_tag>
0700             XSIMD_INLINE batch<T, A> sin(batch<T, A> const& self, Tag = Tag()) noexcept
0701             {
0702                 using batch_type = batch<T, A>;
0703                 const batch_type x = abs(self);
0704                 batch_type xr = constants::nan<batch_type>();
0705                 const batch_type n = detail::trigo_reducer<batch_type, Tag>::reduce(x, xr);
0706                 auto tmp = select(n >= batch_type(2.), batch_type(1.), batch_type(0.));
0707                 auto swap_bit = fma(batch_type(-2.), tmp, n);
0708                 auto sign_bit = bitofsign(self) ^ select(tmp != batch_type(0.), constants::signmask<batch_type>(), batch_type(0.));
0709                 const batch_type z = xr * xr;
0710                 const batch_type se = detail::sin_eval(z, xr);
0711                 const batch_type ce = detail::cos_eval(z);
0712                 const batch_type z1 = select(swap_bit == batch_type(0.), se, ce);
0713                 return z1 ^ sign_bit;
0714             }
0715         }
0716
0717         template <class A, class T>
0718         XSIMD_INLINE batch<T, A> sin(batch<T, A> const& self, requires_arch<common>) noexcept
0719         {
0720             return detail::sin(self);
0721         }
0722
0723         template <class A, class T>
0724         XSIMD_INLINE batch<std::complex<T>, A> sin(batch<std::complex<T>, A> const& z, requires_arch<common>) noexcept
0725         {
0726             return { sin(z.real()) * cosh(z.imag()), cos(z.real()) * sinh(z.imag()) };
0727         }
0728
0729         // sincos
0730         template <class A, class T>
0731         XSIMD_INLINE std::pair<batch<T, A>, batch<T, A>> sincos(batch<T, A> const& self, requires_arch<common>) noexcept
0732         {
0733             using batch_type = batch<T, A>;
0734             const batch_type x = abs(self);
0735             batch_type xr = constants::nan<batch_type>();
0736             const batch_type n = detail::trigo_reducer<batch_type>::reduce(x, xr);
0737             auto tmp = select(n >= batch_type(2.), batch_type(1.), batch_type(0.));
0738             auto swap_bit = fma(batch_type(-2.), tmp, n);
0739             const batch_type z = xr * xr;
0740             const batch_type se = detail::sin_eval(z, xr);
0741             const batch_type ce = detail::cos_eval(z);
0742             auto sin_sign_bit = bitofsign(self) ^ select(tmp != batch_type(0.), constants::signmask<batch_type>(), batch_type(0.));
0743             const batch_type sin_z1 = select(swap_bit == batch_type(0.), se, ce);
0744             auto cos_sign_bit = select((swap_bit ^ tmp) != batch_type(0.), constants::signmask<batch_type>(), batch_type(0.));
0745             const batch_type cos_z1 = select(swap_bit != batch_type(0.), se, ce);
0746             return std::make_pair(sin_z1 ^ sin_sign_bit, cos_z1 ^ cos_sign_bit);
0747         }
0748
0749         template <class A, class T>
0750         XSIMD_INLINE std::pair<batch<std::complex<T>, A>, batch<std::complex<T>, A>>
0751         sincos(batch<std::complex<T>, A> const& z, requires_arch<common>) noexcept
0752         {
0753             using batch_type = batch<std::complex<T>, A>;
0754             using real_batch = typename batch_type::real_batch;
0755             real_batch rcos = cos(z.real());
0756             real_batch rsin = sin(z.real());
0757             real_batch icosh = cosh(z.imag());
0758             real_batch isinh = sinh(z.imag());
0759             return std::make_pair(batch_type(rsin * icosh, rcos * isinh), batch_type(rcos * icosh, -rsin * isinh));
0760         }
0761
0762         // sinh
0763         namespace detail
0764         {
0765             /* origin: boost/simd/arch/common/detail/common/sinh_kernel.hpp */
0766             /*
0767              * ====================================================
0768              * copyright 2016 NumScale SAS
0769              *
0770              * Distributed under the Boost Software License, Version 1.0.
0771              * (See copy at http://boost.org/LICENSE_1_0.txt)
0772              * ====================================================
0773              */
0774             template <class A>
0775             XSIMD_INLINE batch<float, A> sinh_kernel(batch<float, A> const& self) noexcept
0776             {
0777                 using batch_type = batch<float, A>;
0778                 batch_type sqr_self = self * self;
0779                 return detail::horner<batch_type,
0780                                       0x3f800000, // 1.0f
0781                                       0x3e2aaacc, // 1.66667160211E-1f
0782                                       0x3c087bbe, // 8.33028376239E-3f
0783                                       0x39559e2f // 2.03721912945E-4f
0784                                       >(sqr_self)
0785                     * self;
0786             }
0787
0788             template <class A>
0789             XSIMD_INLINE batch<double, A> sinh_kernel(batch<double, A> const& self) noexcept
0790             {
0791                 using batch_type = batch<double, A>;
0792                 batch_type sqrself = self * self;
0793                 return fma(self, (detail::horner<batch_type,
0794                                                  0xc115782bdbf6ab05ull, //  -3.51754964808151394800E5
0795                                                  0xc0c694b8c71d6182ull, //  -1.15614435765005216044E4,
0796                                                  0xc064773a398ff4feull, //  -1.63725857525983828727E2,
0797                                                  0xbfe9435fe8bb3cd6ull //  -7.89474443963537015605E-1
0798                                                  >(sqrself)
0799                                   / detail::horner1<batch_type,
0800                                                     0xc1401a20e4f90044ull, //  -2.11052978884890840399E6
0801                                                     0x40e1a7ba7ed72245ull, //   3.61578279834431989373E4,
0802                                                     0xc0715b6096e96484ull //  -2.77711081420602794433E2,
0803                                                     >(sqrself))
0804                                * sqrself,
0805                            self);
0806             }
0807         }
0808         /* origin: boost/simd/arch/common/simd/function/sinh.hpp */
0809         /*
0810          * ====================================================
0811          * copyright 2016 NumScale SAS
0812          *
0813          * Distributed under the Boost Software License, Version 1.0.
0814          * (See copy at http://boost.org/LICENSE_1_0.txt)
0815          * ====================================================
0816          */
0817         template <class A, class T>
0818         XSIMD_INLINE batch<T, A> sinh(batch<T, A> const& a, requires_arch<common>) noexcept
0819         {
0820             using batch_type = batch<T, A>;
0821             batch_type half(0.5);
0822             batch_type x = abs(a);
0823             auto lt1 = x < batch_type(1.);
0824             batch_type bts = bitofsign(a);
0825             batch_type z(0.);
0826             if (any(lt1))
0827             {
0828                 z = detail::sinh_kernel(x);
0829                 if (all(lt1))
0830                     return z ^ bts;
0831             }
0832             auto test1 = x > (constants::maxlog<batch_type>() - constants::log_2<batch_type>());
0833             batch_type fac = select(test1, half, batch_type(1.));
0834             batch_type tmp = exp(x * fac);
0835             batch_type tmp1 = half * tmp;
0836             batch_type r = select(test1, tmp1 * tmp, tmp1 - half / tmp);
0837             return select(lt1, z, r) ^ bts;
0838         }
0839         template <class A, class T>
0840         XSIMD_INLINE batch<std::complex<T>, A> sinh(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
0841         {
0842             auto x = z.real();
0843             auto y = z.imag();
0844             return { sinh(x) * cos(y), cosh(x) * sin(y) };
0845         }
0846
0847         // tan
0848         template <class A, class T>
0849         XSIMD_INLINE batch<T, A> tan(batch<T, A> const& self, requires_arch<common>) noexcept
0850         {
0851             using batch_type = batch<T, A>;
0852             const batch_type x = abs(self);
0853             batch_type xr = constants::nan<batch_type>();
0854             const batch_type n = detail::trigo_reducer<batch_type>::reduce(x, xr);
0855             auto tmp = select(n >= batch_type(2.), batch_type(1.), batch_type(0.));
0856             auto swap_bit = fma(batch_type(-2.), tmp, n);
0857             auto test = (swap_bit == batch_type(0.));
0858             const batch_type y = detail::tan_eval(xr, test);
0859             return y ^ bitofsign(self);
0860         }
0861         template <class A, class T>
0862         XSIMD_INLINE batch<std::complex<T>, A> tan(batch<std::complex<T>, A> const& z, requires_arch<common>) noexcept
0863         {
0864             using batch_type = batch<std::complex<T>, A>;
0865             using real_batch = typename batch_type::real_batch;
0866             real_batch d = cos(2 * z.real()) + cosh(2 * z.imag());
0867             real_batch wreal = sin(2 * z.real()) / d;
0868             real_batch wimag = sinh(2 * z.imag());
0869 #ifdef __FAST_MATH__
0870             return batch_type(wreal, real_batch(1.)), batch_type(wreal, wimag / d);
0871 #else
0872             batch_type winf(constants::infinity<real_batch>(), constants::infinity<real_batch>());
0873             batch_type wres = select(isinf(wimag), batch_type(wreal, real_batch(1.)), batch_type(wreal, wimag / d));
0874             return select(d == real_batch(0.), winf, wres);
0875 #endif
0876         }
0877
0878         // tanh
0879         namespace detail
0880         {
0881             /* origin: boost/simd/arch/common/detail/common/tanh_kernel.hpp */
0882             /*
0883              * ====================================================
0884              * copyright 2016 NumScale SAS
0885              *
0886              * Distributed under the Boost Software License, Version 1.0.
0887              * (See copy at http://boost.org/LICENSE_1_0.txt)
0888              * ====================================================
0889              */
0890             template <class B>
0891             struct tanh_kernel;
0892
0893             template <class A>
0894             struct tanh_kernel<batch<float, A>>
0895             {
0896                 using batch_type = batch<float, A>;
0897                 static XSIMD_INLINE batch_type tanh(const batch_type& x) noexcept
0898                 {
0899                     batch_type sqrx = x * x;
0900                     return fma(detail::horner<batch_type,
0901                                               0xbeaaaa99, //    -3.33332819422E-1F
0902                                               0x3e088393, //    +1.33314422036E-1F
0903                                               0xbd5c1e2d, //    -5.37397155531E-2F
0904                                               0x3ca9134e, //    +2.06390887954E-2F
0905                                               0xbbbaf0ea //    -5.70498872745E-3F
0906                                               >(sqrx)
0907                                    * sqrx,
0908                                x, x);
0909                 }
0910
0911                 static XSIMD_INLINE batch_type cotanh(const batch_type& x) noexcept
0912                 {
0913                     return batch_type(1.) / tanh(x);
0914                 }
0915             };
0916
0917             template <class A>
0918             struct tanh_kernel<batch<double, A>>
0919             {
0920                 using batch_type = batch<double, A>;
0921                 static XSIMD_INLINE batch_type tanh(const batch_type& x) noexcept
0922                 {
0923                     batch_type sqrx = x * x;
0924                     return fma(sqrx * p(sqrx) / q(sqrx), x, x);
0925                 }
0926
0927                 static XSIMD_INLINE batch_type cotanh(const batch_type& x) noexcept
0928                 {
0929                     batch_type sqrx = x * x;
0930                     batch_type qval = q(sqrx);
0931                     return qval / (x * fma(p(sqrx), sqrx, qval));
0932                 }
0933
0934                 static XSIMD_INLINE batch_type p(const batch_type& x) noexcept
0935                 {
0936                     return detail::horner<batch_type,
0937                                           0xc0993ac030580563, // -1.61468768441708447952E3
0938                                           0xc058d26a0e26682d, // -9.92877231001918586564E1,
0939                                           0xbfeedc5baafd6f4b // -9.64399179425052238628E-1
0940                                           >(x);
0941                 }
0942
0943                 static XSIMD_INLINE batch_type q(const batch_type& x) noexcept
0944                 {
0945                     return detail::horner1<batch_type,
0946                                            0x40b2ec102442040c, //  4.84406305325125486048E3
0947                                            0x40a176fa0e5535fa, //  2.23548839060100448583E3,
0948                                            0x405c33f28a581B86 //  1.12811678491632931402E2,
0949                                            >(x);
0950                 }
0951             };
0952
0953         }
0954         /* origin: boost/simd/arch/common/simd/function/tanh.hpp */
0955         /*
0956          * ====================================================
0957          * copyright 2016 NumScale SAS
0958          *
0959          * Distributed under the Boost Software License, Version 1.0.
0960          * (See copy at http://boost.org/LICENSE_1_0.txt)
0961          * ====================================================
0962          */
0963         template <class A, class T>
0964         XSIMD_INLINE batch<T, A> tanh(batch<T, A> const& self, requires_arch<common>) noexcept
0965         {
0966             using batch_type = batch<T, A>;
0967             batch_type one(1.);
0968             batch_type x = abs(self);
0969             auto test = x < (batch_type(5.) / batch_type(8.));
0970             batch_type bts = bitofsign(self);
0971             batch_type z = one;
0972             if (any(test))
0973             {
0974                 z = detail::tanh_kernel<batch_type>::tanh(x);
0975                 if (all(test))
0976                     return z ^ bts;
0977             }
0978             batch_type r = fma(batch_type(-2.), one / (one + exp(x + x)), one);
0979             return select(test, z, r) ^ bts;
0980         }
0981         template <class A, class T>
0982         XSIMD_INLINE batch<std::complex<T>, A> tanh(const batch<std::complex<T>, A>& z, requires_arch<common>) noexcept
0983         {
0984             using real_batch = typename batch<std::complex<T>, A>::real_batch;
0985             auto x = z.real();
0986             auto y = z.imag();
0987             real_batch two(2);
0988             auto d = cosh(two * x) + cos(two * y);
0989             return { sinh(two * x) / d, sin(two * y) / d };
0990         }
0991
0992     }
0993
0994 }
0995
0996 #endif