BesselFunctionsImpl.h - cart-elc - Source code for CART-ELC

BesselFunctionsImpl.h (69632B)
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2015 Eugene Brevdo <ebrevdo@gmail.com>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_BESSEL_FUNCTIONS_H
     11 #define EIGEN_BESSEL_FUNCTIONS_H
     12 
     13 namespace Eigen {
     14 namespace internal {
     15 
     16 //  Parts of this code are based on the Cephes Math Library.
     17 //
     18 //  Cephes Math Library Release 2.8:  June, 2000
     19 //  Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
     20 //
     21 //  Permission has been kindly provided by the original author
     22 //  to incorporate the Cephes software into the Eigen codebase:
     23 //
     24 //    From: Stephen Moshier
     25 //    To: Eugene Brevdo
     26 //    Subject: Re: Permission to wrap several cephes functions in Eigen
     27 //
     28 //    Hello Eugene,
     29 //
     30 //    Thank you for writing.
     31 //
     32 //    If your licensing is similar to BSD, the formal way that has been
     33 //    handled is simply to add a statement to the effect that you are incorporating
     34 //    the Cephes software by permission of the author.
     35 //
     36 //    Good luck with your project,
     37 //    Steve
     38 
     39 
     40 /****************************************************************************
     41  * Implementation of Bessel function, based on Cephes                       *
     42  ****************************************************************************/
     43 
     44 template <typename Scalar>
     45 struct bessel_i0e_retval {
     46   typedef Scalar type;
     47 };
     48 
     49 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
     50 struct generic_i0e {
     51   EIGEN_DEVICE_FUNC
     52   static EIGEN_STRONG_INLINE T run(const T&) {
     53     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
     54                         THIS_TYPE_IS_NOT_SUPPORTED);
     55     return ScalarType(0);
     56   }
     57 };
     58 
     59 template <typename T>
     60 struct generic_i0e<T, float> {
     61   EIGEN_DEVICE_FUNC
     62   static EIGEN_STRONG_INLINE T run(const T& x) {
     63     /*  i0ef.c
     64      *
     65      *  Modified Bessel function of order zero,
     66      *  exponentially scaled
     67      *
     68      *
     69      *
     70      * SYNOPSIS:
     71      *
     72      * float x, y, i0ef();
     73      *
     74      * y = i0ef( x );
     75      *
     76      *
     77      *
     78      * DESCRIPTION:
     79      *
     80      * Returns exponentially scaled modified Bessel function
     81      * of order zero of the argument.
     82      *
     83      * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
     84      *
     85      *
     86      *
     87      * ACCURACY:
     88      *
     89      *                      Relative error:
     90      * arithmetic   domain     # trials      peak         rms
     91      *    IEEE      0,30        100000      3.7e-7      7.0e-8
     92      * See i0f().
     93      *
     94      */
     95 
     96     const float A[] = {-1.30002500998624804212E-8f, 6.04699502254191894932E-8f,
     97                        -2.67079385394061173391E-7f, 1.11738753912010371815E-6f,
     98                        -4.41673835845875056359E-6f, 1.64484480707288970893E-5f,
     99                        -5.75419501008210370398E-5f, 1.88502885095841655729E-4f,
    100                        -5.76375574538582365885E-4f, 1.63947561694133579842E-3f,
    101                        -4.32430999505057594430E-3f, 1.05464603945949983183E-2f,
    102                        -2.37374148058994688156E-2f, 4.93052842396707084878E-2f,
    103                        -9.49010970480476444210E-2f, 1.71620901522208775349E-1f,
    104                        -3.04682672343198398683E-1f, 6.76795274409476084995E-1f};
    105 
    106     const float B[] = {3.39623202570838634515E-9f, 2.26666899049817806459E-8f,
    107                        2.04891858946906374183E-7f, 2.89137052083475648297E-6f,
    108                        6.88975834691682398426E-5f, 3.36911647825569408990E-3f,
    109                        8.04490411014108831608E-1f};
    110     T y = pabs(x);
    111     T y_le_eight = internal::pchebevl<T, 18>::run(
    112         pmadd(pset1<T>(0.5f), y, pset1<T>(-2.0f)), A);
    113     T y_gt_eight = pmul(
    114         internal::pchebevl<T, 7>::run(
    115             psub(pdiv(pset1<T>(32.0f), y), pset1<T>(2.0f)), B),
    116         prsqrt(y));
    117     // TODO: Perhaps instead check whether all packet elements are in
    118     // [-8, 8] and evaluate a branch based off of that. It's possible
    119     // in practice most elements are in this region.
    120     return pselect(pcmp_le(y, pset1<T>(8.0f)), y_le_eight, y_gt_eight);
    121   }
    122 };
    123 
    124 template <typename T>
    125 struct generic_i0e<T, double> {
    126   EIGEN_DEVICE_FUNC
    127   static EIGEN_STRONG_INLINE T run(const T& x) {
    128     /*  i0e.c
    129      *
    130      *  Modified Bessel function of order zero,
    131      *  exponentially scaled
    132      *
    133      *
    134      *
    135      * SYNOPSIS:
    136      *
    137      * double x, y, i0e();
    138      *
    139      * y = i0e( x );
    140      *
    141      *
    142      *
    143      * DESCRIPTION:
    144      *
    145      * Returns exponentially scaled modified Bessel function
    146      * of order zero of the argument.
    147      *
    148      * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
    149      *
    150      *
    151      *
    152      * ACCURACY:
    153      *
    154      *                      Relative error:
    155      * arithmetic   domain     # trials      peak         rms
    156      *    IEEE      0,30        30000       5.4e-16     1.2e-16
    157      * See i0().
    158      *
    159      */
    160 
    161     const double A[] = {-4.41534164647933937950E-18, 3.33079451882223809783E-17,
    162                         -2.43127984654795469359E-16, 1.71539128555513303061E-15,
    163                         -1.16853328779934516808E-14, 7.67618549860493561688E-14,
    164                         -4.85644678311192946090E-13, 2.95505266312963983461E-12,
    165                         -1.72682629144155570723E-11, 9.67580903537323691224E-11,
    166                         -5.18979560163526290666E-10, 2.65982372468238665035E-9,
    167                         -1.30002500998624804212E-8,  6.04699502254191894932E-8,
    168                         -2.67079385394061173391E-7,  1.11738753912010371815E-6,
    169                         -4.41673835845875056359E-6,  1.64484480707288970893E-5,
    170                         -5.75419501008210370398E-5,  1.88502885095841655729E-4,
    171                         -5.76375574538582365885E-4,  1.63947561694133579842E-3,
    172                         -4.32430999505057594430E-3,  1.05464603945949983183E-2,
    173                         -2.37374148058994688156E-2,  4.93052842396707084878E-2,
    174                         -9.49010970480476444210E-2,  1.71620901522208775349E-1,
    175                         -3.04682672343198398683E-1,  6.76795274409476084995E-1};
    176     const double B[] = {
    177         -7.23318048787475395456E-18, -4.83050448594418207126E-18,
    178         4.46562142029675999901E-17,  3.46122286769746109310E-17,
    179         -2.82762398051658348494E-16, -3.42548561967721913462E-16,
    180         1.77256013305652638360E-15,  3.81168066935262242075E-15,
    181         -9.55484669882830764870E-15, -4.15056934728722208663E-14,
    182         1.54008621752140982691E-14,  3.85277838274214270114E-13,
    183         7.18012445138366623367E-13,  -1.79417853150680611778E-12,
    184         -1.32158118404477131188E-11, -3.14991652796324136454E-11,
    185         1.18891471078464383424E-11,  4.94060238822496958910E-10,
    186         3.39623202570838634515E-9,   2.26666899049817806459E-8,
    187         2.04891858946906374183E-7,   2.89137052083475648297E-6,
    188         6.88975834691682398426E-5,   3.36911647825569408990E-3,
    189         8.04490411014108831608E-1};
    190     T y = pabs(x);
    191     T y_le_eight = internal::pchebevl<T, 30>::run(
    192         pmadd(pset1<T>(0.5), y, pset1<T>(-2.0)), A);
    193     T y_gt_eight = pmul(
    194         internal::pchebevl<T, 25>::run(
    195             psub(pdiv(pset1<T>(32.0), y), pset1<T>(2.0)), B),
    196         prsqrt(y));
    197     // TODO: Perhaps instead check whether all packet elements are in
    198     // [-8, 8] and evaluate a branch based off of that. It's possible
    199     // in practice most elements are in this region.
    200     return pselect(pcmp_le(y, pset1<T>(8.0)), y_le_eight, y_gt_eight);
    201   }
    202 };
    203 
    204 template <typename T>
    205 struct bessel_i0e_impl {
    206   EIGEN_DEVICE_FUNC
    207   static EIGEN_STRONG_INLINE T run(const T x) {
    208     return generic_i0e<T>::run(x);
    209   }
    210 };
    211 
    212 template <typename Scalar>
    213 struct bessel_i0_retval {
    214   typedef Scalar type;
    215 };
    216 
    217 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
    218 struct generic_i0 {
    219   EIGEN_DEVICE_FUNC
    220   static EIGEN_STRONG_INLINE T run(const T& x) {
    221     return pmul(
    222         pexp(pabs(x)),
    223         generic_i0e<T, ScalarType>::run(x));
    224   }
    225 };
    226 
    227 template <typename T>
    228 struct bessel_i0_impl {
    229   EIGEN_DEVICE_FUNC
    230   static EIGEN_STRONG_INLINE T run(const T x) {
    231     return generic_i0<T>::run(x);
    232   }
    233 };
    234 
    235 template <typename Scalar>
    236 struct bessel_i1e_retval {
    237   typedef Scalar type;
    238 };
    239 
    240 template <typename T, typename ScalarType = typename unpacket_traits<T>::type >
    241 struct generic_i1e {
    242   EIGEN_DEVICE_FUNC
    243   static EIGEN_STRONG_INLINE T run(const T&) {
    244     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
    245                         THIS_TYPE_IS_NOT_SUPPORTED);
    246     return ScalarType(0);
    247   }
    248 };
    249 
    250 template <typename T>
    251 struct generic_i1e<T, float> {
    252   EIGEN_DEVICE_FUNC
    253   static EIGEN_STRONG_INLINE T run(const T& x) {
    254     /* i1ef.c
    255      *
    256      *  Modified Bessel function of order one,
    257      *  exponentially scaled
    258      *
    259      *
    260      *
    261      * SYNOPSIS:
    262      *
    263      * float x, y, i1ef();
    264      *
    265      * y = i1ef( x );
    266      *
    267      *
    268      *
    269      * DESCRIPTION:
    270      *
    271      * Returns exponentially scaled modified Bessel function
    272      * of order one of the argument.
    273      *
    274      * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
    275      *
    276      *
    277      *
    278      * ACCURACY:
    279      *
    280      *                      Relative error:
    281      * arithmetic   domain     # trials      peak         rms
    282      *    IEEE      0, 30       30000       1.5e-6      1.5e-7
    283      * See i1().
    284      *
    285      */
    286     const float A[] = {9.38153738649577178388E-9f, -4.44505912879632808065E-8f,
    287                        2.00329475355213526229E-7f, -8.56872026469545474066E-7f,
    288                        3.47025130813767847674E-6f, -1.32731636560394358279E-5f,
    289                        4.78156510755005422638E-5f, -1.61760815825896745588E-4f,
    290                        5.12285956168575772895E-4f, -1.51357245063125314899E-3f,
    291                        4.15642294431288815669E-3f, -1.05640848946261981558E-2f,
    292                        2.47264490306265168283E-2f, -5.29459812080949914269E-2f,
    293                        1.02643658689847095384E-1f, -1.76416518357834055153E-1f,
    294                        2.52587186443633654823E-1f};
    295 
    296     const float B[] = {-3.83538038596423702205E-9f, -2.63146884688951950684E-8f,
    297                        -2.51223623787020892529E-7f, -3.88256480887769039346E-6f,
    298                        -1.10588938762623716291E-4f, -9.76109749136146840777E-3f,
    299                        7.78576235018280120474E-1f};
    300 
    301 
    302     T y = pabs(x);
    303     T y_le_eight = pmul(y, internal::pchebevl<T, 17>::run(
    304         pmadd(pset1<T>(0.5f), y, pset1<T>(-2.0f)), A));
    305     T y_gt_eight = pmul(
    306         internal::pchebevl<T, 7>::run(
    307             psub(pdiv(pset1<T>(32.0f), y),
    308                  pset1<T>(2.0f)), B),
    309         prsqrt(y));
    310     // TODO: Perhaps instead check whether all packet elements are in
    311     // [-8, 8] and evaluate a branch based off of that. It's possible
    312     // in practice most elements are in this region.
    313     y = pselect(pcmp_le(y, pset1<T>(8.0f)), y_le_eight, y_gt_eight);
    314     return pselect(pcmp_lt(x, pset1<T>(0.0f)), pnegate(y), y);
    315   }
    316 };
    317 
    318 template <typename T>
    319 struct generic_i1e<T, double> {
    320   EIGEN_DEVICE_FUNC
    321   static EIGEN_STRONG_INLINE T run(const T& x) {
    322     /*  i1e.c
    323      *
    324      *  Modified Bessel function of order one,
    325      *  exponentially scaled
    326      *
    327      *
    328      *
    329      * SYNOPSIS:
    330      *
    331      * double x, y, i1e();
    332      *
    333      * y = i1e( x );
    334      *
    335      *
    336      *
    337      * DESCRIPTION:
    338      *
    339      * Returns exponentially scaled modified Bessel function
    340      * of order one of the argument.
    341      *
    342      * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
    343      *
    344      *
    345      *
    346      * ACCURACY:
    347      *
    348      *                      Relative error:
    349      * arithmetic   domain     # trials      peak         rms
    350      *    IEEE      0, 30       30000       2.0e-15     2.0e-16
    351      * See i1().
    352      *
    353      */
    354     const double A[] = {2.77791411276104639959E-18, -2.11142121435816608115E-17,
    355                         1.55363195773620046921E-16, -1.10559694773538630805E-15,
    356                         7.60068429473540693410E-15, -5.04218550472791168711E-14,
    357                         3.22379336594557470981E-13, -1.98397439776494371520E-12,
    358                         1.17361862988909016308E-11, -6.66348972350202774223E-11,
    359                         3.62559028155211703701E-10, -1.88724975172282928790E-9,
    360                         9.38153738649577178388E-9,  -4.44505912879632808065E-8,
    361                         2.00329475355213526229E-7,  -8.56872026469545474066E-7,
    362                         3.47025130813767847674E-6,  -1.32731636560394358279E-5,
    363                         4.78156510755005422638E-5,  -1.61760815825896745588E-4,
    364                         5.12285956168575772895E-4,  -1.51357245063125314899E-3,
    365                         4.15642294431288815669E-3,  -1.05640848946261981558E-2,
    366                         2.47264490306265168283E-2,  -5.29459812080949914269E-2,
    367                         1.02643658689847095384E-1,  -1.76416518357834055153E-1,
    368                         2.52587186443633654823E-1};
    369     const double B[] = {
    370         7.51729631084210481353E-18,  4.41434832307170791151E-18,
    371         -4.65030536848935832153E-17, -3.20952592199342395980E-17,
    372         2.96262899764595013876E-16,  3.30820231092092828324E-16,
    373         -1.88035477551078244854E-15, -3.81440307243700780478E-15,
    374         1.04202769841288027642E-14,  4.27244001671195135429E-14,
    375         -2.10154184277266431302E-14, -4.08355111109219731823E-13,
    376         -7.19855177624590851209E-13, 2.03562854414708950722E-12,
    377         1.41258074366137813316E-11,  3.25260358301548823856E-11,
    378         -1.89749581235054123450E-11, -5.58974346219658380687E-10,
    379         -3.83538038596423702205E-9,  -2.63146884688951950684E-8,
    380         -2.51223623787020892529E-7,  -3.88256480887769039346E-6,
    381         -1.10588938762623716291E-4,  -9.76109749136146840777E-3,
    382         7.78576235018280120474E-1};
    383     T y = pabs(x);
    384     T y_le_eight = pmul(y, internal::pchebevl<T, 29>::run(
    385         pmadd(pset1<T>(0.5), y, pset1<T>(-2.0)), A));
    386     T y_gt_eight = pmul(
    387         internal::pchebevl<T, 25>::run(
    388             psub(pdiv(pset1<T>(32.0), y),
    389                  pset1<T>(2.0)), B),
    390         prsqrt(y));
    391     // TODO: Perhaps instead check whether all packet elements are in
    392     // [-8, 8] and evaluate a branch based off of that. It's possible
    393     // in practice most elements are in this region.
    394     y = pselect(pcmp_le(y, pset1<T>(8.0)), y_le_eight, y_gt_eight);
    395     return pselect(pcmp_lt(x, pset1<T>(0.0)), pnegate(y), y);
    396   }
    397 };
    398 
    399 template <typename T>
    400 struct bessel_i1e_impl {
    401   EIGEN_DEVICE_FUNC
    402   static EIGEN_STRONG_INLINE T run(const T x) {
    403     return generic_i1e<T>::run(x);
    404   }
    405 };
    406 
    407 template <typename T>
    408 struct bessel_i1_retval {
    409   typedef T type;
    410 };
    411 
    412 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
    413 struct generic_i1 {
    414   EIGEN_DEVICE_FUNC
    415   static EIGEN_STRONG_INLINE T run(const T& x) {
    416     return pmul(
    417         pexp(pabs(x)),
    418         generic_i1e<T, ScalarType>::run(x));
    419   }
    420 };
    421 
    422 template <typename T>
    423 struct bessel_i1_impl {
    424   EIGEN_DEVICE_FUNC
    425   static EIGEN_STRONG_INLINE T run(const T x) {
    426     return generic_i1<T>::run(x);
    427   }
    428 };
    429 
    430 template <typename T>
    431 struct bessel_k0e_retval {
    432   typedef T type;
    433 };
    434 
    435 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
    436 struct generic_k0e {
    437   EIGEN_DEVICE_FUNC
    438   static EIGEN_STRONG_INLINE T run(const T&) {
    439     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
    440                         THIS_TYPE_IS_NOT_SUPPORTED);
    441     return ScalarType(0);
    442   }
    443 };
    444 
    445 template <typename T>
    446 struct generic_k0e<T, float> {
    447   EIGEN_DEVICE_FUNC
    448   static EIGEN_STRONG_INLINE T run(const T& x) {
    449     /*  k0ef.c
    450      *	Modified Bessel function, third kind, order zero,
    451      *	exponentially scaled
    452      *
    453      *
    454      *
    455      * SYNOPSIS:
    456      *
    457      * float x, y, k0ef();
    458      *
    459      * y = k0ef( x );
    460      *
    461      *
    462      *
    463      * DESCRIPTION:
    464      *
    465      * Returns exponentially scaled modified Bessel function
    466      * of the third kind of order zero of the argument.
    467      *
    468      *
    469      *
    470      * ACCURACY:
    471      *
    472      *                      Relative error:
    473      * arithmetic   domain     # trials      peak         rms
    474      *    IEEE      0, 30       30000       8.1e-7      7.8e-8
    475      * See k0().
    476      *
    477      */
    478 
    479     const float A[] = {1.90451637722020886025E-9f, 2.53479107902614945675E-7f,
    480                        2.28621210311945178607E-5f, 1.26461541144692592338E-3f,
    481                        3.59799365153615016266E-2f, 3.44289899924628486886E-1f,
    482                        -5.35327393233902768720E-1f};
    483 
    484     const float B[] = {-1.69753450938905987466E-9f, 8.57403401741422608519E-9f,
    485                        -4.66048989768794782956E-8f, 2.76681363944501510342E-7f,
    486                        -1.83175552271911948767E-6f, 1.39498137188764993662E-5f,
    487                        -1.28495495816278026384E-4f, 1.56988388573005337491E-3f,
    488                        -3.14481013119645005427E-2f, 2.44030308206595545468E0f};
    489     const T MAXNUM = pset1<T>(NumTraits<float>::infinity());
    490     const T two = pset1<T>(2.0);
    491     T x_le_two = internal::pchebevl<T, 7>::run(
    492         pmadd(x, x, pset1<T>(-2.0)), A);
    493     x_le_two = pmadd(
    494         generic_i0<T, float>::run(x), pnegate(
    495             plog(pmul(pset1<T>(0.5), x))), x_le_two);
    496     x_le_two = pmul(pexp(x), x_le_two);
    497     T x_gt_two = pmul(
    498             internal::pchebevl<T, 10>::run(
    499                 psub(pdiv(pset1<T>(8.0), x), two), B),
    500             prsqrt(x));
    501     return pselect(
    502         pcmp_le(x, pset1<T>(0.0)),
    503         MAXNUM,
    504         pselect(pcmp_le(x, two), x_le_two, x_gt_two));
    505   }
    506 };
    507 
    508 template <typename T>
    509 struct generic_k0e<T, double> {
    510   EIGEN_DEVICE_FUNC
    511   static EIGEN_STRONG_INLINE T run(const T& x) {
    512     /*  k0e.c
    513      *	Modified Bessel function, third kind, order zero,
    514      *	exponentially scaled
    515      *
    516      *
    517      *
    518      * SYNOPSIS:
    519      *
    520      * double x, y, k0e();
    521      *
    522      * y = k0e( x );
    523      *
    524      *
    525      *
    526      * DESCRIPTION:
    527      *
    528      * Returns exponentially scaled modified Bessel function
    529      * of the third kind of order zero of the argument.
    530      *
    531      *
    532      *
    533      * ACCURACY:
    534      *
    535      *                      Relative error:
    536      * arithmetic   domain     # trials      peak         rms
    537      *    IEEE      0, 30       30000       1.4e-15     1.4e-16
    538      * See k0().
    539      *
    540      */
    541 
    542     const double A[] = {
    543       1.37446543561352307156E-16,
    544       4.25981614279661018399E-14,
    545       1.03496952576338420167E-11,
    546       1.90451637722020886025E-9,
    547       2.53479107902614945675E-7,
    548       2.28621210311945178607E-5,
    549       1.26461541144692592338E-3,
    550       3.59799365153615016266E-2,
    551       3.44289899924628486886E-1,
    552       -5.35327393233902768720E-1};
    553     const double B[] = {
    554        5.30043377268626276149E-18, -1.64758043015242134646E-17,
    555        5.21039150503902756861E-17, -1.67823109680541210385E-16,
    556        5.51205597852431940784E-16, -1.84859337734377901440E-15,
    557        6.34007647740507060557E-15, -2.22751332699166985548E-14,
    558        8.03289077536357521100E-14, -2.98009692317273043925E-13,
    559        1.14034058820847496303E-12, -4.51459788337394416547E-12,
    560        1.85594911495471785253E-11, -7.95748924447710747776E-11,
    561        3.57739728140030116597E-10, -1.69753450938905987466E-9,
    562        8.57403401741422608519E-9, -4.66048989768794782956E-8,
    563        2.76681363944501510342E-7, -1.83175552271911948767E-6,
    564        1.39498137188764993662E-5, -1.28495495816278026384E-4,
    565        1.56988388573005337491E-3, -3.14481013119645005427E-2,
    566        2.44030308206595545468E0
    567     };
    568     const T MAXNUM = pset1<T>(NumTraits<double>::infinity());
    569     const T two = pset1<T>(2.0);
    570     T x_le_two = internal::pchebevl<T, 10>::run(
    571         pmadd(x, x, pset1<T>(-2.0)), A);
    572     x_le_two = pmadd(
    573         generic_i0<T, double>::run(x), pmul(
    574             pset1<T>(-1.0), plog(pmul(pset1<T>(0.5), x))), x_le_two);
    575     x_le_two = pmul(pexp(x), x_le_two);
    576     x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
    577     T x_gt_two = pmul(
    578             internal::pchebevl<T, 25>::run(
    579                 psub(pdiv(pset1<T>(8.0), x), two), B),
    580             prsqrt(x));
    581     return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
    582   }
    583 };
    584 
    585 template <typename T>
    586 struct bessel_k0e_impl {
    587   EIGEN_DEVICE_FUNC
    588   static EIGEN_STRONG_INLINE T run(const T x) {
    589     return generic_k0e<T>::run(x);
    590   }
    591 };
    592 
    593 template <typename T>
    594 struct bessel_k0_retval {
    595   typedef T type;
    596 };
    597 
    598 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
    599 struct generic_k0 {
    600   EIGEN_DEVICE_FUNC
    601   static EIGEN_STRONG_INLINE T run(const T&) {
    602     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
    603                         THIS_TYPE_IS_NOT_SUPPORTED);
    604     return ScalarType(0);
    605   }
    606 };
    607 
    608 template <typename T>
    609 struct generic_k0<T, float> {
    610   EIGEN_DEVICE_FUNC
    611   static EIGEN_STRONG_INLINE T run(const T& x) {
    612     /*  k0f.c
    613      *	Modified Bessel function, third kind, order zero
    614      *
    615      *
    616      *
    617      * SYNOPSIS:
    618      *
    619      * float x, y, k0f();
    620      *
    621      * y = k0f( x );
    622      *
    623      *
    624      *
    625      * DESCRIPTION:
    626      *
    627      * Returns modified Bessel function of the third kind
    628      * of order zero of the argument.
    629      *
    630      * The range is partitioned into the two intervals [0,8] and
    631      * (8, infinity).  Chebyshev polynomial expansions are employed
    632      * in each interval.
    633      *
    634      *
    635      *
    636      * ACCURACY:
    637      *
    638      * Tested at 2000 random points between 0 and 8.  Peak absolute
    639      * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
    640      *                      Relative error:
    641      * arithmetic   domain     # trials      peak         rms
    642      *    IEEE      0, 30       30000       7.8e-7      8.5e-8
    643      *
    644      * ERROR MESSAGES:
    645      *
    646      *   message         condition      value returned
    647      *  K0 domain          x <= 0          MAXNUM
    648      *
    649      */
    650 
    651     const float A[] = {1.90451637722020886025E-9f, 2.53479107902614945675E-7f,
    652                        2.28621210311945178607E-5f, 1.26461541144692592338E-3f,
    653                        3.59799365153615016266E-2f, 3.44289899924628486886E-1f,
    654                        -5.35327393233902768720E-1f};
    655 
    656     const float B[] = {-1.69753450938905987466E-9f, 8.57403401741422608519E-9f,
    657                        -4.66048989768794782956E-8f, 2.76681363944501510342E-7f,
    658                        -1.83175552271911948767E-6f, 1.39498137188764993662E-5f,
    659                        -1.28495495816278026384E-4f, 1.56988388573005337491E-3f,
    660                        -3.14481013119645005427E-2f, 2.44030308206595545468E0f};
    661     const T MAXNUM = pset1<T>(NumTraits<float>::infinity());
    662     const T two = pset1<T>(2.0);
    663     T x_le_two = internal::pchebevl<T, 7>::run(
    664         pmadd(x, x, pset1<T>(-2.0)), A);
    665     x_le_two = pmadd(
    666         generic_i0<T, float>::run(x), pnegate(
    667             plog(pmul(pset1<T>(0.5), x))), x_le_two);
    668     x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
    669     T x_gt_two = pmul(
    670         pmul(
    671             pexp(pnegate(x)),
    672             internal::pchebevl<T, 10>::run(
    673                 psub(pdiv(pset1<T>(8.0), x), two), B)),
    674         prsqrt(x));
    675     return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
    676   }
    677 };
    678 
    679 template <typename T>
    680 struct generic_k0<T, double> {
    681   EIGEN_DEVICE_FUNC
    682   static EIGEN_STRONG_INLINE T run(const T& x) {
    683     /*
    684      *
    685      *	Modified Bessel function, third kind, order zero,
    686      *	exponentially scaled
    687      *
    688      *
    689      *
    690      * SYNOPSIS:
    691      *
    692      * double x, y, k0();
    693      *
    694      * y = k0( x );
    695      *
    696      *
    697      *
    698      * DESCRIPTION:
    699      *
    700      * Returns exponentially scaled modified Bessel function
    701      * of the third kind of order zero of the argument.
    702      *
    703      *
    704      *
    705      * ACCURACY:
    706      *
    707      *                      Relative error:
    708      * arithmetic   domain     # trials      peak         rms
    709      *    IEEE      0, 30       30000       1.4e-15     1.4e-16
    710      * See k0().
    711      *
    712      */
    713     const double A[] = {
    714       1.37446543561352307156E-16,
    715       4.25981614279661018399E-14,
    716       1.03496952576338420167E-11,
    717       1.90451637722020886025E-9,
    718       2.53479107902614945675E-7,
    719       2.28621210311945178607E-5,
    720       1.26461541144692592338E-3,
    721       3.59799365153615016266E-2,
    722       3.44289899924628486886E-1,
    723       -5.35327393233902768720E-1};
    724     const double B[] = {
    725        5.30043377268626276149E-18, -1.64758043015242134646E-17,
    726        5.21039150503902756861E-17, -1.67823109680541210385E-16,
    727        5.51205597852431940784E-16, -1.84859337734377901440E-15,
    728        6.34007647740507060557E-15, -2.22751332699166985548E-14,
    729        8.03289077536357521100E-14, -2.98009692317273043925E-13,
    730        1.14034058820847496303E-12, -4.51459788337394416547E-12,
    731        1.85594911495471785253E-11, -7.95748924447710747776E-11,
    732        3.57739728140030116597E-10, -1.69753450938905987466E-9,
    733        8.57403401741422608519E-9, -4.66048989768794782956E-8,
    734        2.76681363944501510342E-7, -1.83175552271911948767E-6,
    735        1.39498137188764993662E-5, -1.28495495816278026384E-4,
    736        1.56988388573005337491E-3, -3.14481013119645005427E-2,
    737        2.44030308206595545468E0
    738     };
    739     const T MAXNUM = pset1<T>(NumTraits<double>::infinity());
    740     const T two = pset1<T>(2.0);
    741     T x_le_two = internal::pchebevl<T, 10>::run(
    742         pmadd(x, x, pset1<T>(-2.0)), A);
    743     x_le_two = pmadd(
    744         generic_i0<T, double>::run(x), pnegate(
    745             plog(pmul(pset1<T>(0.5), x))), x_le_two);
    746     x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
    747     T x_gt_two = pmul(
    748         pmul(
    749             pexp(-x),
    750             internal::pchebevl<T, 25>::run(
    751                 psub(pdiv(pset1<T>(8.0), x), two), B)),
    752         prsqrt(x));
    753     return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
    754   }
    755 };
    756 
    757 template <typename T>
    758 struct bessel_k0_impl {
    759   EIGEN_DEVICE_FUNC
    760   static EIGEN_STRONG_INLINE T run(const T x) {
    761     return generic_k0<T>::run(x);
    762   }
    763 };
    764 
    765 template <typename T>
    766 struct bessel_k1e_retval {
    767   typedef T type;
    768 };
    769 
    770 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
    771 struct generic_k1e {
    772   EIGEN_DEVICE_FUNC
    773   static EIGEN_STRONG_INLINE T run(const T&) {
    774     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
    775                         THIS_TYPE_IS_NOT_SUPPORTED);
    776     return ScalarType(0);
    777   }
    778 };
    779 
    780 template <typename T>
    781 struct generic_k1e<T, float> {
    782   EIGEN_DEVICE_FUNC
    783   static EIGEN_STRONG_INLINE T run(const T& x) {
    784     /* k1ef.c
    785      *
    786      *	Modified Bessel function, third kind, order one,
    787      *	exponentially scaled
    788      *
    789      *
    790      *
    791      * SYNOPSIS:
    792      *
    793      * float x, y, k1ef();
    794      *
    795      * y = k1ef( x );
    796      *
    797      *
    798      *
    799      * DESCRIPTION:
    800      *
    801      * Returns exponentially scaled modified Bessel function
    802      * of the third kind of order one of the argument:
    803      *
    804      *      k1e(x) = exp(x) * k1(x).
    805      *
    806      *
    807      *
    808      * ACCURACY:
    809      *
    810      *                      Relative error:
    811      * arithmetic   domain     # trials      peak         rms
    812      *    IEEE      0, 30       30000       4.9e-7      6.7e-8
    813      * See k1().
    814      *
    815      */
    816 
    817     const float A[] = {-2.21338763073472585583E-8f, -2.43340614156596823496E-6f,
    818                         -1.73028895751305206302E-4f, -6.97572385963986435018E-3f,
    819                         -1.22611180822657148235E-1f, -3.53155960776544875667E-1f,
    820                         1.52530022733894777053E0f};
    821     const float B[] = {2.01504975519703286596E-9f, -1.03457624656780970260E-8f,
    822                        5.74108412545004946722E-8f, -3.50196060308781257119E-7f,
    823                        2.40648494783721712015E-6f, -1.93619797416608296024E-5f,
    824                        1.95215518471351631108E-4f, -2.85781685962277938680E-3f,
    825                        1.03923736576817238437E-1f, 2.72062619048444266945E0f};
    826     const T MAXNUM = pset1<T>(NumTraits<float>::infinity());
    827     const T two = pset1<T>(2.0);
    828     T x_le_two = pdiv(internal::pchebevl<T, 7>::run(
    829         pmadd(x, x, pset1<T>(-2.0)), A), x);
    830     x_le_two = pmadd(
    831         generic_i1<T, float>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two);
    832     x_le_two = pmul(x_le_two, pexp(x));
    833     x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
    834     T x_gt_two = pmul(
    835         internal::pchebevl<T, 10>::run(
    836             psub(pdiv(pset1<T>(8.0), x), two), B),
    837         prsqrt(x));
    838     return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
    839   }
    840 };
    841 
    842 template <typename T>
    843 struct generic_k1e<T, double> {
    844   EIGEN_DEVICE_FUNC
    845   static EIGEN_STRONG_INLINE T run(const T& x) {
    846     /*  k1e.c
    847      *
    848      *	Modified Bessel function, third kind, order one,
    849      *	exponentially scaled
    850      *
    851      *
    852      *
    853      * SYNOPSIS:
    854      *
    855      * double x, y, k1e();
    856      *
    857      * y = k1e( x );
    858      *
    859      *
    860      *
    861      * DESCRIPTION:
    862      *
    863      * Returns exponentially scaled modified Bessel function
    864      * of the third kind of order one of the argument:
    865      *
    866      *      k1e(x) = exp(x) * k1(x).
    867      *
    868      *
    869      *
    870      * ACCURACY:
    871      *
    872      *                      Relative error:
    873      * arithmetic   domain     # trials      peak         rms
    874      *    IEEE      0, 30       30000       7.8e-16     1.2e-16
    875      * See k1().
    876      *
    877      */
    878     const double A[] = {-7.02386347938628759343E-18, -2.42744985051936593393E-15,
    879                         -6.66690169419932900609E-13, -1.41148839263352776110E-10,
    880                         -2.21338763073472585583E-8, -2.43340614156596823496E-6,
    881                         -1.73028895751305206302E-4, -6.97572385963986435018E-3,
    882                         -1.22611180822657148235E-1, -3.53155960776544875667E-1,
    883                         1.52530022733894777053E0};
    884     const double B[] = {-5.75674448366501715755E-18, 1.79405087314755922667E-17,
    885                         -5.68946255844285935196E-17, 1.83809354436663880070E-16,
    886                         -6.05704724837331885336E-16, 2.03870316562433424052E-15,
    887                         -7.01983709041831346144E-15, 2.47715442448130437068E-14,
    888                         -8.97670518232499435011E-14, 3.34841966607842919884E-13,
    889                         -1.28917396095102890680E-12, 5.13963967348173025100E-12,
    890                         -2.12996783842756842877E-11, 9.21831518760500529508E-11,
    891                         -4.19035475934189648750E-10, 2.01504975519703286596E-9,
    892                         -1.03457624656780970260E-8, 5.74108412545004946722E-8,
    893                         -3.50196060308781257119E-7, 2.40648494783721712015E-6,
    894                         -1.93619797416608296024E-5, 1.95215518471351631108E-4,
    895                         -2.85781685962277938680E-3, 1.03923736576817238437E-1,
    896                         2.72062619048444266945E0};
    897     const T MAXNUM = pset1<T>(NumTraits<double>::infinity());
    898     const T two = pset1<T>(2.0);
    899     T x_le_two = pdiv(internal::pchebevl<T, 11>::run(
    900         pmadd(x, x, pset1<T>(-2.0)), A), x);
    901     x_le_two = pmadd(
    902         generic_i1<T, double>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two);
    903     x_le_two = pmul(x_le_two, pexp(x));
    904     x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
    905     T x_gt_two = pmul(
    906         internal::pchebevl<T, 25>::run(
    907             psub(pdiv(pset1<T>(8.0), x), two), B),
    908         prsqrt(x));
    909     return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
    910   }
    911 };
    912 
    913 template <typename T>
    914 struct bessel_k1e_impl {
    915   EIGEN_DEVICE_FUNC
    916   static EIGEN_STRONG_INLINE T run(const T x) {
    917     return generic_k1e<T>::run(x);
    918   }
    919 };
    920 
    921 template <typename T>
    922 struct bessel_k1_retval {
    923   typedef T type;
    924 };
    925 
    926 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
    927 struct generic_k1 {
    928   EIGEN_DEVICE_FUNC
    929   static EIGEN_STRONG_INLINE T run(const T&) {
    930     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
    931                         THIS_TYPE_IS_NOT_SUPPORTED);
    932     return ScalarType(0);
    933   }
    934 };
    935 
    936 template <typename T>
    937 struct generic_k1<T, float> {
    938   EIGEN_DEVICE_FUNC
    939   static EIGEN_STRONG_INLINE T run(const T& x) {
    940     /* k1f.c
    941      *	Modified Bessel function, third kind, order one
    942      *
    943      *
    944      *
    945      * SYNOPSIS:
    946      *
    947      * float x, y, k1f();
    948      *
    949      * y = k1f( x );
    950      *
    951      *
    952      *
    953      * DESCRIPTION:
    954      *
    955      * Computes the modified Bessel function of the third kind
    956      * of order one of the argument.
    957      *
    958      * The range is partitioned into the two intervals [0,2] and
    959      * (2, infinity).  Chebyshev polynomial expansions are employed
    960      * in each interval.
    961      *
    962      *
    963      *
    964      * ACCURACY:
    965      *
    966      *                      Relative error:
    967      * arithmetic   domain     # trials      peak         rms
    968      *    IEEE      0, 30       30000       4.6e-7      7.6e-8
    969      *
    970      * ERROR MESSAGES:
    971      *
    972      *   message         condition      value returned
    973      * k1 domain          x <= 0          MAXNUM
    974      *
    975      */
    976 
    977     const float A[] = {-2.21338763073472585583E-8f, -2.43340614156596823496E-6f,
    978                         -1.73028895751305206302E-4f, -6.97572385963986435018E-3f,
    979                         -1.22611180822657148235E-1f, -3.53155960776544875667E-1f,
    980                         1.52530022733894777053E0f};
    981     const float B[] = {2.01504975519703286596E-9f, -1.03457624656780970260E-8f,
    982                        5.74108412545004946722E-8f, -3.50196060308781257119E-7f,
    983                        2.40648494783721712015E-6f, -1.93619797416608296024E-5f,
    984                        1.95215518471351631108E-4f, -2.85781685962277938680E-3f,
    985                        1.03923736576817238437E-1f, 2.72062619048444266945E0f};
    986     const T MAXNUM = pset1<T>(NumTraits<float>::infinity());
    987     const T two = pset1<T>(2.0);
    988     T x_le_two = pdiv(internal::pchebevl<T, 7>::run(
    989         pmadd(x, x, pset1<T>(-2.0)), A), x);
    990     x_le_two = pmadd(
    991         generic_i1<T, float>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two);
    992     x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
    993     T x_gt_two = pmul(
    994         pexp(pnegate(x)),
    995         pmul(
    996             internal::pchebevl<T, 10>::run(
    997                 psub(pdiv(pset1<T>(8.0), x), two), B),
    998             prsqrt(x)));
    999     return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
   1000   }
   1001 };
   1002 
   1003 template <typename T>
   1004 struct generic_k1<T, double> {
   1005   EIGEN_DEVICE_FUNC
   1006   static EIGEN_STRONG_INLINE T run(const T& x) {
   1007     /*  k1.c
   1008      *	Modified Bessel function, third kind, order one
   1009      *
   1010      *
   1011      *
   1012      * SYNOPSIS:
   1013      *
   1014      * float x, y, k1f();
   1015      *
   1016      * y = k1f( x );
   1017      *
   1018      *
   1019      *
   1020      * DESCRIPTION:
   1021      *
   1022      * Computes the modified Bessel function of the third kind
   1023      * of order one of the argument.
   1024      *
   1025      * The range is partitioned into the two intervals [0,2] and
   1026      * (2, infinity).  Chebyshev polynomial expansions are employed
   1027      * in each interval.
   1028      *
   1029      *
   1030      *
   1031      * ACCURACY:
   1032      *
   1033      *                      Relative error:
   1034      * arithmetic   domain     # trials      peak         rms
   1035      *    IEEE      0, 30       30000       4.6e-7      7.6e-8
   1036      *
   1037      * ERROR MESSAGES:
   1038      *
   1039      *   message         condition      value returned
   1040      * k1 domain          x <= 0          MAXNUM
   1041      *
   1042      */
   1043     const double A[] = {-7.02386347938628759343E-18, -2.42744985051936593393E-15,
   1044                         -6.66690169419932900609E-13, -1.41148839263352776110E-10,
   1045                         -2.21338763073472585583E-8, -2.43340614156596823496E-6,
   1046                         -1.73028895751305206302E-4, -6.97572385963986435018E-3,
   1047                         -1.22611180822657148235E-1, -3.53155960776544875667E-1,
   1048                         1.52530022733894777053E0};
   1049     const double B[] = {-5.75674448366501715755E-18, 1.79405087314755922667E-17,
   1050                         -5.68946255844285935196E-17, 1.83809354436663880070E-16,
   1051                         -6.05704724837331885336E-16, 2.03870316562433424052E-15,
   1052                         -7.01983709041831346144E-15, 2.47715442448130437068E-14,
   1053                         -8.97670518232499435011E-14, 3.34841966607842919884E-13,
   1054                         -1.28917396095102890680E-12, 5.13963967348173025100E-12,
   1055                         -2.12996783842756842877E-11, 9.21831518760500529508E-11,
   1056                         -4.19035475934189648750E-10, 2.01504975519703286596E-9,
   1057                         -1.03457624656780970260E-8, 5.74108412545004946722E-8,
   1058                         -3.50196060308781257119E-7, 2.40648494783721712015E-6,
   1059                         -1.93619797416608296024E-5, 1.95215518471351631108E-4,
   1060                         -2.85781685962277938680E-3, 1.03923736576817238437E-1,
   1061                         2.72062619048444266945E0};
   1062     const T MAXNUM = pset1<T>(NumTraits<double>::infinity());
   1063     const T two = pset1<T>(2.0);
   1064     T x_le_two = pdiv(internal::pchebevl<T, 11>::run(
   1065         pmadd(x, x, pset1<T>(-2.0)), A), x);
   1066     x_le_two = pmadd(
   1067         generic_i1<T, double>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two);
   1068     x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
   1069     T x_gt_two = pmul(
   1070         pexp(-x),
   1071         pmul(
   1072             internal::pchebevl<T, 25>::run(
   1073                 psub(pdiv(pset1<T>(8.0), x), two), B),
   1074             prsqrt(x)));
   1075     return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
   1076   }
   1077 };
   1078 
   1079 template <typename T>
   1080 struct bessel_k1_impl {
   1081   EIGEN_DEVICE_FUNC
   1082   static EIGEN_STRONG_INLINE T run(const T x) {
   1083     return generic_k1<T>::run(x);
   1084   }
   1085 };
   1086 
   1087 template <typename T>
   1088 struct bessel_j0_retval {
   1089   typedef T type;
   1090 };
   1091 
   1092 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
   1093 struct generic_j0 {
   1094   EIGEN_DEVICE_FUNC
   1095   static EIGEN_STRONG_INLINE T run(const T&) {
   1096     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
   1097                         THIS_TYPE_IS_NOT_SUPPORTED);
   1098     return ScalarType(0);
   1099   }
   1100 };
   1101 
   1102 template <typename T>
   1103 struct generic_j0<T, float> {
   1104   EIGEN_DEVICE_FUNC
   1105   static EIGEN_STRONG_INLINE T run(const T& x) {
   1106     /* j0f.c
   1107      *	Bessel function of order zero
   1108      *
   1109      *
   1110      *
   1111      * SYNOPSIS:
   1112      *
   1113      * float x, y, j0f();
   1114      *
   1115      * y = j0f( x );
   1116      *
   1117      *
   1118      *
   1119      * DESCRIPTION:
   1120      *
   1121      * Returns Bessel function of order zero of the argument.
   1122      *
   1123      * The domain is divided into the intervals [0, 2] and
   1124      * (2, infinity). In the first interval the following polynomial
   1125      * approximation is used:
   1126      *
   1127      *
   1128      *        2         2         2
   1129      * (w - r  ) (w - r  ) (w - r  ) P(w)
   1130      *       1         2         3
   1131      *
   1132      *            2
   1133      * where w = x  and the three r's are zeros of the function.
   1134      *
   1135      * In the second interval, the modulus and phase are approximated
   1136      * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
   1137      * and Phase(x) = x + 1/x R(1/x^2) - pi/4.  The function is
   1138      *
   1139      *   j0(x) = Modulus(x) cos( Phase(x) ).
   1140      *
   1141      *
   1142      *
   1143      * ACCURACY:
   1144      *
   1145      *                      Absolute error:
   1146      * arithmetic   domain     # trials      peak         rms
   1147      *    IEEE      0, 2        100000      1.3e-7      3.6e-8
   1148      *    IEEE      2, 32       100000      1.9e-7      5.4e-8
   1149      *
   1150      */
   1151 
   1152     const float JP[] = {-6.068350350393235E-008f, 6.388945720783375E-006f,
   1153                         -3.969646342510940E-004f, 1.332913422519003E-002f,
   1154                         -1.729150680240724E-001f};
   1155     const float MO[] = {-6.838999669318810E-002f, 1.864949361379502E-001f,
   1156                         -2.145007480346739E-001f, 1.197549369473540E-001f,
   1157                         -3.560281861530129E-003f, -4.969382655296620E-002f,
   1158                         -3.355424622293709E-006f, 7.978845717621440E-001f};
   1159     const float PH[] = {3.242077816988247E+001f, -3.630592630518434E+001f,
   1160                         1.756221482109099E+001f, -4.974978466280903E+000f,
   1161                         1.001973420681837E+000f, -1.939906941791308E-001f,
   1162                         6.490598792654666E-002f, -1.249992184872738E-001f};
   1163     const T DR1 =  pset1<T>(5.78318596294678452118f);
   1164     const T NEG_PIO4F = pset1<T>(-0.7853981633974483096f); /* -pi / 4 */
   1165     T y = pabs(x);
   1166     T z = pmul(y, y);
   1167     T y_le_two = pselect(
   1168         pcmp_lt(y, pset1<T>(1.0e-3f)),
   1169         pmadd(z, pset1<T>(-0.25f), pset1<T>(1.0f)),
   1170         pmul(psub(z, DR1), internal::ppolevl<T, 4>::run(z, JP)));
   1171     T q = pdiv(pset1<T>(1.0f), y);
   1172     T w = prsqrt(y);
   1173     T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO));
   1174     w = pmul(q, q);
   1175     T yn = pmadd(q, internal::ppolevl<T, 7>::run(w, PH), NEG_PIO4F);
   1176     T y_gt_two = pmul(p, pcos(padd(yn, y)));
   1177     return pselect(pcmp_le(y, pset1<T>(2.0)), y_le_two, y_gt_two);
   1178   }
   1179 };
   1180 
   1181 template <typename T>
   1182 struct generic_j0<T, double> {
   1183   EIGEN_DEVICE_FUNC
   1184   static EIGEN_STRONG_INLINE T run(const T& x) {
   1185     /*  j0.c
   1186      *	Bessel function of order zero
   1187      *
   1188      *
   1189      *
   1190      * SYNOPSIS:
   1191      *
   1192      * double x, y, j0();
   1193      *
   1194      * y = j0( x );
   1195      *
   1196      *
   1197      *
   1198      * DESCRIPTION:
   1199      *
   1200      * Returns Bessel function of order zero of the argument.
   1201      *
   1202      * The domain is divided into the intervals [0, 5] and
   1203      * (5, infinity). In the first interval the following rational
   1204      * approximation is used:
   1205      *
   1206      *
   1207      *        2         2
   1208      * (w - r  ) (w - r  ) P (w) / Q (w)
   1209      *       1         2    3       8
   1210      *
   1211      *            2
   1212      * where w = x  and the two r's are zeros of the function.
   1213      *
   1214      * In the second interval, the Hankel asymptotic expansion
   1215      * is employed with two rational functions of degree 6/6
   1216      * and 7/7.
   1217      *
   1218      *
   1219      *
   1220      * ACCURACY:
   1221      *
   1222      *                      Absolute error:
   1223      * arithmetic   domain     # trials      peak         rms
   1224      *    DEC       0, 30       10000       4.4e-17     6.3e-18
   1225      *    IEEE      0, 30       60000       4.2e-16     1.1e-16
   1226      *
   1227      */
   1228     const double PP[] = {7.96936729297347051624E-4, 8.28352392107440799803E-2,
   1229                         1.23953371646414299388E0, 5.44725003058768775090E0,
   1230                         8.74716500199817011941E0, 5.30324038235394892183E0,
   1231                         9.99999999999999997821E-1};
   1232     const double PQ[] = {9.24408810558863637013E-4, 8.56288474354474431428E-2,
   1233                          1.25352743901058953537E0, 5.47097740330417105182E0,
   1234                          8.76190883237069594232E0, 5.30605288235394617618E0,
   1235                          1.00000000000000000218E0};
   1236     const double QP[] = {-1.13663838898469149931E-2, -1.28252718670509318512E0,
   1237                          -1.95539544257735972385E1, -9.32060152123768231369E1,
   1238                          -1.77681167980488050595E2, -1.47077505154951170175E2,
   1239                          -5.14105326766599330220E1, -6.05014350600728481186E0};
   1240     const double QQ[] = {1.00000000000000000000E0, 6.43178256118178023184E1,
   1241                          8.56430025976980587198E2, 3.88240183605401609683E3,
   1242                          7.24046774195652478189E3, 5.93072701187316984827E3,
   1243                          2.06209331660327847417E3, 2.42005740240291393179E2};
   1244     const double RP[] = {-4.79443220978201773821E9, 1.95617491946556577543E12,
   1245                          -2.49248344360967716204E14, 9.70862251047306323952E15};
   1246     const double RQ[] = {1.00000000000000000000E0, 4.99563147152651017219E2,
   1247                          1.73785401676374683123E5, 4.84409658339962045305E7,
   1248                          1.11855537045356834862E10, 2.11277520115489217587E12,
   1249                          3.10518229857422583814E14, 3.18121955943204943306E16,
   1250                          1.71086294081043136091E18};
   1251     const T DR1 = pset1<T>(5.78318596294678452118E0);
   1252     const T DR2 = pset1<T>(3.04712623436620863991E1);
   1253     const T SQ2OPI = pset1<T>(7.9788456080286535587989E-1); /* sqrt(2 / pi) */
   1254     const T NEG_PIO4 = pset1<T>(-0.7853981633974483096); /* pi / 4 */
   1255 
   1256     T y = pabs(x);
   1257     T z = pmul(y, y);
   1258     T y_le_five = pselect(
   1259         pcmp_lt(y, pset1<T>(1.0e-5)),
   1260         pmadd(z, pset1<T>(-0.25), pset1<T>(1.0)),
   1261         pmul(pmul(psub(z, DR1), psub(z, DR2)),
   1262              pdiv(internal::ppolevl<T, 3>::run(z, RP),
   1263                   internal::ppolevl<T, 8>::run(z, RQ))));
   1264     T s = pdiv(pset1<T>(25.0), z);
   1265     T p = pdiv(
   1266         internal::ppolevl<T, 6>::run(s, PP),
   1267         internal::ppolevl<T, 6>::run(s, PQ));
   1268     T q = pdiv(
   1269         internal::ppolevl<T, 7>::run(s, QP),
   1270         internal::ppolevl<T, 7>::run(s, QQ));
   1271     T yn = padd(y, NEG_PIO4);
   1272     T w = pdiv(pset1<T>(-5.0), y);
   1273     p = pmadd(p, pcos(yn), pmul(w, pmul(q, psin(yn))));
   1274     T y_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(y)));
   1275     return pselect(pcmp_le(y, pset1<T>(5.0)), y_le_five, y_gt_five);
   1276   }
   1277 };
   1278 
   1279 template <typename T>
   1280 struct bessel_j0_impl {
   1281   EIGEN_DEVICE_FUNC
   1282   static EIGEN_STRONG_INLINE T run(const T x) {
   1283     return generic_j0<T>::run(x);
   1284   }
   1285 };
   1286 
   1287 template <typename T>
   1288 struct bessel_y0_retval {
   1289   typedef T type;
   1290 };
   1291 
   1292 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
   1293 struct generic_y0 {
   1294   EIGEN_DEVICE_FUNC
   1295   static EIGEN_STRONG_INLINE T run(const T&) {
   1296     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
   1297                         THIS_TYPE_IS_NOT_SUPPORTED);
   1298     return ScalarType(0);
   1299   }
   1300 };
   1301 
   1302 template <typename T>
   1303 struct generic_y0<T, float> {
   1304   EIGEN_DEVICE_FUNC
   1305   static EIGEN_STRONG_INLINE T run(const T& x) {
   1306     /* j0f.c
   1307      * 	Bessel function of the second kind, order zero
   1308      *
   1309      *
   1310      *
   1311      * SYNOPSIS:
   1312      *
   1313      * float x, y, y0f();
   1314      *
   1315      * y = y0f( x );
   1316      *
   1317      *
   1318      *
   1319      * DESCRIPTION:
   1320      *
   1321      * Returns Bessel function of the second kind, of order
   1322      * zero, of the argument.
   1323      *
   1324      * The domain is divided into the intervals [0, 2] and
   1325      * (2, infinity). In the first interval a rational approximation
   1326      * R(x) is employed to compute
   1327      *
   1328      *                  2         2         2
   1329      * y0(x)  =  (w - r  ) (w - r  ) (w - r  ) R(x)  +  2/pi ln(x) j0(x).
   1330      *                 1         2         3
   1331      *
   1332      * Thus a call to j0() is required.  The three zeros are removed
   1333      * from R(x) to improve its numerical stability.
   1334      *
   1335      * In the second interval, the modulus and phase are approximated
   1336      * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
   1337      * and Phase(x) = x + 1/x S(1/x^2) - pi/4.  Then the function is
   1338      *
   1339      *   y0(x) = Modulus(x) sin( Phase(x) ).
   1340      *
   1341      *
   1342      *
   1343      *
   1344      * ACCURACY:
   1345      *
   1346      *  Absolute error, when y0(x) < 1; else relative error:
   1347      *
   1348      * arithmetic   domain     # trials      peak         rms
   1349      *    IEEE      0,  2       100000      2.4e-7      3.4e-8
   1350      *    IEEE      2, 32       100000      1.8e-7      5.3e-8
   1351      *
   1352      */
   1353 
   1354     const float YP[] = {9.454583683980369E-008f, -9.413212653797057E-006f,
   1355                         5.344486707214273E-004f, -1.584289289821316E-002f,
   1356                         1.707584643733568E-001f};
   1357     const float MO[] = {-6.838999669318810E-002f, 1.864949361379502E-001f,
   1358                         -2.145007480346739E-001f, 1.197549369473540E-001f,
   1359                         -3.560281861530129E-003f, -4.969382655296620E-002f,
   1360                         -3.355424622293709E-006f, 7.978845717621440E-001f};
   1361     const float PH[] = {3.242077816988247E+001f, -3.630592630518434E+001f,
   1362                         1.756221482109099E+001f, -4.974978466280903E+000f,
   1363                         1.001973420681837E+000f, -1.939906941791308E-001f,
   1364                         6.490598792654666E-002f, -1.249992184872738E-001f};
   1365     const T YZ1 = pset1<T>(0.43221455686510834878f);
   1366     const T TWOOPI =  pset1<T>(0.636619772367581343075535f); /* 2 / pi */
   1367     const T NEG_PIO4F = pset1<T>(-0.7853981633974483096f); /* -pi / 4 */
   1368     const T NEG_MAXNUM = pset1<T>(-NumTraits<float>::infinity());
   1369     T z = pmul(x, x);
   1370     T x_le_two = pmul(TWOOPI, pmul(plog(x), generic_j0<T, float>::run(x)));
   1371     x_le_two = pmadd(
   1372         psub(z, YZ1), internal::ppolevl<T, 4>::run(z, YP), x_le_two);
   1373     x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), NEG_MAXNUM, x_le_two);
   1374     T q = pdiv(pset1<T>(1.0), x);
   1375     T w = prsqrt(x);
   1376     T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO));
   1377     T u = pmul(q, q);
   1378     T xn = pmadd(q, internal::ppolevl<T, 7>::run(u, PH), NEG_PIO4F);
   1379     T x_gt_two = pmul(p, psin(padd(xn, x)));
   1380     return pselect(pcmp_le(x, pset1<T>(2.0)), x_le_two, x_gt_two);
   1381   }
   1382 };
   1383 
   1384 template <typename T>
   1385 struct generic_y0<T, double> {
   1386   EIGEN_DEVICE_FUNC
   1387   static EIGEN_STRONG_INLINE T run(const T& x) {
   1388     /*  j0.c
   1389      *	Bessel function of the second kind, order zero
   1390      *
   1391      *
   1392      *
   1393      * SYNOPSIS:
   1394      *
   1395      * double x, y, y0();
   1396      *
   1397      * y = y0( x );
   1398      *
   1399      *
   1400      *
   1401      * DESCRIPTION:
   1402      *
   1403      * Returns Bessel function of the second kind, of order
   1404      * zero, of the argument.
   1405      *
   1406      * The domain is divided into the intervals [0, 5] and
   1407      * (5, infinity). In the first interval a rational approximation
   1408      * R(x) is employed to compute
   1409      *   y0(x)  = R(x)  +   2 * log(x) * j0(x) / PI.
   1410      * Thus a call to j0() is required.
   1411      *
   1412      * In the second interval, the Hankel asymptotic expansion
   1413      * is employed with two rational functions of degree 6/6
   1414      * and 7/7.
   1415      *
   1416      *
   1417      *
   1418      * ACCURACY:
   1419      *
   1420      *  Absolute error, when y0(x) < 1; else relative error:
   1421      *
   1422      * arithmetic   domain     # trials      peak         rms
   1423      *    DEC       0, 30        9400       7.0e-17     7.9e-18
   1424      *    IEEE      0, 30       30000       1.3e-15     1.6e-16
   1425      *
   1426      */
   1427     const double PP[] = {7.96936729297347051624E-4, 8.28352392107440799803E-2,
   1428                         1.23953371646414299388E0, 5.44725003058768775090E0,
   1429                         8.74716500199817011941E0, 5.30324038235394892183E0,
   1430                         9.99999999999999997821E-1};
   1431     const double PQ[] = {9.24408810558863637013E-4, 8.56288474354474431428E-2,
   1432                          1.25352743901058953537E0, 5.47097740330417105182E0,
   1433                          8.76190883237069594232E0, 5.30605288235394617618E0,
   1434                          1.00000000000000000218E0};
   1435     const double QP[] = {-1.13663838898469149931E-2, -1.28252718670509318512E0,
   1436                          -1.95539544257735972385E1, -9.32060152123768231369E1,
   1437                          -1.77681167980488050595E2, -1.47077505154951170175E2,
   1438                          -5.14105326766599330220E1, -6.05014350600728481186E0};
   1439     const double QQ[] = {1.00000000000000000000E0, 6.43178256118178023184E1,
   1440                          8.56430025976980587198E2, 3.88240183605401609683E3,
   1441                          7.24046774195652478189E3, 5.93072701187316984827E3,
   1442                          2.06209331660327847417E3, 2.42005740240291393179E2};
   1443     const double YP[] = {1.55924367855235737965E4, -1.46639295903971606143E7,
   1444                          5.43526477051876500413E9, -9.82136065717911466409E11,
   1445                          8.75906394395366999549E13, -3.46628303384729719441E15,
   1446                          4.42733268572569800351E16, -1.84950800436986690637E16};
   1447     const double YQ[] = {1.00000000000000000000E0,  1.04128353664259848412E3,
   1448                          6.26107330137134956842E5, 2.68919633393814121987E8,
   1449                          8.64002487103935000337E10, 2.02979612750105546709E13,
   1450                          3.17157752842975028269E15, 2.50596256172653059228E17};
   1451     const T SQ2OPI = pset1<T>(7.9788456080286535587989E-1); /* sqrt(2 / pi) */
   1452     const T TWOOPI =  pset1<T>(0.636619772367581343075535); /* 2 / pi */
   1453     const T NEG_PIO4 = pset1<T>(-0.7853981633974483096); /* -pi / 4 */
   1454     const T NEG_MAXNUM = pset1<T>(-NumTraits<double>::infinity());
   1455 
   1456     T z = pmul(x, x);
   1457     T x_le_five = pdiv(internal::ppolevl<T, 7>::run(z, YP),
   1458                        internal::ppolevl<T, 7>::run(z, YQ));
   1459     x_le_five = pmadd(
   1460         pmul(TWOOPI, plog(x)), generic_j0<T, double>::run(x), x_le_five);
   1461     x_le_five = pselect(pcmp_le(x, pset1<T>(0.0)), NEG_MAXNUM, x_le_five);
   1462     T s = pdiv(pset1<T>(25.0), z);
   1463     T p = pdiv(
   1464         internal::ppolevl<T, 6>::run(s, PP),
   1465         internal::ppolevl<T, 6>::run(s, PQ));
   1466     T q = pdiv(
   1467         internal::ppolevl<T, 7>::run(s, QP),
   1468         internal::ppolevl<T, 7>::run(s, QQ));
   1469     T xn = padd(x, NEG_PIO4);
   1470     T w = pdiv(pset1<T>(5.0), x);
   1471     p = pmadd(p, psin(xn), pmul(w, pmul(q, pcos(xn))));
   1472     T x_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(x)));
   1473     return pselect(pcmp_le(x, pset1<T>(5.0)), x_le_five, x_gt_five);
   1474   }
   1475 };
   1476 
   1477 template <typename T>
   1478 struct bessel_y0_impl {
   1479   EIGEN_DEVICE_FUNC
   1480   static EIGEN_STRONG_INLINE T run(const T x) {
   1481     return generic_y0<T>::run(x);
   1482   }
   1483 };
   1484 
   1485 template <typename T>
   1486 struct bessel_j1_retval {
   1487   typedef T type;
   1488 };
   1489 
   1490 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
   1491 struct generic_j1 {
   1492   EIGEN_DEVICE_FUNC
   1493   static EIGEN_STRONG_INLINE T run(const T&) {
   1494     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
   1495                         THIS_TYPE_IS_NOT_SUPPORTED);
   1496     return ScalarType(0);
   1497   }
   1498 };
   1499 
   1500 template <typename T>
   1501 struct generic_j1<T, float> {
   1502   EIGEN_DEVICE_FUNC
   1503   static EIGEN_STRONG_INLINE T run(const T& x) {
   1504     /* j1f.c
   1505      *	Bessel function of order one
   1506      *
   1507      *
   1508      *
   1509      * SYNOPSIS:
   1510      *
   1511      * float x, y, j1f();
   1512      *
   1513      * y = j1f( x );
   1514      *
   1515      *
   1516      *
   1517      * DESCRIPTION:
   1518      *
   1519      * Returns Bessel function of order one of the argument.
   1520      *
   1521      * The domain is divided into the intervals [0, 2] and
   1522      * (2, infinity). In the first interval a polynomial approximation
   1523      *        2
   1524      * (w - r  ) x P(w)
   1525      *       1
   1526      *                     2
   1527      * is used, where w = x  and r is the first zero of the function.
   1528      *
   1529      * In the second interval, the modulus and phase are approximated
   1530      * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
   1531      * and Phase(x) = x + 1/x R(1/x^2) - 3pi/4.  The function is
   1532      *
   1533      *   j0(x) = Modulus(x) cos( Phase(x) ).
   1534      *
   1535      *
   1536      *
   1537      * ACCURACY:
   1538      *
   1539      *                      Absolute error:
   1540      * arithmetic   domain      # trials      peak       rms
   1541      *    IEEE      0,  2       100000       1.2e-7     2.5e-8
   1542      *    IEEE      2, 32       100000       2.0e-7     5.3e-8
   1543      *
   1544      *
   1545      */
   1546 
   1547     const float JP[] = {-4.878788132172128E-009f, 6.009061827883699E-007f,
   1548                         -4.541343896997497E-005f, 1.937383947804541E-003f,
   1549                         -3.405537384615824E-002f};
   1550     const float MO1[] = {6.913942741265801E-002f, -2.284801500053359E-001f,
   1551                         3.138238455499697E-001f, -2.102302420403875E-001f,
   1552                         5.435364690523026E-003f, 1.493389585089498E-001f,
   1553                         4.976029650847191E-006f, 7.978845453073848E-001f};
   1554     const float PH1[] = {-4.497014141919556E+001f, 5.073465654089319E+001f,
   1555                         -2.485774108720340E+001f, 7.222973196770240E+000f,
   1556                         -1.544842782180211E+000f, 3.503787691653334E-001f,
   1557                         -1.637986776941202E-001f, 3.749989509080821E-001f};
   1558     const T Z1 = pset1<T>(1.46819706421238932572E1f);
   1559     const T NEG_THPIO4F = pset1<T>(-2.35619449019234492885f);    /* -3*pi/4 */
   1560 
   1561     T y = pabs(x);
   1562     T z = pmul(y, y);
   1563     T y_le_two = pmul(
   1564         psub(z, Z1),
   1565         pmul(x, internal::ppolevl<T, 4>::run(z, JP)));
   1566     T q = pdiv(pset1<T>(1.0f), y);
   1567     T w = prsqrt(y);
   1568     T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO1));
   1569     w = pmul(q, q);
   1570     T yn = pmadd(q, internal::ppolevl<T, 7>::run(w, PH1), NEG_THPIO4F);
   1571     T y_gt_two = pmul(p, pcos(padd(yn, y)));
   1572     // j1 is an odd function. This implementation differs from cephes to
   1573     // take this fact in to account. Cephes returns -j1(x) for y > 2 range.
   1574     y_gt_two = pselect(
   1575         pcmp_lt(x, pset1<T>(0.0f)), pnegate(y_gt_two), y_gt_two);
   1576     return pselect(pcmp_le(y, pset1<T>(2.0f)), y_le_two, y_gt_two);
   1577   }
   1578 };
   1579 
   1580 template <typename T>
   1581 struct generic_j1<T, double> {
   1582   EIGEN_DEVICE_FUNC
   1583   static EIGEN_STRONG_INLINE T run(const T& x) {
   1584     /*  j1.c
   1585      *	Bessel function of order one
   1586      *
   1587      *
   1588      *
   1589      * SYNOPSIS:
   1590      *
   1591      * double x, y, j1();
   1592      *
   1593      * y = j1( x );
   1594      *
   1595      *
   1596      *
   1597      * DESCRIPTION:
   1598      *
   1599      * Returns Bessel function of order one of the argument.
   1600      *
   1601      * The domain is divided into the intervals [0, 8] and
   1602      * (8, infinity). In the first interval a 24 term Chebyshev
   1603      * expansion is used. In the second, the asymptotic
   1604      * trigonometric representation is employed using two
   1605      * rational functions of degree 5/5.
   1606      *
   1607      *
   1608      *
   1609      * ACCURACY:
   1610      *
   1611      *                      Absolute error:
   1612      * arithmetic   domain      # trials      peak         rms
   1613      *    DEC       0, 30       10000       4.0e-17     1.1e-17
   1614      *    IEEE      0, 30       30000       2.6e-16     1.1e-16
   1615      *
   1616      */
   1617     const double PP[] = {7.62125616208173112003E-4, 7.31397056940917570436E-2,
   1618                          1.12719608129684925192E0, 5.11207951146807644818E0,
   1619                          8.42404590141772420927E0, 5.21451598682361504063E0,
   1620                          1.00000000000000000254E0};
   1621     const double PQ[] = {5.71323128072548699714E-4, 6.88455908754495404082E-2,
   1622                          1.10514232634061696926E0, 5.07386386128601488557E0,
   1623                          8.39985554327604159757E0, 5.20982848682361821619E0,
   1624                          9.99999999999999997461E-1};
   1625     const double QP[] = {5.10862594750176621635E-2, 4.98213872951233449420E0,
   1626                          7.58238284132545283818E1, 3.66779609360150777800E2,
   1627                          7.10856304998926107277E2, 5.97489612400613639965E2,
   1628                          2.11688757100572135698E2, 2.52070205858023719784E1};
   1629     const double QQ[] = {1.00000000000000000000E0, 7.42373277035675149943E1,
   1630                          1.05644886038262816351E3, 4.98641058337653607651E3,
   1631                          9.56231892404756170795E3, 7.99704160447350683650E3,
   1632                          2.82619278517639096600E3, 3.36093607810698293419E2};
   1633     const double RP[] = {-8.99971225705559398224E8, 4.52228297998194034323E11,
   1634                          -7.27494245221818276015E13, 3.68295732863852883286E15};
   1635     const double RQ[] = {1.00000000000000000000E0, 6.20836478118054335476E2,
   1636                          2.56987256757748830383E5, 8.35146791431949253037E7,
   1637                          2.21511595479792499675E10, 4.74914122079991414898E12,
   1638                          7.84369607876235854894E14, 8.95222336184627338078E16,
   1639                          5.32278620332680085395E18};
   1640     const T Z1 = pset1<T>(1.46819706421238932572E1);
   1641     const T Z2 = pset1<T>(4.92184563216946036703E1);
   1642     const T NEG_THPIO4 = pset1<T>(-2.35619449019234492885);    /* -3*pi/4 */
   1643     const T SQ2OPI = pset1<T>(7.9788456080286535587989E-1); /* sqrt(2 / pi) */
   1644     T y = pabs(x);
   1645     T z = pmul(y, y);
   1646     T y_le_five = pdiv(internal::ppolevl<T, 3>::run(z, RP),
   1647                        internal::ppolevl<T, 8>::run(z, RQ));
   1648     y_le_five = pmul(pmul(pmul(y_le_five, x), psub(z, Z1)), psub(z, Z2));
   1649     T s = pdiv(pset1<T>(25.0), z);
   1650     T p = pdiv(
   1651         internal::ppolevl<T, 6>::run(s, PP),
   1652         internal::ppolevl<T, 6>::run(s, PQ));
   1653     T q = pdiv(
   1654         internal::ppolevl<T, 7>::run(s, QP),
   1655         internal::ppolevl<T, 7>::run(s, QQ));
   1656     T yn = padd(y, NEG_THPIO4);
   1657     T w = pdiv(pset1<T>(-5.0), y);
   1658     p = pmadd(p, pcos(yn), pmul(w, pmul(q, psin(yn))));
   1659     T y_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(y)));
   1660     // j1 is an odd function. This implementation differs from cephes to
   1661     // take this fact in to account. Cephes returns -j1(x) for y > 5 range.
   1662     y_gt_five = pselect(
   1663         pcmp_lt(x, pset1<T>(0.0)), pnegate(y_gt_five), y_gt_five);
   1664     return pselect(pcmp_le(y, pset1<T>(5.0)), y_le_five, y_gt_five);
   1665   }
   1666 };
   1667 
   1668 template <typename T>
   1669 struct bessel_j1_impl {
   1670   EIGEN_DEVICE_FUNC
   1671   static EIGEN_STRONG_INLINE T run(const T x) {
   1672     return generic_j1<T>::run(x);
   1673   }
   1674 };
   1675 
   1676 template <typename T>
   1677 struct bessel_y1_retval {
   1678   typedef T type;
   1679 };
   1680 
   1681 template <typename T, typename ScalarType = typename unpacket_traits<T>::type>
   1682 struct generic_y1 {
   1683   EIGEN_DEVICE_FUNC
   1684   static EIGEN_STRONG_INLINE T run(const T&) {
   1685     EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
   1686                         THIS_TYPE_IS_NOT_SUPPORTED);
   1687     return ScalarType(0);
   1688   }
   1689 };
   1690 
   1691 template <typename T>
   1692 struct generic_y1<T, float> {
   1693   EIGEN_DEVICE_FUNC
   1694   static EIGEN_STRONG_INLINE T run(const T& x) {
   1695     /* j1f.c
   1696      *	Bessel function of second kind of order one
   1697      *
   1698      *
   1699      *
   1700      * SYNOPSIS:
   1701      *
   1702      * double x, y, y1();
   1703      *
   1704      * y = y1( x );
   1705      *
   1706      *
   1707      *
   1708      * DESCRIPTION:
   1709      *
   1710      * Returns Bessel function of the second kind of order one
   1711      * of the argument.
   1712      *
   1713      * The domain is divided into the intervals [0, 2] and
   1714      * (2, infinity). In the first interval a rational approximation
   1715      * R(x) is employed to compute
   1716      *
   1717      *                  2
   1718      * y0(x)  =  (w - r  ) x R(x^2)  +  2/pi (ln(x) j1(x) - 1/x) .
   1719      *                 1
   1720      *
   1721      * Thus a call to j1() is required.
   1722      *
   1723      * In the second interval, the modulus and phase are approximated
   1724      * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
   1725      * and Phase(x) = x + 1/x S(1/x^2) - 3pi/4.  Then the function is
   1726      *
   1727      *   y0(x) = Modulus(x) sin( Phase(x) ).
   1728      *
   1729      *
   1730      *
   1731      *
   1732      * ACCURACY:
   1733      *
   1734      *                      Absolute error:
   1735      * arithmetic   domain      # trials      peak         rms
   1736      *    IEEE      0,  2       100000       2.2e-7     4.6e-8
   1737      *    IEEE      2, 32       100000       1.9e-7     5.3e-8
   1738      *
   1739      * (error criterion relative when |y1| > 1).
   1740      *
   1741      */
   1742 
   1743     const float YP[] = {8.061978323326852E-009f, -9.496460629917016E-007f,
   1744                         6.719543806674249E-005f, -2.641785726447862E-003f,
   1745                         4.202369946500099E-002f};
   1746     const float MO1[] = {6.913942741265801E-002f, -2.284801500053359E-001f,
   1747                         3.138238455499697E-001f, -2.102302420403875E-001f,
   1748                         5.435364690523026E-003f, 1.493389585089498E-001f,
   1749                         4.976029650847191E-006f, 7.978845453073848E-001f};
   1750     const float PH1[] = {-4.497014141919556E+001f, 5.073465654089319E+001f,
   1751                         -2.485774108720340E+001f, 7.222973196770240E+000f,
   1752                         -1.544842782180211E+000f, 3.503787691653334E-001f,
   1753                         -1.637986776941202E-001f, 3.749989509080821E-001f};
   1754     const T YO1 = pset1<T>(4.66539330185668857532f);
   1755     const T NEG_THPIO4F = pset1<T>(-2.35619449019234492885f);    /* -3*pi/4 */
   1756     const T TWOOPI = pset1<T>(0.636619772367581343075535f); /* 2/pi */
   1757     const T NEG_MAXNUM = pset1<T>(-NumTraits<float>::infinity());
   1758 
   1759     T z = pmul(x, x);
   1760     T x_le_two = pmul(psub(z, YO1), internal::ppolevl<T, 4>::run(z, YP));
   1761     x_le_two = pmadd(
   1762        x_le_two, x,
   1763        pmul(TWOOPI, pmadd(
   1764            generic_j1<T, float>::run(x), plog(x),
   1765            pdiv(pset1<T>(-1.0f), x))));
   1766     x_le_two = pselect(pcmp_lt(x, pset1<T>(0.0f)), NEG_MAXNUM, x_le_two);
   1767 
   1768     T q = pdiv(pset1<T>(1.0), x);
   1769     T w = prsqrt(x);
   1770     T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO1));
   1771     w = pmul(q, q);
   1772     T xn = pmadd(q, internal::ppolevl<T, 7>::run(w, PH1), NEG_THPIO4F);
   1773     T x_gt_two = pmul(p, psin(padd(xn, x)));
   1774     return pselect(pcmp_le(x, pset1<T>(2.0)), x_le_two, x_gt_two);
   1775   }
   1776 };
   1777 
   1778 template <typename T>
   1779 struct generic_y1<T, double> {
   1780   EIGEN_DEVICE_FUNC
   1781   static EIGEN_STRONG_INLINE T run(const T& x) {
   1782     /*  j1.c
   1783      *	Bessel function of second kind of order one
   1784      *
   1785      *
   1786      *
   1787      * SYNOPSIS:
   1788      *
   1789      * double x, y, y1();
   1790      *
   1791      * y = y1( x );
   1792      *
   1793      *
   1794      *
   1795      * DESCRIPTION:
   1796      *
   1797      * Returns Bessel function of the second kind of order one
   1798      * of the argument.
   1799      *
   1800      * The domain is divided into the intervals [0, 8] and
   1801      * (8, infinity). In the first interval a 25 term Chebyshev
   1802      * expansion is used, and a call to j1() is required.
   1803      * In the second, the asymptotic trigonometric representation
   1804      * is employed using two rational functions of degree 5/5.
   1805      *
   1806      *
   1807      *
   1808      * ACCURACY:
   1809      *
   1810      *                      Absolute error:
   1811      * arithmetic   domain      # trials      peak         rms
   1812      *    DEC       0, 30       10000       8.6e-17     1.3e-17
   1813      *    IEEE      0, 30       30000       1.0e-15     1.3e-16
   1814      *
   1815      * (error criterion relative when |y1| > 1).
   1816      *
   1817      */
   1818     const double PP[] = {7.62125616208173112003E-4, 7.31397056940917570436E-2,
   1819                          1.12719608129684925192E0, 5.11207951146807644818E0,
   1820                          8.42404590141772420927E0, 5.21451598682361504063E0,
   1821                          1.00000000000000000254E0};
   1822     const double PQ[] = {5.71323128072548699714E-4, 6.88455908754495404082E-2,
   1823                          1.10514232634061696926E0, 5.07386386128601488557E0,
   1824                          8.39985554327604159757E0, 5.20982848682361821619E0,
   1825                          9.99999999999999997461E-1};
   1826     const double QP[] = {5.10862594750176621635E-2, 4.98213872951233449420E0,
   1827                          7.58238284132545283818E1, 3.66779609360150777800E2,
   1828                          7.10856304998926107277E2, 5.97489612400613639965E2,
   1829                          2.11688757100572135698E2, 2.52070205858023719784E1};
   1830     const double QQ[] = {1.00000000000000000000E0, 7.42373277035675149943E1,
   1831                          1.05644886038262816351E3, 4.98641058337653607651E3,
   1832                          9.56231892404756170795E3, 7.99704160447350683650E3,
   1833                          2.82619278517639096600E3, 3.36093607810698293419E2};
   1834     const double YP[] = {1.26320474790178026440E9, -6.47355876379160291031E11,
   1835                          1.14509511541823727583E14, -8.12770255501325109621E15,
   1836                          2.02439475713594898196E17, -7.78877196265950026825E17};
   1837     const double YQ[] = {1.00000000000000000000E0, 5.94301592346128195359E2,
   1838                          2.35564092943068577943E5, 7.34811944459721705660E7,
   1839                          1.87601316108706159478E10, 3.88231277496238566008E12,
   1840                          6.20557727146953693363E14, 6.87141087355300489866E16,
   1841                          3.97270608116560655612E18};
   1842     const T SQ2OPI = pset1<T>(.79788456080286535588);
   1843     const T NEG_THPIO4 = pset1<T>(-2.35619449019234492885);    /* -3*pi/4 */
   1844     const T TWOOPI = pset1<T>(0.636619772367581343075535); /* 2/pi */
   1845     const T NEG_MAXNUM = pset1<T>(-NumTraits<double>::infinity());
   1846 
   1847     T z = pmul(x, x);
   1848     T x_le_five = pdiv(internal::ppolevl<T, 5>::run(z, YP),
   1849                    internal::ppolevl<T, 8>::run(z, YQ));
   1850     x_le_five = pmadd(
   1851         x_le_five, x, pmul(
   1852             TWOOPI, pmadd(generic_j1<T, double>::run(x), plog(x),
   1853                           pdiv(pset1<T>(-1.0), x))));
   1854 
   1855     x_le_five = pselect(pcmp_le(x, pset1<T>(0.0)), NEG_MAXNUM, x_le_five);
   1856     T s = pdiv(pset1<T>(25.0), z);
   1857     T p = pdiv(
   1858         internal::ppolevl<T, 6>::run(s, PP),
   1859         internal::ppolevl<T, 6>::run(s, PQ));
   1860     T q = pdiv(
   1861         internal::ppolevl<T, 7>::run(s, QP),
   1862         internal::ppolevl<T, 7>::run(s, QQ));
   1863     T xn = padd(x, NEG_THPIO4);
   1864     T w = pdiv(pset1<T>(5.0), x);
   1865     p = pmadd(p, psin(xn), pmul(w, pmul(q, pcos(xn))));
   1866     T x_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(x)));
   1867     return pselect(pcmp_le(x, pset1<T>(5.0)), x_le_five, x_gt_five);
   1868   }
   1869 };
   1870 
   1871 template <typename T>
   1872 struct bessel_y1_impl {
   1873   EIGEN_DEVICE_FUNC
   1874   static EIGEN_STRONG_INLINE T run(const T x) {
   1875     return generic_y1<T>::run(x);
   1876   }
   1877 };
   1878 
   1879 }  // end namespace internal
   1880 
   1881 namespace numext {
   1882 
   1883 template <typename Scalar>
   1884 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_i0, Scalar)
   1885     bessel_i0(const Scalar& x) {
   1886   return EIGEN_MATHFUNC_IMPL(bessel_i0, Scalar)::run(x);
   1887 }
   1888 
   1889 template <typename Scalar>
   1890 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_i0e, Scalar)
   1891     bessel_i0e(const Scalar& x) {
   1892   return EIGEN_MATHFUNC_IMPL(bessel_i0e, Scalar)::run(x);
   1893 }
   1894 
   1895 template <typename Scalar>
   1896 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_i1, Scalar)
   1897     bessel_i1(const Scalar& x) {
   1898   return EIGEN_MATHFUNC_IMPL(bessel_i1, Scalar)::run(x);
   1899 }
   1900 
   1901 template <typename Scalar>
   1902 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_i1e, Scalar)
   1903     bessel_i1e(const Scalar& x) {
   1904   return EIGEN_MATHFUNC_IMPL(bessel_i1e, Scalar)::run(x);
   1905 }
   1906 
   1907 template <typename Scalar>
   1908 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_k0, Scalar)
   1909     bessel_k0(const Scalar& x) {
   1910   return EIGEN_MATHFUNC_IMPL(bessel_k0, Scalar)::run(x);
   1911 }
   1912 
   1913 template <typename Scalar>
   1914 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_k0e, Scalar)
   1915     bessel_k0e(const Scalar& x) {
   1916   return EIGEN_MATHFUNC_IMPL(bessel_k0e, Scalar)::run(x);
   1917 }
   1918 
   1919 template <typename Scalar>
   1920 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_k1, Scalar)
   1921     bessel_k1(const Scalar& x) {
   1922   return EIGEN_MATHFUNC_IMPL(bessel_k1, Scalar)::run(x);
   1923 }
   1924 
   1925 template <typename Scalar>
   1926 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_k1e, Scalar)
   1927     bessel_k1e(const Scalar& x) {
   1928   return EIGEN_MATHFUNC_IMPL(bessel_k1e, Scalar)::run(x);
   1929 }
   1930 
   1931 template <typename Scalar>
   1932 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_j0, Scalar)
   1933     bessel_j0(const Scalar& x) {
   1934   return EIGEN_MATHFUNC_IMPL(bessel_j0, Scalar)::run(x);
   1935 }
   1936 
   1937 template <typename Scalar>
   1938 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_y0, Scalar)
   1939     bessel_y0(const Scalar& x) {
   1940   return EIGEN_MATHFUNC_IMPL(bessel_y0, Scalar)::run(x);
   1941 }
   1942 
   1943 template <typename Scalar>
   1944 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_j1, Scalar)
   1945     bessel_j1(const Scalar& x) {
   1946   return EIGEN_MATHFUNC_IMPL(bessel_j1, Scalar)::run(x);
   1947 }
   1948 
   1949 template <typename Scalar>
   1950 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(bessel_y1, Scalar)
   1951     bessel_y1(const Scalar& x) {
   1952   return EIGEN_MATHFUNC_IMPL(bessel_y1, Scalar)::run(x);
   1953 }
   1954 
   1955 }  // end namespace numext
   1956 
   1957 }  // end namespace Eigen
   1958 
   1959 #endif  // EIGEN_BESSEL_FUNCTIONS_H
	cart-elc Source code for CART-ELC
	git clone git://git.laack.co/cart-elc.git
	Log \| Files \| Refs \| README \| LICENSE