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) {
124 #if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__) 126 return ::lgammaf_r(x, &signgam);
134 struct lgamma_impl<double> {
136 static EIGEN_STRONG_INLINE
double run(
double x) {
137 #if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__) 139 return ::lgamma_r(x, &signgam);
151 template <
typename Scalar>
152 struct digamma_retval {
169 template <
typename Scalar>
170 struct digamma_impl_maybe_poly {
172 static EIGEN_STRONG_INLINE Scalar run(
const Scalar) {
173 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
174 THIS_TYPE_IS_NOT_SUPPORTED);
181 struct digamma_impl_maybe_poly<float> {
183 static EIGEN_STRONG_INLINE
float run(
const float s) {
185 -4.16666666666666666667E-3f,
186 3.96825396825396825397E-3f,
187 -8.33333333333333333333E-3f,
188 8.33333333333333333333E-2f
194 return z * cephes::polevl<float, 3>::run(z, A);
200 struct digamma_impl_maybe_poly<double> {
202 static EIGEN_STRONG_INLINE
double run(
const double s) {
204 8.33333333333333333333E-2,
205 -2.10927960927960927961E-2,
206 7.57575757575757575758E-3,
207 -4.16666666666666666667E-3,
208 3.96825396825396825397E-3,
209 -8.33333333333333333333E-3,
210 8.33333333333333333333E-2
216 return z * cephes::polevl<double, 6>::run(z, A);
222 template <
typename Scalar>
223 struct digamma_impl {
225 static Scalar run(Scalar x) {
283 Scalar p, q, nz, s, w, y;
284 bool negative =
false;
286 const Scalar maxnum = NumTraits<Scalar>::infinity();
287 const Scalar m_pi = Scalar(EIGEN_PI);
289 const Scalar zero = Scalar(0);
290 const Scalar one = Scalar(1);
291 const Scalar half = Scalar(0.5);
297 p = numext::floor(q);
310 nz = m_pi / numext::tan(m_pi * nz);
321 while (s < Scalar(10)) {
326 y = digamma_impl_maybe_poly<Scalar>::run(s);
328 y = numext::log(s) - (half / s) - y - w;
330 return (negative) ? y - nz : y;
338 template <
typename Scalar>
341 static EIGEN_STRONG_INLINE Scalar run(
const Scalar) {
342 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
343 THIS_TYPE_IS_NOT_SUPPORTED);
348 template <
typename Scalar>
353 #if EIGEN_HAS_C99_MATH 355 struct erf_impl<float> {
357 static EIGEN_STRONG_INLINE
float run(
float x) { return ::erff(x); }
361 struct erf_impl<double> {
363 static EIGEN_STRONG_INLINE
double run(
double x) { return ::erf(x); }
365 #endif // EIGEN_HAS_C99_MATH 371 template <
typename Scalar>
374 static EIGEN_STRONG_INLINE Scalar run(
const Scalar) {
375 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
376 THIS_TYPE_IS_NOT_SUPPORTED);
381 template <
typename Scalar>
386 #if EIGEN_HAS_C99_MATH 388 struct erfc_impl<float> {
390 static EIGEN_STRONG_INLINE
float run(
const float x) { return ::erfcf(x); }
394 struct erfc_impl<double> {
396 static EIGEN_STRONG_INLINE
double run(
const double x) { return ::erfc(x); }
398 #endif // EIGEN_HAS_C99_MATH 404 template <
typename Scalar>
405 struct igammac_retval {
410 template <
typename Scalar>
411 struct cephes_helper {
413 static EIGEN_STRONG_INLINE Scalar machep() { assert(
false &&
"machep not supported for this type");
return 0.0; }
415 static EIGEN_STRONG_INLINE Scalar big() { assert(
false &&
"big not supported for this type");
return 0.0; }
417 static EIGEN_STRONG_INLINE Scalar biginv() { assert(
false &&
"biginv not supported for this type");
return 0.0; }
421 struct cephes_helper<float> {
423 static EIGEN_STRONG_INLINE
float machep() {
424 return NumTraits<float>::epsilon() / 2;
427 static EIGEN_STRONG_INLINE
float big() {
429 return 1.0f / (NumTraits<float>::epsilon() / 2);
432 static EIGEN_STRONG_INLINE
float biginv() {
439 struct cephes_helper<double> {
441 static EIGEN_STRONG_INLINE
double machep() {
442 return NumTraits<double>::epsilon() / 2;
445 static EIGEN_STRONG_INLINE
double big() {
446 return 1.0 / NumTraits<double>::epsilon();
449 static EIGEN_STRONG_INLINE
double biginv() {
451 return NumTraits<double>::epsilon();
455 #if !EIGEN_HAS_C99_MATH 457 template <
typename Scalar>
458 struct igammac_impl {
460 static Scalar run(Scalar a, Scalar x) {
461 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
462 THIS_TYPE_IS_NOT_SUPPORTED);
469 template <
typename Scalar>
struct igamma_impl;
471 template <
typename Scalar>
472 struct igammac_impl {
474 static Scalar run(Scalar a, Scalar x) {
529 const Scalar zero = 0;
530 const Scalar one = 1;
531 const Scalar nan = NumTraits<Scalar>::quiet_NaN();
533 if ((x < zero) || (a <= zero)) {
538 if ((x < one) || (x < a)) {
546 return (one - igamma_impl<Scalar>::Impl(a, x));
554 friend struct igamma_impl<Scalar>;
563 EIGEN_DEVICE_FUNC
static Scalar Impl(Scalar a, Scalar x) {
564 const Scalar zero = 0;
565 const Scalar one = 1;
566 const Scalar two = 2;
567 const Scalar machep = cephes_helper<Scalar>::machep();
568 const Scalar maxlog = numext::log(NumTraits<Scalar>::highest());
569 const Scalar big = cephes_helper<Scalar>::big();
570 const Scalar biginv = cephes_helper<Scalar>::biginv();
571 const Scalar inf = NumTraits<Scalar>::infinity();
573 Scalar ans, ax, c, yc, r, t, y, z;
574 Scalar pk, pkm1, pkm2, qk, qkm1, qkm2;
576 if (x == inf)
return zero;
579 ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
583 ax = numext::exp(ax);
600 pk = pkm1 * z - pkm2 * yc;
601 qk = qkm1 * z - qkm2 * yc;
604 t = numext::abs((ans - r) / r);
613 if (numext::abs(pk) > big) {
628 #endif // EIGEN_HAS_C99_MATH 634 template <
typename Scalar>
635 struct igamma_retval {
639 #if !EIGEN_HAS_C99_MATH 641 template <
typename Scalar>
644 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) {
645 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
646 THIS_TYPE_IS_NOT_SUPPORTED);
653 template <
typename Scalar>
656 static Scalar run(Scalar a, Scalar x) {
717 const Scalar zero = 0;
718 const Scalar one = 1;
719 const Scalar nan = NumTraits<Scalar>::quiet_NaN();
721 if (x == zero)
return zero;
723 if ((x < zero) || (a <= zero)) {
727 if ((x > one) && (x > a)) {
735 return (one - igammac_impl<Scalar>::Impl(a, x));
743 friend struct igammac_impl<Scalar>;
752 EIGEN_DEVICE_FUNC
static Scalar Impl(Scalar a, Scalar x) {
753 const Scalar zero = 0;
754 const Scalar one = 1;
755 const Scalar machep = cephes_helper<Scalar>::machep();
756 const Scalar maxlog = numext::log(NumTraits<Scalar>::highest());
758 Scalar ans, ax, c, r;
761 ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
766 ax = numext::exp(ax);
777 if (c/ans <= machep) {
782 return (ans * ax / a);
786 #endif // EIGEN_HAS_C99_MATH 792 template <
typename Scalar>
797 template <
typename Scalar>
798 struct zeta_impl_series {
800 static EIGEN_STRONG_INLINE Scalar run(
const Scalar) {
801 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
802 THIS_TYPE_IS_NOT_SUPPORTED);
808 struct zeta_impl_series<float> {
810 static EIGEN_STRONG_INLINE
bool run(
float& a,
float& b,
float& s,
const float x,
const float machep) {
816 b = numext::pow( a, -x );
818 if( numext::abs(b/s) < machep )
828 struct zeta_impl_series<double> {
830 static EIGEN_STRONG_INLINE
bool run(
double& a,
double& b,
double& s,
const double x,
const double machep) {
832 while( (i < 9) || (a <= 9.0) )
836 b = numext::pow( a, -x );
838 if( numext::abs(b/s) < machep )
847 template <
typename Scalar>
850 static Scalar run(Scalar x, Scalar q) {
913 Scalar p, r, a, b, k, s, t, w;
921 Scalar(-1.8924375803183791606e9),
922 Scalar(7.47242496e10),
923 Scalar(-2.950130727918164224e12),
924 Scalar(1.1646782814350067249e14),
925 Scalar(-4.5979787224074726105e15),
926 Scalar(1.8152105401943546773e17),
927 Scalar(-7.1661652561756670113e18)
930 const Scalar maxnum = NumTraits<Scalar>::infinity();
931 const Scalar zero = 0.0, half = 0.5, one = 1.0;
932 const Scalar machep = cephes_helper<Scalar>::machep();
933 const Scalar nan = NumTraits<Scalar>::quiet_NaN();
945 if(q == numext::floor(q))
950 r = numext::floor(p);
960 s = numext::pow( q, -x );
964 if (zeta_impl_series<Scalar>::run(a, b, s, x, machep)) {
973 for( i=0; i<12; i++ )
979 t = numext::abs(t/s);
996 template <
typename Scalar>
997 struct polygamma_retval {
1001 #if !EIGEN_HAS_C99_MATH 1003 template <
typename Scalar>
1004 struct polygamma_impl {
1006 static EIGEN_STRONG_INLINE Scalar run(Scalar n, Scalar x) {
1007 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
1008 THIS_TYPE_IS_NOT_SUPPORTED);
1015 template <
typename Scalar>
1016 struct polygamma_impl {
1018 static Scalar run(Scalar n, Scalar x) {
1019 Scalar zero = 0.0, one = 1.0;
1020 Scalar nplus = n + one;
1021 const Scalar nan = NumTraits<Scalar>::quiet_NaN();
1024 if (numext::floor(n) != n) {
1028 else if (n == zero) {
1029 return digamma_impl<Scalar>::run(x);
1033 Scalar factorial = numext::exp(lgamma_impl<Scalar>::run(nplus));
1034 return numext::pow(-one, nplus) * factorial * zeta_impl<Scalar>::run(nplus, x);
1039 #endif // EIGEN_HAS_C99_MATH 1045 template <
typename Scalar>
1046 struct betainc_retval {
1047 typedef Scalar type;
1050 #if !EIGEN_HAS_C99_MATH 1052 template <
typename Scalar>
1053 struct betainc_impl {
1055 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x) {
1056 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
1057 THIS_TYPE_IS_NOT_SUPPORTED);
1064 template <
typename Scalar>
1065 struct betainc_impl {
1067 static EIGEN_STRONG_INLINE Scalar run(Scalar, Scalar, Scalar) {
1137 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value ==
false),
1138 THIS_TYPE_IS_NOT_SUPPORTED);
1146 template <
typename Scalar>
1147 struct incbeta_cfe {
1149 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x,
bool small_branch) {
1150 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, float>::value ||
1151 internal::is_same<Scalar, double>::value),
1152 THIS_TYPE_IS_NOT_SUPPORTED);
1153 const Scalar big = cephes_helper<Scalar>::big();
1154 const Scalar machep = cephes_helper<Scalar>::machep();
1155 const Scalar biginv = cephes_helper<Scalar>::biginv();
1157 const Scalar zero = 0;
1158 const Scalar one = 1;
1159 const Scalar two = 2;
1161 Scalar xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
1162 Scalar k1, k2, k3, k4, k5, k6, k7, k8, k26update;
1166 const int num_iters = (internal::is_same<Scalar, float>::value) ? 100 : 300;
1167 const Scalar thresh =
1168 (internal::is_same<Scalar, float>::value) ? machep : Scalar(3) * machep;
1169 Scalar r = (internal::is_same<Scalar, float>::value) ? zero : one;
1202 xk = -(x * k1 * k2) / (k3 * k4);
1203 pk = pkm1 + pkm2 * xk;
1204 qk = qkm1 + qkm2 * xk;
1210 xk = (x * k5 * k6) / (k7 * k8);
1211 pk = pkm1 + pkm2 * xk;
1212 qk = qkm1 + qkm2 * xk;
1220 if (numext::abs(ans - r) < numext::abs(r) * thresh) {
1235 if ((numext::abs(qk) + numext::abs(pk)) > big) {
1241 if ((numext::abs(qk) < biginv) || (numext::abs(pk) < biginv)) {
1247 }
while (++n < num_iters);
1254 template <
typename Scalar>
1255 struct betainc_helper {};
1258 struct betainc_helper<float> {
1260 EIGEN_DEVICE_FUNC
static EIGEN_STRONG_INLINE
float incbsa(
float aa,
float bb,
1262 float ans, a, b, t, x, onemx;
1263 bool reversed_a_b =
false;
1268 if (xx > (aa / (aa + bb))) {
1269 reversed_a_b =
true;
1283 if (numext::abs(b * x / a) < 0.3f) {
1284 t = betainc_helper<float>::incbps(a, b, x);
1285 if (reversed_a_b) t = 1.0f - t;
1290 ans = x * (a + b - 2.0f) / (a - 1.0f);
1292 ans = incbeta_cfe<float>::run(a, b, x,
true );
1293 t = b * numext::log(t);
1295 ans = incbeta_cfe<float>::run(a, b, x,
false );
1296 t = (b - 1.0f) * numext::log(t);
1299 t += a * numext::log(x) + lgamma_impl<float>::run(a + b) -
1300 lgamma_impl<float>::run(a) - lgamma_impl<float>::run(b);
1301 t += numext::log(ans / a);
1304 if (reversed_a_b) t = 1.0f - t;
1309 static EIGEN_STRONG_INLINE
float incbps(
float a,
float b,
float x) {
1311 const float machep = cephes_helper<float>::machep();
1313 y = a * numext::log(x) + (b - 1.0f) * numext::log1p(-x) - numext::log(a);
1314 y -= lgamma_impl<float>::run(a) + lgamma_impl<float>::run(b);
1315 y += lgamma_impl<float>::run(a + b);
1328 }
while (numext::abs(u) > machep);
1330 return numext::exp(y) * (1.0f + s);
1335 struct betainc_impl<float> {
1337 static float run(
float a,
float b,
float x) {
1338 const float nan = NumTraits<float>::quiet_NaN();
1341 if (a <= 0.0f)
return nan;
1342 if (b <= 0.0f)
return nan;
1343 if ((x <= 0.0f) || (x >= 1.0f)) {
1344 if (x == 0.0f)
return 0.0f;
1345 if (x == 1.0f)
return 1.0f;
1352 ans = betainc_helper<float>::incbsa(a + 1.0f, b, x);
1353 t = a * numext::log(x) + b * numext::log1p(-x) +
1354 lgamma_impl<float>::run(a + b) - lgamma_impl<float>::run(a + 1.0f) -
1355 lgamma_impl<float>::run(b);
1356 return (ans + numext::exp(t));
1358 return betainc_helper<float>::incbsa(a, b, x);
1364 struct betainc_helper<double> {
1366 static EIGEN_STRONG_INLINE
double incbps(
double a,
double b,
double x) {
1367 const double machep = cephes_helper<double>::machep();
1369 double s, t, u, v, n, t1, z, ai;
1379 while (numext::abs(v) > z) {
1380 u = (n - b) * x / n;
1389 u = a * numext::log(x);
1397 t = lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) -
1398 lgamma_impl<double>::run(b) + u + numext::log(s);
1399 return s = numext::exp(t);
1404 struct betainc_impl<double> {
1406 static double run(
double aa,
double bb,
double xx) {
1407 const double nan = NumTraits<double>::quiet_NaN();
1408 const double machep = cephes_helper<double>::machep();
1411 double a, b, t, x, xc, w, y;
1412 bool reversed_a_b =
false;
1414 if (aa <= 0.0 || bb <= 0.0) {
1418 if ((xx <= 0.0) || (xx >= 1.0)) {
1419 if (xx == 0.0)
return (0.0);
1420 if (xx == 1.0)
return (1.0);
1425 if ((bb * xx) <= 1.0 && xx <= 0.95) {
1426 return betainc_helper<double>::incbps(aa, bb, xx);
1432 if (xx > (aa / (aa + bb))) {
1433 reversed_a_b =
true;
1445 if (reversed_a_b && (b * x) <= 1.0 && x <= 0.95) {
1446 t = betainc_helper<double>::incbps(a, b, x);
1456 y = x * (a + b - 2.0) - (a - 1.0);
1458 w = incbeta_cfe<double>::run(a, b, x,
true );
1460 w = incbeta_cfe<double>::run(a, b, x,
false ) / xc;
1467 y = a * numext::log(x);
1468 t = b * numext::log(xc);
1481 y += t + lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) -
1482 lgamma_impl<double>::run(b);
1483 y += numext::log(w / a);
1500 #endif // EIGEN_HAS_C99_MATH 1506 template <
typename Scalar>
1507 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(lgamma, Scalar)
1508 lgamma(
const Scalar& x) {
1509 return EIGEN_MATHFUNC_IMPL(lgamma, Scalar)::run(x);
1512 template <
typename Scalar>
1513 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(digamma, Scalar)
1514 digamma(
const Scalar& x) {
1515 return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x);
1518 template <
typename Scalar>
1519 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
zeta, Scalar)
1520 zeta(
const Scalar& x,
const Scalar& q) {
1521 return EIGEN_MATHFUNC_IMPL(
zeta, Scalar)::run(x, q);
1524 template <
typename Scalar>
1525 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
polygamma, Scalar)
1526 polygamma(
const Scalar& n,
const Scalar& x) {
1527 return EIGEN_MATHFUNC_IMPL(
polygamma, Scalar)::run(n, x);
1530 template <
typename Scalar>
1531 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(erf, Scalar)
1532 erf(
const Scalar& x) {
1533 return EIGEN_MATHFUNC_IMPL(erf, Scalar)::run(x);
1536 template <
typename Scalar>
1537 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar)
1538 erfc(
const Scalar& x) {
1539 return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x);
1542 template <
typename Scalar>
1543 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
igamma, Scalar)
1544 igamma(
const Scalar& a,
const Scalar& x) {
1545 return EIGEN_MATHFUNC_IMPL(
igamma, Scalar)::run(a, x);
1548 template <
typename Scalar>
1549 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
igammac, Scalar)
1550 igammac(
const Scalar& a,
const Scalar& x) {
1551 return EIGEN_MATHFUNC_IMPL(
igammac, Scalar)::run(a, x);
1554 template <
typename Scalar>
1555 EIGEN_DEVICE_FUNC
inline EIGEN_MATHFUNC_RETVAL(
betainc, Scalar)
1556 betainc(
const Scalar& a,
const Scalar& b,
const Scalar& x) {
1557 return EIGEN_MATHFUNC_IMPL(
betainc, Scalar)::run(a, b, x);
1565 #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