![]() |
Eigen-unsupported
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2015 Eugene Brevdo <ebrevdo@gmail.com> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_SPECIAL_FUNCTIONS_H 00011 #define EIGEN_SPECIAL_FUNCTIONS_H 00012 00013 namespace Eigen { 00014 namespace internal { 00015 00016 // Parts of this code are based on the Cephes Math Library. 00017 // 00018 // Cephes Math Library Release 2.8: June, 2000 00019 // Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier 00020 // 00021 // Permission has been kindly provided by the original author 00022 // to incorporate the Cephes software into the Eigen codebase: 00023 // 00024 // From: Stephen Moshier 00025 // To: Eugene Brevdo 00026 // Subject: Re: Permission to wrap several cephes functions in Eigen 00027 // 00028 // Hello Eugene, 00029 // 00030 // Thank you for writing. 00031 // 00032 // If your licensing is similar to BSD, the formal way that has been 00033 // handled is simply to add a statement to the effect that you are incorporating 00034 // the Cephes software by permission of the author. 00035 // 00036 // Good luck with your project, 00037 // Steve 00038 00039 namespace cephes { 00040 00041 /* polevl (modified for Eigen) 00042 * 00043 * Evaluate polynomial 00044 * 00045 * 00046 * 00047 * SYNOPSIS: 00048 * 00049 * int N; 00050 * Scalar x, y, coef[N+1]; 00051 * 00052 * y = polevl<decltype(x), N>( x, coef); 00053 * 00054 * 00055 * 00056 * DESCRIPTION: 00057 * 00058 * Evaluates polynomial of degree N: 00059 * 00060 * 2 N 00061 * y = C + C x + C x +...+ C x 00062 * 0 1 2 N 00063 * 00064 * Coefficients are stored in reverse order: 00065 * 00066 * coef[0] = C , ..., coef[N] = C . 00067 * N 0 00068 * 00069 * The function p1evl() assumes that coef[N] = 1.0 and is 00070 * omitted from the array. Its calling arguments are 00071 * otherwise the same as polevl(). 00072 * 00073 * 00074 * The Eigen implementation is templatized. For best speed, store 00075 * coef as a const array (constexpr), e.g. 00076 * 00077 * const double coef[] = {1.0, 2.0, 3.0, ...}; 00078 * 00079 */ 00080 template <typename Scalar, int N> 00081 struct polevl { 00082 EIGEN_DEVICE_FUNC 00083 static EIGEN_STRONG_INLINE Scalar run(const Scalar x, const Scalar coef[]) { 00084 EIGEN_STATIC_ASSERT((N > 0), YOU_MADE_A_PROGRAMMING_MISTAKE); 00085 00086 return polevl<Scalar, N - 1>::run(x, coef) * x + coef[N]; 00087 } 00088 }; 00089 00090 template <typename Scalar> 00091 struct polevl<Scalar, 0> { 00092 EIGEN_DEVICE_FUNC 00093 static EIGEN_STRONG_INLINE Scalar run(const Scalar, const Scalar coef[]) { 00094 return coef[0]; 00095 } 00096 }; 00097 00098 } // end namespace cephes 00099 00100 /**************************************************************************** 00101 * Implementation of lgamma, requires C++11/C99 * 00102 ****************************************************************************/ 00103 00104 template <typename Scalar> 00105 struct lgamma_impl { 00106 EIGEN_DEVICE_FUNC 00107 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 00108 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 00109 THIS_TYPE_IS_NOT_SUPPORTED); 00110 return Scalar(0); 00111 } 00112 }; 00113 00114 template <typename Scalar> 00115 struct lgamma_retval { 00116 typedef Scalar type; 00117 }; 00118 00119 #if EIGEN_HAS_C99_MATH 00120 template <> 00121 struct lgamma_impl<float> { 00122 EIGEN_DEVICE_FUNC 00123 static EIGEN_STRONG_INLINE float run(float x) { 00124 #if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__) 00125 int signgam; 00126 return ::lgammaf_r(x, &signgam); 00127 #else 00128 return ::lgammaf(x); 00129 #endif 00130 } 00131 }; 00132 00133 template <> 00134 struct lgamma_impl<double> { 00135 EIGEN_DEVICE_FUNC 00136 static EIGEN_STRONG_INLINE double run(double x) { 00137 #if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__) 00138 int signgam; 00139 return ::lgamma_r(x, &signgam); 00140 #else 00141 return ::lgamma(x); 00142 #endif 00143 } 00144 }; 00145 #endif 00146 00147 /**************************************************************************** 00148 * Implementation of digamma (psi), based on Cephes * 00149 ****************************************************************************/ 00150 00151 template <typename Scalar> 00152 struct digamma_retval { 00153 typedef Scalar type; 00154 }; 00155 00156 /* 00157 * 00158 * Polynomial evaluation helper for the Psi (digamma) function. 00159 * 00160 * digamma_impl_maybe_poly::run(s) evaluates the asymptotic Psi expansion for 00161 * input Scalar s, assuming s is above 10.0. 00162 * 00163 * If s is above a certain threshold for the given Scalar type, zero 00164 * is returned. Otherwise the polynomial is evaluated with enough 00165 * coefficients for results matching Scalar machine precision. 00166 * 00167 * 00168 */ 00169 template <typename Scalar> 00170 struct digamma_impl_maybe_poly { 00171 EIGEN_DEVICE_FUNC 00172 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 00173 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 00174 THIS_TYPE_IS_NOT_SUPPORTED); 00175 return Scalar(0); 00176 } 00177 }; 00178 00179 00180 template <> 00181 struct digamma_impl_maybe_poly<float> { 00182 EIGEN_DEVICE_FUNC 00183 static EIGEN_STRONG_INLINE float run(const float s) { 00184 const float A[] = { 00185 -4.16666666666666666667E-3f, 00186 3.96825396825396825397E-3f, 00187 -8.33333333333333333333E-3f, 00188 8.33333333333333333333E-2f 00189 }; 00190 00191 float z; 00192 if (s < 1.0e8f) { 00193 z = 1.0f / (s * s); 00194 return z * cephes::polevl<float, 3>::run(z, A); 00195 } else return 0.0f; 00196 } 00197 }; 00198 00199 template <> 00200 struct digamma_impl_maybe_poly<double> { 00201 EIGEN_DEVICE_FUNC 00202 static EIGEN_STRONG_INLINE double run(const double s) { 00203 const double A[] = { 00204 8.33333333333333333333E-2, 00205 -2.10927960927960927961E-2, 00206 7.57575757575757575758E-3, 00207 -4.16666666666666666667E-3, 00208 3.96825396825396825397E-3, 00209 -8.33333333333333333333E-3, 00210 8.33333333333333333333E-2 00211 }; 00212 00213 double z; 00214 if (s < 1.0e17) { 00215 z = 1.0 / (s * s); 00216 return z * cephes::polevl<double, 6>::run(z, A); 00217 } 00218 else return 0.0; 00219 } 00220 }; 00221 00222 template <typename Scalar> 00223 struct digamma_impl { 00224 EIGEN_DEVICE_FUNC 00225 static Scalar run(Scalar x) { 00226 /* 00227 * 00228 * Psi (digamma) function (modified for Eigen) 00229 * 00230 * 00231 * SYNOPSIS: 00232 * 00233 * double x, y, psi(); 00234 * 00235 * y = psi( x ); 00236 * 00237 * 00238 * DESCRIPTION: 00239 * 00240 * d - 00241 * psi(x) = -- ln | (x) 00242 * dx 00243 * 00244 * is the logarithmic derivative of the gamma function. 00245 * For integer x, 00246 * n-1 00247 * - 00248 * psi(n) = -EUL + > 1/k. 00249 * - 00250 * k=1 00251 * 00252 * If x is negative, it is transformed to a positive argument by the 00253 * reflection formula psi(1-x) = psi(x) + pi cot(pi x). 00254 * For general positive x, the argument is made greater than 10 00255 * using the recurrence psi(x+1) = psi(x) + 1/x. 00256 * Then the following asymptotic expansion is applied: 00257 * 00258 * inf. B 00259 * - 2k 00260 * psi(x) = log(x) - 1/2x - > ------- 00261 * - 2k 00262 * k=1 2k x 00263 * 00264 * where the B2k are Bernoulli numbers. 00265 * 00266 * ACCURACY (float): 00267 * Relative error (except absolute when |psi| < 1): 00268 * arithmetic domain # trials peak rms 00269 * IEEE 0,30 30000 1.3e-15 1.4e-16 00270 * IEEE -30,0 40000 1.5e-15 2.2e-16 00271 * 00272 * ACCURACY (double): 00273 * Absolute error, relative when |psi| > 1 : 00274 * arithmetic domain # trials peak rms 00275 * IEEE -33,0 30000 8.2e-7 1.2e-7 00276 * IEEE 0,33 100000 7.3e-7 7.7e-8 00277 * 00278 * ERROR MESSAGES: 00279 * message condition value returned 00280 * psi singularity x integer <=0 INFINITY 00281 */ 00282 00283 Scalar p, q, nz, s, w, y; 00284 bool negative = false; 00285 00286 const Scalar maxnum = NumTraits<Scalar>::infinity(); 00287 const Scalar m_pi = Scalar(EIGEN_PI); 00288 00289 const Scalar zero = Scalar(0); 00290 const Scalar one = Scalar(1); 00291 const Scalar half = Scalar(0.5); 00292 nz = zero; 00293 00294 if (x <= zero) { 00295 negative = true; 00296 q = x; 00297 p = numext::floor(q); 00298 if (p == q) { 00299 return maxnum; 00300 } 00301 /* Remove the zeros of tan(m_pi x) 00302 * by subtracting the nearest integer from x 00303 */ 00304 nz = q - p; 00305 if (nz != half) { 00306 if (nz > half) { 00307 p += one; 00308 nz = q - p; 00309 } 00310 nz = m_pi / numext::tan(m_pi * nz); 00311 } 00312 else { 00313 nz = zero; 00314 } 00315 x = one - x; 00316 } 00317 00318 /* use the recurrence psi(x+1) = psi(x) + 1/x. */ 00319 s = x; 00320 w = zero; 00321 while (s < Scalar(10)) { 00322 w += one / s; 00323 s += one; 00324 } 00325 00326 y = digamma_impl_maybe_poly<Scalar>::run(s); 00327 00328 y = numext::log(s) - (half / s) - y - w; 00329 00330 return (negative) ? y - nz : y; 00331 } 00332 }; 00333 00334 /**************************************************************************** 00335 * Implementation of erf, requires C++11/C99 * 00336 ****************************************************************************/ 00337 00338 template <typename Scalar> 00339 struct erf_impl { 00340 EIGEN_DEVICE_FUNC 00341 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 00342 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 00343 THIS_TYPE_IS_NOT_SUPPORTED); 00344 return Scalar(0); 00345 } 00346 }; 00347 00348 template <typename Scalar> 00349 struct erf_retval { 00350 typedef Scalar type; 00351 }; 00352 00353 #if EIGEN_HAS_C99_MATH 00354 template <> 00355 struct erf_impl<float> { 00356 EIGEN_DEVICE_FUNC 00357 static EIGEN_STRONG_INLINE float run(float x) { return ::erff(x); } 00358 }; 00359 00360 template <> 00361 struct erf_impl<double> { 00362 EIGEN_DEVICE_FUNC 00363 static EIGEN_STRONG_INLINE double run(double x) { return ::erf(x); } 00364 }; 00365 #endif // EIGEN_HAS_C99_MATH 00366 00367 /*************************************************************************** 00368 * Implementation of erfc, requires C++11/C99 * 00369 ****************************************************************************/ 00370 00371 template <typename Scalar> 00372 struct erfc_impl { 00373 EIGEN_DEVICE_FUNC 00374 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 00375 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 00376 THIS_TYPE_IS_NOT_SUPPORTED); 00377 return Scalar(0); 00378 } 00379 }; 00380 00381 template <typename Scalar> 00382 struct erfc_retval { 00383 typedef Scalar type; 00384 }; 00385 00386 #if EIGEN_HAS_C99_MATH 00387 template <> 00388 struct erfc_impl<float> { 00389 EIGEN_DEVICE_FUNC 00390 static EIGEN_STRONG_INLINE float run(const float x) { return ::erfcf(x); } 00391 }; 00392 00393 template <> 00394 struct erfc_impl<double> { 00395 EIGEN_DEVICE_FUNC 00396 static EIGEN_STRONG_INLINE double run(const double x) { return ::erfc(x); } 00397 }; 00398 #endif // EIGEN_HAS_C99_MATH 00399 00400 /************************************************************************************************************** 00401 * Implementation of igammac (complemented incomplete gamma integral), based on Cephes but requires C++11/C99 * 00402 **************************************************************************************************************/ 00403 00404 template <typename Scalar> 00405 struct igammac_retval { 00406 typedef Scalar type; 00407 }; 00408 00409 // NOTE: cephes_helper is also used to implement zeta 00410 template <typename Scalar> 00411 struct cephes_helper { 00412 EIGEN_DEVICE_FUNC 00413 static EIGEN_STRONG_INLINE Scalar machep() { assert(false && "machep not supported for this type"); return 0.0; } 00414 EIGEN_DEVICE_FUNC 00415 static EIGEN_STRONG_INLINE Scalar big() { assert(false && "big not supported for this type"); return 0.0; } 00416 EIGEN_DEVICE_FUNC 00417 static EIGEN_STRONG_INLINE Scalar biginv() { assert(false && "biginv not supported for this type"); return 0.0; } 00418 }; 00419 00420 template <> 00421 struct cephes_helper<float> { 00422 EIGEN_DEVICE_FUNC 00423 static EIGEN_STRONG_INLINE float machep() { 00424 return NumTraits<float>::epsilon() / 2; // 1.0 - machep == 1.0 00425 } 00426 EIGEN_DEVICE_FUNC 00427 static EIGEN_STRONG_INLINE float big() { 00428 // use epsneg (1.0 - epsneg == 1.0) 00429 return 1.0f / (NumTraits<float>::epsilon() / 2); 00430 } 00431 EIGEN_DEVICE_FUNC 00432 static EIGEN_STRONG_INLINE float biginv() { 00433 // epsneg 00434 return machep(); 00435 } 00436 }; 00437 00438 template <> 00439 struct cephes_helper<double> { 00440 EIGEN_DEVICE_FUNC 00441 static EIGEN_STRONG_INLINE double machep() { 00442 return NumTraits<double>::epsilon() / 2; // 1.0 - machep == 1.0 00443 } 00444 EIGEN_DEVICE_FUNC 00445 static EIGEN_STRONG_INLINE double big() { 00446 return 1.0 / NumTraits<double>::epsilon(); 00447 } 00448 EIGEN_DEVICE_FUNC 00449 static EIGEN_STRONG_INLINE double biginv() { 00450 // inverse of eps 00451 return NumTraits<double>::epsilon(); 00452 } 00453 }; 00454 00455 #if !EIGEN_HAS_C99_MATH 00456 00457 template <typename Scalar> 00458 struct igammac_impl { 00459 EIGEN_DEVICE_FUNC 00460 static Scalar run(Scalar a, Scalar x) { 00461 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 00462 THIS_TYPE_IS_NOT_SUPPORTED); 00463 return Scalar(0); 00464 } 00465 }; 00466 00467 #else 00468 00469 template <typename Scalar> struct igamma_impl; // predeclare igamma_impl 00470 00471 template <typename Scalar> 00472 struct igammac_impl { 00473 EIGEN_DEVICE_FUNC 00474 static Scalar run(Scalar a, Scalar x) { 00475 /* igamc() 00476 * 00477 * Incomplete gamma integral (modified for Eigen) 00478 * 00479 * 00480 * 00481 * SYNOPSIS: 00482 * 00483 * double a, x, y, igamc(); 00484 * 00485 * y = igamc( a, x ); 00486 * 00487 * DESCRIPTION: 00488 * 00489 * The function is defined by 00490 * 00491 * 00492 * igamc(a,x) = 1 - igam(a,x) 00493 * 00494 * inf. 00495 * - 00496 * 1 | | -t a-1 00497 * = ----- | e t dt. 00498 * - | | 00499 * | (a) - 00500 * x 00501 * 00502 * 00503 * In this implementation both arguments must be positive. 00504 * The integral is evaluated by either a power series or 00505 * continued fraction expansion, depending on the relative 00506 * values of a and x. 00507 * 00508 * ACCURACY (float): 00509 * 00510 * Relative error: 00511 * arithmetic domain # trials peak rms 00512 * IEEE 0,30 30000 7.8e-6 5.9e-7 00513 * 00514 * 00515 * ACCURACY (double): 00516 * 00517 * Tested at random a, x. 00518 * a x Relative error: 00519 * arithmetic domain domain # trials peak rms 00520 * IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15 00521 * IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15 00522 * 00523 */ 00524 /* 00525 Cephes Math Library Release 2.2: June, 1992 00526 Copyright 1985, 1987, 1992 by Stephen L. Moshier 00527 Direct inquiries to 30 Frost Street, Cambridge, MA 02140 00528 */ 00529 const Scalar zero = 0; 00530 const Scalar one = 1; 00531 const Scalar nan = NumTraits<Scalar>::quiet_NaN(); 00532 00533 if ((x < zero) || (a <= zero)) { 00534 // domain error 00535 return nan; 00536 } 00537 00538 if ((x < one) || (x < a)) { 00539 /* The checks above ensure that we meet the preconditions for 00540 * igamma_impl::Impl(), so call it, rather than igamma_impl::Run(). 00541 * Calling Run() would also work, but in that case the compiler may not be 00542 * able to prove that igammac_impl::Run and igamma_impl::Run are not 00543 * mutually recursive. This leads to worse code, particularly on 00544 * platforms like nvptx, where recursion is allowed only begrudgingly. 00545 */ 00546 return (one - igamma_impl<Scalar>::Impl(a, x)); 00547 } 00548 00549 return Impl(a, x); 00550 } 00551 00552 private: 00553 /* igamma_impl calls igammac_impl::Impl. */ 00554 friend struct igamma_impl<Scalar>; 00555 00556 /* Actually computes igamc(a, x). 00557 * 00558 * Preconditions: 00559 * a > 0 00560 * x >= 1 00561 * x >= a 00562 */ 00563 EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) { 00564 const Scalar zero = 0; 00565 const Scalar one = 1; 00566 const Scalar two = 2; 00567 const Scalar machep = cephes_helper<Scalar>::machep(); 00568 const Scalar maxlog = numext::log(NumTraits<Scalar>::highest()); 00569 const Scalar big = cephes_helper<Scalar>::big(); 00570 const Scalar biginv = cephes_helper<Scalar>::biginv(); 00571 const Scalar inf = NumTraits<Scalar>::infinity(); 00572 00573 Scalar ans, ax, c, yc, r, t, y, z; 00574 Scalar pk, pkm1, pkm2, qk, qkm1, qkm2; 00575 00576 if (x == inf) return zero; // std::isinf crashes on CUDA 00577 00578 /* Compute x**a * exp(-x) / gamma(a) */ 00579 ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); 00580 if (ax < -maxlog) { // underflow 00581 return zero; 00582 } 00583 ax = numext::exp(ax); 00584 00585 // continued fraction 00586 y = one - a; 00587 z = x + y + one; 00588 c = zero; 00589 pkm2 = one; 00590 qkm2 = x; 00591 pkm1 = x + one; 00592 qkm1 = z * x; 00593 ans = pkm1 / qkm1; 00594 00595 while (true) { 00596 c += one; 00597 y += one; 00598 z += two; 00599 yc = y * c; 00600 pk = pkm1 * z - pkm2 * yc; 00601 qk = qkm1 * z - qkm2 * yc; 00602 if (qk != zero) { 00603 r = pk / qk; 00604 t = numext::abs((ans - r) / r); 00605 ans = r; 00606 } else { 00607 t = one; 00608 } 00609 pkm2 = pkm1; 00610 pkm1 = pk; 00611 qkm2 = qkm1; 00612 qkm1 = qk; 00613 if (numext::abs(pk) > big) { 00614 pkm2 *= biginv; 00615 pkm1 *= biginv; 00616 qkm2 *= biginv; 00617 qkm1 *= biginv; 00618 } 00619 if (t <= machep) { 00620 break; 00621 } 00622 } 00623 00624 return (ans * ax); 00625 } 00626 }; 00627 00628 #endif // EIGEN_HAS_C99_MATH 00629 00630 /************************************************************************************************ 00631 * Implementation of igamma (incomplete gamma integral), based on Cephes but requires C++11/C99 * 00632 ************************************************************************************************/ 00633 00634 template <typename Scalar> 00635 struct igamma_retval { 00636 typedef Scalar type; 00637 }; 00638 00639 #if !EIGEN_HAS_C99_MATH 00640 00641 template <typename Scalar> 00642 struct igamma_impl { 00643 EIGEN_DEVICE_FUNC 00644 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) { 00645 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 00646 THIS_TYPE_IS_NOT_SUPPORTED); 00647 return Scalar(0); 00648 } 00649 }; 00650 00651 #else 00652 00653 template <typename Scalar> 00654 struct igamma_impl { 00655 EIGEN_DEVICE_FUNC 00656 static Scalar run(Scalar a, Scalar x) { 00657 /* igam() 00658 * Incomplete gamma integral 00659 * 00660 * 00661 * 00662 * SYNOPSIS: 00663 * 00664 * double a, x, y, igam(); 00665 * 00666 * y = igam( a, x ); 00667 * 00668 * DESCRIPTION: 00669 * 00670 * The function is defined by 00671 * 00672 * x 00673 * - 00674 * 1 | | -t a-1 00675 * igam(a,x) = ----- | e t dt. 00676 * - | | 00677 * | (a) - 00678 * 0 00679 * 00680 * 00681 * In this implementation both arguments must be positive. 00682 * The integral is evaluated by either a power series or 00683 * continued fraction expansion, depending on the relative 00684 * values of a and x. 00685 * 00686 * ACCURACY (double): 00687 * 00688 * Relative error: 00689 * arithmetic domain # trials peak rms 00690 * IEEE 0,30 200000 3.6e-14 2.9e-15 00691 * IEEE 0,100 300000 9.9e-14 1.5e-14 00692 * 00693 * 00694 * ACCURACY (float): 00695 * 00696 * Relative error: 00697 * arithmetic domain # trials peak rms 00698 * IEEE 0,30 20000 7.8e-6 5.9e-7 00699 * 00700 */ 00701 /* 00702 Cephes Math Library Release 2.2: June, 1992 00703 Copyright 1985, 1987, 1992 by Stephen L. Moshier 00704 Direct inquiries to 30 Frost Street, Cambridge, MA 02140 00705 */ 00706 00707 00708 /* left tail of incomplete gamma function: 00709 * 00710 * inf. k 00711 * a -x - x 00712 * x e > ---------- 00713 * - - 00714 * k=0 | (a+k+1) 00715 * 00716 */ 00717 const Scalar zero = 0; 00718 const Scalar one = 1; 00719 const Scalar nan = NumTraits<Scalar>::quiet_NaN(); 00720 00721 if (x == zero) return zero; 00722 00723 if ((x < zero) || (a <= zero)) { // domain error 00724 return nan; 00725 } 00726 00727 if ((x > one) && (x > a)) { 00728 /* The checks above ensure that we meet the preconditions for 00729 * igammac_impl::Impl(), so call it, rather than igammac_impl::Run(). 00730 * Calling Run() would also work, but in that case the compiler may not be 00731 * able to prove that igammac_impl::Run and igamma_impl::Run are not 00732 * mutually recursive. This leads to worse code, particularly on 00733 * platforms like nvptx, where recursion is allowed only begrudgingly. 00734 */ 00735 return (one - igammac_impl<Scalar>::Impl(a, x)); 00736 } 00737 00738 return Impl(a, x); 00739 } 00740 00741 private: 00742 /* igammac_impl calls igamma_impl::Impl. */ 00743 friend struct igammac_impl<Scalar>; 00744 00745 /* Actually computes igam(a, x). 00746 * 00747 * Preconditions: 00748 * x > 0 00749 * a > 0 00750 * !(x > 1 && x > a) 00751 */ 00752 EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) { 00753 const Scalar zero = 0; 00754 const Scalar one = 1; 00755 const Scalar machep = cephes_helper<Scalar>::machep(); 00756 const Scalar maxlog = numext::log(NumTraits<Scalar>::highest()); 00757 00758 Scalar ans, ax, c, r; 00759 00760 /* Compute x**a * exp(-x) / gamma(a) */ 00761 ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); 00762 if (ax < -maxlog) { 00763 // underflow 00764 return zero; 00765 } 00766 ax = numext::exp(ax); 00767 00768 /* power series */ 00769 r = a; 00770 c = one; 00771 ans = one; 00772 00773 while (true) { 00774 r += one; 00775 c *= x/r; 00776 ans += c; 00777 if (c/ans <= machep) { 00778 break; 00779 } 00780 } 00781 00782 return (ans * ax / a); 00783 } 00784 }; 00785 00786 #endif // EIGEN_HAS_C99_MATH 00787 00788 /***************************************************************************** 00789 * Implementation of Riemann zeta function of two arguments, based on Cephes * 00790 *****************************************************************************/ 00791 00792 template <typename Scalar> 00793 struct zeta_retval { 00794 typedef Scalar type; 00795 }; 00796 00797 template <typename Scalar> 00798 struct zeta_impl_series { 00799 EIGEN_DEVICE_FUNC 00800 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 00801 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 00802 THIS_TYPE_IS_NOT_SUPPORTED); 00803 return Scalar(0); 00804 } 00805 }; 00806 00807 template <> 00808 struct zeta_impl_series<float> { 00809 EIGEN_DEVICE_FUNC 00810 static EIGEN_STRONG_INLINE bool run(float& a, float& b, float& s, const float x, const float machep) { 00811 int i = 0; 00812 while(i < 9) 00813 { 00814 i += 1; 00815 a += 1.0f; 00816 b = numext::pow( a, -x ); 00817 s += b; 00818 if( numext::abs(b/s) < machep ) 00819 return true; 00820 } 00821 00822 //Return whether we are done 00823 return false; 00824 } 00825 }; 00826 00827 template <> 00828 struct zeta_impl_series<double> { 00829 EIGEN_DEVICE_FUNC 00830 static EIGEN_STRONG_INLINE bool run(double& a, double& b, double& s, const double x, const double machep) { 00831 int i = 0; 00832 while( (i < 9) || (a <= 9.0) ) 00833 { 00834 i += 1; 00835 a += 1.0; 00836 b = numext::pow( a, -x ); 00837 s += b; 00838 if( numext::abs(b/s) < machep ) 00839 return true; 00840 } 00841 00842 //Return whether we are done 00843 return false; 00844 } 00845 }; 00846 00847 template <typename Scalar> 00848 struct zeta_impl { 00849 EIGEN_DEVICE_FUNC 00850 static Scalar run(Scalar x, Scalar q) { 00851 /* zeta.c 00852 * 00853 * Riemann zeta function of two arguments 00854 * 00855 * 00856 * 00857 * SYNOPSIS: 00858 * 00859 * double x, q, y, zeta(); 00860 * 00861 * y = zeta( x, q ); 00862 * 00863 * 00864 * 00865 * DESCRIPTION: 00866 * 00867 * 00868 * 00869 * inf. 00870 * - -x 00871 * zeta(x,q) = > (k+q) 00872 * - 00873 * k=0 00874 * 00875 * where x > 1 and q is not a negative integer or zero. 00876 * The Euler-Maclaurin summation formula is used to obtain 00877 * the expansion 00878 * 00879 * n 00880 * - -x 00881 * zeta(x,q) = > (k+q) 00882 * - 00883 * k=1 00884 * 00885 * 1-x inf. B x(x+1)...(x+2j) 00886 * (n+q) 1 - 2j 00887 * + --------- - ------- + > -------------------- 00888 * x-1 x - x+2j+1 00889 * 2(n+q) j=1 (2j)! (n+q) 00890 * 00891 * where the B2j are Bernoulli numbers. Note that (see zetac.c) 00892 * zeta(x,1) = zetac(x) + 1. 00893 * 00894 * 00895 * 00896 * ACCURACY: 00897 * 00898 * Relative error for single precision: 00899 * arithmetic domain # trials peak rms 00900 * IEEE 0,25 10000 6.9e-7 1.0e-7 00901 * 00902 * Large arguments may produce underflow in powf(), in which 00903 * case the results are inaccurate. 00904 * 00905 * REFERENCE: 00906 * 00907 * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, 00908 * Series, and Products, p. 1073; Academic Press, 1980. 00909 * 00910 */ 00911 00912 int i; 00913 Scalar p, r, a, b, k, s, t, w; 00914 00915 const Scalar A[] = { 00916 Scalar(12.0), 00917 Scalar(-720.0), 00918 Scalar(30240.0), 00919 Scalar(-1209600.0), 00920 Scalar(47900160.0), 00921 Scalar(-1.8924375803183791606e9), /*1.307674368e12/691*/ 00922 Scalar(7.47242496e10), 00923 Scalar(-2.950130727918164224e12), /*1.067062284288e16/3617*/ 00924 Scalar(1.1646782814350067249e14), /*5.109094217170944e18/43867*/ 00925 Scalar(-4.5979787224074726105e15), /*8.028576626982912e20/174611*/ 00926 Scalar(1.8152105401943546773e17), /*1.5511210043330985984e23/854513*/ 00927 Scalar(-7.1661652561756670113e18) /*1.6938241367317436694528e27/236364091*/ 00928 }; 00929 00930 const Scalar maxnum = NumTraits<Scalar>::infinity(); 00931 const Scalar zero = 0.0, half = 0.5, one = 1.0; 00932 const Scalar machep = cephes_helper<Scalar>::machep(); 00933 const Scalar nan = NumTraits<Scalar>::quiet_NaN(); 00934 00935 if( x == one ) 00936 return maxnum; 00937 00938 if( x < one ) 00939 { 00940 return nan; 00941 } 00942 00943 if( q <= zero ) 00944 { 00945 if(q == numext::floor(q)) 00946 { 00947 return maxnum; 00948 } 00949 p = x; 00950 r = numext::floor(p); 00951 if (p != r) 00952 return nan; 00953 } 00954 00955 /* Permit negative q but continue sum until n+q > +9 . 00956 * This case should be handled by a reflection formula. 00957 * If q<0 and x is an integer, there is a relation to 00958 * the polygamma function. 00959 */ 00960 s = numext::pow( q, -x ); 00961 a = q; 00962 b = zero; 00963 // Run the summation in a helper function that is specific to the floating precision 00964 if (zeta_impl_series<Scalar>::run(a, b, s, x, machep)) { 00965 return s; 00966 } 00967 00968 w = a; 00969 s += b*w/(x-one); 00970 s -= half * b; 00971 a = one; 00972 k = zero; 00973 for( i=0; i<12; i++ ) 00974 { 00975 a *= x + k; 00976 b /= w; 00977 t = a*b/A[i]; 00978 s = s + t; 00979 t = numext::abs(t/s); 00980 if( t < machep ) { 00981 break; 00982 } 00983 k += one; 00984 a *= x + k; 00985 b /= w; 00986 k += one; 00987 } 00988 return s; 00989 } 00990 }; 00991 00992 /**************************************************************************** 00993 * Implementation of polygamma function, requires C++11/C99 * 00994 ****************************************************************************/ 00995 00996 template <typename Scalar> 00997 struct polygamma_retval { 00998 typedef Scalar type; 00999 }; 01000 01001 #if !EIGEN_HAS_C99_MATH 01002 01003 template <typename Scalar> 01004 struct polygamma_impl { 01005 EIGEN_DEVICE_FUNC 01006 static EIGEN_STRONG_INLINE Scalar run(Scalar n, Scalar x) { 01007 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 01008 THIS_TYPE_IS_NOT_SUPPORTED); 01009 return Scalar(0); 01010 } 01011 }; 01012 01013 #else 01014 01015 template <typename Scalar> 01016 struct polygamma_impl { 01017 EIGEN_DEVICE_FUNC 01018 static Scalar run(Scalar n, Scalar x) { 01019 Scalar zero = 0.0, one = 1.0; 01020 Scalar nplus = n + one; 01021 const Scalar nan = NumTraits<Scalar>::quiet_NaN(); 01022 01023 // Check that n is an integer 01024 if (numext::floor(n) != n) { 01025 return nan; 01026 } 01027 // Just return the digamma function for n = 1 01028 else if (n == zero) { 01029 return digamma_impl<Scalar>::run(x); 01030 } 01031 // Use the same implementation as scipy 01032 else { 01033 Scalar factorial = numext::exp(lgamma_impl<Scalar>::run(nplus)); 01034 return numext::pow(-one, nplus) * factorial * zeta_impl<Scalar>::run(nplus, x); 01035 } 01036 } 01037 }; 01038 01039 #endif // EIGEN_HAS_C99_MATH 01040 01041 /************************************************************************************************ 01042 * Implementation of betainc (incomplete beta integral), based on Cephes but requires C++11/C99 * 01043 ************************************************************************************************/ 01044 01045 template <typename Scalar> 01046 struct betainc_retval { 01047 typedef Scalar type; 01048 }; 01049 01050 #if !EIGEN_HAS_C99_MATH 01051 01052 template <typename Scalar> 01053 struct betainc_impl { 01054 EIGEN_DEVICE_FUNC 01055 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x) { 01056 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 01057 THIS_TYPE_IS_NOT_SUPPORTED); 01058 return Scalar(0); 01059 } 01060 }; 01061 01062 #else 01063 01064 template <typename Scalar> 01065 struct betainc_impl { 01066 EIGEN_DEVICE_FUNC 01067 static EIGEN_STRONG_INLINE Scalar run(Scalar, Scalar, Scalar) { 01068 /* betaincf.c 01069 * 01070 * Incomplete beta integral 01071 * 01072 * 01073 * SYNOPSIS: 01074 * 01075 * float a, b, x, y, betaincf(); 01076 * 01077 * y = betaincf( a, b, x ); 01078 * 01079 * 01080 * DESCRIPTION: 01081 * 01082 * Returns incomplete beta integral of the arguments, evaluated 01083 * from zero to x. The function is defined as 01084 * 01085 * x 01086 * - - 01087 * | (a+b) | | a-1 b-1 01088 * ----------- | t (1-t) dt. 01089 * - - | | 01090 * | (a) | (b) - 01091 * 0 01092 * 01093 * The domain of definition is 0 <= x <= 1. In this 01094 * implementation a and b are restricted to positive values. 01095 * The integral from x to 1 may be obtained by the symmetry 01096 * relation 01097 * 01098 * 1 - betainc( a, b, x ) = betainc( b, a, 1-x ). 01099 * 01100 * The integral is evaluated by a continued fraction expansion. 01101 * If a < 1, the function calls itself recursively after a 01102 * transformation to increase a to a+1. 01103 * 01104 * ACCURACY (float): 01105 * 01106 * Tested at random points (a,b,x) with a and b in the indicated 01107 * interval and x between 0 and 1. 01108 * 01109 * arithmetic domain # trials peak rms 01110 * Relative error: 01111 * IEEE 0,30 10000 3.7e-5 5.1e-6 01112 * IEEE 0,100 10000 1.7e-4 2.5e-5 01113 * The useful domain for relative error is limited by underflow 01114 * of the single precision exponential function. 01115 * Absolute error: 01116 * IEEE 0,30 100000 2.2e-5 9.6e-7 01117 * IEEE 0,100 10000 6.5e-5 3.7e-6 01118 * 01119 * Larger errors may occur for extreme ratios of a and b. 01120 * 01121 * ACCURACY (double): 01122 * arithmetic domain # trials peak rms 01123 * IEEE 0,5 10000 6.9e-15 4.5e-16 01124 * IEEE 0,85 250000 2.2e-13 1.7e-14 01125 * IEEE 0,1000 30000 5.3e-12 6.3e-13 01126 * IEEE 0,10000 250000 9.3e-11 7.1e-12 01127 * IEEE 0,100000 10000 8.7e-10 4.8e-11 01128 * Outputs smaller than the IEEE gradual underflow threshold 01129 * were excluded from these statistics. 01130 * 01131 * ERROR MESSAGES: 01132 * message condition value returned 01133 * incbet domain x<0, x>1 nan 01134 * incbet underflow nan 01135 */ 01136 01137 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 01138 THIS_TYPE_IS_NOT_SUPPORTED); 01139 return Scalar(0); 01140 } 01141 }; 01142 01143 /* Continued fraction expansion #1 for incomplete beta integral (small_branch = True) 01144 * Continued fraction expansion #2 for incomplete beta integral (small_branch = False) 01145 */ 01146 template <typename Scalar> 01147 struct incbeta_cfe { 01148 EIGEN_DEVICE_FUNC 01149 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x, bool small_branch) { 01150 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, float>::value || 01151 internal::is_same<Scalar, double>::value), 01152 THIS_TYPE_IS_NOT_SUPPORTED); 01153 const Scalar big = cephes_helper<Scalar>::big(); 01154 const Scalar machep = cephes_helper<Scalar>::machep(); 01155 const Scalar biginv = cephes_helper<Scalar>::biginv(); 01156 01157 const Scalar zero = 0; 01158 const Scalar one = 1; 01159 const Scalar two = 2; 01160 01161 Scalar xk, pk, pkm1, pkm2, qk, qkm1, qkm2; 01162 Scalar k1, k2, k3, k4, k5, k6, k7, k8, k26update; 01163 Scalar ans; 01164 int n; 01165 01166 const int num_iters = (internal::is_same<Scalar, float>::value) ? 100 : 300; 01167 const Scalar thresh = 01168 (internal::is_same<Scalar, float>::value) ? machep : Scalar(3) * machep; 01169 Scalar r = (internal::is_same<Scalar, float>::value) ? zero : one; 01170 01171 if (small_branch) { 01172 k1 = a; 01173 k2 = a + b; 01174 k3 = a; 01175 k4 = a + one; 01176 k5 = one; 01177 k6 = b - one; 01178 k7 = k4; 01179 k8 = a + two; 01180 k26update = one; 01181 } else { 01182 k1 = a; 01183 k2 = b - one; 01184 k3 = a; 01185 k4 = a + one; 01186 k5 = one; 01187 k6 = a + b; 01188 k7 = a + one; 01189 k8 = a + two; 01190 k26update = -one; 01191 x = x / (one - x); 01192 } 01193 01194 pkm2 = zero; 01195 qkm2 = one; 01196 pkm1 = one; 01197 qkm1 = one; 01198 ans = one; 01199 n = 0; 01200 01201 do { 01202 xk = -(x * k1 * k2) / (k3 * k4); 01203 pk = pkm1 + pkm2 * xk; 01204 qk = qkm1 + qkm2 * xk; 01205 pkm2 = pkm1; 01206 pkm1 = pk; 01207 qkm2 = qkm1; 01208 qkm1 = qk; 01209 01210 xk = (x * k5 * k6) / (k7 * k8); 01211 pk = pkm1 + pkm2 * xk; 01212 qk = qkm1 + qkm2 * xk; 01213 pkm2 = pkm1; 01214 pkm1 = pk; 01215 qkm2 = qkm1; 01216 qkm1 = qk; 01217 01218 if (qk != zero) { 01219 r = pk / qk; 01220 if (numext::abs(ans - r) < numext::abs(r) * thresh) { 01221 return r; 01222 } 01223 ans = r; 01224 } 01225 01226 k1 += one; 01227 k2 += k26update; 01228 k3 += two; 01229 k4 += two; 01230 k5 += one; 01231 k6 -= k26update; 01232 k7 += two; 01233 k8 += two; 01234 01235 if ((numext::abs(qk) + numext::abs(pk)) > big) { 01236 pkm2 *= biginv; 01237 pkm1 *= biginv; 01238 qkm2 *= biginv; 01239 qkm1 *= biginv; 01240 } 01241 if ((numext::abs(qk) < biginv) || (numext::abs(pk) < biginv)) { 01242 pkm2 *= big; 01243 pkm1 *= big; 01244 qkm2 *= big; 01245 qkm1 *= big; 01246 } 01247 } while (++n < num_iters); 01248 01249 return ans; 01250 } 01251 }; 01252 01253 /* Helper functions depending on the Scalar type */ 01254 template <typename Scalar> 01255 struct betainc_helper {}; 01256 01257 template <> 01258 struct betainc_helper<float> { 01259 /* Core implementation, assumes a large (> 1.0) */ 01260 EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE float incbsa(float aa, float bb, 01261 float xx) { 01262 float ans, a, b, t, x, onemx; 01263 bool reversed_a_b = false; 01264 01265 onemx = 1.0f - xx; 01266 01267 /* see if x is greater than the mean */ 01268 if (xx > (aa / (aa + bb))) { 01269 reversed_a_b = true; 01270 a = bb; 01271 b = aa; 01272 t = xx; 01273 x = onemx; 01274 } else { 01275 a = aa; 01276 b = bb; 01277 t = onemx; 01278 x = xx; 01279 } 01280 01281 /* Choose expansion for optimal convergence */ 01282 if (b > 10.0f) { 01283 if (numext::abs(b * x / a) < 0.3f) { 01284 t = betainc_helper<float>::incbps(a, b, x); 01285 if (reversed_a_b) t = 1.0f - t; 01286 return t; 01287 } 01288 } 01289 01290 ans = x * (a + b - 2.0f) / (a - 1.0f); 01291 if (ans < 1.0f) { 01292 ans = incbeta_cfe<float>::run(a, b, x, true /* small_branch */); 01293 t = b * numext::log(t); 01294 } else { 01295 ans = incbeta_cfe<float>::run(a, b, x, false /* small_branch */); 01296 t = (b - 1.0f) * numext::log(t); 01297 } 01298 01299 t += a * numext::log(x) + lgamma_impl<float>::run(a + b) - 01300 lgamma_impl<float>::run(a) - lgamma_impl<float>::run(b); 01301 t += numext::log(ans / a); 01302 t = numext::exp(t); 01303 01304 if (reversed_a_b) t = 1.0f - t; 01305 return t; 01306 } 01307 01308 EIGEN_DEVICE_FUNC 01309 static EIGEN_STRONG_INLINE float incbps(float a, float b, float x) { 01310 float t, u, y, s; 01311 const float machep = cephes_helper<float>::machep(); 01312 01313 y = a * numext::log(x) + (b - 1.0f) * numext::log1p(-x) - numext::log(a); 01314 y -= lgamma_impl<float>::run(a) + lgamma_impl<float>::run(b); 01315 y += lgamma_impl<float>::run(a + b); 01316 01317 t = x / (1.0f - x); 01318 s = 0.0f; 01319 u = 1.0f; 01320 do { 01321 b -= 1.0f; 01322 if (b == 0.0f) { 01323 break; 01324 } 01325 a += 1.0f; 01326 u *= t * b / a; 01327 s += u; 01328 } while (numext::abs(u) > machep); 01329 01330 return numext::exp(y) * (1.0f + s); 01331 } 01332 }; 01333 01334 template <> 01335 struct betainc_impl<float> { 01336 EIGEN_DEVICE_FUNC 01337 static float run(float a, float b, float x) { 01338 const float nan = NumTraits<float>::quiet_NaN(); 01339 float ans, t; 01340 01341 if (a <= 0.0f) return nan; 01342 if (b <= 0.0f) return nan; 01343 if ((x <= 0.0f) || (x >= 1.0f)) { 01344 if (x == 0.0f) return 0.0f; 01345 if (x == 1.0f) return 1.0f; 01346 // mtherr("betaincf", DOMAIN); 01347 return nan; 01348 } 01349 01350 /* transformation for small aa */ 01351 if (a <= 1.0f) { 01352 ans = betainc_helper<float>::incbsa(a + 1.0f, b, x); 01353 t = a * numext::log(x) + b * numext::log1p(-x) + 01354 lgamma_impl<float>::run(a + b) - lgamma_impl<float>::run(a + 1.0f) - 01355 lgamma_impl<float>::run(b); 01356 return (ans + numext::exp(t)); 01357 } else { 01358 return betainc_helper<float>::incbsa(a, b, x); 01359 } 01360 } 01361 }; 01362 01363 template <> 01364 struct betainc_helper<double> { 01365 EIGEN_DEVICE_FUNC 01366 static EIGEN_STRONG_INLINE double incbps(double a, double b, double x) { 01367 const double machep = cephes_helper<double>::machep(); 01368 01369 double s, t, u, v, n, t1, z, ai; 01370 01371 ai = 1.0 / a; 01372 u = (1.0 - b) * x; 01373 v = u / (a + 1.0); 01374 t1 = v; 01375 t = u; 01376 n = 2.0; 01377 s = 0.0; 01378 z = machep * ai; 01379 while (numext::abs(v) > z) { 01380 u = (n - b) * x / n; 01381 t *= u; 01382 v = t / (a + n); 01383 s += v; 01384 n += 1.0; 01385 } 01386 s += t1; 01387 s += ai; 01388 01389 u = a * numext::log(x); 01390 // TODO: gamma() is not directly implemented in Eigen. 01391 /* 01392 if ((a + b) < maxgam && numext::abs(u) < maxlog) { 01393 t = gamma(a + b) / (gamma(a) * gamma(b)); 01394 s = s * t * pow(x, a); 01395 } else { 01396 */ 01397 t = lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) - 01398 lgamma_impl<double>::run(b) + u + numext::log(s); 01399 return s = numext::exp(t); 01400 } 01401 }; 01402 01403 template <> 01404 struct betainc_impl<double> { 01405 EIGEN_DEVICE_FUNC 01406 static double run(double aa, double bb, double xx) { 01407 const double nan = NumTraits<double>::quiet_NaN(); 01408 const double machep = cephes_helper<double>::machep(); 01409 // const double maxgam = 171.624376956302725; 01410 01411 double a, b, t, x, xc, w, y; 01412 bool reversed_a_b = false; 01413 01414 if (aa <= 0.0 || bb <= 0.0) { 01415 return nan; // goto domerr; 01416 } 01417 01418 if ((xx <= 0.0) || (xx >= 1.0)) { 01419 if (xx == 0.0) return (0.0); 01420 if (xx == 1.0) return (1.0); 01421 // mtherr("incbet", DOMAIN); 01422 return nan; 01423 } 01424 01425 if ((bb * xx) <= 1.0 && xx <= 0.95) { 01426 return betainc_helper<double>::incbps(aa, bb, xx); 01427 } 01428 01429 w = 1.0 - xx; 01430 01431 /* Reverse a and b if x is greater than the mean. */ 01432 if (xx > (aa / (aa + bb))) { 01433 reversed_a_b = true; 01434 a = bb; 01435 b = aa; 01436 xc = xx; 01437 x = w; 01438 } else { 01439 a = aa; 01440 b = bb; 01441 xc = w; 01442 x = xx; 01443 } 01444 01445 if (reversed_a_b && (b * x) <= 1.0 && x <= 0.95) { 01446 t = betainc_helper<double>::incbps(a, b, x); 01447 if (t <= machep) { 01448 t = 1.0 - machep; 01449 } else { 01450 t = 1.0 - t; 01451 } 01452 return t; 01453 } 01454 01455 /* Choose expansion for better convergence. */ 01456 y = x * (a + b - 2.0) - (a - 1.0); 01457 if (y < 0.0) { 01458 w = incbeta_cfe<double>::run(a, b, x, true /* small_branch */); 01459 } else { 01460 w = incbeta_cfe<double>::run(a, b, x, false /* small_branch */) / xc; 01461 } 01462 01463 /* Multiply w by the factor 01464 a b _ _ _ 01465 x (1-x) | (a+b) / ( a | (a) | (b) ) . */ 01466 01467 y = a * numext::log(x); 01468 t = b * numext::log(xc); 01469 // TODO: gamma is not directly implemented in Eigen. 01470 /* 01471 if ((a + b) < maxgam && numext::abs(y) < maxlog && numext::abs(t) < maxlog) 01472 { 01473 t = pow(xc, b); 01474 t *= pow(x, a); 01475 t /= a; 01476 t *= w; 01477 t *= gamma(a + b) / (gamma(a) * gamma(b)); 01478 } else { 01479 */ 01480 /* Resort to logarithms. */ 01481 y += t + lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) - 01482 lgamma_impl<double>::run(b); 01483 y += numext::log(w / a); 01484 t = numext::exp(y); 01485 01486 /* } */ 01487 // done: 01488 01489 if (reversed_a_b) { 01490 if (t <= machep) { 01491 t = 1.0 - machep; 01492 } else { 01493 t = 1.0 - t; 01494 } 01495 } 01496 return t; 01497 } 01498 }; 01499 01500 #endif // EIGEN_HAS_C99_MATH 01501 01502 } // end namespace internal 01503 01504 namespace numext { 01505 01506 template <typename Scalar> 01507 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(lgamma, Scalar) 01508 lgamma(const Scalar& x) { 01509 return EIGEN_MATHFUNC_IMPL(lgamma, Scalar)::run(x); 01510 } 01511 01512 template <typename Scalar> 01513 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(digamma, Scalar) 01514 digamma(const Scalar& x) { 01515 return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x); 01516 } 01517 01518 template <typename Scalar> 01519 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(zeta, Scalar) 01520 zeta(const Scalar& x, const Scalar& q) { 01521 return EIGEN_MATHFUNC_IMPL(zeta, Scalar)::run(x, q); 01522 } 01523 01524 template <typename Scalar> 01525 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(polygamma, Scalar) 01526 polygamma(const Scalar& n, const Scalar& x) { 01527 return EIGEN_MATHFUNC_IMPL(polygamma, Scalar)::run(n, x); 01528 } 01529 01530 template <typename Scalar> 01531 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erf, Scalar) 01532 erf(const Scalar& x) { 01533 return EIGEN_MATHFUNC_IMPL(erf, Scalar)::run(x); 01534 } 01535 01536 template <typename Scalar> 01537 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar) 01538 erfc(const Scalar& x) { 01539 return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x); 01540 } 01541 01542 template <typename Scalar> 01543 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma, Scalar) 01544 igamma(const Scalar& a, const Scalar& x) { 01545 return EIGEN_MATHFUNC_IMPL(igamma, Scalar)::run(a, x); 01546 } 01547 01548 template <typename Scalar> 01549 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igammac, Scalar) 01550 igammac(const Scalar& a, const Scalar& x) { 01551 return EIGEN_MATHFUNC_IMPL(igammac, Scalar)::run(a, x); 01552 } 01553 01554 template <typename Scalar> 01555 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(betainc, Scalar) 01556 betainc(const Scalar& a, const Scalar& b, const Scalar& x) { 01557 return EIGEN_MATHFUNC_IMPL(betainc, Scalar)::run(a, b, x); 01558 } 01559 01560 } // end namespace numext 01561 01562 01563 } // end namespace Eigen 01564 01565 #endif // EIGEN_SPECIAL_FUNCTIONS_H