src/SpecialFunctions/BesselFunctionsFunctors.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2016 Eugene Brevdo <ebrevdo@gmail.com>
0005 // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_BESSELFUNCTIONS_FUNCTORS_H
0012 #define EIGEN_BESSELFUNCTIONS_FUNCTORS_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017
0018 /** \internal
0019  * \brief Template functor to compute the modified Bessel function of the first
0020  * kind of order zero.
0021  * \sa class CwiseUnaryOp, Cwise::bessel_i0()
0022  */
0023 template <typename Scalar>
0024 struct scalar_bessel_i0_op {
0025   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i0_op)
0026   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0027     using numext::bessel_i0;
0028     return bessel_i0(x);
0029   }
0030   typedef typename packet_traits<Scalar>::type Packet;
0031   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0032     return internal::pbessel_i0(x);
0033   }
0034 };
0035 template <typename Scalar>
0036 struct functor_traits<scalar_bessel_i0_op<Scalar> > {
0037   enum {
0038     // On average, a Chebyshev polynomial of order N=20 is computed.
0039     // The cost is N multiplications and 2N additions. We also add
0040     // the cost of an additional exp over i0e.
0041     Cost = 28 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
0042     PacketAccess = packet_traits<Scalar>::HasBessel
0043   };
0044 };
0045
0046 /** \internal
0047  * \brief Template functor to compute the exponentially scaled modified Bessel
0048  * function of the first kind of order zero
0049  * \sa class CwiseUnaryOp, Cwise::bessel_i0e()
0050  */
0051 template <typename Scalar>
0052 struct scalar_bessel_i0e_op {
0053   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i0e_op)
0054   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0055     using numext::bessel_i0e;
0056     return bessel_i0e(x);
0057   }
0058   typedef typename packet_traits<Scalar>::type Packet;
0059   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0060     return internal::pbessel_i0e(x);
0061   }
0062 };
0063 template <typename Scalar>
0064 struct functor_traits<scalar_bessel_i0e_op<Scalar> > {
0065   enum {
0066     // On average, a Chebyshev polynomial of order N=20 is computed.
0067     // The cost is N multiplications and 2N additions.
0068     Cost = 20 * NumTraits<Scalar>::MulCost + 40 * NumTraits<Scalar>::AddCost,
0069     PacketAccess = packet_traits<Scalar>::HasBessel
0070   };
0071 };
0072
0073 /** \internal
0074  * \brief Template functor to compute the modified Bessel function of the first
0075  * kind of order one
0076  * \sa class CwiseUnaryOp, Cwise::bessel_i1()
0077  */
0078 template <typename Scalar>
0079 struct scalar_bessel_i1_op {
0080   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i1_op)
0081   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0082     using numext::bessel_i1;
0083     return bessel_i1(x);
0084   }
0085   typedef typename packet_traits<Scalar>::type Packet;
0086   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0087     return internal::pbessel_i1(x);
0088   }
0089 };
0090 template <typename Scalar>
0091 struct functor_traits<scalar_bessel_i1_op<Scalar> > {
0092   enum {
0093     // On average, a Chebyshev polynomial of order N=20 is computed.
0094     // The cost is N multiplications and 2N additions. We also add
0095     // the cost of an additional exp over i1e.
0096     Cost = 28 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
0097     PacketAccess = packet_traits<Scalar>::HasBessel
0098   };
0099 };
0100
0101 /** \internal
0102  * \brief Template functor to compute the exponentially scaled modified Bessel
0103  * function of the first kind of order zero
0104  * \sa class CwiseUnaryOp, Cwise::bessel_i1e()
0105  */
0106 template <typename Scalar>
0107 struct scalar_bessel_i1e_op {
0108   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i1e_op)
0109   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0110     using numext::bessel_i1e;
0111     return bessel_i1e(x);
0112   }
0113   typedef typename packet_traits<Scalar>::type Packet;
0114   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0115     return internal::pbessel_i1e(x);
0116   }
0117 };
0118 template <typename Scalar>
0119 struct functor_traits<scalar_bessel_i1e_op<Scalar> > {
0120   enum {
0121     // On average, a Chebyshev polynomial of order N=20 is computed.
0122     // The cost is N multiplications and 2N additions.
0123     Cost = 20 * NumTraits<Scalar>::MulCost + 40 * NumTraits<Scalar>::AddCost,
0124     PacketAccess = packet_traits<Scalar>::HasBessel
0125   };
0126 };
0127
0128 /** \internal
0129  * \brief Template functor to compute the Bessel function of the second kind of
0130  * order zero
0131  * \sa class CwiseUnaryOp, Cwise::bessel_j0()
0132  */
0133 template <typename Scalar>
0134 struct scalar_bessel_j0_op {
0135   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_j0_op)
0136   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0137     using numext::bessel_j0;
0138     return bessel_j0(x);
0139   }
0140   typedef typename packet_traits<Scalar>::type Packet;
0141   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0142     return internal::pbessel_j0(x);
0143   }
0144 };
0145 template <typename Scalar>
0146 struct functor_traits<scalar_bessel_j0_op<Scalar> > {
0147   enum {
0148     // 6 polynomial of order ~N=8 is computed.
0149     // The cost is N multiplications and N additions each, along with a
0150     // sine, cosine and rsqrt cost.
0151     Cost = 63 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
0152     PacketAccess = packet_traits<Scalar>::HasBessel
0153   };
0154 };
0155
0156 /** \internal
0157  * \brief Template functor to compute the Bessel function of the second kind of
0158  * order zero
0159  * \sa class CwiseUnaryOp, Cwise::bessel_y0()
0160  */
0161 template <typename Scalar>
0162 struct scalar_bessel_y0_op {
0163   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_y0_op)
0164   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0165     using numext::bessel_y0;
0166     return bessel_y0(x);
0167   }
0168   typedef typename packet_traits<Scalar>::type Packet;
0169   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0170     return internal::pbessel_y0(x);
0171   }
0172 };
0173 template <typename Scalar>
0174 struct functor_traits<scalar_bessel_y0_op<Scalar> > {
0175   enum {
0176     // 6 polynomial of order ~N=8 is computed.
0177     // The cost is N multiplications and N additions each, along with a
0178     // sine, cosine, rsqrt and j0 cost.
0179     Cost = 126 * NumTraits<Scalar>::MulCost + 96 * NumTraits<Scalar>::AddCost,
0180     PacketAccess = packet_traits<Scalar>::HasBessel
0181   };
0182 };
0183
0184 /** \internal
0185  * \brief Template functor to compute the Bessel function of the first kind of
0186  * order one
0187  * \sa class CwiseUnaryOp, Cwise::bessel_j1()
0188  */
0189 template <typename Scalar>
0190 struct scalar_bessel_j1_op {
0191   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_j1_op)
0192   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0193     using numext::bessel_j1;
0194     return bessel_j1(x);
0195   }
0196   typedef typename packet_traits<Scalar>::type Packet;
0197   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0198     return internal::pbessel_j1(x);
0199   }
0200 };
0201 template <typename Scalar>
0202 struct functor_traits<scalar_bessel_j1_op<Scalar> > {
0203   enum {
0204     // 6 polynomial of order ~N=8 is computed.
0205     // The cost is N multiplications and N additions each, along with a
0206     // sine, cosine and rsqrt cost.
0207     Cost = 63 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
0208     PacketAccess = packet_traits<Scalar>::HasBessel
0209   };
0210 };
0211
0212 /** \internal
0213  * \brief Template functor to compute the Bessel function of the second kind of
0214  * order one
0215  * \sa class CwiseUnaryOp, Cwise::bessel_j1e()
0216  */
0217 template <typename Scalar>
0218 struct scalar_bessel_y1_op {
0219   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_y1_op)
0220   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0221     using numext::bessel_y1;
0222     return bessel_y1(x);
0223   }
0224   typedef typename packet_traits<Scalar>::type Packet;
0225   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0226     return internal::pbessel_y1(x);
0227   }
0228 };
0229 template <typename Scalar>
0230 struct functor_traits<scalar_bessel_y1_op<Scalar> > {
0231   enum {
0232     // 6 polynomial of order ~N=8 is computed.
0233     // The cost is N multiplications and N additions each, along with a
0234     // sine, cosine, rsqrt and j1 cost.
0235     Cost = 126 * NumTraits<Scalar>::MulCost + 96 * NumTraits<Scalar>::AddCost,
0236     PacketAccess = packet_traits<Scalar>::HasBessel
0237   };
0238 };
0239
0240 /** \internal
0241  * \brief Template functor to compute the modified Bessel function of the second
0242  * kind of order zero
0243  * \sa class CwiseUnaryOp, Cwise::bessel_k0()
0244  */
0245 template <typename Scalar>
0246 struct scalar_bessel_k0_op {
0247   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k0_op)
0248   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0249     using numext::bessel_k0;
0250     return bessel_k0(x);
0251   }
0252   typedef typename packet_traits<Scalar>::type Packet;
0253   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0254     return internal::pbessel_k0(x);
0255   }
0256 };
0257 template <typename Scalar>
0258 struct functor_traits<scalar_bessel_k0_op<Scalar> > {
0259   enum {
0260     // On average, a Chebyshev polynomial of order N=10 is computed.
0261     // The cost is N multiplications and 2N additions. In addition we compute
0262     // i0, a log, exp and prsqrt and sin and cos.
0263     Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
0264     PacketAccess = packet_traits<Scalar>::HasBessel
0265   };
0266 };
0267
0268 /** \internal
0269  * \brief Template functor to compute the exponentially scaled modified Bessel
0270  * function of the second kind of order zero
0271  * \sa class CwiseUnaryOp, Cwise::bessel_k0e()
0272  */
0273 template <typename Scalar>
0274 struct scalar_bessel_k0e_op {
0275   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k0e_op)
0276   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0277     using numext::bessel_k0e;
0278     return bessel_k0e(x);
0279   }
0280   typedef typename packet_traits<Scalar>::type Packet;
0281   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0282     return internal::pbessel_k0e(x);
0283   }
0284 };
0285 template <typename Scalar>
0286 struct functor_traits<scalar_bessel_k0e_op<Scalar> > {
0287   enum {
0288     // On average, a Chebyshev polynomial of order N=10 is computed.
0289     // The cost is N multiplications and 2N additions. In addition we compute
0290     // i0, a log, exp and prsqrt and sin and cos.
0291     Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
0292     PacketAccess = packet_traits<Scalar>::HasBessel
0293   };
0294 };
0295
0296 /** \internal
0297  * \brief Template functor to compute the modified Bessel function of the
0298  * second kind of order one
0299  * \sa class CwiseUnaryOp, Cwise::bessel_k1()
0300  */
0301 template <typename Scalar>
0302 struct scalar_bessel_k1_op {
0303   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k1_op)
0304   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0305     using numext::bessel_k1;
0306     return bessel_k1(x);
0307   }
0308   typedef typename packet_traits<Scalar>::type Packet;
0309   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0310     return internal::pbessel_k1(x);
0311   }
0312 };
0313 template <typename Scalar>
0314 struct functor_traits<scalar_bessel_k1_op<Scalar> > {
0315   enum {
0316     // On average, a Chebyshev polynomial of order N=10 is computed.
0317     // The cost is N multiplications and 2N additions. In addition we compute
0318     // i1, a log, exp and prsqrt and sin and cos.
0319     Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
0320     PacketAccess = packet_traits<Scalar>::HasBessel
0321   };
0322 };
0323
0324 /** \internal
0325  * \brief Template functor to compute the exponentially scaled modified Bessel
0326  * function of the second kind of order one
0327  * \sa class CwiseUnaryOp, Cwise::bessel_k1e()
0328  */
0329 template <typename Scalar>
0330 struct scalar_bessel_k1e_op {
0331   EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k1e_op)
0332   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
0333     using numext::bessel_k1e;
0334     return bessel_k1e(x);
0335   }
0336   typedef typename packet_traits<Scalar>::type Packet;
0337   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0338     return internal::pbessel_k1e(x);
0339   }
0340 };
0341 template <typename Scalar>
0342 struct functor_traits<scalar_bessel_k1e_op<Scalar> > {
0343   enum {
0344     // On average, a Chebyshev polynomial of order N=10 is computed.
0345     // The cost is N multiplications and 2N additions. In addition we compute
0346     // i1, a log, exp and prsqrt and sin and cos.
0347     Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
0348     PacketAccess = packet_traits<Scalar>::HasBessel
0349   };
0350 };
0351
0352
0353 } // end namespace internal
0354
0355 } // end namespace Eigen
0356
0357 #endif // EIGEN_BESSELFUNCTIONS_FUNCTORS_H