SpecialFunctionsImpl.h
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
 All Classes Functions Variables Typedefs Enumerator