10 #ifndef EIGEN_SPECIAL_FUNCTIONS_H 11 #define EIGEN_SPECIAL_FUNCTIONS_H 80 template <
typename Scalar,
int N>
83 static EIGEN_STRONG_INLINE Scalar run(
const Scalar x,
const Scalar coef[]) {
84 EIGEN_STATIC_ASSERT((N > 0), YOU_MADE_A_PROGRAMMING_MISTAKE);
86 return polevl<Scalar, N - 1>::run(x, coef) * x + coef[N];
90 template <
typename Scalar>
91 struct polevl<Scalar, 0> {
93 static EIGEN_STRONG_INLINE Scalar run(
const Scalar,
const Scalar coef[]) {
104 template <
typename Scalar>
107 static EIGEN_STRONG_INLINE Scalar run(
const Scalar) {
108 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
109 THIS_TYPE_IS_NOT_SUPPORTED);
114 template <
typename Scalar>
115 struct lgamma_retval {
119 #if EIGEN_HAS_C99_MATH 121 struct lgamma_impl<float> {
123 static EIGEN_STRONG_INLINE
float run(
float x) { return ::lgammaf(x); }
127 struct lgamma_impl<double> {
129 static EIGEN_STRONG_INLINE
double run(
double x) { return ::lgamma(x); }
137 template <
typename Scalar>
138 struct digamma_retval {
155 template <
typename Scalar>
156 struct digamma_impl_maybe_poly {
158 static EIGEN_STRONG_INLINE Scalar run(
const Scalar) {
159 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
160 THIS_TYPE_IS_NOT_SUPPORTED);
167 struct digamma_impl_maybe_poly<float> {
169 static EIGEN_STRONG_INLINE
float run(
const float s) {
171 -4.16666666666666666667E-3f,
172 3.96825396825396825397E-3f,
173 -8.33333333333333333333E-3f,
174 8.33333333333333333333E-2f
180 return z * cephes::polevl<float, 3>::run(z, A);
186 struct digamma_impl_maybe_poly<double> {
188 static EIGEN_STRONG_INLINE
double run(
const double s) {
190 8.33333333333333333333E-2,
191 -2.10927960927960927961E-2,
192 7.57575757575757575758E-3,
193 -4.16666666666666666667E-3,
194 3.96825396825396825397E-3,
195 -8.33333333333333333333E-3,
196 8.33333333333333333333E-2
202 return z * cephes::polevl<double, 6>::run(z, A);
208 template <
typename Scalar>
209 struct digamma_impl {
211 static Scalar run(Scalar x) {
269 Scalar p, q, nz, s, w, y;
270 bool negative =
false;
272 const Scalar maxnum = NumTraits<Scalar>::infinity();
273 const Scalar m_pi = Scalar(EIGEN_PI);
275 const Scalar zero = Scalar(0);
276 const Scalar one = Scalar(1);
277 const Scalar half = Scalar(0.5);
283 p = numext::floor(q);
296 nz = m_pi / numext::tan(m_pi * nz);
307 while (s < Scalar(10)) {
312 y = digamma_impl_maybe_poly<Scalar>::run(s);
314 y = numext::log(s) - (half / s) - y - w;
316 return (negative) ? y - nz : y;
324 template <
typename Scalar>
327 static EIGEN_STRONG_INLINE Scalar run(
const Scalar) {
328 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
329 THIS_TYPE_IS_NOT_SUPPORTED);
334 template <
typename Scalar>
339 #if EIGEN_HAS_C99_MATH 341 struct erf_impl<float> {
343 static EIGEN_STRONG_INLINE
float run(
float x) { return ::erff(x); }
347 struct erf_impl<double> {
349 static EIGEN_STRONG_INLINE
double run(
double x) { return ::erf(x); }
351 #endif // EIGEN_HAS_C99_MATH 357 template <
typename Scalar>
360 static EIGEN_STRONG_INLINE Scalar run(
const Scalar) {
361 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
362 THIS_TYPE_IS_NOT_SUPPORTED);
367 template <
typename Scalar>
372 #if EIGEN_HAS_C99_MATH 374 struct erfc_impl<float> {
376 static EIGEN_STRONG_INLINE
float run(
const float x) { return ::erfcf(x); }
380 struct erfc_impl<double> {
382 static EIGEN_STRONG_INLINE
double run(
const double x) { return ::erfc(x); }
384 #endif // EIGEN_HAS_C99_MATH 390 template <
typename Scalar>
391 struct igammac_retval {
396 template <
typename Scalar>
397 struct cephes_helper {
399 static EIGEN_STRONG_INLINE Scalar machep() { assert(
false &&
"machep not supported for this type");
return 0.0; }
401 static EIGEN_STRONG_INLINE Scalar big() { assert(
false &&
"big not supported for this type");
return 0.0; }
403 static EIGEN_STRONG_INLINE Scalar biginv() { assert(
false &&
"biginv not supported for this type");
return 0.0; }
407 struct cephes_helper<float> {
409 static EIGEN_STRONG_INLINE
float machep() {
410 return NumTraits<float>::epsilon() / 2;
413 static EIGEN_STRONG_INLINE
float big() {
415 return 1.0f / (NumTraits<float>::epsilon() / 2);
418 static EIGEN_STRONG_INLINE
float biginv() {
425 struct cephes_helper<double> {
427 static EIGEN_STRONG_INLINE
double machep() {
428 return NumTraits<double>::epsilon() / 2;
431 static EIGEN_STRONG_INLINE
double big() {
432 return 1.0 / NumTraits<double>::epsilon();
435 static EIGEN_STRONG_INLINE
double biginv() {
437 return NumTraits<double>::epsilon();
441 #if !EIGEN_HAS_C99_MATH 443 template <
typename Scalar>
444 struct igammac_impl {
446 static Scalar run(Scalar a, Scalar x) {
447 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
448 THIS_TYPE_IS_NOT_SUPPORTED);
455 template <
typename Scalar>
struct igamma_impl;
457 template <
typename Scalar>
458 struct igammac_impl {
460 static Scalar run(Scalar a, Scalar x) {
515 const Scalar zero = 0;
516 const Scalar one = 1;
517 const Scalar nan = NumTraits<Scalar>::quiet_NaN();
519 if ((x < zero) || (a <= zero)) {
524 if ((x < one) || (x < a)) {
532 return (one - igamma_impl<Scalar>::Impl(a, x));
540 friend struct igamma_impl<Scalar>;
549 EIGEN_DEVICE_FUNC
static Scalar Impl(Scalar a, Scalar x) {
550 const Scalar zero = 0;
551 const Scalar one = 1;
552 const Scalar two = 2;
553 const Scalar machep = cephes_helper<Scalar>::machep();
554 const Scalar maxlog = numext::log(NumTraits<Scalar>::highest());
555 const Scalar big = cephes_helper<Scalar>::big();
556 const Scalar biginv = cephes_helper<Scalar>::biginv();
557 const Scalar inf = NumTraits<Scalar>::infinity();
559 Scalar ans, ax, c, yc, r, t, y, z;
560 Scalar pk, pkm1, pkm2, qk, qkm1, qkm2;
562 if (x == inf)
return zero;
565 ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
569 ax = numext::exp(ax);
586 pk = pkm1 * z - pkm2 * yc;
587 qk = qkm1 * z - qkm2 * yc;
590 t = numext::abs((ans - r) / r);
599 if (numext::abs(pk) > big) {
614 #endif // EIGEN_HAS_C99_MATH 620 template <
typename Scalar>
621 struct igamma_retval {
625 #if !EIGEN_HAS_C99_MATH 627 template <
typename Scalar>
630 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) {
631 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
632 THIS_TYPE_IS_NOT_SUPPORTED);
639 template <
typename Scalar>
642 static Scalar run(Scalar a, Scalar x) {
703 const Scalar zero = 0;
704 const Scalar one = 1;
705 const Scalar nan = NumTraits<Scalar>::quiet_NaN();
707 if (x == zero)
return zero;
709 if ((x < zero) || (a <= zero)) {
713 if ((x > one) && (x > a)) {
721 return (one - igammac_impl<Scalar>::Impl(a, x));
729 friend struct igammac_impl<Scalar>;
738 EIGEN_DEVICE_FUNC
static Scalar Impl(Scalar a, Scalar x) {
739 const Scalar zero = 0;
740 const Scalar one = 1;
741 const Scalar machep = cephes_helper<Scalar>::machep();
742 const Scalar maxlog = numext::log(NumTraits<Scalar>::highest());
744 Scalar ans, ax, c, r;
747 ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
752 ax = numext::exp(ax);
763 if (c/ans <= machep) {
768 return (ans * ax / a);
772 #endif // EIGEN_HAS_C99_MATH 778 template <
typename Scalar>
783 template <
typename Scalar>
784 struct zeta_impl_series {
786 static EIGEN_STRONG_INLINE Scalar run(
const Scalar) {
787 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
788 THIS_TYPE_IS_NOT_SUPPORTED);
794 struct zeta_impl_series<float> {
796 static EIGEN_STRONG_INLINE
bool run(
float& a,
float& b,
float& s,
const float x,
const float machep) {
802 b = numext::pow( a, -x );
804 if( numext::abs(b/s) < machep )
814 struct zeta_impl_series<double> {
816 static EIGEN_STRONG_INLINE
bool run(
double& a,
double& b,
double& s,
const double x,
const double machep) {
818 while( (i < 9) || (a <= 9.0) )
822 b = numext::pow( a, -x );
824 if( numext::abs(b/s) < machep )
833 template <
typename Scalar>
836 static Scalar run(Scalar x, Scalar q) {
899 Scalar p, r, a, b, k, s, t, w;
907 Scalar(-1.8924375803183791606e9),
908 Scalar(7.47242496e10),
909 Scalar(-2.950130727918164224e12),
910 Scalar(1.1646782814350067249e14),
911 Scalar(-4.5979787224074726105e15),
912 Scalar(1.8152105401943546773e17),
913 Scalar(-7.1661652561756670113e18)
916 const Scalar maxnum = NumTraits<Scalar>::infinity();
917 const Scalar zero = 0.0, half = 0.5, one = 1.0;
918 const Scalar machep = cephes_helper<Scalar>::machep();
919 const Scalar nan = NumTraits<Scalar>::quiet_NaN();
931 if(q == numext::floor(q))
936 r = numext::floor(p);
946 s = numext::pow( q, -x );
950 if (zeta_impl_series<Scalar>::run(a, b, s, x, machep)) {
959 for( i=0; i<12; i++ )
965 t = numext::abs(t/s);
982 template <
typename Scalar>
983 struct polygamma_retval {
987 #if !EIGEN_HAS_C99_MATH 989 template <
typename Scalar>
990 struct polygamma_impl {
992 static EIGEN_STRONG_INLINE Scalar run(Scalar n, Scalar x) {
993 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
994 THIS_TYPE_IS_NOT_SUPPORTED);
1001 template <
typename Scalar>
1002 struct polygamma_impl {
1004 static Scalar run(Scalar n, Scalar x) {
1005 Scalar zero = 0.0, one = 1.0;
1006 Scalar nplus = n + one;
1007 const Scalar nan = NumTraits<Scalar>::quiet_NaN();
1010 if (numext::floor(n) != n) {
1014 else if (n == zero) {
1015 return digamma_impl<Scalar>::run(x);
1019 Scalar factorial = numext::exp(lgamma_impl<Scalar>::run(nplus));
1020 return numext::pow(-one, nplus) * factorial * zeta_impl<Scalar>::run(nplus, x);
1025 #endif // EIGEN_HAS_C99_MATH 1031 template <
typename Scalar>
1032 struct betainc_retval {
1033 typedef Scalar type;
1036 #if !EIGEN_HAS_C99_MATH 1038 template <
typename Scalar>
1039 struct betainc_impl {
1041 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x) {
1042 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
1043 THIS_TYPE_IS_NOT_SUPPORTED);
1050 template <
typename Scalar>
1051 struct betainc_impl {
1053 static EIGEN_STRONG_INLINE Scalar run(Scalar, Scalar, Scalar) {
1123 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
1124 THIS_TYPE_IS_NOT_SUPPORTED);
1132 template <
typename Scalar>
1133 struct incbeta_cfe {
1135 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x,
bool small_branch) {
1136 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, float>::value ||
1137 internal::is_same<Scalar, double>::value),
1138 THIS_TYPE_IS_NOT_SUPPORTED);
1139 const Scalar big = cephes_helper<Scalar>::big();
1140 const Scalar machep = cephes_helper<Scalar>::machep();
1141 const Scalar biginv = cephes_helper<Scalar>::biginv();
1143 const Scalar zero = 0;
1144 const Scalar one = 1;
1145 const Scalar two = 2;
1147 Scalar xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
1148 Scalar k1, k2, k3, k4, k5, k6, k7, k8, k26update;
1152 const int num_iters = (internal::is_same<Scalar, float>::value) ? 100 : 300;
1153 const Scalar thresh =
1154 (internal::is_same<Scalar, float>::value) ? machep : Scalar(3) * machep;
1155 Scalar r = (internal::is_same<Scalar, float>::value) ? zero : one;
1188 xk = -(x * k1 * k2) / (k3 * k4);
1189 pk = pkm1 + pkm2 * xk;
1190 qk = qkm1 + qkm2 * xk;
1196 xk = (x * k5 * k6) / (k7 * k8);
1197 pk = pkm1 + pkm2 * xk;
1198 qk = qkm1 + qkm2 * xk;
1206 if (numext::abs(ans - r) < numext::abs(r) * thresh) {
1221 if ((numext::abs(qk) + numext::abs(pk)) > big) {
1227 if ((numext::abs(qk) < biginv) || (numext::abs(pk) < biginv)) {
1233 }
while (++n < num_iters);
1240 template <
typename Scalar>
1241 struct betainc_helper {};
1244 struct betainc_helper<float> {
1246 EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE
float incbsa(
float aa,
float bb,
1248 float ans, a, b, t, x, onemx;
1249 bool reversed_a_b =
false;
1254 if (xx > (aa / (aa + bb))) {
1255 reversed_a_b =
true;
1269 if (numext::abs(b * x / a) < 0.3f) {
1270 t = betainc_helper<float>::incbps(a, b, x);
1271 if (reversed_a_b) t = 1.0f - t;
1276 ans = x * (a + b - 2.0f) / (a - 1.0f);
1278 ans = incbeta_cfe<float>::run(a, b, x,
true );
1279 t = b * numext::log(t);
1281 ans = incbeta_cfe<float>::run(a, b, x,
false );
1282 t = (b - 1.0f) * numext::log(t);
1285 t += a * numext::log(x) + lgamma_impl<float>::run(a + b) -
1286 lgamma_impl<float>::run(a) - lgamma_impl<float>::run(b);
1287 t += numext::log(ans / a);
1290 if (reversed_a_b) t = 1.0f - t;
1295 static EIGEN_STRONG_INLINE
float incbps(
float a,
float b,
float x) {
1297 const float machep = cephes_helper<float>::machep();
1299 y = a * numext::log(x) + (b - 1.0f) * numext::log1p(-x) - numext::log(a);
1300 y -= lgamma_impl<float>::run(a) + lgamma_impl<float>::run(b);
1301 y += lgamma_impl<float>::run(a + b);
1314 }
while (numext::abs(u) > machep);
1316 return numext::exp(y) * (1.0f + s);
1321 struct betainc_impl<float> {
1323 static float run(
float a,
float b,
float x) {
1324 const float nan = NumTraits<float>::quiet_NaN();
1327 if (a <= 0.0f)
return nan;
1328 if (b <= 0.0f)
return nan;
1329 if ((x <= 0.0f) || (x >= 1.0f)) {
1330 if (x == 0.0f)
return 0.0f;
1331 if (x == 1.0f)
return 1.0f;
1338 ans = betainc_helper<float>::incbsa(a + 1.0f, b, x);
1339 t = a * numext::log(x) + b * numext::log1p(-x) +
1340 lgamma_impl<float>::run(a + b) - lgamma_impl<float>::run(a + 1.0f) -
1341 lgamma_impl<float>::run(b);
1342 return (ans + numext::exp(t));
1344 return betainc_helper<float>::incbsa(a, b, x);
1350 struct betainc_helper<double> {
1352 static EIGEN_STRONG_INLINE
double incbps(
double a,
double b,
double x) {
1353 const double machep = cephes_helper<double>::machep();
1355 double s, t, u, v, n, t1, z, ai;
1365 while (numext::abs(v) > z) {
1366 u = (n - b) * x / n;
1375 u = a * numext::log(x);
1383 t = lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) -
1384 lgamma_impl<double>::run(b) + u + numext::log(s);
1385 return s = numext::exp(t);
1390 struct betainc_impl<double> {
1392 static double run(
double aa,
double bb,
double xx) {
1393 const double nan = NumTraits<double>::quiet_NaN();
1394 const double machep = cephes_helper<double>::machep();
1397 double a, b, t, x, xc, w, y;
1398 bool reversed_a_b =
false;
1400 if (aa <= 0.0 || bb <= 0.0) {
1404 if ((xx <= 0.0) || (xx >= 1.0)) {
1405 if (xx == 0.0)
return (0.0);
1406 if (xx == 1.0)
return (1.0);
1411 if ((bb * xx) <= 1.0 && xx <= 0.95) {
1412 return betainc_helper<double>::incbps(aa, bb, xx);
1418 if (xx > (aa / (aa + bb))) {
1419 reversed_a_b =
true;
1431 if (reversed_a_b && (b * x) <= 1.0 && x <= 0.95) {
1432 t = betainc_helper<double>::incbps(a, b, x);
1442 y = x * (a + b - 2.0) - (a - 1.0);
1444 w = incbeta_cfe<double>::run(a, b, x,
true );
1446 w = incbeta_cfe<double>::run(a, b, x,
false ) / xc;
1453 y = a * numext::log(x);
1454 t = b * numext::log(xc);
1467 y += t + lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) -
1468 lgamma_impl<double>::run(b);
1469 y += numext::log(w / a);
1486 #endif // EIGEN_HAS_C99_MATH 1492 template <
typename Scalar>
1493 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(lgamma, Scalar)
1494 lgamma(const Scalar& x) {
1495 return EIGEN_MATHFUNC_IMPL(lgamma, Scalar)::run(x);
1498 template <
typename Scalar>
1499 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(digamma, Scalar)
1500 digamma(const Scalar& x) {
1501 return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x);
1504 template <
typename Scalar>
1505 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
zeta, Scalar)
1506 zeta(const Scalar& x, const Scalar& q) {
1507 return EIGEN_MATHFUNC_IMPL(
zeta, Scalar)::run(x, q);
1510 template <
typename Scalar>
1511 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
polygamma, Scalar)
1512 polygamma(const Scalar& n, const Scalar& x) {
1513 return EIGEN_MATHFUNC_IMPL(
polygamma, Scalar)::run(n, x);
1516 template <
typename Scalar>
1517 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(erf, Scalar)
1518 erf(const Scalar& x) {
1519 return EIGEN_MATHFUNC_IMPL(erf, Scalar)::run(x);
1522 template <
typename Scalar>
1523 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar)
1524 erfc(const Scalar& x) {
1525 return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x);
1528 template <
typename Scalar>
1529 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
igamma, Scalar)
1530 igamma(const Scalar& a, const Scalar& x) {
1531 return EIGEN_MATHFUNC_IMPL(
igamma, Scalar)::run(a, x);
1534 template <
typename Scalar>
1535 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
igammac, Scalar)
1536 igammac(const Scalar& a, const Scalar& x) {
1537 return EIGEN_MATHFUNC_IMPL(
igammac, Scalar)::run(a, x);
1540 template <
typename Scalar>
1541 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
betainc, Scalar)
1542 betainc(const Scalar& a, const Scalar& b, const Scalar& x) {
1543 return EIGEN_MATHFUNC_IMPL(
betainc, Scalar)::run(a, b, x);
1551 #endif // EIGEN_SPECIAL_FUNCTIONS_H Namespace containing all symbols from the Eigen library.
Definition: AdolcForward:45
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_igammac_op< typename Derived::Scalar >, const Derived, const ExponentDerived > igammac(const Eigen::ArrayBase< Derived > &a, const Eigen::ArrayBase< ExponentDerived > &x)
Definition: SpecialFunctionsArrayAPI.h:48
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_igamma_op< typename Derived::Scalar >, const Derived, const ExponentDerived > igamma(const Eigen::ArrayBase< Derived > &a, const Eigen::ArrayBase< ExponentDerived > &x)
Definition: SpecialFunctionsArrayAPI.h:28
const TensorCwiseTernaryOp< internal::scalar_betainc_op< typename XDerived::Scalar >, const ADerived, const BDerived, const XDerived > betainc(const ADerived &a, const BDerived &b, const XDerived &x)
Definition: TensorGlobalFunctions.h:24
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:114
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_polygamma_op< typename DerivedX::Scalar >, const DerivedN, const DerivedX > polygamma(const Eigen::ArrayBase< DerivedN > &n, const Eigen::ArrayBase< DerivedX > &x)
Definition: SpecialFunctionsArrayAPI.h:70