1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2015 Eugene Brevdo <ebrevdo@gmail.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_SPECIAL_FUNCTIONS_H 11 #define EIGEN_SPECIAL_FUNCTIONS_H 12 13 namespace Eigen { 14 namespace internal { 15 16 // Parts of this code are based on the Cephes Math Library. 17 // 18 // Cephes Math Library Release 2.8: June, 2000 19 // Copyright 1984, 1987, 1992, 2000 by Stephen L. Moshier 20 // 21 // Permission has been kindly provided by the original author 22 // to incorporate the Cephes software into the Eigen codebase: 23 // 24 // From: Stephen Moshier 25 // To: Eugene Brevdo 26 // Subject: Re: Permission to wrap several cephes functions in Eigen 27 // 28 // Hello Eugene, 29 // 30 // Thank you for writing. 31 // 32 // If your licensing is similar to BSD, the formal way that has been 33 // handled is simply to add a statement to the effect that you are incorporating 34 // the Cephes software by permission of the author. 35 // 36 // Good luck with your project, 37 // Steve 38 39 namespace cephes { 40 41 /* polevl (modified for Eigen) 42 * 43 * Evaluate polynomial 44 * 45 * 46 * 47 * SYNOPSIS: 48 * 49 * int N; 50 * Scalar x, y, coef[N+1]; 51 * 52 * y = polevl<decltype(x), N>( x, coef); 53 * 54 * 55 * 56 * DESCRIPTION: 57 * 58 * Evaluates polynomial of degree N: 59 * 60 * 2 N 61 * y = C + C x + C x +...+ C x 62 * 0 1 2 N 63 * 64 * Coefficients are stored in reverse order: 65 * 66 * coef[0] = C , ..., coef[N] = C . 67 * N 0 68 * 69 * The function p1evl() assumes that coef[N] = 1.0 and is 70 * omitted from the array. Its calling arguments are 71 * otherwise the same as polevl(). 72 * 73 * 74 * The Eigen implementation is templatized. For best speed, store 75 * coef as a const array (constexpr), e.g. 76 * 77 * const double coef[] = {1.0, 2.0, 3.0, ...}; 78 * 79 */ 80 template <typename Scalar, int N> 81 struct polevl { 82 EIGEN_DEVICE_FUNC runpolevl83 static EIGEN_STRONG_INLINE Scalar run(const Scalar x, const Scalar coef[]) { 84 EIGEN_STATIC_ASSERT((N > 0), YOU_MADE_A_PROGRAMMING_MISTAKE); 85 86 return polevl<Scalar, N - 1>::run(x, coef) * x + coef[N]; 87 } 88 }; 89 90 template <typename Scalar> 91 struct polevl<Scalar, 0> { 92 EIGEN_DEVICE_FUNC 93 static EIGEN_STRONG_INLINE Scalar run(const Scalar, const Scalar coef[]) { 94 return coef[0]; 95 } 96 }; 97 98 } // end namespace cephes 99 100 /**************************************************************************** 101 * Implementation of lgamma, requires C++11/C99 * 102 ****************************************************************************/ 103 104 template <typename Scalar> 105 struct lgamma_impl { 106 EIGEN_DEVICE_FUNC 107 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 108 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 109 THIS_TYPE_IS_NOT_SUPPORTED); 110 return Scalar(0); 111 } 112 }; 113 114 template <typename Scalar> 115 struct lgamma_retval { 116 typedef Scalar type; 117 }; 118 119 #if EIGEN_HAS_C99_MATH 120 template <> 121 struct lgamma_impl<float> { 122 EIGEN_DEVICE_FUNC 123 static EIGEN_STRONG_INLINE float run(float x) { 124 #if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__) 125 int signgam; 126 return ::lgammaf_r(x, &signgam); 127 #else 128 return ::lgammaf(x); 129 #endif 130 } 131 }; 132 133 template <> 134 struct lgamma_impl<double> { 135 EIGEN_DEVICE_FUNC 136 static EIGEN_STRONG_INLINE double run(double x) { 137 #if !defined(__CUDA_ARCH__) && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE)) && !defined(__APPLE__) 138 int signgam; 139 return ::lgamma_r(x, &signgam); 140 #else 141 return ::lgamma(x); 142 #endif 143 } 144 }; 145 #endif 146 147 /**************************************************************************** 148 * Implementation of digamma (psi), based on Cephes * 149 ****************************************************************************/ 150 151 template <typename Scalar> 152 struct digamma_retval { 153 typedef Scalar type; 154 }; 155 156 /* 157 * 158 * Polynomial evaluation helper for the Psi (digamma) function. 159 * 160 * digamma_impl_maybe_poly::run(s) evaluates the asymptotic Psi expansion for 161 * input Scalar s, assuming s is above 10.0. 162 * 163 * If s is above a certain threshold for the given Scalar type, zero 164 * is returned. Otherwise the polynomial is evaluated with enough 165 * coefficients for results matching Scalar machine precision. 166 * 167 * 168 */ 169 template <typename Scalar> 170 struct digamma_impl_maybe_poly { 171 EIGEN_DEVICE_FUNC 172 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 173 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 174 THIS_TYPE_IS_NOT_SUPPORTED); 175 return Scalar(0); 176 } 177 }; 178 179 180 template <> 181 struct digamma_impl_maybe_poly<float> { 182 EIGEN_DEVICE_FUNC 183 static EIGEN_STRONG_INLINE float run(const float s) { 184 const float A[] = { 185 -4.16666666666666666667E-3f, 186 3.96825396825396825397E-3f, 187 -8.33333333333333333333E-3f, 188 8.33333333333333333333E-2f 189 }; 190 191 float z; 192 if (s < 1.0e8f) { 193 z = 1.0f / (s * s); 194 return z * cephes::polevl<float, 3>::run(z, A); 195 } else return 0.0f; 196 } 197 }; 198 199 template <> 200 struct digamma_impl_maybe_poly<double> { 201 EIGEN_DEVICE_FUNC 202 static EIGEN_STRONG_INLINE double run(const double s) { 203 const double A[] = { 204 8.33333333333333333333E-2, 205 -2.10927960927960927961E-2, 206 7.57575757575757575758E-3, 207 -4.16666666666666666667E-3, 208 3.96825396825396825397E-3, 209 -8.33333333333333333333E-3, 210 8.33333333333333333333E-2 211 }; 212 213 double z; 214 if (s < 1.0e17) { 215 z = 1.0 / (s * s); 216 return z * cephes::polevl<double, 6>::run(z, A); 217 } 218 else return 0.0; 219 } 220 }; 221 222 template <typename Scalar> 223 struct digamma_impl { 224 EIGEN_DEVICE_FUNC 225 static Scalar run(Scalar x) { 226 /* 227 * 228 * Psi (digamma) function (modified for Eigen) 229 * 230 * 231 * SYNOPSIS: 232 * 233 * double x, y, psi(); 234 * 235 * y = psi( x ); 236 * 237 * 238 * DESCRIPTION: 239 * 240 * d - 241 * psi(x) = -- ln | (x) 242 * dx 243 * 244 * is the logarithmic derivative of the gamma function. 245 * For integer x, 246 * n-1 247 * - 248 * psi(n) = -EUL + > 1/k. 249 * - 250 * k=1 251 * 252 * If x is negative, it is transformed to a positive argument by the 253 * reflection formula psi(1-x) = psi(x) + pi cot(pi x). 254 * For general positive x, the argument is made greater than 10 255 * using the recurrence psi(x+1) = psi(x) + 1/x. 256 * Then the following asymptotic expansion is applied: 257 * 258 * inf. B 259 * - 2k 260 * psi(x) = log(x) - 1/2x - > ------- 261 * - 2k 262 * k=1 2k x 263 * 264 * where the B2k are Bernoulli numbers. 265 * 266 * ACCURACY (float): 267 * Relative error (except absolute when |psi| < 1): 268 * arithmetic domain # trials peak rms 269 * IEEE 0,30 30000 1.3e-15 1.4e-16 270 * IEEE -30,0 40000 1.5e-15 2.2e-16 271 * 272 * ACCURACY (double): 273 * Absolute error, relative when |psi| > 1 : 274 * arithmetic domain # trials peak rms 275 * IEEE -33,0 30000 8.2e-7 1.2e-7 276 * IEEE 0,33 100000 7.3e-7 7.7e-8 277 * 278 * ERROR MESSAGES: 279 * message condition value returned 280 * psi singularity x integer <=0 INFINITY 281 */ 282 283 Scalar p, q, nz, s, w, y; 284 bool negative = false; 285 286 const Scalar maxnum = NumTraits<Scalar>::infinity(); 287 const Scalar m_pi = Scalar(EIGEN_PI); 288 289 const Scalar zero = Scalar(0); 290 const Scalar one = Scalar(1); 291 const Scalar half = Scalar(0.5); 292 nz = zero; 293 294 if (x <= zero) { 295 negative = true; 296 q = x; 297 p = numext::floor(q); 298 if (p == q) { 299 return maxnum; 300 } 301 /* Remove the zeros of tan(m_pi x) 302 * by subtracting the nearest integer from x 303 */ 304 nz = q - p; 305 if (nz != half) { 306 if (nz > half) { 307 p += one; 308 nz = q - p; 309 } 310 nz = m_pi / numext::tan(m_pi * nz); 311 } 312 else { 313 nz = zero; 314 } 315 x = one - x; 316 } 317 318 /* use the recurrence psi(x+1) = psi(x) + 1/x. */ 319 s = x; 320 w = zero; 321 while (s < Scalar(10)) { 322 w += one / s; 323 s += one; 324 } 325 326 y = digamma_impl_maybe_poly<Scalar>::run(s); 327 328 y = numext::log(s) - (half / s) - y - w; 329 330 return (negative) ? y - nz : y; 331 } 332 }; 333 334 /**************************************************************************** 335 * Implementation of erf, requires C++11/C99 * 336 ****************************************************************************/ 337 338 template <typename Scalar> 339 struct erf_impl { 340 EIGEN_DEVICE_FUNC 341 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 342 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 343 THIS_TYPE_IS_NOT_SUPPORTED); 344 return Scalar(0); 345 } 346 }; 347 348 template <typename Scalar> 349 struct erf_retval { 350 typedef Scalar type; 351 }; 352 353 #if EIGEN_HAS_C99_MATH 354 template <> 355 struct erf_impl<float> { 356 EIGEN_DEVICE_FUNC 357 static EIGEN_STRONG_INLINE float run(float x) { return ::erff(x); } 358 }; 359 360 template <> 361 struct erf_impl<double> { 362 EIGEN_DEVICE_FUNC 363 static EIGEN_STRONG_INLINE double run(double x) { return ::erf(x); } 364 }; 365 #endif // EIGEN_HAS_C99_MATH 366 367 /*************************************************************************** 368 * Implementation of erfc, requires C++11/C99 * 369 ****************************************************************************/ 370 371 template <typename Scalar> 372 struct erfc_impl { 373 EIGEN_DEVICE_FUNC 374 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 375 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 376 THIS_TYPE_IS_NOT_SUPPORTED); 377 return Scalar(0); 378 } 379 }; 380 381 template <typename Scalar> 382 struct erfc_retval { 383 typedef Scalar type; 384 }; 385 386 #if EIGEN_HAS_C99_MATH 387 template <> 388 struct erfc_impl<float> { 389 EIGEN_DEVICE_FUNC 390 static EIGEN_STRONG_INLINE float run(const float x) { return ::erfcf(x); } 391 }; 392 393 template <> 394 struct erfc_impl<double> { 395 EIGEN_DEVICE_FUNC 396 static EIGEN_STRONG_INLINE double run(const double x) { return ::erfc(x); } 397 }; 398 #endif // EIGEN_HAS_C99_MATH 399 400 /************************************************************************************************************** 401 * Implementation of igammac (complemented incomplete gamma integral), based on Cephes but requires C++11/C99 * 402 **************************************************************************************************************/ 403 404 template <typename Scalar> 405 struct igammac_retval { 406 typedef Scalar type; 407 }; 408 409 // NOTE: cephes_helper is also used to implement zeta 410 template <typename Scalar> 411 struct cephes_helper { 412 EIGEN_DEVICE_FUNC 413 static EIGEN_STRONG_INLINE Scalar machep() { assert(false && "machep not supported for this type"); return 0.0; } 414 EIGEN_DEVICE_FUNC 415 static EIGEN_STRONG_INLINE Scalar big() { assert(false && "big not supported for this type"); return 0.0; } 416 EIGEN_DEVICE_FUNC 417 static EIGEN_STRONG_INLINE Scalar biginv() { assert(false && "biginv not supported for this type"); return 0.0; } 418 }; 419 420 template <> 421 struct cephes_helper<float> { 422 EIGEN_DEVICE_FUNC 423 static EIGEN_STRONG_INLINE float machep() { 424 return NumTraits<float>::epsilon() / 2; // 1.0 - machep == 1.0 425 } 426 EIGEN_DEVICE_FUNC 427 static EIGEN_STRONG_INLINE float big() { 428 // use epsneg (1.0 - epsneg == 1.0) 429 return 1.0f / (NumTraits<float>::epsilon() / 2); 430 } 431 EIGEN_DEVICE_FUNC 432 static EIGEN_STRONG_INLINE float biginv() { 433 // epsneg 434 return machep(); 435 } 436 }; 437 438 template <> 439 struct cephes_helper<double> { 440 EIGEN_DEVICE_FUNC 441 static EIGEN_STRONG_INLINE double machep() { 442 return NumTraits<double>::epsilon() / 2; // 1.0 - machep == 1.0 443 } 444 EIGEN_DEVICE_FUNC 445 static EIGEN_STRONG_INLINE double big() { 446 return 1.0 / NumTraits<double>::epsilon(); 447 } 448 EIGEN_DEVICE_FUNC 449 static EIGEN_STRONG_INLINE double biginv() { 450 // inverse of eps 451 return NumTraits<double>::epsilon(); 452 } 453 }; 454 455 #if !EIGEN_HAS_C99_MATH 456 457 template <typename Scalar> 458 struct igammac_impl { 459 EIGEN_DEVICE_FUNC 460 static Scalar run(Scalar a, Scalar x) { 461 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 462 THIS_TYPE_IS_NOT_SUPPORTED); 463 return Scalar(0); 464 } 465 }; 466 467 #else 468 469 template <typename Scalar> struct igamma_impl; // predeclare igamma_impl 470 471 template <typename Scalar> 472 struct igammac_impl { 473 EIGEN_DEVICE_FUNC 474 static Scalar run(Scalar a, Scalar x) { 475 /* igamc() 476 * 477 * Incomplete gamma integral (modified for Eigen) 478 * 479 * 480 * 481 * SYNOPSIS: 482 * 483 * double a, x, y, igamc(); 484 * 485 * y = igamc( a, x ); 486 * 487 * DESCRIPTION: 488 * 489 * The function is defined by 490 * 491 * 492 * igamc(a,x) = 1 - igam(a,x) 493 * 494 * inf. 495 * - 496 * 1 | | -t a-1 497 * = ----- | e t dt. 498 * - | | 499 * | (a) - 500 * x 501 * 502 * 503 * In this implementation both arguments must be positive. 504 * The integral is evaluated by either a power series or 505 * continued fraction expansion, depending on the relative 506 * values of a and x. 507 * 508 * ACCURACY (float): 509 * 510 * Relative error: 511 * arithmetic domain # trials peak rms 512 * IEEE 0,30 30000 7.8e-6 5.9e-7 513 * 514 * 515 * ACCURACY (double): 516 * 517 * Tested at random a, x. 518 * a x Relative error: 519 * arithmetic domain domain # trials peak rms 520 * IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15 521 * IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15 522 * 523 */ 524 /* 525 Cephes Math Library Release 2.2: June, 1992 526 Copyright 1985, 1987, 1992 by Stephen L. Moshier 527 Direct inquiries to 30 Frost Street, Cambridge, MA 02140 528 */ 529 const Scalar zero = 0; 530 const Scalar one = 1; 531 const Scalar nan = NumTraits<Scalar>::quiet_NaN(); 532 533 if ((x < zero) || (a <= zero)) { 534 // domain error 535 return nan; 536 } 537 538 if ((x < one) || (x < a)) { 539 /* The checks above ensure that we meet the preconditions for 540 * igamma_impl::Impl(), so call it, rather than igamma_impl::Run(). 541 * Calling Run() would also work, but in that case the compiler may not be 542 * able to prove that igammac_impl::Run and igamma_impl::Run are not 543 * mutually recursive. This leads to worse code, particularly on 544 * platforms like nvptx, where recursion is allowed only begrudgingly. 545 */ 546 return (one - igamma_impl<Scalar>::Impl(a, x)); 547 } 548 549 return Impl(a, x); 550 } 551 552 private: 553 /* igamma_impl calls igammac_impl::Impl. */ 554 friend struct igamma_impl<Scalar>; 555 556 /* Actually computes igamc(a, x). 557 * 558 * Preconditions: 559 * a > 0 560 * x >= 1 561 * x >= a 562 */ 563 EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) { 564 const Scalar zero = 0; 565 const Scalar one = 1; 566 const Scalar two = 2; 567 const Scalar machep = cephes_helper<Scalar>::machep(); 568 const Scalar maxlog = numext::log(NumTraits<Scalar>::highest()); 569 const Scalar big = cephes_helper<Scalar>::big(); 570 const Scalar biginv = cephes_helper<Scalar>::biginv(); 571 const Scalar inf = NumTraits<Scalar>::infinity(); 572 573 Scalar ans, ax, c, yc, r, t, y, z; 574 Scalar pk, pkm1, pkm2, qk, qkm1, qkm2; 575 576 if (x == inf) return zero; // std::isinf crashes on CUDA 577 578 /* Compute x**a * exp(-x) / gamma(a) */ 579 ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); 580 if (ax < -maxlog) { // underflow 581 return zero; 582 } 583 ax = numext::exp(ax); 584 585 // continued fraction 586 y = one - a; 587 z = x + y + one; 588 c = zero; 589 pkm2 = one; 590 qkm2 = x; 591 pkm1 = x + one; 592 qkm1 = z * x; 593 ans = pkm1 / qkm1; 594 595 while (true) { 596 c += one; 597 y += one; 598 z += two; 599 yc = y * c; 600 pk = pkm1 * z - pkm2 * yc; 601 qk = qkm1 * z - qkm2 * yc; 602 if (qk != zero) { 603 r = pk / qk; 604 t = numext::abs((ans - r) / r); 605 ans = r; 606 } else { 607 t = one; 608 } 609 pkm2 = pkm1; 610 pkm1 = pk; 611 qkm2 = qkm1; 612 qkm1 = qk; 613 if (numext::abs(pk) > big) { 614 pkm2 *= biginv; 615 pkm1 *= biginv; 616 qkm2 *= biginv; 617 qkm1 *= biginv; 618 } 619 if (t <= machep) { 620 break; 621 } 622 } 623 624 return (ans * ax); 625 } 626 }; 627 628 #endif // EIGEN_HAS_C99_MATH 629 630 /************************************************************************************************ 631 * Implementation of igamma (incomplete gamma integral), based on Cephes but requires C++11/C99 * 632 ************************************************************************************************/ 633 634 template <typename Scalar> 635 struct igamma_retval { 636 typedef Scalar type; 637 }; 638 639 #if !EIGEN_HAS_C99_MATH 640 641 template <typename Scalar> 642 struct igamma_impl { 643 EIGEN_DEVICE_FUNC 644 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x) { 645 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 646 THIS_TYPE_IS_NOT_SUPPORTED); 647 return Scalar(0); 648 } 649 }; 650 651 #else 652 653 template <typename Scalar> 654 struct igamma_impl { 655 EIGEN_DEVICE_FUNC 656 static Scalar run(Scalar a, Scalar x) { 657 /* igam() 658 * Incomplete gamma integral 659 * 660 * 661 * 662 * SYNOPSIS: 663 * 664 * double a, x, y, igam(); 665 * 666 * y = igam( a, x ); 667 * 668 * DESCRIPTION: 669 * 670 * The function is defined by 671 * 672 * x 673 * - 674 * 1 | | -t a-1 675 * igam(a,x) = ----- | e t dt. 676 * - | | 677 * | (a) - 678 * 0 679 * 680 * 681 * In this implementation both arguments must be positive. 682 * The integral is evaluated by either a power series or 683 * continued fraction expansion, depending on the relative 684 * values of a and x. 685 * 686 * ACCURACY (double): 687 * 688 * Relative error: 689 * arithmetic domain # trials peak rms 690 * IEEE 0,30 200000 3.6e-14 2.9e-15 691 * IEEE 0,100 300000 9.9e-14 1.5e-14 692 * 693 * 694 * ACCURACY (float): 695 * 696 * Relative error: 697 * arithmetic domain # trials peak rms 698 * IEEE 0,30 20000 7.8e-6 5.9e-7 699 * 700 */ 701 /* 702 Cephes Math Library Release 2.2: June, 1992 703 Copyright 1985, 1987, 1992 by Stephen L. Moshier 704 Direct inquiries to 30 Frost Street, Cambridge, MA 02140 705 */ 706 707 708 /* left tail of incomplete gamma function: 709 * 710 * inf. k 711 * a -x - x 712 * x e > ---------- 713 * - - 714 * k=0 | (a+k+1) 715 * 716 */ 717 const Scalar zero = 0; 718 const Scalar one = 1; 719 const Scalar nan = NumTraits<Scalar>::quiet_NaN(); 720 721 if (x == zero) return zero; 722 723 if ((x < zero) || (a <= zero)) { // domain error 724 return nan; 725 } 726 727 if ((x > one) && (x > a)) { 728 /* The checks above ensure that we meet the preconditions for 729 * igammac_impl::Impl(), so call it, rather than igammac_impl::Run(). 730 * Calling Run() would also work, but in that case the compiler may not be 731 * able to prove that igammac_impl::Run and igamma_impl::Run are not 732 * mutually recursive. This leads to worse code, particularly on 733 * platforms like nvptx, where recursion is allowed only begrudgingly. 734 */ 735 return (one - igammac_impl<Scalar>::Impl(a, x)); 736 } 737 738 return Impl(a, x); 739 } 740 741 private: 742 /* igammac_impl calls igamma_impl::Impl. */ 743 friend struct igammac_impl<Scalar>; 744 745 /* Actually computes igam(a, x). 746 * 747 * Preconditions: 748 * x > 0 749 * a > 0 750 * !(x > 1 && x > a) 751 */ 752 EIGEN_DEVICE_FUNC static Scalar Impl(Scalar a, Scalar x) { 753 const Scalar zero = 0; 754 const Scalar one = 1; 755 const Scalar machep = cephes_helper<Scalar>::machep(); 756 const Scalar maxlog = numext::log(NumTraits<Scalar>::highest()); 757 758 Scalar ans, ax, c, r; 759 760 /* Compute x**a * exp(-x) / gamma(a) */ 761 ax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a); 762 if (ax < -maxlog) { 763 // underflow 764 return zero; 765 } 766 ax = numext::exp(ax); 767 768 /* power series */ 769 r = a; 770 c = one; 771 ans = one; 772 773 while (true) { 774 r += one; 775 c *= x/r; 776 ans += c; 777 if (c/ans <= machep) { 778 break; 779 } 780 } 781 782 return (ans * ax / a); 783 } 784 }; 785 786 #endif // EIGEN_HAS_C99_MATH 787 788 /***************************************************************************** 789 * Implementation of Riemann zeta function of two arguments, based on Cephes * 790 *****************************************************************************/ 791 792 template <typename Scalar> 793 struct zeta_retval { 794 typedef Scalar type; 795 }; 796 797 template <typename Scalar> 798 struct zeta_impl_series { 799 EIGEN_DEVICE_FUNC 800 static EIGEN_STRONG_INLINE Scalar run(const Scalar) { 801 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 802 THIS_TYPE_IS_NOT_SUPPORTED); 803 return Scalar(0); 804 } 805 }; 806 807 template <> 808 struct zeta_impl_series<float> { 809 EIGEN_DEVICE_FUNC 810 static EIGEN_STRONG_INLINE bool run(float& a, float& b, float& s, const float x, const float machep) { 811 int i = 0; 812 while(i < 9) 813 { 814 i += 1; 815 a += 1.0f; 816 b = numext::pow( a, -x ); 817 s += b; 818 if( numext::abs(b/s) < machep ) 819 return true; 820 } 821 822 //Return whether we are done 823 return false; 824 } 825 }; 826 827 template <> 828 struct zeta_impl_series<double> { 829 EIGEN_DEVICE_FUNC 830 static EIGEN_STRONG_INLINE bool run(double& a, double& b, double& s, const double x, const double machep) { 831 int i = 0; 832 while( (i < 9) || (a <= 9.0) ) 833 { 834 i += 1; 835 a += 1.0; 836 b = numext::pow( a, -x ); 837 s += b; 838 if( numext::abs(b/s) < machep ) 839 return true; 840 } 841 842 //Return whether we are done 843 return false; 844 } 845 }; 846 847 template <typename Scalar> 848 struct zeta_impl { 849 EIGEN_DEVICE_FUNC 850 static Scalar run(Scalar x, Scalar q) { 851 /* zeta.c 852 * 853 * Riemann zeta function of two arguments 854 * 855 * 856 * 857 * SYNOPSIS: 858 * 859 * double x, q, y, zeta(); 860 * 861 * y = zeta( x, q ); 862 * 863 * 864 * 865 * DESCRIPTION: 866 * 867 * 868 * 869 * inf. 870 * - -x 871 * zeta(x,q) = > (k+q) 872 * - 873 * k=0 874 * 875 * where x > 1 and q is not a negative integer or zero. 876 * The Euler-Maclaurin summation formula is used to obtain 877 * the expansion 878 * 879 * n 880 * - -x 881 * zeta(x,q) = > (k+q) 882 * - 883 * k=1 884 * 885 * 1-x inf. B x(x+1)...(x+2j) 886 * (n+q) 1 - 2j 887 * + --------- - ------- + > -------------------- 888 * x-1 x - x+2j+1 889 * 2(n+q) j=1 (2j)! (n+q) 890 * 891 * where the B2j are Bernoulli numbers. Note that (see zetac.c) 892 * zeta(x,1) = zetac(x) + 1. 893 * 894 * 895 * 896 * ACCURACY: 897 * 898 * Relative error for single precision: 899 * arithmetic domain # trials peak rms 900 * IEEE 0,25 10000 6.9e-7 1.0e-7 901 * 902 * Large arguments may produce underflow in powf(), in which 903 * case the results are inaccurate. 904 * 905 * REFERENCE: 906 * 907 * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, 908 * Series, and Products, p. 1073; Academic Press, 1980. 909 * 910 */ 911 912 int i; 913 Scalar p, r, a, b, k, s, t, w; 914 915 const Scalar A[] = { 916 Scalar(12.0), 917 Scalar(-720.0), 918 Scalar(30240.0), 919 Scalar(-1209600.0), 920 Scalar(47900160.0), 921 Scalar(-1.8924375803183791606e9), /*1.307674368e12/691*/ 922 Scalar(7.47242496e10), 923 Scalar(-2.950130727918164224e12), /*1.067062284288e16/3617*/ 924 Scalar(1.1646782814350067249e14), /*5.109094217170944e18/43867*/ 925 Scalar(-4.5979787224074726105e15), /*8.028576626982912e20/174611*/ 926 Scalar(1.8152105401943546773e17), /*1.5511210043330985984e23/854513*/ 927 Scalar(-7.1661652561756670113e18) /*1.6938241367317436694528e27/236364091*/ 928 }; 929 930 const Scalar maxnum = NumTraits<Scalar>::infinity(); 931 const Scalar zero = 0.0, half = 0.5, one = 1.0; 932 const Scalar machep = cephes_helper<Scalar>::machep(); 933 const Scalar nan = NumTraits<Scalar>::quiet_NaN(); 934 935 if( x == one ) 936 return maxnum; 937 938 if( x < one ) 939 { 940 return nan; 941 } 942 943 if( q <= zero ) 944 { 945 if(q == numext::floor(q)) 946 { 947 return maxnum; 948 } 949 p = x; 950 r = numext::floor(p); 951 if (p != r) 952 return nan; 953 } 954 955 /* Permit negative q but continue sum until n+q > +9 . 956 * This case should be handled by a reflection formula. 957 * If q<0 and x is an integer, there is a relation to 958 * the polygamma function. 959 */ 960 s = numext::pow( q, -x ); 961 a = q; 962 b = zero; 963 // Run the summation in a helper function that is specific to the floating precision 964 if (zeta_impl_series<Scalar>::run(a, b, s, x, machep)) { 965 return s; 966 } 967 968 w = a; 969 s += b*w/(x-one); 970 s -= half * b; 971 a = one; 972 k = zero; 973 for( i=0; i<12; i++ ) 974 { 975 a *= x + k; 976 b /= w; 977 t = a*b/A[i]; 978 s = s + t; 979 t = numext::abs(t/s); 980 if( t < machep ) { 981 break; 982 } 983 k += one; 984 a *= x + k; 985 b /= w; 986 k += one; 987 } 988 return s; 989 } 990 }; 991 992 /**************************************************************************** 993 * Implementation of polygamma function, requires C++11/C99 * 994 ****************************************************************************/ 995 996 template <typename Scalar> 997 struct polygamma_retval { 998 typedef Scalar type; 999 }; 1000 1001 #if !EIGEN_HAS_C99_MATH 1002 1003 template <typename Scalar> 1004 struct polygamma_impl { 1005 EIGEN_DEVICE_FUNC 1006 static EIGEN_STRONG_INLINE Scalar run(Scalar n, Scalar x) { 1007 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 1008 THIS_TYPE_IS_NOT_SUPPORTED); 1009 return Scalar(0); 1010 } 1011 }; 1012 1013 #else 1014 1015 template <typename Scalar> 1016 struct polygamma_impl { 1017 EIGEN_DEVICE_FUNC 1018 static Scalar run(Scalar n, Scalar x) { 1019 Scalar zero = 0.0, one = 1.0; 1020 Scalar nplus = n + one; 1021 const Scalar nan = NumTraits<Scalar>::quiet_NaN(); 1022 1023 // Check that n is an integer 1024 if (numext::floor(n) != n) { 1025 return nan; 1026 } 1027 // Just return the digamma function for n = 1 1028 else if (n == zero) { 1029 return digamma_impl<Scalar>::run(x); 1030 } 1031 // Use the same implementation as scipy 1032 else { 1033 Scalar factorial = numext::exp(lgamma_impl<Scalar>::run(nplus)); 1034 return numext::pow(-one, nplus) * factorial * zeta_impl<Scalar>::run(nplus, x); 1035 } 1036 } 1037 }; 1038 1039 #endif // EIGEN_HAS_C99_MATH 1040 1041 /************************************************************************************************ 1042 * Implementation of betainc (incomplete beta integral), based on Cephes but requires C++11/C99 * 1043 ************************************************************************************************/ 1044 1045 template <typename Scalar> 1046 struct betainc_retval { 1047 typedef Scalar type; 1048 }; 1049 1050 #if !EIGEN_HAS_C99_MATH 1051 1052 template <typename Scalar> 1053 struct betainc_impl { 1054 EIGEN_DEVICE_FUNC 1055 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x) { 1056 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 1057 THIS_TYPE_IS_NOT_SUPPORTED); 1058 return Scalar(0); 1059 } 1060 }; 1061 1062 #else 1063 1064 template <typename Scalar> 1065 struct betainc_impl { 1066 EIGEN_DEVICE_FUNC 1067 static EIGEN_STRONG_INLINE Scalar run(Scalar, Scalar, Scalar) { 1068 /* betaincf.c 1069 * 1070 * Incomplete beta integral 1071 * 1072 * 1073 * SYNOPSIS: 1074 * 1075 * float a, b, x, y, betaincf(); 1076 * 1077 * y = betaincf( a, b, x ); 1078 * 1079 * 1080 * DESCRIPTION: 1081 * 1082 * Returns incomplete beta integral of the arguments, evaluated 1083 * from zero to x. The function is defined as 1084 * 1085 * x 1086 * - - 1087 * | (a+b) | | a-1 b-1 1088 * ----------- | t (1-t) dt. 1089 * - - | | 1090 * | (a) | (b) - 1091 * 0 1092 * 1093 * The domain of definition is 0 <= x <= 1. In this 1094 * implementation a and b are restricted to positive values. 1095 * The integral from x to 1 may be obtained by the symmetry 1096 * relation 1097 * 1098 * 1 - betainc( a, b, x ) = betainc( b, a, 1-x ). 1099 * 1100 * The integral is evaluated by a continued fraction expansion. 1101 * If a < 1, the function calls itself recursively after a 1102 * transformation to increase a to a+1. 1103 * 1104 * ACCURACY (float): 1105 * 1106 * Tested at random points (a,b,x) with a and b in the indicated 1107 * interval and x between 0 and 1. 1108 * 1109 * arithmetic domain # trials peak rms 1110 * Relative error: 1111 * IEEE 0,30 10000 3.7e-5 5.1e-6 1112 * IEEE 0,100 10000 1.7e-4 2.5e-5 1113 * The useful domain for relative error is limited by underflow 1114 * of the single precision exponential function. 1115 * Absolute error: 1116 * IEEE 0,30 100000 2.2e-5 9.6e-7 1117 * IEEE 0,100 10000 6.5e-5 3.7e-6 1118 * 1119 * Larger errors may occur for extreme ratios of a and b. 1120 * 1121 * ACCURACY (double): 1122 * arithmetic domain # trials peak rms 1123 * IEEE 0,5 10000 6.9e-15 4.5e-16 1124 * IEEE 0,85 250000 2.2e-13 1.7e-14 1125 * IEEE 0,1000 30000 5.3e-12 6.3e-13 1126 * IEEE 0,10000 250000 9.3e-11 7.1e-12 1127 * IEEE 0,100000 10000 8.7e-10 4.8e-11 1128 * Outputs smaller than the IEEE gradual underflow threshold 1129 * were excluded from these statistics. 1130 * 1131 * ERROR MESSAGES: 1132 * message condition value returned 1133 * incbet domain x<0, x>1 nan 1134 * incbet underflow nan 1135 */ 1136 1137 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, Scalar>::value == false), 1138 THIS_TYPE_IS_NOT_SUPPORTED); 1139 return Scalar(0); 1140 } 1141 }; 1142 1143 /* Continued fraction expansion #1 for incomplete beta integral (small_branch = True) 1144 * Continued fraction expansion #2 for incomplete beta integral (small_branch = False) 1145 */ 1146 template <typename Scalar> 1147 struct incbeta_cfe { 1148 EIGEN_DEVICE_FUNC 1149 static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x, bool small_branch) { 1150 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, float>::value || 1151 internal::is_same<Scalar, double>::value), 1152 THIS_TYPE_IS_NOT_SUPPORTED); 1153 const Scalar big = cephes_helper<Scalar>::big(); 1154 const Scalar machep = cephes_helper<Scalar>::machep(); 1155 const Scalar biginv = cephes_helper<Scalar>::biginv(); 1156 1157 const Scalar zero = 0; 1158 const Scalar one = 1; 1159 const Scalar two = 2; 1160 1161 Scalar xk, pk, pkm1, pkm2, qk, qkm1, qkm2; 1162 Scalar k1, k2, k3, k4, k5, k6, k7, k8, k26update; 1163 Scalar ans; 1164 int n; 1165 1166 const int num_iters = (internal::is_same<Scalar, float>::value) ? 100 : 300; 1167 const Scalar thresh = 1168 (internal::is_same<Scalar, float>::value) ? machep : Scalar(3) * machep; 1169 Scalar r = (internal::is_same<Scalar, float>::value) ? zero : one; 1170 1171 if (small_branch) { 1172 k1 = a; 1173 k2 = a + b; 1174 k3 = a; 1175 k4 = a + one; 1176 k5 = one; 1177 k6 = b - one; 1178 k7 = k4; 1179 k8 = a + two; 1180 k26update = one; 1181 } else { 1182 k1 = a; 1183 k2 = b - one; 1184 k3 = a; 1185 k4 = a + one; 1186 k5 = one; 1187 k6 = a + b; 1188 k7 = a + one; 1189 k8 = a + two; 1190 k26update = -one; 1191 x = x / (one - x); 1192 } 1193 1194 pkm2 = zero; 1195 qkm2 = one; 1196 pkm1 = one; 1197 qkm1 = one; 1198 ans = one; 1199 n = 0; 1200 1201 do { 1202 xk = -(x * k1 * k2) / (k3 * k4); 1203 pk = pkm1 + pkm2 * xk; 1204 qk = qkm1 + qkm2 * xk; 1205 pkm2 = pkm1; 1206 pkm1 = pk; 1207 qkm2 = qkm1; 1208 qkm1 = qk; 1209 1210 xk = (x * k5 * k6) / (k7 * k8); 1211 pk = pkm1 + pkm2 * xk; 1212 qk = qkm1 + qkm2 * xk; 1213 pkm2 = pkm1; 1214 pkm1 = pk; 1215 qkm2 = qkm1; 1216 qkm1 = qk; 1217 1218 if (qk != zero) { 1219 r = pk / qk; 1220 if (numext::abs(ans - r) < numext::abs(r) * thresh) { 1221 return r; 1222 } 1223 ans = r; 1224 } 1225 1226 k1 += one; 1227 k2 += k26update; 1228 k3 += two; 1229 k4 += two; 1230 k5 += one; 1231 k6 -= k26update; 1232 k7 += two; 1233 k8 += two; 1234 1235 if ((numext::abs(qk) + numext::abs(pk)) > big) { 1236 pkm2 *= biginv; 1237 pkm1 *= biginv; 1238 qkm2 *= biginv; 1239 qkm1 *= biginv; 1240 } 1241 if ((numext::abs(qk) < biginv) || (numext::abs(pk) < biginv)) { 1242 pkm2 *= big; 1243 pkm1 *= big; 1244 qkm2 *= big; 1245 qkm1 *= big; 1246 } 1247 } while (++n < num_iters); 1248 1249 return ans; 1250 } 1251 }; 1252 1253 /* Helper functions depending on the Scalar type */ 1254 template <typename Scalar> 1255 struct betainc_helper {}; 1256 1257 template <> 1258 struct betainc_helper<float> { 1259 /* Core implementation, assumes a large (> 1.0) */ 1260 EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE float incbsa(float aa, float bb, 1261 float xx) { 1262 float ans, a, b, t, x, onemx; 1263 bool reversed_a_b = false; 1264 1265 onemx = 1.0f - xx; 1266 1267 /* see if x is greater than the mean */ 1268 if (xx > (aa / (aa + bb))) { 1269 reversed_a_b = true; 1270 a = bb; 1271 b = aa; 1272 t = xx; 1273 x = onemx; 1274 } else { 1275 a = aa; 1276 b = bb; 1277 t = onemx; 1278 x = xx; 1279 } 1280 1281 /* Choose expansion for optimal convergence */ 1282 if (b > 10.0f) { 1283 if (numext::abs(b * x / a) < 0.3f) { 1284 t = betainc_helper<float>::incbps(a, b, x); 1285 if (reversed_a_b) t = 1.0f - t; 1286 return t; 1287 } 1288 } 1289 1290 ans = x * (a + b - 2.0f) / (a - 1.0f); 1291 if (ans < 1.0f) { 1292 ans = incbeta_cfe<float>::run(a, b, x, true /* small_branch */); 1293 t = b * numext::log(t); 1294 } else { 1295 ans = incbeta_cfe<float>::run(a, b, x, false /* small_branch */); 1296 t = (b - 1.0f) * numext::log(t); 1297 } 1298 1299 t += a * numext::log(x) + lgamma_impl<float>::run(a + b) - 1300 lgamma_impl<float>::run(a) - lgamma_impl<float>::run(b); 1301 t += numext::log(ans / a); 1302 t = numext::exp(t); 1303 1304 if (reversed_a_b) t = 1.0f - t; 1305 return t; 1306 } 1307 1308 EIGEN_DEVICE_FUNC 1309 static EIGEN_STRONG_INLINE float incbps(float a, float b, float x) { 1310 float t, u, y, s; 1311 const float machep = cephes_helper<float>::machep(); 1312 1313 y = a * numext::log(x) + (b - 1.0f) * numext::log1p(-x) - numext::log(a); 1314 y -= lgamma_impl<float>::run(a) + lgamma_impl<float>::run(b); 1315 y += lgamma_impl<float>::run(a + b); 1316 1317 t = x / (1.0f - x); 1318 s = 0.0f; 1319 u = 1.0f; 1320 do { 1321 b -= 1.0f; 1322 if (b == 0.0f) { 1323 break; 1324 } 1325 a += 1.0f; 1326 u *= t * b / a; 1327 s += u; 1328 } while (numext::abs(u) > machep); 1329 1330 return numext::exp(y) * (1.0f + s); 1331 } 1332 }; 1333 1334 template <> 1335 struct betainc_impl<float> { 1336 EIGEN_DEVICE_FUNC 1337 static float run(float a, float b, float x) { 1338 const float nan = NumTraits<float>::quiet_NaN(); 1339 float ans, t; 1340 1341 if (a <= 0.0f) return nan; 1342 if (b <= 0.0f) return nan; 1343 if ((x <= 0.0f) || (x >= 1.0f)) { 1344 if (x == 0.0f) return 0.0f; 1345 if (x == 1.0f) return 1.0f; 1346 // mtherr("betaincf", DOMAIN); 1347 return nan; 1348 } 1349 1350 /* transformation for small aa */ 1351 if (a <= 1.0f) { 1352 ans = betainc_helper<float>::incbsa(a + 1.0f, b, x); 1353 t = a * numext::log(x) + b * numext::log1p(-x) + 1354 lgamma_impl<float>::run(a + b) - lgamma_impl<float>::run(a + 1.0f) - 1355 lgamma_impl<float>::run(b); 1356 return (ans + numext::exp(t)); 1357 } else { 1358 return betainc_helper<float>::incbsa(a, b, x); 1359 } 1360 } 1361 }; 1362 1363 template <> 1364 struct betainc_helper<double> { 1365 EIGEN_DEVICE_FUNC 1366 static EIGEN_STRONG_INLINE double incbps(double a, double b, double x) { 1367 const double machep = cephes_helper<double>::machep(); 1368 1369 double s, t, u, v, n, t1, z, ai; 1370 1371 ai = 1.0 / a; 1372 u = (1.0 - b) * x; 1373 v = u / (a + 1.0); 1374 t1 = v; 1375 t = u; 1376 n = 2.0; 1377 s = 0.0; 1378 z = machep * ai; 1379 while (numext::abs(v) > z) { 1380 u = (n - b) * x / n; 1381 t *= u; 1382 v = t / (a + n); 1383 s += v; 1384 n += 1.0; 1385 } 1386 s += t1; 1387 s += ai; 1388 1389 u = a * numext::log(x); 1390 // TODO: gamma() is not directly implemented in Eigen. 1391 /* 1392 if ((a + b) < maxgam && numext::abs(u) < maxlog) { 1393 t = gamma(a + b) / (gamma(a) * gamma(b)); 1394 s = s * t * pow(x, a); 1395 } else { 1396 */ 1397 t = lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) - 1398 lgamma_impl<double>::run(b) + u + numext::log(s); 1399 return s = numext::exp(t); 1400 } 1401 }; 1402 1403 template <> 1404 struct betainc_impl<double> { 1405 EIGEN_DEVICE_FUNC 1406 static double run(double aa, double bb, double xx) { 1407 const double nan = NumTraits<double>::quiet_NaN(); 1408 const double machep = cephes_helper<double>::machep(); 1409 // const double maxgam = 171.624376956302725; 1410 1411 double a, b, t, x, xc, w, y; 1412 bool reversed_a_b = false; 1413 1414 if (aa <= 0.0 || bb <= 0.0) { 1415 return nan; // goto domerr; 1416 } 1417 1418 if ((xx <= 0.0) || (xx >= 1.0)) { 1419 if (xx == 0.0) return (0.0); 1420 if (xx == 1.0) return (1.0); 1421 // mtherr("incbet", DOMAIN); 1422 return nan; 1423 } 1424 1425 if ((bb * xx) <= 1.0 && xx <= 0.95) { 1426 return betainc_helper<double>::incbps(aa, bb, xx); 1427 } 1428 1429 w = 1.0 - xx; 1430 1431 /* Reverse a and b if x is greater than the mean. */ 1432 if (xx > (aa / (aa + bb))) { 1433 reversed_a_b = true; 1434 a = bb; 1435 b = aa; 1436 xc = xx; 1437 x = w; 1438 } else { 1439 a = aa; 1440 b = bb; 1441 xc = w; 1442 x = xx; 1443 } 1444 1445 if (reversed_a_b && (b * x) <= 1.0 && x <= 0.95) { 1446 t = betainc_helper<double>::incbps(a, b, x); 1447 if (t <= machep) { 1448 t = 1.0 - machep; 1449 } else { 1450 t = 1.0 - t; 1451 } 1452 return t; 1453 } 1454 1455 /* Choose expansion for better convergence. */ 1456 y = x * (a + b - 2.0) - (a - 1.0); 1457 if (y < 0.0) { 1458 w = incbeta_cfe<double>::run(a, b, x, true /* small_branch */); 1459 } else { 1460 w = incbeta_cfe<double>::run(a, b, x, false /* small_branch */) / xc; 1461 } 1462 1463 /* Multiply w by the factor 1464 a b _ _ _ 1465 x (1-x) | (a+b) / ( a | (a) | (b) ) . */ 1466 1467 y = a * numext::log(x); 1468 t = b * numext::log(xc); 1469 // TODO: gamma is not directly implemented in Eigen. 1470 /* 1471 if ((a + b) < maxgam && numext::abs(y) < maxlog && numext::abs(t) < maxlog) 1472 { 1473 t = pow(xc, b); 1474 t *= pow(x, a); 1475 t /= a; 1476 t *= w; 1477 t *= gamma(a + b) / (gamma(a) * gamma(b)); 1478 } else { 1479 */ 1480 /* Resort to logarithms. */ 1481 y += t + lgamma_impl<double>::run(a + b) - lgamma_impl<double>::run(a) - 1482 lgamma_impl<double>::run(b); 1483 y += numext::log(w / a); 1484 t = numext::exp(y); 1485 1486 /* } */ 1487 // done: 1488 1489 if (reversed_a_b) { 1490 if (t <= machep) { 1491 t = 1.0 - machep; 1492 } else { 1493 t = 1.0 - t; 1494 } 1495 } 1496 return t; 1497 } 1498 }; 1499 1500 #endif // EIGEN_HAS_C99_MATH 1501 1502 } // end namespace internal 1503 1504 namespace numext { 1505 1506 template <typename Scalar> 1507 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(lgamma, Scalar) 1508 lgamma(const Scalar& x) { 1509 return EIGEN_MATHFUNC_IMPL(lgamma, Scalar)::run(x); 1510 } 1511 1512 template <typename Scalar> 1513 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(digamma, Scalar) 1514 digamma(const Scalar& x) { 1515 return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x); 1516 } 1517 1518 template <typename Scalar> 1519 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(zeta, Scalar) 1520 zeta(const Scalar& x, const Scalar& q) { 1521 return EIGEN_MATHFUNC_IMPL(zeta, Scalar)::run(x, q); 1522 } 1523 1524 template <typename Scalar> 1525 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(polygamma, Scalar) 1526 polygamma(const Scalar& n, const Scalar& x) { 1527 return EIGEN_MATHFUNC_IMPL(polygamma, Scalar)::run(n, x); 1528 } 1529 1530 template <typename Scalar> 1531 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erf, Scalar) 1532 erf(const Scalar& x) { 1533 return EIGEN_MATHFUNC_IMPL(erf, Scalar)::run(x); 1534 } 1535 1536 template <typename Scalar> 1537 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(erfc, Scalar) 1538 erfc(const Scalar& x) { 1539 return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x); 1540 } 1541 1542 template <typename Scalar> 1543 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igamma, Scalar) 1544 igamma(const Scalar& a, const Scalar& x) { 1545 return EIGEN_MATHFUNC_IMPL(igamma, Scalar)::run(a, x); 1546 } 1547 1548 template <typename Scalar> 1549 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(igammac, Scalar) 1550 igammac(const Scalar& a, const Scalar& x) { 1551 return EIGEN_MATHFUNC_IMPL(igammac, Scalar)::run(a, x); 1552 } 1553 1554 template <typename Scalar> 1555 EIGEN_DEVICE_FUNC inline EIGEN_MATHFUNC_RETVAL(betainc, Scalar) 1556 betainc(const Scalar& a, const Scalar& b, const Scalar& x) { 1557 return EIGEN_MATHFUNC_IMPL(betainc, Scalar)::run(a, b, x); 1558 } 1559 1560 } // end namespace numext 1561 1562 1563 } // end namespace Eigen 1564 1565 #endif // EIGEN_SPECIAL_FUNCTIONS_H 1566