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