cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

BesselFunctionsFunctors.h (12641B)


      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2016 Eugene Brevdo <ebrevdo@gmail.com>
      5 // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
      6 //
      7 // This Source Code Form is subject to the terms of the Mozilla
      8 // Public License v. 2.0. If a copy of the MPL was not distributed
      9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
     10 
     11 #ifndef EIGEN_BESSELFUNCTIONS_FUNCTORS_H
     12 #define EIGEN_BESSELFUNCTIONS_FUNCTORS_H
     13 
     14 namespace Eigen {
     15 
     16 namespace internal {
     17 
     18 /** \internal
     19  * \brief Template functor to compute the modified Bessel function of the first
     20  * kind of order zero.
     21  * \sa class CwiseUnaryOp, Cwise::bessel_i0()
     22  */
     23 template <typename Scalar>
     24 struct scalar_bessel_i0_op {
     25   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i0_op)
     26   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
     27     using numext::bessel_i0;
     28     return bessel_i0(x);
     29   }
     30   typedef typename packet_traits<Scalar>::type Packet;
     31   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
     32     return internal::pbessel_i0(x);
     33   }
     34 };
     35 template <typename Scalar>
     36 struct functor_traits<scalar_bessel_i0_op<Scalar> > {
     37   enum {
     38     // On average, a Chebyshev polynomial of order N=20 is computed.
     39     // The cost is N multiplications and 2N additions. We also add
     40     // the cost of an additional exp over i0e.
     41     Cost = 28 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
     42     PacketAccess = packet_traits<Scalar>::HasBessel
     43   };
     44 };
     45 
     46 /** \internal
     47  * \brief Template functor to compute the exponentially scaled modified Bessel
     48  * function of the first kind of order zero
     49  * \sa class CwiseUnaryOp, Cwise::bessel_i0e()
     50  */
     51 template <typename Scalar>
     52 struct scalar_bessel_i0e_op {
     53   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i0e_op)
     54   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
     55     using numext::bessel_i0e;
     56     return bessel_i0e(x);
     57   }
     58   typedef typename packet_traits<Scalar>::type Packet;
     59   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
     60     return internal::pbessel_i0e(x);
     61   }
     62 };
     63 template <typename Scalar>
     64 struct functor_traits<scalar_bessel_i0e_op<Scalar> > {
     65   enum {
     66     // On average, a Chebyshev polynomial of order N=20 is computed.
     67     // The cost is N multiplications and 2N additions.
     68     Cost = 20 * NumTraits<Scalar>::MulCost + 40 * NumTraits<Scalar>::AddCost,
     69     PacketAccess = packet_traits<Scalar>::HasBessel
     70   };
     71 };
     72 
     73 /** \internal
     74  * \brief Template functor to compute the modified Bessel function of the first
     75  * kind of order one
     76  * \sa class CwiseUnaryOp, Cwise::bessel_i1()
     77  */
     78 template <typename Scalar>
     79 struct scalar_bessel_i1_op {
     80   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i1_op)
     81   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
     82     using numext::bessel_i1;
     83     return bessel_i1(x);
     84   }
     85   typedef typename packet_traits<Scalar>::type Packet;
     86   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
     87     return internal::pbessel_i1(x);
     88   }
     89 };
     90 template <typename Scalar>
     91 struct functor_traits<scalar_bessel_i1_op<Scalar> > {
     92   enum {
     93     // On average, a Chebyshev polynomial of order N=20 is computed.
     94     // The cost is N multiplications and 2N additions. We also add
     95     // the cost of an additional exp over i1e.
     96     Cost = 28 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
     97     PacketAccess = packet_traits<Scalar>::HasBessel
     98   };
     99 };
    100 
    101 /** \internal
    102  * \brief Template functor to compute the exponentially scaled modified Bessel
    103  * function of the first kind of order zero
    104  * \sa class CwiseUnaryOp, Cwise::bessel_i1e()
    105  */
    106 template <typename Scalar>
    107 struct scalar_bessel_i1e_op {
    108   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i1e_op)
    109   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    110     using numext::bessel_i1e;
    111     return bessel_i1e(x);
    112   }
    113   typedef typename packet_traits<Scalar>::type Packet;
    114   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    115     return internal::pbessel_i1e(x);
    116   }
    117 };
    118 template <typename Scalar>
    119 struct functor_traits<scalar_bessel_i1e_op<Scalar> > {
    120   enum {
    121     // On average, a Chebyshev polynomial of order N=20 is computed.
    122     // The cost is N multiplications and 2N additions.
    123     Cost = 20 * NumTraits<Scalar>::MulCost + 40 * NumTraits<Scalar>::AddCost,
    124     PacketAccess = packet_traits<Scalar>::HasBessel
    125   };
    126 };
    127 
    128 /** \internal
    129  * \brief Template functor to compute the Bessel function of the second kind of
    130  * order zero
    131  * \sa class CwiseUnaryOp, Cwise::bessel_j0()
    132  */
    133 template <typename Scalar>
    134 struct scalar_bessel_j0_op {
    135   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_j0_op)
    136   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    137     using numext::bessel_j0;
    138     return bessel_j0(x);
    139   }
    140   typedef typename packet_traits<Scalar>::type Packet;
    141   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    142     return internal::pbessel_j0(x);
    143   }
    144 };
    145 template <typename Scalar>
    146 struct functor_traits<scalar_bessel_j0_op<Scalar> > {
    147   enum {
    148     // 6 polynomial of order ~N=8 is computed.
    149     // The cost is N multiplications and N additions each, along with a
    150     // sine, cosine and rsqrt cost.
    151     Cost = 63 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
    152     PacketAccess = packet_traits<Scalar>::HasBessel
    153   };
    154 };
    155 
    156 /** \internal
    157  * \brief Template functor to compute the Bessel function of the second kind of
    158  * order zero
    159  * \sa class CwiseUnaryOp, Cwise::bessel_y0()
    160  */
    161 template <typename Scalar>
    162 struct scalar_bessel_y0_op {
    163   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_y0_op)
    164   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    165     using numext::bessel_y0;
    166     return bessel_y0(x);
    167   }
    168   typedef typename packet_traits<Scalar>::type Packet;
    169   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    170     return internal::pbessel_y0(x);
    171   }
    172 };
    173 template <typename Scalar>
    174 struct functor_traits<scalar_bessel_y0_op<Scalar> > {
    175   enum {
    176     // 6 polynomial of order ~N=8 is computed.
    177     // The cost is N multiplications and N additions each, along with a
    178     // sine, cosine, rsqrt and j0 cost.
    179     Cost = 126 * NumTraits<Scalar>::MulCost + 96 * NumTraits<Scalar>::AddCost,
    180     PacketAccess = packet_traits<Scalar>::HasBessel
    181   };
    182 };
    183 
    184 /** \internal
    185  * \brief Template functor to compute the Bessel function of the first kind of
    186  * order one
    187  * \sa class CwiseUnaryOp, Cwise::bessel_j1()
    188  */
    189 template <typename Scalar>
    190 struct scalar_bessel_j1_op {
    191   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_j1_op)
    192   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    193     using numext::bessel_j1;
    194     return bessel_j1(x);
    195   }
    196   typedef typename packet_traits<Scalar>::type Packet;
    197   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    198     return internal::pbessel_j1(x);
    199   }
    200 };
    201 template <typename Scalar>
    202 struct functor_traits<scalar_bessel_j1_op<Scalar> > {
    203   enum {
    204     // 6 polynomial of order ~N=8 is computed.
    205     // The cost is N multiplications and N additions each, along with a
    206     // sine, cosine and rsqrt cost.
    207     Cost = 63 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
    208     PacketAccess = packet_traits<Scalar>::HasBessel
    209   };
    210 };
    211 
    212 /** \internal
    213  * \brief Template functor to compute the Bessel function of the second kind of
    214  * order one
    215  * \sa class CwiseUnaryOp, Cwise::bessel_j1e()
    216  */
    217 template <typename Scalar>
    218 struct scalar_bessel_y1_op {
    219   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_y1_op)
    220   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    221     using numext::bessel_y1;
    222     return bessel_y1(x);
    223   }
    224   typedef typename packet_traits<Scalar>::type Packet;
    225   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    226     return internal::pbessel_y1(x);
    227   }
    228 };
    229 template <typename Scalar>
    230 struct functor_traits<scalar_bessel_y1_op<Scalar> > {
    231   enum {
    232     // 6 polynomial of order ~N=8 is computed.
    233     // The cost is N multiplications and N additions each, along with a
    234     // sine, cosine, rsqrt and j1 cost.
    235     Cost = 126 * NumTraits<Scalar>::MulCost + 96 * NumTraits<Scalar>::AddCost,
    236     PacketAccess = packet_traits<Scalar>::HasBessel
    237   };
    238 };
    239 
    240 /** \internal
    241  * \brief Template functor to compute the modified Bessel function of the second
    242  * kind of order zero
    243  * \sa class CwiseUnaryOp, Cwise::bessel_k0()
    244  */
    245 template <typename Scalar>
    246 struct scalar_bessel_k0_op {
    247   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k0_op)
    248   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    249     using numext::bessel_k0;
    250     return bessel_k0(x);
    251   }
    252   typedef typename packet_traits<Scalar>::type Packet;
    253   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    254     return internal::pbessel_k0(x);
    255   }
    256 };
    257 template <typename Scalar>
    258 struct functor_traits<scalar_bessel_k0_op<Scalar> > {
    259   enum {
    260     // On average, a Chebyshev polynomial of order N=10 is computed.
    261     // The cost is N multiplications and 2N additions. In addition we compute
    262     // i0, a log, exp and prsqrt and sin and cos.
    263     Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
    264     PacketAccess = packet_traits<Scalar>::HasBessel
    265   };
    266 };
    267 
    268 /** \internal
    269  * \brief Template functor to compute the exponentially scaled modified Bessel
    270  * function of the second kind of order zero
    271  * \sa class CwiseUnaryOp, Cwise::bessel_k0e()
    272  */
    273 template <typename Scalar>
    274 struct scalar_bessel_k0e_op {
    275   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k0e_op)
    276   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    277     using numext::bessel_k0e;
    278     return bessel_k0e(x);
    279   }
    280   typedef typename packet_traits<Scalar>::type Packet;
    281   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    282     return internal::pbessel_k0e(x);
    283   }
    284 };
    285 template <typename Scalar>
    286 struct functor_traits<scalar_bessel_k0e_op<Scalar> > {
    287   enum {
    288     // On average, a Chebyshev polynomial of order N=10 is computed.
    289     // The cost is N multiplications and 2N additions. In addition we compute
    290     // i0, a log, exp and prsqrt and sin and cos.
    291     Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
    292     PacketAccess = packet_traits<Scalar>::HasBessel
    293   };
    294 };
    295 
    296 /** \internal
    297  * \brief Template functor to compute the modified Bessel function of the
    298  * second kind of order one
    299  * \sa class CwiseUnaryOp, Cwise::bessel_k1()
    300  */
    301 template <typename Scalar>
    302 struct scalar_bessel_k1_op {
    303   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k1_op)
    304   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    305     using numext::bessel_k1;
    306     return bessel_k1(x);
    307   }
    308   typedef typename packet_traits<Scalar>::type Packet;
    309   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    310     return internal::pbessel_k1(x);
    311   }
    312 };
    313 template <typename Scalar>
    314 struct functor_traits<scalar_bessel_k1_op<Scalar> > {
    315   enum {
    316     // On average, a Chebyshev polynomial of order N=10 is computed.
    317     // The cost is N multiplications and 2N additions. In addition we compute
    318     // i1, a log, exp and prsqrt and sin and cos.
    319     Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
    320     PacketAccess = packet_traits<Scalar>::HasBessel
    321   };
    322 };
    323 
    324 /** \internal
    325  * \brief Template functor to compute the exponentially scaled modified Bessel
    326  * function of the second kind of order one
    327  * \sa class CwiseUnaryOp, Cwise::bessel_k1e()
    328  */
    329 template <typename Scalar>
    330 struct scalar_bessel_k1e_op {
    331   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k1e_op)
    332   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
    333     using numext::bessel_k1e;
    334     return bessel_k1e(x);
    335   }
    336   typedef typename packet_traits<Scalar>::type Packet;
    337   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    338     return internal::pbessel_k1e(x);
    339   }
    340 };
    341 template <typename Scalar>
    342 struct functor_traits<scalar_bessel_k1e_op<Scalar> > {
    343   enum {
    344     // On average, a Chebyshev polynomial of order N=10 is computed.
    345     // The cost is N multiplications and 2N additions. In addition we compute
    346     // i1, a log, exp and prsqrt and sin and cos.
    347     Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
    348     PacketAccess = packet_traits<Scalar>::HasBessel
    349   };
    350 };
    351 
    352 
    353 } // end namespace internal
    354 
    355 } // end namespace Eigen
    356 
    357 #endif // EIGEN_BESSELFUNCTIONS_FUNCTORS_H