1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2016 Eugene Brevdo <[email protected]> 5 // Copyright (C) 2016 Gael Guennebaud <[email protected]> 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 { EIGEN_EMPTY_STRUCT_CTORscalar_bessel_i0_op25 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; packetOpscalar_bessel_i0_op31 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 358