src/SpecialFunctions/SpecialFunctionsFunctors.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2016 Eugene Brevdo <ebrevdo@gmail.com>
0005 // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_SPECIALFUNCTIONS_FUNCTORS_H
0012 #define EIGEN_SPECIALFUNCTIONS_FUNCTORS_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017
0018
0019 /** \internal
0020   * \brief Template functor to compute the incomplete gamma function igamma(a, x)
0021   *
0022   * \sa class CwiseBinaryOp, Cwise::igamma
0023   */
0024 template<typename Scalar> struct scalar_igamma_op : binary_op_base<Scalar,Scalar>
0025 {
0026   EIGEN_EMPTY_STRUCT_CTOR(scalar_igamma_op)
0027   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& x) const {
0028     using numext::igamma; return igamma(a, x);
0029   }
0030   template<typename Packet>
0031   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const {
0032     return internal::pigamma(a, x);
0033   }
0034 };
0035 template<typename Scalar>
0036 struct functor_traits<scalar_igamma_op<Scalar> > {
0037   enum {
0038     // Guesstimate
0039     Cost = 20 * NumTraits<Scalar>::MulCost + 10 * NumTraits<Scalar>::AddCost,
0040     PacketAccess = packet_traits<Scalar>::HasIGamma
0041   };
0042 };
0043
0044 /** \internal
0045   * \brief Template functor to compute the derivative of the incomplete gamma
0046   * function igamma_der_a(a, x)
0047   *
0048   * \sa class CwiseBinaryOp, Cwise::igamma_der_a
0049   */
0050 template <typename Scalar>
0051 struct scalar_igamma_der_a_op {
0052   EIGEN_EMPTY_STRUCT_CTOR(scalar_igamma_der_a_op)
0053   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& a, const Scalar& x) const {
0054     using numext::igamma_der_a;
0055     return igamma_der_a(a, x);
0056   }
0057   template <typename Packet>
0058   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const {
0059     return internal::pigamma_der_a(a, x);
0060   }
0061 };
0062 template <typename Scalar>
0063 struct functor_traits<scalar_igamma_der_a_op<Scalar> > {
0064   enum {
0065     // 2x the cost of igamma
0066     Cost = 40 * NumTraits<Scalar>::MulCost + 20 * NumTraits<Scalar>::AddCost,
0067     PacketAccess = packet_traits<Scalar>::HasIGammaDerA
0068   };
0069 };
0070
0071 /** \internal
0072   * \brief Template functor to compute the derivative of the sample
0073   * of a Gamma(alpha, 1) random variable with respect to the parameter alpha
0074   * gamma_sample_der_alpha(alpha, sample)
0075   *
0076   * \sa class CwiseBinaryOp, Cwise::gamma_sample_der_alpha
0077   */
0078 template <typename Scalar>
0079 struct scalar_gamma_sample_der_alpha_op {
0080   EIGEN_EMPTY_STRUCT_CTOR(scalar_gamma_sample_der_alpha_op)
0081   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator()(const Scalar& alpha, const Scalar& sample) const {
0082     using numext::gamma_sample_der_alpha;
0083     return gamma_sample_der_alpha(alpha, sample);
0084   }
0085   template <typename Packet>
0086   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& alpha, const Packet& sample) const {
0087     return internal::pgamma_sample_der_alpha(alpha, sample);
0088   }
0089 };
0090 template <typename Scalar>
0091 struct functor_traits<scalar_gamma_sample_der_alpha_op<Scalar> > {
0092   enum {
0093     // 2x the cost of igamma, minus the lgamma cost (the lgamma cancels out)
0094     Cost = 30 * NumTraits<Scalar>::MulCost + 15 * NumTraits<Scalar>::AddCost,
0095     PacketAccess = packet_traits<Scalar>::HasGammaSampleDerAlpha
0096   };
0097 };
0098
0099 /** \internal
0100   * \brief Template functor to compute the complementary incomplete gamma function igammac(a, x)
0101   *
0102   * \sa class CwiseBinaryOp, Cwise::igammac
0103   */
0104 template<typename Scalar> struct scalar_igammac_op : binary_op_base<Scalar,Scalar>
0105 {
0106   EIGEN_EMPTY_STRUCT_CTOR(scalar_igammac_op)
0107   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& x) const {
0108     using numext::igammac; return igammac(a, x);
0109   }
0110   template<typename Packet>
0111   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& x) const
0112   {
0113     return internal::pigammac(a, x);
0114   }
0115 };
0116 template<typename Scalar>
0117 struct functor_traits<scalar_igammac_op<Scalar> > {
0118   enum {
0119     // Guesstimate
0120     Cost = 20 * NumTraits<Scalar>::MulCost + 10 * NumTraits<Scalar>::AddCost,
0121     PacketAccess = packet_traits<Scalar>::HasIGammac
0122   };
0123 };
0124
0125
0126 /** \internal
0127   * \brief Template functor to compute the incomplete beta integral betainc(a, b, x)
0128   *
0129   */
0130 template<typename Scalar> struct scalar_betainc_op {
0131   EIGEN_EMPTY_STRUCT_CTOR(scalar_betainc_op)
0132   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& x, const Scalar& a, const Scalar& b) const {
0133     using numext::betainc; return betainc(x, a, b);
0134   }
0135   template<typename Packet>
0136   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& x, const Packet& a, const Packet& b) const
0137   {
0138     return internal::pbetainc(x, a, b);
0139   }
0140 };
0141 template<typename Scalar>
0142 struct functor_traits<scalar_betainc_op<Scalar> > {
0143   enum {
0144     // Guesstimate
0145     Cost = 400 * NumTraits<Scalar>::MulCost + 400 * NumTraits<Scalar>::AddCost,
0146     PacketAccess = packet_traits<Scalar>::HasBetaInc
0147   };
0148 };
0149
0150
0151 /** \internal
0152  * \brief Template functor to compute the natural log of the absolute
0153  * value of Gamma of a scalar
0154  * \sa class CwiseUnaryOp, Cwise::lgamma()
0155  */
0156 template<typename Scalar> struct scalar_lgamma_op {
0157   EIGEN_EMPTY_STRUCT_CTOR(scalar_lgamma_op)
0158   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const {
0159     using numext::lgamma; return lgamma(a);
0160   }
0161   typedef typename packet_traits<Scalar>::type Packet;
0162   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& a) const { return internal::plgamma(a); }
0163 };
0164 template<typename Scalar>
0165 struct functor_traits<scalar_lgamma_op<Scalar> >
0166 {
0167   enum {
0168     // Guesstimate
0169     Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
0170     PacketAccess = packet_traits<Scalar>::HasLGamma
0171   };
0172 };
0173
0174 /** \internal
0175  * \brief Template functor to compute psi, the derivative of lgamma of a scalar.
0176  * \sa class CwiseUnaryOp, Cwise::digamma()
0177  */
0178 template<typename Scalar> struct scalar_digamma_op {
0179   EIGEN_EMPTY_STRUCT_CTOR(scalar_digamma_op)
0180   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const {
0181     using numext::digamma; return digamma(a);
0182   }
0183   typedef typename packet_traits<Scalar>::type Packet;
0184   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& a) const { return internal::pdigamma(a); }
0185 };
0186 template<typename Scalar>
0187 struct functor_traits<scalar_digamma_op<Scalar> >
0188 {
0189   enum {
0190     // Guesstimate
0191     Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
0192     PacketAccess = packet_traits<Scalar>::HasDiGamma
0193   };
0194 };
0195
0196 /** \internal
0197  * \brief Template functor to compute the Riemann Zeta function of two arguments.
0198  * \sa class CwiseUnaryOp, Cwise::zeta()
0199  */
0200 template<typename Scalar> struct scalar_zeta_op {
0201     EIGEN_EMPTY_STRUCT_CTOR(scalar_zeta_op)
0202     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& x, const Scalar& q) const {
0203         using numext::zeta; return zeta(x, q);
0204     }
0205     typedef typename packet_traits<Scalar>::type Packet;
0206     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x, const Packet& q) const { return internal::pzeta(x, q); }
0207 };
0208 template<typename Scalar>
0209 struct functor_traits<scalar_zeta_op<Scalar> >
0210 {
0211     enum {
0212         // Guesstimate
0213         Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
0214         PacketAccess = packet_traits<Scalar>::HasZeta
0215     };
0216 };
0217
0218 /** \internal
0219  * \brief Template functor to compute the polygamma function.
0220  * \sa class CwiseUnaryOp, Cwise::polygamma()
0221  */
0222 template<typename Scalar> struct scalar_polygamma_op {
0223     EIGEN_EMPTY_STRUCT_CTOR(scalar_polygamma_op)
0224     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& n, const Scalar& x) const {
0225         using numext::polygamma; return polygamma(n, x);
0226     }
0227     typedef typename packet_traits<Scalar>::type Packet;
0228     EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& n, const Packet& x) const { return internal::ppolygamma(n, x); }
0229 };
0230 template<typename Scalar>
0231 struct functor_traits<scalar_polygamma_op<Scalar> >
0232 {
0233     enum {
0234         // Guesstimate
0235         Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
0236         PacketAccess = packet_traits<Scalar>::HasPolygamma
0237     };
0238 };
0239
0240 /** \internal
0241  * \brief Template functor to compute the error function of a scalar
0242  * \sa class CwiseUnaryOp, ArrayBase::erf()
0243  */
0244 template<typename Scalar> struct scalar_erf_op {
0245   EIGEN_EMPTY_STRUCT_CTOR(scalar_erf_op)
0246   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar
0247   operator()(const Scalar& a) const {
0248     return numext::erf(a);
0249   }
0250   template <typename Packet>
0251   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
0252     return perf(x);
0253   }
0254 };
0255 template <typename Scalar>
0256 struct functor_traits<scalar_erf_op<Scalar> > {
0257   enum {
0258     PacketAccess = packet_traits<Scalar>::HasErf,
0259     Cost =
0260         (PacketAccess
0261 #ifdef EIGEN_VECTORIZE_FMA
0262              // TODO(rmlarsen): Move the FMA cost model to a central location.
0263              // Haswell can issue 2 add/mul/madd per cycle.
0264              // 10 pmadd, 2 pmul, 1 div, 2 other
0265              ? (2 * NumTraits<Scalar>::AddCost +
0266                 7 * NumTraits<Scalar>::MulCost +
0267                 scalar_div_cost<Scalar, packet_traits<Scalar>::HasDiv>::value)
0268 #else
0269              ? (12 * NumTraits<Scalar>::AddCost +
0270                 12 * NumTraits<Scalar>::MulCost +
0271                 scalar_div_cost<Scalar, packet_traits<Scalar>::HasDiv>::value)
0272 #endif
0273              // Assume for simplicity that this is as expensive as an exp().
0274              : (functor_traits<scalar_exp_op<Scalar> >::Cost))
0275   };
0276 };
0277
0278 /** \internal
0279  * \brief Template functor to compute the Complementary Error Function
0280  * of a scalar
0281  * \sa class CwiseUnaryOp, Cwise::erfc()
0282  */
0283 template<typename Scalar> struct scalar_erfc_op {
0284   EIGEN_EMPTY_STRUCT_CTOR(scalar_erfc_op)
0285   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const {
0286     using numext::erfc; return erfc(a);
0287   }
0288   typedef typename packet_traits<Scalar>::type Packet;
0289   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& a) const { return internal::perfc(a); }
0290 };
0291 template<typename Scalar>
0292 struct functor_traits<scalar_erfc_op<Scalar> >
0293 {
0294   enum {
0295     // Guesstimate
0296     Cost = 10 * NumTraits<Scalar>::MulCost + 5 * NumTraits<Scalar>::AddCost,
0297     PacketAccess = packet_traits<Scalar>::HasErfc
0298   };
0299 };
0300
0301 /** \internal
0302  * \brief Template functor to compute the Inverse of the normal distribution
0303  * function of a scalar
0304  * \sa class CwiseUnaryOp, Cwise::ndtri()
0305  */
0306 template<typename Scalar> struct scalar_ndtri_op {
0307   EIGEN_EMPTY_STRUCT_CTOR(scalar_ndtri_op)
0308   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const {
0309     using numext::ndtri; return ndtri(a);
0310   }
0311   typedef typename packet_traits<Scalar>::type Packet;
0312   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& a) const { return internal::pndtri(a); }
0313 };
0314 template<typename Scalar>
0315 struct functor_traits<scalar_ndtri_op<Scalar> >
0316 {
0317   enum {
0318     // On average, We are evaluating rational functions with degree N=9 in the
0319     // numerator and denominator. This results in 2*N additions and 2*N
0320     // multiplications.
0321     Cost = 18 * NumTraits<Scalar>::MulCost + 18 * NumTraits<Scalar>::AddCost,
0322     PacketAccess = packet_traits<Scalar>::HasNdtri
0323   };
0324 };
0325
0326 } // end namespace internal
0327
0328 } // end namespace Eigen
0329
0330 #endif // EIGEN_SPECIALFUNCTIONS_FUNCTORS_H