aboutsummaryrefslogtreecommitdiffhomepage
path: root/unsupported/Eigen/src
diff options
context:
space:
mode:
authorGravatar Srinivas Vasudevan <srvasude@google.com>2019-09-14 12:16:47 -0400
committerGravatar Srinivas Vasudevan <srvasude@google.com>2019-09-14 12:16:47 -0400
commit6e215cf109073da9ffb5b491171613b8db24fd9d (patch)
tree1c171abbf72628ed0dbe37574e8d07c7953b4816 /unsupported/Eigen/src
parentfacdec5aa7d947d5462c9dbaefa7a50c4cabff3b (diff)
Add Bessel functions to SpecialFunctions.
- Split SpecialFunctions files in to a separate BesselFunctions file. In particular add: - Modified bessel functions of the second kind k0, k1, k0e, k1e - Bessel functions of the first kind j0, j1 - Bessel functions of the second kind y0, y1
Diffstat (limited to 'unsupported/Eigen/src')
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/BesselFunctionsArrayAPI.h286
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/BesselFunctionsFunctors.h357
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/BesselFunctionsHalf.h66
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h1959
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/BesselFunctionsPacketMath.h130
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h45
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h54
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h8
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h361
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h18
-rw-r--r--unsupported/Eigen/src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h130
11 files changed, 2929 insertions, 485 deletions
diff --git a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsArrayAPI.h b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsArrayAPI.h
new file mode 100644
index 000000000..8f96c2ae7
--- /dev/null
+++ b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsArrayAPI.h
@@ -0,0 +1,286 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+
+#ifndef EIGEN_BESSELFUNCTIONS_ARRAYAPI_H
+#define EIGEN_BESSELFUNCTIONS_ARRAYAPI_H
+
+namespace Eigen {
+
+/** \returns an expression of the coefficient-wise i0(\a x) to the given
+ * arrays.
+ *
+ * It returns the modified Bessel function of the first kind of order zero.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of i0(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::i0()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_i0_op<typename Derived::Scalar>, const Derived>
+i0(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_i0_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise i0e(\a x) to the given
+ * arrays.
+ *
+ * It returns the exponentially scaled modified Bessel
+ * function of the first kind of order zero.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of i0e(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::i0e()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_i0e_op<typename Derived::Scalar>, const Derived>
+i0e(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_i0e_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise i1(\a x) to the given
+ * arrays.
+ *
+ * It returns the modified Bessel function of the first kind of order one.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of i1(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::i1()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_i1_op<typename Derived::Scalar>, const Derived>
+i1(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_i1_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise i1e(\a x) to the given
+ * arrays.
+ *
+ * It returns the exponentially scaled modified Bessel
+ * function of the first kind of order one.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of i1e(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::i1e()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_i1e_op<typename Derived::Scalar>, const Derived>
+i1e(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_i1e_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise k0(\a x) to the given
+ * arrays.
+ *
+ * It returns the modified Bessel function of the second kind of order zero.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of k0(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::k0()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_k0_op<typename Derived::Scalar>, const Derived>
+k0(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_k0_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise k0e(\a x) to the given
+ * arrays.
+ *
+ * It returns the exponentially scaled modified Bessel
+ * function of the second kind of order zero.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of k0e(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::k0e()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_k0e_op<typename Derived::Scalar>, const Derived>
+k0e(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_k0e_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise k1(\a x) to the given
+ * arrays.
+ *
+ * It returns the modified Bessel function of the second kind of order one.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of k1(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::k1()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_k1_op<typename Derived::Scalar>, const Derived>
+k1(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_k1_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise k1e(\a x) to the given
+ * arrays.
+ *
+ * It returns the exponentially scaled modified Bessel
+ * function of the second kind of order one.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of k1e(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::k1e()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_k1e_op<typename Derived::Scalar>, const Derived>
+k1e(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_k1e_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise j0(\a x) to the given
+ * arrays.
+ *
+ * It returns the Bessel function of the first kind of order zero.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of j0(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::j0()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_j0_op<typename Derived::Scalar>, const Derived>
+j0(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_j0_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise y0(\a x) to the given
+ * arrays.
+ *
+ * It returns the Bessel function of the second kind of order zero.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of y0(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::y0()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_y0_op<typename Derived::Scalar>, const Derived>
+y0(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_y0_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise j1(\a x) to the given
+ * arrays.
+ *
+ * It returns the modified Bessel function of the first kind of order one.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of j1(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::j1()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_j1_op<typename Derived::Scalar>, const Derived>
+j1(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_j1_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+/** \returns an expression of the coefficient-wise y1(\a x) to the given
+ * arrays.
+ *
+ * It returns the Bessel function of the second kind of order one.
+ *
+ * \param x is the argument
+ *
+ * \note This function supports only float and double scalar types. To support
+ * other scalar types, the user has to provide implementations of y1(T) for
+ * any scalar type T to be supported.
+ *
+ * \sa ArrayBase::y1()
+ */
+template <typename Derived>
+EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_y1_op<typename Derived::Scalar>, const Derived>
+y1(const Eigen::ArrayBase<Derived>& x) {
+ return Eigen::CwiseUnaryOp<
+ Eigen::internal::scalar_bessel_y1_op<typename Derived::Scalar>,
+ const Derived>(x.derived());
+}
+
+} // end namespace Eigen
+
+#endif // EIGEN_BESSELFUNCTIONS_ARRAYAPI_H
diff --git a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsFunctors.h b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsFunctors.h
new file mode 100644
index 000000000..e57d5042b
--- /dev/null
+++ b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsFunctors.h
@@ -0,0 +1,357 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2016 Eugene Brevdo <ebrevdo@gmail.com>
+// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_BESSELFUNCTIONS_FUNCTORS_H
+#define EIGEN_BESSELFUNCTIONS_FUNCTORS_H
+
+namespace Eigen {
+
+namespace internal {
+
+/** \internal
+ * \brief Template functor to compute the modified Bessel function of the first
+ * kind of order zero.
+ * \sa class CwiseUnaryOp, Cwise::i0()
+ */
+template <typename Scalar>
+struct scalar_bessel_i0_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i0_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::i0;
+ return i0(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pi0(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_i0_op<Scalar> > {
+ enum {
+ // On average, a Chebyshev polynomial of order N=20 is computed.
+ // The cost is N multiplications and 2N additions. We also add
+ // the cost of an additional exp over i0e.
+ Cost = 28 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the exponentially scaled modified Bessel
+ * function of the first kind of order zero
+ * \sa class CwiseUnaryOp, Cwise::i0e()
+ */
+template <typename Scalar>
+struct scalar_bessel_i0e_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i0e_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::i0e;
+ return i0e(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pi0e(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_i0e_op<Scalar> > {
+ enum {
+ // On average, a Chebyshev polynomial of order N=20 is computed.
+ // The cost is N multiplications and 2N additions.
+ Cost = 20 * NumTraits<Scalar>::MulCost + 40 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the modified Bessel function of the first
+ * kind of order one
+ * \sa class CwiseUnaryOp, Cwise::i1()
+ */
+template <typename Scalar>
+struct scalar_bessel_i1_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i1_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::i1;
+ return i1(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pi1(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_i1_op<Scalar> > {
+ enum {
+ // On average, a Chebyshev polynomial of order N=20 is computed.
+ // The cost is N multiplications and 2N additions. We also add
+ // the cost of an additional exp over i1e.
+ Cost = 28 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the exponentially scaled modified Bessel
+ * function of the first kind of order zero
+ * \sa class CwiseUnaryOp, Cwise::i1e()
+ */
+template <typename Scalar>
+struct scalar_bessel_i1e_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_i1e_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::i1e;
+ return i1e(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pi1e(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_i1e_op<Scalar> > {
+ enum {
+ // On average, a Chebyshev polynomial of order N=20 is computed.
+ // The cost is N multiplications and 2N additions.
+ Cost = 20 * NumTraits<Scalar>::MulCost + 40 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the Bessel function of the second kind of
+ * order zero
+ * \sa class CwiseUnaryOp, Cwise::j0()
+ */
+template <typename Scalar>
+struct scalar_bessel_j0_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_j0_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::j0;
+ return j0(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pj0(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_j0_op<Scalar> > {
+ enum {
+ // 6 polynomial of order ~N=8 is computed.
+ // The cost is N multiplications and N additions each, along with a
+ // sine, cosine and rsqrt cost.
+ Cost = 63 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the Bessel function of the second kind of
+ * order zero
+ * \sa class CwiseUnaryOp, Cwise::y0()
+ */
+template <typename Scalar>
+struct scalar_bessel_y0_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_y0_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::y0;
+ return y0(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::py0(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_y0_op<Scalar> > {
+ enum {
+ // 6 polynomial of order ~N=8 is computed.
+ // The cost is N multiplications and N additions each, along with a
+ // sine, cosine, rsqrt and j0 cost.
+ Cost = 126 * NumTraits<Scalar>::MulCost + 96 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the Bessel function of the first kind of
+ * order one
+ * \sa class CwiseUnaryOp, Cwise::j1()
+ */
+template <typename Scalar>
+struct scalar_bessel_j1_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_j1_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::j1;
+ return j1(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pj1(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_j1_op<Scalar> > {
+ enum {
+ // 6 polynomial of order ~N=8 is computed.
+ // The cost is N multiplications and N additions each, along with a
+ // sine, cosine and rsqrt cost.
+ Cost = 63 * NumTraits<Scalar>::MulCost + 48 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the Bessel function of the second kind of
+ * order one
+ * \sa class CwiseUnaryOp, Cwise::j1e()
+ */
+template <typename Scalar>
+struct scalar_bessel_y1_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_y1_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::y1;
+ return y1(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::py1(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_y1_op<Scalar> > {
+ enum {
+ // 6 polynomial of order ~N=8 is computed.
+ // The cost is N multiplications and N additions each, along with a
+ // sine, cosine, rsqrt and j1 cost.
+ Cost = 126 * NumTraits<Scalar>::MulCost + 96 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the modified Bessel function of the second
+ * kind of order zero
+ * \sa class CwiseUnaryOp, Cwise::k0()
+ */
+template <typename Scalar>
+struct scalar_bessel_k0_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k0_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::k0;
+ return k0(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pk0(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_k0_op<Scalar> > {
+ enum {
+ // On average, a Chebyshev polynomial of order N=10 is computed.
+ // The cost is N multiplications and 2N additions. In addition we compute
+ // i0, a log, exp and prsqrt and sin and cos.
+ Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the exponentially scaled modified Bessel
+ * function of the second kind of order zero
+ * \sa class CwiseUnaryOp, Cwise::k0e()
+ */
+template <typename Scalar>
+struct scalar_bessel_k0e_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k0e_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::k0e;
+ return k0e(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pk0e(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_k0e_op<Scalar> > {
+ enum {
+ // On average, a Chebyshev polynomial of order N=10 is computed.
+ // The cost is N multiplications and 2N additions. In addition we compute
+ // i0, a log, exp and prsqrt and sin and cos.
+ Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the modified Bessel function of the
+ * second kind of order one
+ * \sa class CwiseUnaryOp, Cwise::k1()
+ */
+template <typename Scalar>
+struct scalar_bessel_k1_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k1_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::k1;
+ return k1(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pk1(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_k1_op<Scalar> > {
+ enum {
+ // On average, a Chebyshev polynomial of order N=10 is computed.
+ // The cost is N multiplications and 2N additions. In addition we compute
+ // i1, a log, exp and prsqrt and sin and cos.
+ Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+/** \internal
+ * \brief Template functor to compute the exponentially scaled modified Bessel
+ * function of the second kind of order one
+ * \sa class CwiseUnaryOp, Cwise::k1e()
+ */
+template <typename Scalar>
+struct scalar_bessel_k1e_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_bessel_k1e_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
+ using numext::k1e;
+ return k1e(x);
+ }
+ typedef typename packet_traits<Scalar>::type Packet;
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
+ return internal::pk1e(x);
+ }
+};
+template <typename Scalar>
+struct functor_traits<scalar_bessel_k1e_op<Scalar> > {
+ enum {
+ // On average, a Chebyshev polynomial of order N=10 is computed.
+ // The cost is N multiplications and 2N additions. In addition we compute
+ // i1, a log, exp and prsqrt and sin and cos.
+ Cost = 68 * NumTraits<Scalar>::MulCost + 88 * NumTraits<Scalar>::AddCost,
+ PacketAccess = packet_traits<Scalar>::HasBessel
+ };
+};
+
+
+} // end namespace internal
+
+} // end namespace Eigen
+
+#endif // EIGEN_BESSELFUNCTIONS_FUNCTORS_H
diff --git a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsHalf.h b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsHalf.h
new file mode 100644
index 000000000..9aad7ac96
--- /dev/null
+++ b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsHalf.h
@@ -0,0 +1,66 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_BESSELFUNCTIONS_HALF_H
+#define EIGEN_BESSELFUNCTIONS_HALF_H
+
+namespace Eigen {
+namespace numext {
+
+#if EIGEN_HAS_C99_MATH
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half i0(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::i0(static_cast<float>(x)));
+}
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half i0e(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::i0e(static_cast<float>(x)));
+}
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half i1(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::i1(static_cast<float>(x)));
+}
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half i1e(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::i1e(static_cast<float>(x)));
+}
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half j0(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::j0(static_cast<float>(x)));
+}
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half j1(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::j1(static_cast<float>(x)));
+}
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half y0(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::y0(static_cast<float>(x)));
+}
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half y1(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::y1(static_cast<float>(x)));
+}
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half k0(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::k0(static_cast<float>(x)));
+}
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half k0e(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::k0e(static_cast<float>(x)));
+}
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half k1(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::k1(static_cast<float>(x)));
+}
+template <>
+EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half k1e(const Eigen::half& x) {
+ return Eigen::half(Eigen::numext::k1e(static_cast<float>(x)));
+}
+#endif
+
+} // end namespace numext
+} // end namespace Eigen
+
+#endif // EIGEN_BESSELFUNCTIONS_HALF_H
diff --git a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h
new file mode 100644
index 000000000..b279687c2
--- /dev/null
+++ b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsImpl.h
@@ -0,0 +1,1959 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2015 Eugene Brevdo <ebrevdo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_BESSEL_FUNCTIONS_H
+#define EIGEN_BESSEL_FUNCTIONS_H
+
+namespace Eigen {
+namespace internal {
+
+// Parts of this code are based on the Cephes Math Library.
+//
+// Cephes Math Library Release 2.8: June, 2000
+// Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier
+//
+// Permission has been kindly provided by the original author
+// to incorporate the Cephes software into the Eigen codebase:
+//
+// From: Stephen Moshier
+// To: Eugene Brevdo
+// Subject: Re: Permission to wrap several cephes functions in Eigen
+//
+// Hello Eugene,
+//
+// Thank you for writing.
+//
+// If your licensing is similar to BSD, the formal way that has been
+// handled is simply to add a statement to the effect that you are incorporating
+// the Cephes software by permission of the author.
+//
+// Good luck with your project,
+// Steve
+
+
+/****************************************************************************
+ * Implementation of Bessel function, based on Cephes *
+ ****************************************************************************/
+
+template <typename Scalar>
+struct i0e_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_i0e {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_i0e<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* i0ef.c
+ *
+ * Modified Bessel function of order zero,
+ * exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, i0ef();
+ *
+ * y = i0ef( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of order zero of the argument.
+ *
+ * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0,30 100000 3.7e-7 7.0e-8
+ * See i0f().
+ *
+ */
+
+ const float A[] = {-1.30002500998624804212E-8f, 6.04699502254191894932E-8f,
+ -2.67079385394061173391E-7f, 1.11738753912010371815E-6f,
+ -4.41673835845875056359E-6f, 1.64484480707288970893E-5f,
+ -5.75419501008210370398E-5f, 1.88502885095841655729E-4f,
+ -5.76375574538582365885E-4f, 1.63947561694133579842E-3f,
+ -4.32430999505057594430E-3f, 1.05464603945949983183E-2f,
+ -2.37374148058994688156E-2f, 4.93052842396707084878E-2f,
+ -9.49010970480476444210E-2f, 1.71620901522208775349E-1f,
+ -3.04682672343198398683E-1f, 6.76795274409476084995E-1f};
+
+ const float B[] = {3.39623202570838634515E-9f, 2.26666899049817806459E-8f,
+ 2.04891858946906374183E-7f, 2.89137052083475648297E-6f,
+ 6.88975834691682398426E-5f, 3.36911647825569408990E-3f,
+ 8.04490411014108831608E-1f};
+ T y = pabs(x);
+ T y_le_eight = internal::pchebevl<T, 18>::run(
+ pmadd(pset1<T>(0.5f), y, pset1<T>(-2.0f)), A);
+ T y_gt_eight = pmul(
+ internal::pchebevl<T, 7>::run(
+ psub(pdiv(pset1<T>(32.0f), y), pset1<T>(2.0f)), B),
+ prsqrt(y));
+ // TODO: Perhaps instead check whether all packet elements are in
+ // [-8, 8] and evaluate a branch based off of that. It's possible
+ // in practice most elements are in this region.
+ return pselect(pcmp_le(y, pset1<T>(8.0f)), y_le_eight, y_gt_eight);
+ }
+};
+
+template <typename T>
+struct generic_i0e<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* i0e.c
+ *
+ * Modified Bessel function of order zero,
+ * exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, i0e();
+ *
+ * y = i0e( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of order zero of the argument.
+ *
+ * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0,30 30000 5.4e-16 1.2e-16
+ * See i0().
+ *
+ */
+
+ const double A[] = {-4.41534164647933937950E-18, 3.33079451882223809783E-17,
+ -2.43127984654795469359E-16, 1.71539128555513303061E-15,
+ -1.16853328779934516808E-14, 7.67618549860493561688E-14,
+ -4.85644678311192946090E-13, 2.95505266312963983461E-12,
+ -1.72682629144155570723E-11, 9.67580903537323691224E-11,
+ -5.18979560163526290666E-10, 2.65982372468238665035E-9,
+ -1.30002500998624804212E-8, 6.04699502254191894932E-8,
+ -2.67079385394061173391E-7, 1.11738753912010371815E-6,
+ -4.41673835845875056359E-6, 1.64484480707288970893E-5,
+ -5.75419501008210370398E-5, 1.88502885095841655729E-4,
+ -5.76375574538582365885E-4, 1.63947561694133579842E-3,
+ -4.32430999505057594430E-3, 1.05464603945949983183E-2,
+ -2.37374148058994688156E-2, 4.93052842396707084878E-2,
+ -9.49010970480476444210E-2, 1.71620901522208775349E-1,
+ -3.04682672343198398683E-1, 6.76795274409476084995E-1};
+ const double B[] = {
+ -7.23318048787475395456E-18, -4.83050448594418207126E-18,
+ 4.46562142029675999901E-17, 3.46122286769746109310E-17,
+ -2.82762398051658348494E-16, -3.42548561967721913462E-16,
+ 1.77256013305652638360E-15, 3.81168066935262242075E-15,
+ -9.55484669882830764870E-15, -4.15056934728722208663E-14,
+ 1.54008621752140982691E-14, 3.85277838274214270114E-13,
+ 7.18012445138366623367E-13, -1.79417853150680611778E-12,
+ -1.32158118404477131188E-11, -3.14991652796324136454E-11,
+ 1.18891471078464383424E-11, 4.94060238822496958910E-10,
+ 3.39623202570838634515E-9, 2.26666899049817806459E-8,
+ 2.04891858946906374183E-7, 2.89137052083475648297E-6,
+ 6.88975834691682398426E-5, 3.36911647825569408990E-3,
+ 8.04490411014108831608E-1};
+ T y = pabs(x);
+ T y_le_eight = internal::pchebevl<T, 30>::run(
+ pmadd(pset1<T>(0.5), y, pset1<T>(-2.0)), A);
+ T y_gt_eight = pmul(
+ internal::pchebevl<T, 25>::run(
+ psub(pdiv(pset1<T>(32.0), y), pset1<T>(2.0)), B),
+ prsqrt(y));
+ // TODO: Perhaps instead check whether all packet elements are in
+ // [-8, 8] and evaluate a branch based off of that. It's possible
+ // in practice most elements are in this region.
+ return pselect(pcmp_le(y, pset1<T>(8.0)), y_le_eight, y_gt_eight);
+ }
+};
+
+template <typename Scalar>
+struct i0e_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_i0e<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct i0_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_i0 {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ return pmul(
+ pexp(pabs(x)),
+ generic_i0e<T, ScalarType>::run(x));
+ }
+};
+
+template <typename Scalar>
+struct i0_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_i0<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct i1e_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_i1e {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_i1e<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* i1ef.c
+ *
+ * Modified Bessel function of order one,
+ * exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, i1ef();
+ *
+ * y = i1ef( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of order one of the argument.
+ *
+ * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 1.5e-6 1.5e-7
+ * See i1().
+ *
+ */
+ const float A[] = {9.38153738649577178388E-9f, -4.44505912879632808065E-8f,
+ 2.00329475355213526229E-7f, -8.56872026469545474066E-7f,
+ 3.47025130813767847674E-6f, -1.32731636560394358279E-5f,
+ 4.78156510755005422638E-5f, -1.61760815825896745588E-4f,
+ 5.12285956168575772895E-4f, -1.51357245063125314899E-3f,
+ 4.15642294431288815669E-3f, -1.05640848946261981558E-2f,
+ 2.47264490306265168283E-2f, -5.29459812080949914269E-2f,
+ 1.02643658689847095384E-1f, -1.76416518357834055153E-1f,
+ 2.52587186443633654823E-1f};
+
+ const float B[] = {-3.83538038596423702205E-9f, -2.63146884688951950684E-8f,
+ -2.51223623787020892529E-7f, -3.88256480887769039346E-6f,
+ -1.10588938762623716291E-4f, -9.76109749136146840777E-3f,
+ 7.78576235018280120474E-1f};
+
+
+ T y = pabs(x);
+ T y_le_eight = pmul(y, internal::pchebevl<T, 17>::run(
+ pmadd(pset1<T>(0.5f), y, pset1<T>(-2.0f)), A));
+ T y_gt_eight = pmul(
+ internal::pchebevl<T, 7>::run(
+ psub(pdiv(pset1<T>(32.0f), y),
+ pset1<T>(2.0f)), B),
+ prsqrt(y));
+ // TODO: Perhaps instead check whether all packet elements are in
+ // [-8, 8] and evaluate a branch based off of that. It's possible
+ // in practice most elements are in this region.
+ y = pselect(pcmp_le(y, pset1<T>(8.0f)), y_le_eight, y_gt_eight);
+ return pselect(pcmp_lt(x, pset1<T>(0.0f)), -y, y);
+ }
+};
+
+template <typename T>
+struct generic_i1e<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* i1e.c
+ *
+ * Modified Bessel function of order one,
+ * exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, i1e();
+ *
+ * y = i1e( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of order one of the argument.
+ *
+ * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 2.0e-15 2.0e-16
+ * See i1().
+ *
+ */
+ const double A[] = {2.77791411276104639959E-18, -2.11142121435816608115E-17,
+ 1.55363195773620046921E-16, -1.10559694773538630805E-15,
+ 7.60068429473540693410E-15, -5.04218550472791168711E-14,
+ 3.22379336594557470981E-13, -1.98397439776494371520E-12,
+ 1.17361862988909016308E-11, -6.66348972350202774223E-11,
+ 3.62559028155211703701E-10, -1.88724975172282928790E-9,
+ 9.38153738649577178388E-9, -4.44505912879632808065E-8,
+ 2.00329475355213526229E-7, -8.56872026469545474066E-7,
+ 3.47025130813767847674E-6, -1.32731636560394358279E-5,
+ 4.78156510755005422638E-5, -1.61760815825896745588E-4,
+ 5.12285956168575772895E-4, -1.51357245063125314899E-3,
+ 4.15642294431288815669E-3, -1.05640848946261981558E-2,
+ 2.47264490306265168283E-2, -5.29459812080949914269E-2,
+ 1.02643658689847095384E-1, -1.76416518357834055153E-1,
+ 2.52587186443633654823E-1};
+ const double B[] = {
+ 7.51729631084210481353E-18, 4.41434832307170791151E-18,
+ -4.65030536848935832153E-17, -3.20952592199342395980E-17,
+ 2.96262899764595013876E-16, 3.30820231092092828324E-16,
+ -1.88035477551078244854E-15, -3.81440307243700780478E-15,
+ 1.04202769841288027642E-14, 4.27244001671195135429E-14,
+ -2.10154184277266431302E-14, -4.08355111109219731823E-13,
+ -7.19855177624590851209E-13, 2.03562854414708950722E-12,
+ 1.41258074366137813316E-11, 3.25260358301548823856E-11,
+ -1.89749581235054123450E-11, -5.58974346219658380687E-10,
+ -3.83538038596423702205E-9, -2.63146884688951950684E-8,
+ -2.51223623787020892529E-7, -3.88256480887769039346E-6,
+ -1.10588938762623716291E-4, -9.76109749136146840777E-3,
+ 7.78576235018280120474E-1};
+ T y = pabs(x);
+ T y_le_eight = pmul(y, internal::pchebevl<T, 29>::run(
+ pmadd(pset1<T>(0.5), y, pset1<T>(-2.0)), A));
+ T y_gt_eight = pmul(
+ internal::pchebevl<T, 25>::run(
+ psub(pdiv(pset1<T>(32.0), y),
+ pset1<T>(2.0)), B),
+ prsqrt(y));
+ // TODO: Perhaps instead check whether all packet elements are in
+ // [-8, 8] and evaluate a branch based off of that. It's possible
+ // in practice most elements are in this region.
+ y = pselect(pcmp_le(y, pset1<T>(8.0)), y_le_eight, y_gt_eight);
+ return pselect(pcmp_lt(x, pset1<T>(0.0f)), -y, y);
+ }
+};
+
+template <typename Scalar>
+struct i1e_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_i1e<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct i1_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_i1 {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ return pmul(
+ pexp(pabs(x)),
+ generic_i1e<T, ScalarType>::run(x));
+ }
+};
+
+template <typename Scalar>
+struct i1_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_i1<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct k0e_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_k0e {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_k0e<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* k0ef.c
+ * Modified Bessel function, third kind, order zero,
+ * exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, k0ef();
+ *
+ * y = k0ef( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of the third kind of order zero of the argument.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 8.1e-7 7.8e-8
+ * See k0().
+ *
+ */
+
+ const float A[] = {1.90451637722020886025E-9f, 2.53479107902614945675E-7f,
+ 2.28621210311945178607E-5f, 1.26461541144692592338E-3f,
+ 3.59799365153615016266E-2f, 3.44289899924628486886E-1f,
+ -5.35327393233902768720E-1f};
+
+ const float B[] = {-1.69753450938905987466E-9f, 8.57403401741422608519E-9f,
+ -4.66048989768794782956E-8f, 2.76681363944501510342E-7f,
+ -1.83175552271911948767E-6f, 1.39498137188764993662E-5f,
+ -1.28495495816278026384E-4f, 1.56988388573005337491E-3f,
+ -3.14481013119645005427E-2f, 2.44030308206595545468E0f};
+ const T MAXNUM = pset1<T>(NumTraits<float>::infinity());
+ const T two = pset1<T>(2.0);
+ T x_le_two = internal::pchebevl<T, 7>::run(
+ pmadd(x, x, pset1<T>(-2.0)), A);
+ x_le_two = pmadd(
+ generic_i0<T, float>::run(x), pmul(
+ pset1<T>(-1.0), plog(pmul(pset1<T>(0.5), x))), x_le_two);
+ x_le_two = pmul(pexp(x), x_le_two);
+ T x_gt_two = pmul(
+ internal::pchebevl<T, 10>::run(
+ psub(pdiv(pset1<T>(8.0), x), two), B),
+ prsqrt(x));
+ return pselect(
+ pcmp_le(x, pset1<T>(0.0)),
+ MAXNUM,
+ pselect(pcmp_le(x, two), x_le_two, x_gt_two));
+ }
+};
+
+template <typename T>
+struct generic_k0e<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* k0e.c
+ * Modified Bessel function, third kind, order zero,
+ * exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, k0e();
+ *
+ * y = k0e( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of the third kind of order zero of the argument.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 1.4e-15 1.4e-16
+ * See k0().
+ *
+ */
+
+ const double A[] = {
+ 1.37446543561352307156E-16,
+ 4.25981614279661018399E-14,
+ 1.03496952576338420167E-11,
+ 1.90451637722020886025E-9,
+ 2.53479107902614945675E-7,
+ 2.28621210311945178607E-5,
+ 1.26461541144692592338E-3,
+ 3.59799365153615016266E-2,
+ 3.44289899924628486886E-1,
+ -5.35327393233902768720E-1};
+ const double B[] = {
+ 5.30043377268626276149E-18, -1.64758043015242134646E-17,
+ 5.21039150503902756861E-17, -1.67823109680541210385E-16,
+ 5.51205597852431940784E-16, -1.84859337734377901440E-15,
+ 6.34007647740507060557E-15, -2.22751332699166985548E-14,
+ 8.03289077536357521100E-14, -2.98009692317273043925E-13,
+ 1.14034058820847496303E-12, -4.51459788337394416547E-12,
+ 1.85594911495471785253E-11, -7.95748924447710747776E-11,
+ 3.57739728140030116597E-10, -1.69753450938905987466E-9,
+ 8.57403401741422608519E-9, -4.66048989768794782956E-8,
+ 2.76681363944501510342E-7, -1.83175552271911948767E-6,
+ 1.39498137188764993662E-5, -1.28495495816278026384E-4,
+ 1.56988388573005337491E-3, -3.14481013119645005427E-2,
+ 2.44030308206595545468E0
+ };
+ const T MAXNUM = pset1<T>(NumTraits<double>::infinity());
+ const T two = pset1<T>(2.0);
+ T x_le_two = internal::pchebevl<T, 10>::run(
+ pmadd(x, x, pset1<T>(-2.0)), A);
+ x_le_two = pmadd(
+ generic_i0<T, double>::run(x), pmul(
+ pset1<T>(-1.0), plog(pmul(pset1<T>(0.5), x))), x_le_two);
+ x_le_two = pmul(pexp(x), x_le_two);
+ x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
+ T x_gt_two = pmul(
+ internal::pchebevl<T, 25>::run(
+ psub(pdiv(pset1<T>(8.0), x), two), B),
+ prsqrt(x));
+ return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
+ }
+};
+
+template <typename Scalar>
+struct k0e_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_k0e<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct k0_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_k0 {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_k0<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* k0f.c
+ * Modified Bessel function, third kind, order zero
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, k0f();
+ *
+ * y = k0f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns modified Bessel function of the third kind
+ * of order zero of the argument.
+ *
+ * The range is partitioned into the two intervals [0,8] and
+ * (8, infinity). Chebyshev polynomial expansions are employed
+ * in each interval.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Tested at 2000 random points between 0 and 8. Peak absolute
+ * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 7.8e-7 8.5e-8
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * K0 domain x <= 0 MAXNUM
+ *
+ */
+
+ const float A[] = {1.90451637722020886025E-9f, 2.53479107902614945675E-7f,
+ 2.28621210311945178607E-5f, 1.26461541144692592338E-3f,
+ 3.59799365153615016266E-2f, 3.44289899924628486886E-1f,
+ -5.35327393233902768720E-1f};
+
+ const float B[] = {-1.69753450938905987466E-9f, 8.57403401741422608519E-9f,
+ -4.66048989768794782956E-8f, 2.76681363944501510342E-7f,
+ -1.83175552271911948767E-6f, 1.39498137188764993662E-5f,
+ -1.28495495816278026384E-4f, 1.56988388573005337491E-3f,
+ -3.14481013119645005427E-2f, 2.44030308206595545468E0f};
+ const T MAXNUM = pset1<T>(NumTraits<float>::infinity());
+ const T two = pset1<T>(2.0);
+ T x_le_two = internal::pchebevl<T, 7>::run(
+ pmadd(x, x, pset1<T>(-2.0)), A);
+ x_le_two = pmadd(
+ generic_i0<T, float>::run(x), pmul(
+ pset1<T>(-1.0), plog(pmul(pset1<T>(0.5), x))), x_le_two);
+ x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
+ T x_gt_two = pmul(
+ pmul(
+ pexp(-x),
+ internal::pchebevl<T, 10>::run(
+ psub(pdiv(pset1<T>(8.0), x), two), B)),
+ prsqrt(x));
+ return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
+ }
+};
+
+template <typename T>
+struct generic_k0<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /*
+ *
+ * Modified Bessel function, third kind, order zero,
+ * exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, k0();
+ *
+ * y = k0( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of the third kind of order zero of the argument.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 1.4e-15 1.4e-16
+ * See k0().
+ *
+ */
+ const double A[] = {
+ 1.37446543561352307156E-16,
+ 4.25981614279661018399E-14,
+ 1.03496952576338420167E-11,
+ 1.90451637722020886025E-9,
+ 2.53479107902614945675E-7,
+ 2.28621210311945178607E-5,
+ 1.26461541144692592338E-3,
+ 3.59799365153615016266E-2,
+ 3.44289899924628486886E-1,
+ -5.35327393233902768720E-1};
+ const double B[] = {
+ 5.30043377268626276149E-18, -1.64758043015242134646E-17,
+ 5.21039150503902756861E-17, -1.67823109680541210385E-16,
+ 5.51205597852431940784E-16, -1.84859337734377901440E-15,
+ 6.34007647740507060557E-15, -2.22751332699166985548E-14,
+ 8.03289077536357521100E-14, -2.98009692317273043925E-13,
+ 1.14034058820847496303E-12, -4.51459788337394416547E-12,
+ 1.85594911495471785253E-11, -7.95748924447710747776E-11,
+ 3.57739728140030116597E-10, -1.69753450938905987466E-9,
+ 8.57403401741422608519E-9, -4.66048989768794782956E-8,
+ 2.76681363944501510342E-7, -1.83175552271911948767E-6,
+ 1.39498137188764993662E-5, -1.28495495816278026384E-4,
+ 1.56988388573005337491E-3, -3.14481013119645005427E-2,
+ 2.44030308206595545468E0
+ };
+ const T MAXNUM = pset1<T>(NumTraits<double>::infinity());
+ const T two = pset1<T>(2.0);
+ T x_le_two = internal::pchebevl<T, 10>::run(
+ pmadd(x, x, pset1<T>(-2.0)), A);
+ x_le_two = pmadd(
+ generic_i0<T, double>::run(x), pmul(
+ pset1<T>(-1.0), plog(pmul(pset1<T>(0.5), x))), x_le_two);
+ x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
+ T x_gt_two = pmul(
+ pmul(
+ pexp(-x),
+ internal::pchebevl<T, 25>::run(
+ psub(pdiv(pset1<T>(8.0), x), two), B)),
+ prsqrt(x));
+ return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
+ }
+};
+
+template <typename Scalar>
+struct k0_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_k0<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct k1e_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_k1e {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_k1e<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* k1ef.c
+ *
+ * Modified Bessel function, third kind, order one,
+ * exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, k1ef();
+ *
+ * y = k1ef( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of the third kind of order one of the argument:
+ *
+ * k1e(x) = exp(x) * k1(x).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 4.9e-7 6.7e-8
+ * See k1().
+ *
+ */
+
+ const float A[] = {-2.21338763073472585583E-8f, -2.43340614156596823496E-6f,
+ -1.73028895751305206302E-4f, -6.97572385963986435018E-3f,
+ -1.22611180822657148235E-1f, -3.53155960776544875667E-1f,
+ 1.52530022733894777053E0f};
+ const float B[] = {2.01504975519703286596E-9f, -1.03457624656780970260E-8f,
+ 5.74108412545004946722E-8f, -3.50196060308781257119E-7f,
+ 2.40648494783721712015E-6f, -1.93619797416608296024E-5f,
+ 1.95215518471351631108E-4f, -2.85781685962277938680E-3f,
+ 1.03923736576817238437E-1f, 2.72062619048444266945E0f};
+ const T MAXNUM = pset1<T>(NumTraits<float>::infinity());
+ const T two = pset1<T>(2.0);
+ T x_le_two = pdiv(internal::pchebevl<T, 7>::run(
+ pmadd(x, x, pset1<T>(-2.0)), A), x);
+ x_le_two = pmadd(
+ generic_i1<T, float>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two);
+ x_le_two = pmul(x_le_two, pexp(x));
+ x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
+ T x_gt_two = pmul(
+ internal::pchebevl<T, 10>::run(
+ psub(pdiv(pset1<T>(8.0), x), two), B),
+ prsqrt(x));
+ return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
+ }
+};
+
+template <typename T>
+struct generic_k1e<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* k1e.c
+ *
+ * Modified Bessel function, third kind, order one,
+ * exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, k1e();
+ *
+ * y = k1e( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of the third kind of order one of the argument:
+ *
+ * k1e(x) = exp(x) * k1(x).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 7.8e-16 1.2e-16
+ * See k1().
+ *
+ */
+ const double A[] = {-7.02386347938628759343E-18, -2.42744985051936593393E-15,
+ -6.66690169419932900609E-13, -1.41148839263352776110E-10,
+ -2.21338763073472585583E-8, -2.43340614156596823496E-6,
+ -1.73028895751305206302E-4, -6.97572385963986435018E-3,
+ -1.22611180822657148235E-1, -3.53155960776544875667E-1,
+ 1.52530022733894777053E0};
+ const double B[] = {-5.75674448366501715755E-18, 1.79405087314755922667E-17,
+ -5.68946255844285935196E-17, 1.83809354436663880070E-16,
+ -6.05704724837331885336E-16, 2.03870316562433424052E-15,
+ -7.01983709041831346144E-15, 2.47715442448130437068E-14,
+ -8.97670518232499435011E-14, 3.34841966607842919884E-13,
+ -1.28917396095102890680E-12, 5.13963967348173025100E-12,
+ -2.12996783842756842877E-11, 9.21831518760500529508E-11,
+ -4.19035475934189648750E-10, 2.01504975519703286596E-9,
+ -1.03457624656780970260E-8, 5.74108412545004946722E-8,
+ -3.50196060308781257119E-7, 2.40648494783721712015E-6,
+ -1.93619797416608296024E-5, 1.95215518471351631108E-4,
+ -2.85781685962277938680E-3, 1.03923736576817238437E-1,
+ 2.72062619048444266945E0};
+ const T MAXNUM = pset1<T>(NumTraits<double>::infinity());
+ const T two = pset1<T>(2.0);
+ T x_le_two = pdiv(internal::pchebevl<T, 11>::run(
+ pmadd(x, x, pset1<T>(-2.0)), A), x);
+ x_le_two = pmadd(
+ generic_i1<T, double>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two);
+ x_le_two = pmul(x_le_two, pexp(x));
+ x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
+ T x_gt_two = pmul(
+ internal::pchebevl<T, 25>::run(
+ psub(pdiv(pset1<T>(8.0), x), two), B),
+ prsqrt(x));
+ return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
+ }
+};
+
+template <typename Scalar>
+struct k1e_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_k1e<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct k1_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_k1 {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_k1<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* k1f.c
+ * Modified Bessel function, third kind, order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, k1f();
+ *
+ * y = k1f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes the modified Bessel function of the third kind
+ * of order one of the argument.
+ *
+ * The range is partitioned into the two intervals [0,2] and
+ * (2, infinity). Chebyshev polynomial expansions are employed
+ * in each interval.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 4.6e-7 7.6e-8
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * k1 domain x <= 0 MAXNUM
+ *
+ */
+
+ const float A[] = {-2.21338763073472585583E-8f, -2.43340614156596823496E-6f,
+ -1.73028895751305206302E-4f, -6.97572385963986435018E-3f,
+ -1.22611180822657148235E-1f, -3.53155960776544875667E-1f,
+ 1.52530022733894777053E0f};
+ const float B[] = {2.01504975519703286596E-9f, -1.03457624656780970260E-8f,
+ 5.74108412545004946722E-8f, -3.50196060308781257119E-7f,
+ 2.40648494783721712015E-6f, -1.93619797416608296024E-5f,
+ 1.95215518471351631108E-4f, -2.85781685962277938680E-3f,
+ 1.03923736576817238437E-1f, 2.72062619048444266945E0f};
+ const T MAXNUM = pset1<T>(NumTraits<float>::infinity());
+ const T two = pset1<T>(2.0);
+ T x_le_two = pdiv(internal::pchebevl<T, 7>::run(
+ pmadd(x, x, pset1<T>(-2.0)), A), x);
+ x_le_two = pmadd(
+ generic_i1<T, float>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two);
+ x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
+ T x_gt_two = pmul(
+ pexp(-x),
+ pmul(
+ internal::pchebevl<T, 10>::run(
+ psub(pdiv(pset1<T>(8.0), x), two), B),
+ prsqrt(x)));
+ return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
+ }
+};
+
+template <typename T>
+struct generic_k1<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* k1.c
+ * Modified Bessel function, third kind, order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, k1f();
+ *
+ * y = k1f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes the modified Bessel function of the third kind
+ * of order one of the argument.
+ *
+ * The range is partitioned into the two intervals [0,2] and
+ * (2, infinity). Chebyshev polynomial expansions are employed
+ * in each interval.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 30 30000 4.6e-7 7.6e-8
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * k1 domain x <= 0 MAXNUM
+ *
+ */
+ const double A[] = {-7.02386347938628759343E-18, -2.42744985051936593393E-15,
+ -6.66690169419932900609E-13, -1.41148839263352776110E-10,
+ -2.21338763073472585583E-8, -2.43340614156596823496E-6,
+ -1.73028895751305206302E-4, -6.97572385963986435018E-3,
+ -1.22611180822657148235E-1, -3.53155960776544875667E-1,
+ 1.52530022733894777053E0};
+ const double B[] = {-5.75674448366501715755E-18, 1.79405087314755922667E-17,
+ -5.68946255844285935196E-17, 1.83809354436663880070E-16,
+ -6.05704724837331885336E-16, 2.03870316562433424052E-15,
+ -7.01983709041831346144E-15, 2.47715442448130437068E-14,
+ -8.97670518232499435011E-14, 3.34841966607842919884E-13,
+ -1.28917396095102890680E-12, 5.13963967348173025100E-12,
+ -2.12996783842756842877E-11, 9.21831518760500529508E-11,
+ -4.19035475934189648750E-10, 2.01504975519703286596E-9,
+ -1.03457624656780970260E-8, 5.74108412545004946722E-8,
+ -3.50196060308781257119E-7, 2.40648494783721712015E-6,
+ -1.93619797416608296024E-5, 1.95215518471351631108E-4,
+ -2.85781685962277938680E-3, 1.03923736576817238437E-1,
+ 2.72062619048444266945E0};
+ const T MAXNUM = pset1<T>(NumTraits<double>::infinity());
+ const T two = pset1<T>(2.0);
+ T x_le_two = pdiv(internal::pchebevl<T, 11>::run(
+ pmadd(x, x, pset1<T>(-2.0)), A), x);
+ x_le_two = pmadd(
+ generic_i1<T, double>::run(x), plog(pmul(pset1<T>(0.5), x)), x_le_two);
+ x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), MAXNUM, x_le_two);
+ T x_gt_two = pmul(
+ pexp(-x),
+ pmul(
+ internal::pchebevl<T, 25>::run(
+ psub(pdiv(pset1<T>(8.0), x), two), B),
+ prsqrt(x)));
+ return pselect(pcmp_le(x, two), x_le_two, x_gt_two);
+ }
+};
+
+template <typename Scalar>
+struct k1_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_k1<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct j0_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_j0 {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_j0<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* j0f.c
+ * Bessel function of order zero
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, j0f();
+ *
+ * y = j0f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of order zero of the argument.
+ *
+ * The domain is divided into the intervals [0, 2] and
+ * (2, infinity). In the first interval the following polynomial
+ * approximation is used:
+ *
+ *
+ * 2 2 2
+ * (w - r ) (w - r ) (w - r ) P(w)
+ * 1 2 3
+ *
+ * 2
+ * where w = x and the three r's are zeros of the function.
+ *
+ * In the second interval, the modulus and phase are approximated
+ * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
+ * and Phase(x) = x + 1/x R(1/x^2) - pi/4. The function is
+ *
+ * j0(x) = Modulus(x) cos( Phase(x) ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 2 100000 1.3e-7 3.6e-8
+ * IEEE 2, 32 100000 1.9e-7 5.4e-8
+ *
+ */
+
+ const float JP[] = {-6.068350350393235E-008f, 6.388945720783375E-006f,
+ -3.969646342510940E-004f, 1.332913422519003E-002f,
+ -1.729150680240724E-001f};
+ const float MO[] = {-6.838999669318810E-002f, 1.864949361379502E-001f,
+ -2.145007480346739E-001f, 1.197549369473540E-001f,
+ -3.560281861530129E-003f, -4.969382655296620E-002f,
+ -3.355424622293709E-006f, 7.978845717621440E-001f};
+ const float PH[] = {3.242077816988247E+001f, -3.630592630518434E+001f,
+ 1.756221482109099E+001f, -4.974978466280903E+000f,
+ 1.001973420681837E+000f, -1.939906941791308E-001f,
+ 6.490598792654666E-002f, -1.249992184872738E-001f};
+ const T DR1 = pset1<T>(5.78318596294678452118f);
+ const T NEG_PIO4F = pset1<T>(-0.7853981633974483096f); /* -pi / 4 */
+ T y = pabs(x);
+ T z = pmul(y, y);
+ T y_le_two = pselect(
+ pcmp_lt(y, pset1<T>(1.0e-3f)),
+ pmadd(z, pset1<T>(-0.25f), pset1<T>(1.0f)),
+ pmul(psub(z, DR1), internal::ppolevl<T, 4>::run(z, JP)));
+ T q = pdiv(pset1<T>(1.0f), y);
+ T w = prsqrt(y);
+ T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO));
+ w = pmul(q, q);
+ T yn = pmadd(q, internal::ppolevl<T, 7>::run(w, PH), NEG_PIO4F);
+ T y_gt_two = pmul(p, pcos(padd(yn, y)));
+ return pselect(pcmp_le(y, pset1<T>(2.0)), y_le_two, y_gt_two);
+ }
+};
+
+template <typename T>
+struct generic_j0<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* j0.c
+ * Bessel function of order zero
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, j0();
+ *
+ * y = j0( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of order zero of the argument.
+ *
+ * The domain is divided into the intervals [0, 5] and
+ * (5, infinity). In the first interval the following rational
+ * approximation is used:
+ *
+ *
+ * 2 2
+ * (w - r ) (w - r ) P (w) / Q (w)
+ * 1 2 3 8
+ *
+ * 2
+ * where w = x and the two r's are zeros of the function.
+ *
+ * In the second interval, the Hankel asymptotic expansion
+ * is employed with two rational functions of degree 6/6
+ * and 7/7.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error:
+ * arithmetic domain # trials peak rms
+ * DEC 0, 30 10000 4.4e-17 6.3e-18
+ * IEEE 0, 30 60000 4.2e-16 1.1e-16
+ *
+ */
+ const double PP[] = {7.96936729297347051624E-4, 8.28352392107440799803E-2,
+ 1.23953371646414299388E0, 5.44725003058768775090E0,
+ 8.74716500199817011941E0, 5.30324038235394892183E0,
+ 9.99999999999999997821E-1};
+ const double PQ[] = {9.24408810558863637013E-4, 8.56288474354474431428E-2,
+ 1.25352743901058953537E0, 5.47097740330417105182E0,
+ 8.76190883237069594232E0, 5.30605288235394617618E0,
+ 1.00000000000000000218E0};
+ const double QP[] = {-1.13663838898469149931E-2, -1.28252718670509318512E0,
+ -1.95539544257735972385E1, -9.32060152123768231369E1,
+ -1.77681167980488050595E2, -1.47077505154951170175E2,
+ -5.14105326766599330220E1, -6.05014350600728481186E0};
+ const double QQ[] = {1.00000000000000000000E0, 6.43178256118178023184E1,
+ 8.56430025976980587198E2, 3.88240183605401609683E3,
+ 7.24046774195652478189E3, 5.93072701187316984827E3,
+ 2.06209331660327847417E3, 2.42005740240291393179E2};
+ const double RP[] = {-4.79443220978201773821E9, 1.95617491946556577543E12,
+ -2.49248344360967716204E14, 9.70862251047306323952E15};
+ const double RQ[] = {1.00000000000000000000E0, 4.99563147152651017219E2,
+ 1.73785401676374683123E5, 4.84409658339962045305E7,
+ 1.11855537045356834862E10, 2.11277520115489217587E12,
+ 3.10518229857422583814E14, 3.18121955943204943306E16,
+ 1.71086294081043136091E18};
+ const T DR1 = pset1<T>(5.78318596294678452118E0);
+ const T DR2 = pset1<T>(3.04712623436620863991E1);
+ const T SQ2OPI = pset1<T>(7.9788456080286535587989E-1); /* sqrt(2 / pi) */
+ const T NEG_PIO4 = pset1<T>(-0.7853981633974483096); /* pi / 4 */
+
+ T y = pabs(x);
+ T z = pmul(y, y);
+ T y_le_five = pselect(
+ pcmp_lt(y, pset1<T>(1.0e-5)),
+ pmadd(z, pset1<T>(-0.25), pset1<T>(1.0)),
+ pmul(pmul(psub(z, DR1), psub(z, DR2)),
+ pdiv(internal::ppolevl<T, 3>::run(z, RP),
+ internal::ppolevl<T, 8>::run(z, RQ))));
+ T s = pdiv(pset1<T>(25.0), z);
+ T p = pdiv(
+ internal::ppolevl<T, 6>::run(s, PP),
+ internal::ppolevl<T, 6>::run(s, PQ));
+ T q = pdiv(
+ internal::ppolevl<T, 7>::run(s, QP),
+ internal::ppolevl<T, 7>::run(s, QQ));
+ T yn = padd(y, NEG_PIO4);
+ T w = pdiv(pset1<T>(-5.0), y);
+ p = pmadd(p, pcos(yn), pmul(w, pmul(q, psin(yn))));
+ T y_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(y)));
+ return pselect(pcmp_le(y, pset1<T>(5.0)), y_le_five, y_gt_five);
+ }
+};
+
+template <typename Scalar>
+struct j0_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_j0<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct y0_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_y0 {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_y0<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* j0f.c
+ * Bessel function of the second kind, order zero
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, y0f();
+ *
+ * y = y0f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of the second kind, of order
+ * zero, of the argument.
+ *
+ * The domain is divided into the intervals [0, 2] and
+ * (2, infinity). In the first interval a rational approximation
+ * R(x) is employed to compute
+ *
+ * 2 2 2
+ * y0(x) = (w - r ) (w - r ) (w - r ) R(x) + 2/pi ln(x) j0(x).
+ * 1 2 3
+ *
+ * Thus a call to j0() is required. The three zeros are removed
+ * from R(x) to improve its numerical stability.
+ *
+ * In the second interval, the modulus and phase are approximated
+ * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
+ * and Phase(x) = x + 1/x S(1/x^2) - pi/4. Then the function is
+ *
+ * y0(x) = Modulus(x) sin( Phase(x) ).
+ *
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error, when y0(x) < 1; else relative error:
+ *
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 2 100000 2.4e-7 3.4e-8
+ * IEEE 2, 32 100000 1.8e-7 5.3e-8
+ *
+ */
+
+ const float YP[] = {9.454583683980369E-008f, -9.413212653797057E-006f,
+ 5.344486707214273E-004f, -1.584289289821316E-002f,
+ 1.707584643733568E-001f};
+ const float MO[] = {-6.838999669318810E-002f, 1.864949361379502E-001f,
+ -2.145007480346739E-001f, 1.197549369473540E-001f,
+ -3.560281861530129E-003f, -4.969382655296620E-002f,
+ -3.355424622293709E-006f, 7.978845717621440E-001f};
+ const float PH[] = {3.242077816988247E+001f, -3.630592630518434E+001f,
+ 1.756221482109099E+001f, -4.974978466280903E+000f,
+ 1.001973420681837E+000f, -1.939906941791308E-001f,
+ 6.490598792654666E-002f, -1.249992184872738E-001f};
+ const T YZ1 = pset1<T>(0.43221455686510834878f);
+ const T TWOOPI = pset1<T>(0.636619772367581343075535f); /* 2 / pi */
+ const T NEG_PIO4F = pset1<T>(-0.7853981633974483096f); /* -pi / 4 */
+ const T NEG_MAXNUM = pset1<T>(-NumTraits<float>::infinity());
+ T z = pmul(x, x);
+ T x_le_two = pmul(TWOOPI, pmul(plog(x), generic_j0<T, float>::run(x)));
+ x_le_two = pmadd(
+ psub(z, YZ1), internal::ppolevl<T, 4>::run(z, YP), x_le_two);
+ x_le_two = pselect(pcmp_le(x, pset1<T>(0.0)), NEG_MAXNUM, x_le_two);
+ T q = pdiv(pset1<T>(1.0), x);
+ T w = prsqrt(x);
+ T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO));
+ T u = pmul(q, q);
+ T xn = pmadd(q, internal::ppolevl<T, 7>::run(u, PH), NEG_PIO4F);
+ T x_gt_two = pmul(p, psin(padd(xn, x)));
+ return pselect(pcmp_le(x, pset1<T>(2.0)), x_le_two, x_gt_two);
+ }
+};
+
+template <typename T>
+struct generic_y0<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* j0.c
+ * Bessel function of the second kind, order zero
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, y0();
+ *
+ * y = y0( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of the second kind, of order
+ * zero, of the argument.
+ *
+ * The domain is divided into the intervals [0, 5] and
+ * (5, infinity). In the first interval a rational approximation
+ * R(x) is employed to compute
+ * y0(x) = R(x) + 2 * log(x) * j0(x) / PI.
+ * Thus a call to j0() is required.
+ *
+ * In the second interval, the Hankel asymptotic expansion
+ * is employed with two rational functions of degree 6/6
+ * and 7/7.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error, when y0(x) < 1; else relative error:
+ *
+ * arithmetic domain # trials peak rms
+ * DEC 0, 30 9400 7.0e-17 7.9e-18
+ * IEEE 0, 30 30000 1.3e-15 1.6e-16
+ *
+ */
+ const double PP[] = {7.96936729297347051624E-4, 8.28352392107440799803E-2,
+ 1.23953371646414299388E0, 5.44725003058768775090E0,
+ 8.74716500199817011941E0, 5.30324038235394892183E0,
+ 9.99999999999999997821E-1};
+ const double PQ[] = {9.24408810558863637013E-4, 8.56288474354474431428E-2,
+ 1.25352743901058953537E0, 5.47097740330417105182E0,
+ 8.76190883237069594232E0, 5.30605288235394617618E0,
+ 1.00000000000000000218E0};
+ const double QP[] = {-1.13663838898469149931E-2, -1.28252718670509318512E0,
+ -1.95539544257735972385E1, -9.32060152123768231369E1,
+ -1.77681167980488050595E2, -1.47077505154951170175E2,
+ -5.14105326766599330220E1, -6.05014350600728481186E0};
+ const double QQ[] = {1.00000000000000000000E0, 6.43178256118178023184E1,
+ 8.56430025976980587198E2, 3.88240183605401609683E3,
+ 7.24046774195652478189E3, 5.93072701187316984827E3,
+ 2.06209331660327847417E3, 2.42005740240291393179E2};
+ const double YP[] = {1.55924367855235737965E4, -1.46639295903971606143E7,
+ 5.43526477051876500413E9, -9.82136065717911466409E11,
+ 8.75906394395366999549E13, -3.46628303384729719441E15,
+ 4.42733268572569800351E16, -1.84950800436986690637E16};
+ const double YQ[] = {1.00000000000000000000E0, 1.04128353664259848412E3,
+ 6.26107330137134956842E5, 2.68919633393814121987E8,
+ 8.64002487103935000337E10, 2.02979612750105546709E13,
+ 3.17157752842975028269E15, 2.50596256172653059228E17};
+ const T SQ2OPI = pset1<T>(7.9788456080286535587989E-1); /* sqrt(2 / pi) */
+ const T TWOOPI = pset1<T>(0.636619772367581343075535); /* 2 / pi */
+ const T NEG_PIO4 = pset1<T>(-0.7853981633974483096); /* -pi / 4 */
+ const T NEG_MAXNUM = pset1<T>(-NumTraits<double>::infinity());
+
+ T z = pmul(x, x);
+ T x_le_five = pdiv(internal::ppolevl<T, 7>::run(z, YP),
+ internal::ppolevl<T, 7>::run(z, YQ));
+ x_le_five = pmadd(
+ pmul(TWOOPI, plog(x)), generic_j0<T, double>::run(x), x_le_five);
+ x_le_five = pselect(pcmp_le(x, pset1<T>(0.0)), NEG_MAXNUM, x_le_five);
+ T s = pdiv(pset1<T>(25.0), z);
+ T p = pdiv(
+ internal::ppolevl<T, 6>::run(s, PP),
+ internal::ppolevl<T, 6>::run(s, PQ));
+ T q = pdiv(
+ internal::ppolevl<T, 7>::run(s, QP),
+ internal::ppolevl<T, 7>::run(s, QQ));
+ T xn = padd(x, NEG_PIO4);
+ T w = pdiv(pset1<T>(5.0), x);
+ p = pmadd(p, psin(xn), pmul(w, pmul(q, pcos(xn))));
+ T x_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(x)));
+ return pselect(pcmp_le(x, pset1<T>(5.0)), x_le_five, x_gt_five);
+ }
+};
+
+template <typename Scalar>
+struct y0_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_y0<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct j1_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_j1 {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_j1<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* j1f.c
+ * Bessel function of order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, j1f();
+ *
+ * y = j1f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of order one of the argument.
+ *
+ * The domain is divided into the intervals [0, 2] and
+ * (2, infinity). In the first interval a polynomial approximation
+ * 2
+ * (w - r ) x P(w)
+ * 1
+ * 2
+ * is used, where w = x and r is the first zero of the function.
+ *
+ * In the second interval, the modulus and phase are approximated
+ * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
+ * and Phase(x) = x + 1/x R(1/x^2) - 3pi/4. The function is
+ *
+ * j0(x) = Modulus(x) cos( Phase(x) ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 2 100000 1.2e-7 2.5e-8
+ * IEEE 2, 32 100000 2.0e-7 5.3e-8
+ *
+ *
+ */
+
+ const float JP[] = {-4.878788132172128E-009f, 6.009061827883699E-007f,
+ -4.541343896997497E-005f, 1.937383947804541E-003f,
+ -3.405537384615824E-002f};
+ const float MO1[] = {6.913942741265801E-002f, -2.284801500053359E-001f,
+ 3.138238455499697E-001f, -2.102302420403875E-001f,
+ 5.435364690523026E-003f, 1.493389585089498E-001f,
+ 4.976029650847191E-006f, 7.978845453073848E-001f};
+ const float PH1[] = {-4.497014141919556E+001f, 5.073465654089319E+001f,
+ -2.485774108720340E+001f, 7.222973196770240E+000f,
+ -1.544842782180211E+000f, 3.503787691653334E-001f,
+ -1.637986776941202E-001f, 3.749989509080821E-001f};
+ const T Z1 = pset1<T>(1.46819706421238932572E1f);
+ const T NEG_THPIO4F = pset1<T>(-2.35619449019234492885f); /* -3*pi/4 */
+
+ T y = pabs(x);
+ T z = pmul(y, y);
+ T y_le_two = pmul(
+ psub(z, Z1),
+ pmul(x, internal::ppolevl<T, 4>::run(z, JP)));
+ T q = pdiv(pset1<T>(1.0f), y);
+ T w = prsqrt(y);
+ T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO1));
+ w = pmul(q, q);
+ T yn = pmadd(q, internal::ppolevl<T, 7>::run(w, PH1), NEG_THPIO4F);
+ T y_gt_two = pmul(p, pcos(padd(yn, y)));
+ // j1 is an odd function. This implementation differs from cephes to
+ // take this fact in to account. Cephes returns -j1(x) for y > 2 range.
+ y_gt_two = pselect(
+ pcmp_lt(x, pset1<T>(0.0f)), pmul(pset1<T>(-1.0f), y_gt_two), y_gt_two);
+ return pselect(pcmp_le(y, pset1<T>(2.0f)), y_le_two, y_gt_two);
+ }
+};
+
+template <typename T>
+struct generic_j1<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* j1.c
+ * Bessel function of order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, j1();
+ *
+ * y = j1( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of order one of the argument.
+ *
+ * The domain is divided into the intervals [0, 8] and
+ * (8, infinity). In the first interval a 24 term Chebyshev
+ * expansion is used. In the second, the asymptotic
+ * trigonometric representation is employed using two
+ * rational functions of degree 5/5.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error:
+ * arithmetic domain # trials peak rms
+ * DEC 0, 30 10000 4.0e-17 1.1e-17
+ * IEEE 0, 30 30000 2.6e-16 1.1e-16
+ *
+ */
+ const double PP[] = {7.62125616208173112003E-4, 7.31397056940917570436E-2,
+ 1.12719608129684925192E0, 5.11207951146807644818E0,
+ 8.42404590141772420927E0, 5.21451598682361504063E0,
+ 1.00000000000000000254E0};
+ const double PQ[] = {5.71323128072548699714E-4, 6.88455908754495404082E-2,
+ 1.10514232634061696926E0, 5.07386386128601488557E0,
+ 8.39985554327604159757E0, 5.20982848682361821619E0,
+ 9.99999999999999997461E-1};
+ const double QP[] = {5.10862594750176621635E-2, 4.98213872951233449420E0,
+ 7.58238284132545283818E1, 3.66779609360150777800E2,
+ 7.10856304998926107277E2, 5.97489612400613639965E2,
+ 2.11688757100572135698E2, 2.52070205858023719784E1};
+ const double QQ[] = {1.00000000000000000000E0, 7.42373277035675149943E1,
+ 1.05644886038262816351E3, 4.98641058337653607651E3,
+ 9.56231892404756170795E3, 7.99704160447350683650E3,
+ 2.82619278517639096600E3, 3.36093607810698293419E2};
+ const double RP[] = {-8.99971225705559398224E8, 4.52228297998194034323E11,
+ -7.27494245221818276015E13, 3.68295732863852883286E15};
+ const double RQ[] = {1.00000000000000000000E0, 6.20836478118054335476E2,
+ 2.56987256757748830383E5, 8.35146791431949253037E7,
+ 2.21511595479792499675E10, 4.74914122079991414898E12,
+ 7.84369607876235854894E14, 8.95222336184627338078E16,
+ 5.32278620332680085395E18};
+ const T Z1 = pset1<T>(1.46819706421238932572E1);
+ const T Z2 = pset1<T>(4.92184563216946036703E1);
+ const T NEG_THPIO4 = pset1<T>(-2.35619449019234492885); /* -3*pi/4 */
+ const T SQ2OPI = pset1<T>(7.9788456080286535587989E-1); /* sqrt(2 / pi) */
+ T y = pabs(x);
+ T z = pmul(y, y);
+ T y_le_five = pdiv(internal::ppolevl<T, 3>::run(z, RP),
+ internal::ppolevl<T, 8>::run(z, RQ));
+ y_le_five = pmul(pmul(pmul(y_le_five, x), psub(z, Z1)), psub(z, Z2));
+ T s = pdiv(pset1<T>(25.0), z);
+ T p = pdiv(
+ internal::ppolevl<T, 6>::run(s, PP),
+ internal::ppolevl<T, 6>::run(s, PQ));
+ T q = pdiv(
+ internal::ppolevl<T, 7>::run(s, QP),
+ internal::ppolevl<T, 7>::run(s, QQ));
+ T yn = padd(y, NEG_THPIO4);
+ T w = pdiv(pset1<T>(-5.0), y);
+ p = pmadd(p, pcos(yn), pmul(w, pmul(q, psin(yn))));
+ T y_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(y)));
+ // j1 is an odd function. This implementation differs from cephes to
+ // take this fact in to account. Cephes returns -j1(x) for y > 5 range.
+ y_gt_five = pselect(
+ pcmp_lt(x, pset1<T>(0.0f)), pmul(pset1<T>(-1.0), y_gt_five), y_gt_five);
+ return pselect(pcmp_le(y, pset1<T>(5.0)), y_le_five, y_gt_five);
+ }
+};
+
+template <typename Scalar>
+struct j1_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_j1<Scalar, Scalar>::run(x);
+ }
+};
+
+template <typename Scalar>
+struct y1_retval {
+ typedef Scalar type;
+};
+
+template <typename T, typename ScalarType>
+struct generic_y1 {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T&) {
+ EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
+ THIS_TYPE_IS_NOT_SUPPORTED);
+ return ScalarType(0);
+ }
+};
+
+template <typename T>
+struct generic_y1<T, float> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* j1f.c
+ * Bessel function of second kind of order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, y1();
+ *
+ * y = y1( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of the second kind of order one
+ * of the argument.
+ *
+ * The domain is divided into the intervals [0, 2] and
+ * (2, infinity). In the first interval a rational approximation
+ * R(x) is employed to compute
+ *
+ * 2
+ * y0(x) = (w - r ) x R(x^2) + 2/pi (ln(x) j1(x) - 1/x) .
+ * 1
+ *
+ * Thus a call to j1() is required.
+ *
+ * In the second interval, the modulus and phase are approximated
+ * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
+ * and Phase(x) = x + 1/x S(1/x^2) - 3pi/4. Then the function is
+ *
+ * y0(x) = Modulus(x) sin( Phase(x) ).
+ *
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 2 100000 2.2e-7 4.6e-8
+ * IEEE 2, 32 100000 1.9e-7 5.3e-8
+ *
+ * (error criterion relative when |y1| > 1).
+ *
+ */
+
+ const float YP[] = {8.061978323326852E-009f, -9.496460629917016E-007f,
+ 6.719543806674249E-005f, -2.641785726447862E-003f,
+ 4.202369946500099E-002f};
+ const float MO1[] = {6.913942741265801E-002f, -2.284801500053359E-001f,
+ 3.138238455499697E-001f, -2.102302420403875E-001f,
+ 5.435364690523026E-003f, 1.493389585089498E-001f,
+ 4.976029650847191E-006f, 7.978845453073848E-001f};
+ const float PH1[] = {-4.497014141919556E+001f, 5.073465654089319E+001f,
+ -2.485774108720340E+001f, 7.222973196770240E+000f,
+ -1.544842782180211E+000f, 3.503787691653334E-001f,
+ -1.637986776941202E-001f, 3.749989509080821E-001f};
+ const T YO1 = pset1<T>(4.66539330185668857532f);
+ const T NEG_THPIO4F = pset1<T>(-2.35619449019234492885f); /* -3*pi/4 */
+ const T TWOOPI = pset1<T>(0.636619772367581343075535f); /* 2/pi */
+ const T NEG_MAXNUM = pset1<T>(-NumTraits<float>::infinity());
+
+ T z = pmul(x, x);
+ T x_le_two = pmul(psub(z, YO1), internal::ppolevl<T, 4>::run(z, YP));
+ x_le_two = pmadd(
+ x_le_two, x,
+ pmul(TWOOPI, pmadd(
+ generic_j1<T, float>::run(x), plog(x),
+ pdiv(pset1<T>(-1.0f), x))));
+ x_le_two = pselect(pcmp_lt(x, pset1<T>(0.0f)), NEG_MAXNUM, x_le_two);
+
+ T q = pdiv(pset1<T>(1.0), x);
+ T w = prsqrt(x);
+ T p = pmul(w, internal::ppolevl<T, 7>::run(q, MO1));
+ w = pmul(q, q);
+ T xn = pmadd(q, internal::ppolevl<T, 7>::run(w, PH1), NEG_THPIO4F);
+ T x_gt_two = pmul(p, psin(padd(xn, x)));
+ return pselect(pcmp_le(x, pset1<T>(2.0)), x_le_two, x_gt_two);
+ }
+};
+
+template <typename T>
+struct generic_y1<T, double> {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE T run(const T& x) {
+ /* j1.c
+ * Bessel function of second kind of order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, y1();
+ *
+ * y = y1( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of the second kind of order one
+ * of the argument.
+ *
+ * The domain is divided into the intervals [0, 8] and
+ * (8, infinity). In the first interval a 25 term Chebyshev
+ * expansion is used, and a call to j1() is required.
+ * In the second, the asymptotic trigonometric representation
+ * is employed using two rational functions of degree 5/5.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error:
+ * arithmetic domain # trials peak rms
+ * DEC 0, 30 10000 8.6e-17 1.3e-17
+ * IEEE 0, 30 30000 1.0e-15 1.3e-16
+ *
+ * (error criterion relative when |y1| > 1).
+ *
+ */
+ const double PP[] = {7.62125616208173112003E-4, 7.31397056940917570436E-2,
+ 1.12719608129684925192E0, 5.11207951146807644818E0,
+ 8.42404590141772420927E0, 5.21451598682361504063E0,
+ 1.00000000000000000254E0};
+ const double PQ[] = {5.71323128072548699714E-4, 6.88455908754495404082E-2,
+ 1.10514232634061696926E0, 5.07386386128601488557E0,
+ 8.39985554327604159757E0, 5.20982848682361821619E0,
+ 9.99999999999999997461E-1};
+ const double QP[] = {5.10862594750176621635E-2, 4.98213872951233449420E0,
+ 7.58238284132545283818E1, 3.66779609360150777800E2,
+ 7.10856304998926107277E2, 5.97489612400613639965E2,
+ 2.11688757100572135698E2, 2.52070205858023719784E1};
+ const double QQ[] = {1.00000000000000000000E0, 7.42373277035675149943E1,
+ 1.05644886038262816351E3, 4.98641058337653607651E3,
+ 9.56231892404756170795E3, 7.99704160447350683650E3,
+ 2.82619278517639096600E3, 3.36093607810698293419E2};
+ const double YP[] = {1.26320474790178026440E9, -6.47355876379160291031E11,
+ 1.14509511541823727583E14, -8.12770255501325109621E15,
+ 2.02439475713594898196E17, -7.78877196265950026825E17};
+ const double YQ[] = {1.00000000000000000000E0, 5.94301592346128195359E2,
+ 2.35564092943068577943E5, 7.34811944459721705660E7,
+ 1.87601316108706159478E10, 3.88231277496238566008E12,
+ 6.20557727146953693363E14, 6.87141087355300489866E16,
+ 3.97270608116560655612E18};
+ const T SQ2OPI = pset1<T>(.79788456080286535588);
+ const T NEG_THPIO4 = pset1<T>(-2.35619449019234492885); /* -3*pi/4 */
+ const T TWOOPI = pset1<T>(0.636619772367581343075535); /* 2/pi */
+ const T NEG_MAXNUM = pset1<T>(-NumTraits<double>::infinity());
+
+ T z = pmul(x, x);
+ T x_le_five = pdiv(internal::ppolevl<T, 5>::run(z, YP),
+ internal::ppolevl<T, 8>::run(z, YQ));
+ x_le_five = pmadd(
+ x_le_five, x, pmul(
+ TWOOPI, pmadd(generic_j1<T, double>::run(x), plog(x),
+ pdiv(pset1<T>(-1.0), x))));
+
+ x_le_five = pselect(pcmp_le(x, pset1<T>(0.0)), NEG_MAXNUM, x_le_five);
+ T s = pdiv(pset1<T>(25.0), z);
+ T p = pdiv(
+ internal::ppolevl<T, 6>::run(s, PP),
+ internal::ppolevl<T, 6>::run(s, PQ));
+ T q = pdiv(
+ internal::ppolevl<T, 7>::run(s, QP),
+ internal::ppolevl<T, 7>::run(s, QQ));
+ T xn = padd(x, NEG_THPIO4);
+ T w = pdiv(pset1<T>(5.0), x);
+ p = pmadd(p, psin(xn), pmul(w, pmul(q, pcos(xn))));
+ T x_gt_five = pmul(p, pmul(SQ2OPI, prsqrt(x)));
+ return pselect(pcmp_le(x, pset1<T>(5.0)), x_le_five, x_gt_five);
+ }
+};
+
+template <typename Scalar>
+struct y1_impl {
+ EIGEN_DEVICE_FUNC
+ static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
+ return generic_y1<Scalar, Scalar>::run(x);
+ }
+};
+
+} // end namespace internal
+
+namespace numext {
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(i0, Scalar)
+ i0(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(i0, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(i0e, Scalar)
+ i0e(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(i0e, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(i1, Scalar)
+ i1(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(i1, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(i1e, Scalar)
+ i1e(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(i1e, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(k0, Scalar)
+ k0(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(k0, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(k0e, Scalar)
+ k0e(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(k0e, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(k1, Scalar)
+ k1(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(k1, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(k1e, Scalar)
+ k1e(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(k1e, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(j0, Scalar)
+ j0(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(j0, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(y0, Scalar)
+ y0(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(y0, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(j1, Scalar)
+ j1(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(j1, Scalar)::run(x);
+}
+
+template <typename Scalar>
+EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(y1, Scalar)
+ y1(const Scalar& x) {
+ return EIGEN_MATHFUNC_IMPL(y1, Scalar)::run(x);
+}
+
+} // end namespace numext
+
+} // end namespace Eigen
+
+#endif // EIGEN_BESSEL_FUNCTIONS_H
diff --git a/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsPacketMath.h b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsPacketMath.h
new file mode 100644
index 000000000..70eaad5cd
--- /dev/null
+++ b/unsupported/Eigen/src/SpecialFunctions/BesselFunctionsPacketMath.h
@@ -0,0 +1,130 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_BESSELFUNCTIONS_PACKETMATH_H
+#define EIGEN_BESSELFUNCTIONS_PACKETMATH_H
+
+namespace Eigen {
+
+namespace internal {
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order zero i0(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pi0(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_i0; return generic_i0<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order zero i0e(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pi0e(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_i0e; return generic_i0e<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order one i1(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pi1(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_i1; return generic_i1<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order one i1e(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pi1e(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_i1e; return generic_i1e<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order zero j0(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pj0(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_j0; return generic_j0<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order zero j1(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pj1(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_j1; return generic_j1<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order one y0(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet py0(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_y0; return generic_y0<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order one y1(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet py1(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_y1; return generic_y1<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order zero k0(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pk0(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_k0; return generic_k0<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order zero k0e(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pk0e(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_k0e; return generic_k0e<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order one k1e(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pk1(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_k1; return generic_k1<Packet, ScalarType>::run(x);
+}
+
+/** \internal \returns the exponentially scaled modified Bessel function of
+ * order one k1e(\a a) (coeff-wise) */
+template <typename Packet>
+EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
+Packet pk1e(const Packet& x) {
+ typedef typename unpacket_traits<Packet>::type ScalarType;
+ using internal::generic_k1e; return generic_k1e<Packet, ScalarType>::run(x);
+}
+
+} // end namespace internal
+
+} // end namespace Eigen
+
+#endif // EIGEN_BESSELFUNCTIONS_PACKETMATH_H
+
diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h
index 617401e9d..691ff4d03 100644
--- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h
+++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsArrayAPI.h
@@ -161,51 +161,6 @@ zeta(const Eigen::ArrayBase<DerivedX>& x, const Eigen::ArrayBase<DerivedQ>& q)
);
}
-/** \returns an expression of the coefficient-wise i0e(\a x) to the given
- * arrays.
- *
- * It returns the exponentially scaled modified Bessel
- * function of order zero.
- *
- * \param x is the argument
- *
- * \note This function supports only float and double scalar types. To support
- * other scalar types, the user has to provide implementations of i0e(T) for
- * any scalar type T to be supported.
- *
- * \sa ArrayBase::i0e()
- */
-template <typename Derived>
-EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
- Eigen::internal::scalar_i0e_op<typename Derived::Scalar>, const Derived>
-i0e(const Eigen::ArrayBase<Derived>& x) {
- return Eigen::CwiseUnaryOp<
- Eigen::internal::scalar_i0e_op<typename Derived::Scalar>,
- const Derived>(x.derived());
-}
-
-/** \returns an expression of the coefficient-wise i1e(\a x) to the given
- * arrays.
- *
- * It returns the exponentially scaled modified Bessel
- * function of order one.
- *
- * \param x is the argument
- *
- * \note This function supports only float and double scalar types. To support
- * other scalar types, the user has to provide implementations of i1e(T) for
- * any scalar type T to be supported.
- *
- * \sa ArrayBase::i1e()
- */
-template <typename Derived>
-EIGEN_STRONG_INLINE const Eigen::CwiseUnaryOp<
- Eigen::internal::scalar_i1e_op<typename Derived::Scalar>, const Derived>
-i1e(const Eigen::ArrayBase<Derived>& x) {
- return Eigen::CwiseUnaryOp<
- Eigen::internal::scalar_i1e_op<typename Derived::Scalar>,
- const Derived>(x.derived());
-}
} // end namespace Eigen
diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h
index 13a72a3ee..a4287c31f 100644
--- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h
+++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsFunctors.h
@@ -308,60 +308,6 @@ struct functor_traits<scalar_ndtri_op<Scalar> >
};
};
-/** \internal
- * \brief Template functor to compute the exponentially scaled modified Bessel
- * function of order zero
- * \sa class CwiseUnaryOp, Cwise::i0e()
- */
-template <typename Scalar>
-struct scalar_i0e_op {
- EIGEN_EMPTY_STRUCT_CTOR(scalar_i0e_op)
- EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
- using numext::i0e;
- return i0e(x);
- }
- typedef typename packet_traits<Scalar>::type Packet;
- EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
- return internal::pi0e(x);
- }
-};
-template <typename Scalar>
-struct functor_traits<scalar_i0e_op<Scalar> > {
- enum {
- // On average, a Chebyshev polynomial of order N=20 is computed.
- // The cost is N multiplications and 2N additions.
- Cost = 20 * NumTraits<Scalar>::MulCost + 40 * NumTraits<Scalar>::AddCost,
- PacketAccess = packet_traits<Scalar>::HasI0e
- };
-};
-
-/** \internal
- * \brief Template functor to compute the exponentially scaled modified Bessel
- * function of order zero
- * \sa class CwiseUnaryOp, Cwise::i1e()
- */
-template <typename Scalar>
-struct scalar_i1e_op {
- EIGEN_EMPTY_STRUCT_CTOR(scalar_i1e_op)
- EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& x) const {
- using numext::i1e;
- return i1e(x);
- }
- typedef typename packet_traits<Scalar>::type Packet;
- EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
- return internal::pi1e(x);
- }
-};
-template <typename Scalar>
-struct functor_traits<scalar_i1e_op<Scalar> > {
- enum {
- // On average, a Chebyshev polynomial of order N=20 is computed.
- // The cost is N multiplications and 2N additions.
- Cost = 20 * NumTraits<Scalar>::MulCost + 40 * NumTraits<Scalar>::AddCost,
- PacketAccess = packet_traits<Scalar>::HasI1e
- };
-};
-
} // end namespace internal
} // end namespace Eigen
diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h
index 538db2afa..2a3a53168 100644
--- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h
+++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsHalf.h
@@ -50,14 +50,6 @@ template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half igammac(const Eigen
template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half betainc(const Eigen::half& a, const Eigen::half& b, const Eigen::half& x) {
return Eigen::half(Eigen::numext::betainc(static_cast<float>(a), static_cast<float>(b), static_cast<float>(x)));
}
-template <>
-EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half i0e(const Eigen::half& x) {
- return Eigen::half(Eigen::numext::i0e(static_cast<float>(x)));
-}
-template <>
-EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC Eigen::half i1e(const Eigen::half& x) {
- return Eigen::half(Eigen::numext::i1e(static_cast<float>(x)));
-}
#endif
} // end namespace numext
diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
index 7c6d32049..ea00bd96e 100644
--- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
+++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsImpl.h
@@ -1757,7 +1757,7 @@ struct betainc_helper<double> {
if ((a + b) < maxgam && numext::abs(u) < maxlog) {
t = gamma(a + b) / (gamma(a) * gamma(b));
s = s * t * pow(x, a);
- } else {
+ }
*/
t = lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) -
lgamma_impl<double>::run(b) + u + numext::log(s);
@@ -1864,351 +1864,6 @@ struct betainc_impl<double> {
#endif // EIGEN_HAS_C99_MATH
-/****************************************************************************
- * Implementation of Bessel function, based on Cephes *
- ****************************************************************************/
-
-template <typename Scalar>
-struct i0e_retval {
- typedef Scalar type;
-};
-
-template <typename T, typename ScalarType>
-struct generic_i0e {
- EIGEN_DEVICE_FUNC
- static EIGEN_STRONG_INLINE T run(const T&) {
- EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
- THIS_TYPE_IS_NOT_SUPPORTED);
- return ScalarType(0);
- }
-};
-
-template <typename T>
-struct generic_i0e<T, float> {
- EIGEN_DEVICE_FUNC
- static EIGEN_STRONG_INLINE T run(const T& x) {
- /* i0ef.c
- *
- * Modified Bessel function of order zero,
- * exponentially scaled
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, i0ef();
- *
- * y = i0ef( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns exponentially scaled modified Bessel function
- * of order zero of the argument.
- *
- * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
- *
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * IEEE 0,30 100000 3.7e-7 7.0e-8
- * See i0f().
- *
- */
-
- const float A[] = {-1.30002500998624804212E-8f, 6.04699502254191894932E-8f,
- -2.67079385394061173391E-7f, 1.11738753912010371815E-6f,
- -4.41673835845875056359E-6f, 1.64484480707288970893E-5f,
- -5.75419501008210370398E-5f, 1.88502885095841655729E-4f,
- -5.76375574538582365885E-4f, 1.63947561694133579842E-3f,
- -4.32430999505057594430E-3f, 1.05464603945949983183E-2f,
- -2.37374148058994688156E-2f, 4.93052842396707084878E-2f,
- -9.49010970480476444210E-2f, 1.71620901522208775349E-1f,
- -3.04682672343198398683E-1f, 6.76795274409476084995E-1f};
-
- const float B[] = {3.39623202570838634515E-9f, 2.26666899049817806459E-8f,
- 2.04891858946906374183E-7f, 2.89137052083475648297E-6f,
- 6.88975834691682398426E-5f, 3.36911647825569408990E-3f,
- 8.04490411014108831608E-1f};
- T y = pabs(x);
- T y_le_eight = internal::pchebevl<T, 18>::run(
- pmadd(pset1<T>(0.5f), y, pset1<T>(-2.0f)), A);
- T y_gt_eight = pdiv(
- internal::pchebevl<T, 7>::run(
- psub(pdiv(pset1<T>(32.0f), y), pset1<T>(2.0f)), B),
- psqrt(y));
- // TODO: Perhaps instead check whether all packet elements are in
- // [-8, 8] and evaluate a branch based off of that. It's possible
- // in practice most elements are in this region.
- return pselect(pcmp_le(y, pset1<T>(8.0f)), y_le_eight, y_gt_eight);
- }
-};
-
-template <typename T>
-struct generic_i0e<T, double> {
- EIGEN_DEVICE_FUNC
- static EIGEN_STRONG_INLINE T run(const T& x) {
- /* i0e.c
- *
- * Modified Bessel function of order zero,
- * exponentially scaled
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, i0e();
- *
- * y = i0e( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns exponentially scaled modified Bessel function
- * of order zero of the argument.
- *
- * The function is defined as i0e(x) = exp(-|x|) j0( ix ).
- *
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * IEEE 0,30 30000 5.4e-16 1.2e-16
- * See i0().
- *
- */
-
- const double A[] = {-4.41534164647933937950E-18, 3.33079451882223809783E-17,
- -2.43127984654795469359E-16, 1.71539128555513303061E-15,
- -1.16853328779934516808E-14, 7.67618549860493561688E-14,
- -4.85644678311192946090E-13, 2.95505266312963983461E-12,
- -1.72682629144155570723E-11, 9.67580903537323691224E-11,
- -5.18979560163526290666E-10, 2.65982372468238665035E-9,
- -1.30002500998624804212E-8, 6.04699502254191894932E-8,
- -2.67079385394061173391E-7, 1.11738753912010371815E-6,
- -4.41673835845875056359E-6, 1.64484480707288970893E-5,
- -5.75419501008210370398E-5, 1.88502885095841655729E-4,
- -5.76375574538582365885E-4, 1.63947561694133579842E-3,
- -4.32430999505057594430E-3, 1.05464603945949983183E-2,
- -2.37374148058994688156E-2, 4.93052842396707084878E-2,
- -9.49010970480476444210E-2, 1.71620901522208775349E-1,
- -3.04682672343198398683E-1, 6.76795274409476084995E-1};
- const double B[] = {
- -7.23318048787475395456E-18, -4.83050448594418207126E-18,
- 4.46562142029675999901E-17, 3.46122286769746109310E-17,
- -2.82762398051658348494E-16, -3.42548561967721913462E-16,
- 1.77256013305652638360E-15, 3.81168066935262242075E-15,
- -9.55484669882830764870E-15, -4.15056934728722208663E-14,
- 1.54008621752140982691E-14, 3.85277838274214270114E-13,
- 7.18012445138366623367E-13, -1.79417853150680611778E-12,
- -1.32158118404477131188E-11, -3.14991652796324136454E-11,
- 1.18891471078464383424E-11, 4.94060238822496958910E-10,
- 3.39623202570838634515E-9, 2.26666899049817806459E-8,
- 2.04891858946906374183E-7, 2.89137052083475648297E-6,
- 6.88975834691682398426E-5, 3.36911647825569408990E-3,
- 8.04490411014108831608E-1};
- T y = pabs(x);
- T y_le_eight = internal::pchebevl<T, 30>::run(
- pmadd(pset1<T>(0.5), y, pset1<T>(-2.0)), A);
- T y_gt_eight = pdiv(
- internal::pchebevl<T, 25>::run(
- psub(pdiv(pset1<T>(32.0), y), pset1<T>(2.0)), B),
- psqrt(y));
- // TODO: Perhaps instead check whether all packet elements are in
- // [-8, 8] and evaluate a branch based off of that. It's possible
- // in practice most elements are in this region.
- return pselect(pcmp_le(y, pset1<T>(8.0)), y_le_eight, y_gt_eight);
- }
-};
-
-template <typename Scalar>
-struct i0e_impl {
- EIGEN_DEVICE_FUNC
- static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
- return generic_i0e<Scalar, Scalar>::run(x);
- }
-};
-
-
-template <typename Scalar>
-struct i1e_retval {
- typedef Scalar type;
-};
-
-template <typename T, typename ScalarType>
-struct generic_i1e {
- EIGEN_DEVICE_FUNC
- static EIGEN_STRONG_INLINE T run(const T&) {
- EIGEN_STATIC_ASSERT((internal::is_same<T, T>::value == false),
- THIS_TYPE_IS_NOT_SUPPORTED);
- return ScalarType(0);
- }
-};
-
-template <typename T>
-struct generic_i1e<T, float> {
- EIGEN_DEVICE_FUNC
- static EIGEN_STRONG_INLINE T run(const T& x) {
- /* i1ef.c
- *
- * Modified Bessel function of order one,
- * exponentially scaled
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, y, i1ef();
- *
- * y = i1ef( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns exponentially scaled modified Bessel function
- * of order one of the argument.
- *
- * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
- *
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * IEEE 0, 30 30000 1.5e-6 1.5e-7
- * See i1().
- *
- */
- const float A[] = {9.38153738649577178388E-9f, -4.44505912879632808065E-8f,
- 2.00329475355213526229E-7f, -8.56872026469545474066E-7f,
- 3.47025130813767847674E-6f, -1.32731636560394358279E-5f,
- 4.78156510755005422638E-5f, -1.61760815825896745588E-4f,
- 5.12285956168575772895E-4f, -1.51357245063125314899E-3f,
- 4.15642294431288815669E-3f, -1.05640848946261981558E-2f,
- 2.47264490306265168283E-2f, -5.29459812080949914269E-2f,
- 1.02643658689847095384E-1f, -1.76416518357834055153E-1f,
- 2.52587186443633654823E-1f};
-
- const float B[] = {-3.83538038596423702205E-9f, -2.63146884688951950684E-8f,
- -2.51223623787020892529E-7f, -3.88256480887769039346E-6f,
- -1.10588938762623716291E-4f, -9.76109749136146840777E-3f,
- 7.78576235018280120474E-1f};
-
-
- T y = pabs(x);
- T y_le_eight = pmul(y, internal::pchebevl<T, 17>::run(
- pmadd(pset1<T>(0.5f), y, pset1<T>(-2.0f)), A));
- T y_gt_eight = pdiv(
- internal::pchebevl<T, 7>::run(
- psub(pdiv(pset1<T>(32.0f), y),
- pset1<T>(2.0f)), B),
- psqrt(y));
- // TODO: Perhaps instead check whether all packet elements are in
- // [-8, 8] and evaluate a branch based off of that. It's possible
- // in practice most elements are in this region.
- y = pselect(pcmp_le(y, pset1<T>(8.0f)), y_le_eight, y_gt_eight);
- return pselect(pcmp_lt(x, pset1<T>(0.0f)), -y, y);
- }
-};
-
-template <typename T>
-struct generic_i1e<T, double> {
- EIGEN_DEVICE_FUNC
- static EIGEN_STRONG_INLINE T run(const T& x) {
- /* i1e.c
- *
- * Modified Bessel function of order one,
- * exponentially scaled
- *
- *
- *
- * SYNOPSIS:
- *
- * double x, y, i1e();
- *
- * y = i1e( x );
- *
- *
- *
- * DESCRIPTION:
- *
- * Returns exponentially scaled modified Bessel function
- * of order one of the argument.
- *
- * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
- *
- *
- *
- * ACCURACY:
- *
- * Relative error:
- * arithmetic domain # trials peak rms
- * IEEE 0, 30 30000 2.0e-15 2.0e-16
- * See i1().
- *
- */
- const double A[] = {2.77791411276104639959E-18, -2.11142121435816608115E-17,
- 1.55363195773620046921E-16, -1.10559694773538630805E-15,
- 7.60068429473540693410E-15, -5.04218550472791168711E-14,
- 3.22379336594557470981E-13, -1.98397439776494371520E-12,
- 1.17361862988909016308E-11, -6.66348972350202774223E-11,
- 3.62559028155211703701E-10, -1.88724975172282928790E-9,
- 9.38153738649577178388E-9, -4.44505912879632808065E-8,
- 2.00329475355213526229E-7, -8.56872026469545474066E-7,
- 3.47025130813767847674E-6, -1.32731636560394358279E-5,
- 4.78156510755005422638E-5, -1.61760815825896745588E-4,
- 5.12285956168575772895E-4, -1.51357245063125314899E-3,
- 4.15642294431288815669E-3, -1.05640848946261981558E-2,
- 2.47264490306265168283E-2, -5.29459812080949914269E-2,
- 1.02643658689847095384E-1, -1.76416518357834055153E-1,
- 2.52587186443633654823E-1};
- const double B[] = {
- 7.51729631084210481353E-18, 4.41434832307170791151E-18,
- -4.65030536848935832153E-17, -3.20952592199342395980E-17,
- 2.96262899764595013876E-16, 3.30820231092092828324E-16,
- -1.88035477551078244854E-15, -3.81440307243700780478E-15,
- 1.04202769841288027642E-14, 4.27244001671195135429E-14,
- -2.10154184277266431302E-14, -4.08355111109219731823E-13,
- -7.19855177624590851209E-13, 2.03562854414708950722E-12,
- 1.41258074366137813316E-11, 3.25260358301548823856E-11,
- -1.89749581235054123450E-11, -5.58974346219658380687E-10,
- -3.83538038596423702205E-9, -2.63146884688951950684E-8,
- -2.51223623787020892529E-7, -3.88256480887769039346E-6,
- -1.10588938762623716291E-4, -9.76109749136146840777E-3,
- 7.78576235018280120474E-1};
- T y = pabs(x);
- T y_le_eight = pmul(y, internal::pchebevl<T, 29>::run(
- pmadd(pset1<T>(0.5), y, pset1<T>(-2.0)), A));
- T y_gt_eight = pdiv(
- internal::pchebevl<T, 25>::run(
- psub(pdiv(pset1<T>(32.0), y),
- pset1<T>(2.0)), B),
- psqrt(y));
- // TODO: Perhaps instead check whether all packet elements are in
- // [-8, 8] and evaluate a branch based off of that. It's possible
- // in practice most elements are in this region.
- y = pselect(pcmp_le(y, pset1<T>(8.0)), y_le_eight, y_gt_eight);
- return pselect(pcmp_lt(x, pset1<T>(0.0f)), -y, y);
- }
-};
-
-template <typename Scalar>
-struct i1e_impl {
- EIGEN_DEVICE_FUNC
- static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
- return generic_i1e<Scalar, Scalar>::run(x);
- }
-};
-
} // end namespace internal
namespace numext {
@@ -2285,21 +1940,7 @@ EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(betainc, Scalar)
return EIGEN_MATHFUNC_IMPL(betainc, Scalar)::run(a, b, x);
}
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(i0e, Scalar)
- i0e(const Scalar& x) {
- return EIGEN_MATHFUNC_IMPL(i0e, Scalar)::run(x);
-}
-
-template <typename Scalar>
-EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(i1e, Scalar)
- i1e(const Scalar& x) {
- return EIGEN_MATHFUNC_IMPL(i1e, Scalar)::run(x);
-}
-
} // end namespace numext
-
-
} // end namespace Eigen
#endif // EIGEN_SPECIAL_FUNCTIONS_H
diff --git a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h
index 21908e512..577015690 100644
--- a/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h
+++ b/unsupported/Eigen/src/SpecialFunctions/SpecialFunctionsPacketMath.h
@@ -72,24 +72,6 @@ Packet pigammac(const Packet& a, const Packet& x) { using numext::igammac; retur
template<typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
Packet pbetainc(const Packet& a, const Packet& b,const Packet& x) { using numext::betainc; return betainc(a, b, x); }
-/** \internal \returns the exponentially scaled modified Bessel function of
- * order zero i0e(\a a) (coeff-wise) */
-template <typename Packet>
-EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
-Packet pi0e(const Packet& x) {
- typedef typename unpacket_traits<Packet>::type ScalarType;
- using internal::generic_i0e; return generic_i0e<Packet, ScalarType>::run(x);
-}
-
-/** \internal \returns the exponentially scaled modified Bessel function of
- * order one i1e(\a a) (coeff-wise) */
-template <typename Packet>
-EIGEN_DEVICE_FUNC EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
-Packet pi1e(const Packet& x) {
- typedef typename unpacket_traits<Packet>::type ScalarType;
- using internal::generic_i1e; return generic_i1e<Packet, ScalarType>::run(x);
-}
-
} // end namespace internal
} // end namespace Eigen
diff --git a/unsupported/Eigen/src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h b/unsupported/Eigen/src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h
index c831edc17..b886e278c 100644
--- a/unsupported/Eigen/src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h
+++ b/unsupported/Eigen/src/SpecialFunctions/arch/GPU/GpuSpecialFunctions.h
@@ -218,6 +218,19 @@ pi0e<double2>(const double2& x) {
}
template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pi0<float4>(const float4& x) {
+ using numext::i0;
+ return make_float4(i0(x.x), i0(x.y), i0(x.z), i0(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+pi0<double2>(const double2& x) {
+ using numext::i0;
+ return make_double2(i0(x.x), i0(x.y));
+}
+
+template <>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pi1e<float4>(const float4& x) {
using numext::i1e;
return make_float4(i1e(x.x), i1e(x.y), i1e(x.z), i1e(x.w));
@@ -230,6 +243,123 @@ pi1e<double2>(const double2& x) {
return make_double2(i1e(x.x), i1e(x.y));
}
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pi1<float4>(const float4& x) {
+ using numext::i1;
+ return make_float4(i1(x.x), i1(x.y), i1(x.z), i1(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+pi1<double2>(const double2& x) {
+ using numext::i1;
+ return make_double2(i1(x.x), i1(x.y));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pk0e<float4>(const float4& x) {
+ using numext::k0e;
+ return make_float4(k0e(x.x), k0e(x.y), k0e(x.z), k0e(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+pk0e<double2>(const double2& x) {
+ using numext::k0e;
+ return make_double2(k0e(x.x), k0e(x.y));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pk0<float4>(const float4& x) {
+ using numext::k0;
+ return make_float4(k0(x.x), k0(x.y), k0(x.z), k0(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+pk0<double2>(const double2& x) {
+ using numext::k0;
+ return make_double2(k0(x.x), k0(x.y));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pk1e<float4>(const float4& x) {
+ using numext::k1e;
+ return make_float4(k1e(x.x), k1e(x.y), k1e(x.z), k1e(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+pk1e<double2>(const double2& x) {
+ using numext::k1e;
+ return make_double2(k1e(x.x), k1e(x.y));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pk1<float4>(const float4& x) {
+ using numext::k1;
+ return make_float4(k1(x.x), k1(x.y), k1(x.z), k1(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+pk1<double2>(const double2& x) {
+ using numext::k1;
+ return make_double2(k1(x.x), k1(x.y));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pj0<float4>(const float4& x) {
+ using numext::j0;
+ return make_float4(j0(x.x), j0(x.y), j0(x.z), j0(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+pj0<double2>(const double2& x) {
+ using numext::j0;
+ return make_double2(j0(x.x), j0(x.y));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 pj1<float4>(const float4& x) {
+ using numext::j1;
+ return make_float4(j1(x.x), j1(x.y), j1(x.z), j1(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+pj1<double2>(const double2& x) {
+ using numext::j1;
+ return make_double2(j1(x.x), j1(x.y));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 py0<float4>(const float4& x) {
+ using numext::y0;
+ return make_float4(y0(x.x), y0(x.y), y0(x.z), y0(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+py0<double2>(const double2& x) {
+ using numext::y0;
+ return make_double2(y0(x.x), y0(x.y));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float4 py1<float4>(const float4& x) {
+ using numext::y1;
+ return make_float4(y1(x.x), y1(x.y), y1(x.z), y1(x.w));
+}
+
+template <>
+EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double2
+py1<double2>(const double2& x) {
+ using numext::y1;
+ return make_double2(y1(x.x), y1(x.y));
+}
+
#endif
} // end namespace internal
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-15 16:07:32 +0000