SpecialFunctionsImpl.h
Go to the documentation of this file.
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 
40 /****************************************************************************
41  * Implementation of lgamma, requires C++11/C99 *
42  ****************************************************************************/
43 
44 template <typename Scalar>
45 struct lgamma_impl {
47  static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
49  THIS_TYPE_IS_NOT_SUPPORTED);
50  return Scalar(0);
51  }
52 };
53 
54 template <typename Scalar>
55 struct lgamma_retval {
56  typedef Scalar type;
57 };
58 
59 #if EIGEN_HAS_C99_MATH
60 // Since glibc 2.19
61 #if defined(__GLIBC__) && ((__GLIBC__>=2 && __GLIBC_MINOR__ >= 19) || __GLIBC__>2) \
62  && (defined(_DEFAULT_SOURCE) || defined(_BSD_SOURCE) || defined(_SVID_SOURCE))
63 #define EIGEN_HAS_LGAMMA_R
64 #endif
65 
66 // Glibc versions before 2.19
67 #if defined(__GLIBC__) && ((__GLIBC__==2 && __GLIBC_MINOR__ < 19) || __GLIBC__<2) \
68  && (defined(_BSD_SOURCE) || defined(_SVID_SOURCE))
69 #define EIGEN_HAS_LGAMMA_R
70 #endif
71 
72 template <>
73 struct lgamma_impl<float> {
75  static EIGEN_STRONG_INLINE float run(float x) {
76 #if !defined(EIGEN_GPU_COMPILE_PHASE) && defined (EIGEN_HAS_LGAMMA_R) && !defined(__APPLE__)
77  int dummy;
78  return ::lgammaf_r(x, &dummy);
79 #elif defined(SYCL_DEVICE_ONLY)
80  return cl::sycl::lgamma(x);
81 #else
82  return ::lgammaf(x);
83 #endif
84  }
85 };
86 
87 template <>
88 struct lgamma_impl<double> {
90  static EIGEN_STRONG_INLINE double run(double x) {
91 #if !defined(EIGEN_GPU_COMPILE_PHASE) && defined(EIGEN_HAS_LGAMMA_R) && !defined(__APPLE__)
92  int dummy;
93  return ::lgamma_r(x, &dummy);
94 #elif defined(SYCL_DEVICE_ONLY)
95  return cl::sycl::lgamma(x);
96 #else
98 #endif
99  }
100 };
101 
102 #undef EIGEN_HAS_LGAMMA_R
103 #endif
104 
105 /****************************************************************************
106  * Implementation of digamma (psi), based on Cephes *
107  ****************************************************************************/
108 
109 template <typename Scalar>
111  typedef Scalar type;
112 };
113 
114 /*
115  *
116  * Polynomial evaluation helper for the Psi (digamma) function.
117  *
118  * digamma_impl_maybe_poly::run(s) evaluates the asymptotic Psi expansion for
119  * input Scalar s, assuming s is above 10.0.
120  *
121  * If s is above a certain threshold for the given Scalar type, zero
122  * is returned. Otherwise the polynomial is evaluated with enough
123  * coefficients for results matching Scalar machine precision.
124  *
125  *
126  */
127 template <typename Scalar>
132  THIS_TYPE_IS_NOT_SUPPORTED);
133  return Scalar(0);
134  }
135 };
136 
137 
138 template <>
141  static EIGEN_STRONG_INLINE float run(const float s) {
142  const float A[] = {
143  -4.16666666666666666667E-3f,
144  3.96825396825396825397E-3f,
145  -8.33333333333333333333E-3f,
146  8.33333333333333333333E-2f
147  };
148 
149  float z;
150  if (s < 1.0e8f) {
151  z = 1.0f / (s * s);
153  } else return 0.0f;
154  }
155 };
156 
157 template <>
158 struct digamma_impl_maybe_poly<double> {
160  static EIGEN_STRONG_INLINE double run(const double s) {
161  const double A[] = {
162  8.33333333333333333333E-2,
163  -2.10927960927960927961E-2,
164  7.57575757575757575758E-3,
165  -4.16666666666666666667E-3,
166  3.96825396825396825397E-3,
167  -8.33333333333333333333E-3,
168  8.33333333333333333333E-2
169  };
170 
171  double z;
172  if (s < 1.0e17) {
173  z = 1.0 / (s * s);
175  }
176  else return 0.0;
177  }
178 };
179 
180 template <typename Scalar>
181 struct digamma_impl {
183  static Scalar run(Scalar x) {
184  /*
185  *
186  * Psi (digamma) function (modified for Eigen)
187  *
188  *
189  * SYNOPSIS:
190  *
191  * double x, y, psi();
192  *
193  * y = psi( x );
194  *
195  *
196  * DESCRIPTION:
197  *
198  * d -
199  * psi(x) = -- ln | (x)
200  * dx
201  *
202  * is the logarithmic derivative of the gamma function.
203  * For integer x,
204  * n-1
205  * -
206  * psi(n) = -EUL + > 1/k.
207  * -
208  * k=1
209  *
210  * If x is negative, it is transformed to a positive argument by the
211  * reflection formula psi(1-x) = psi(x) + pi cot(pi x).
212  * For general positive x, the argument is made greater than 10
213  * using the recurrence psi(x+1) = psi(x) + 1/x.
214  * Then the following asymptotic expansion is applied:
215  *
216  * inf. B
217  * - 2k
218  * psi(x) = log(x) - 1/2x - > -------
219  * - 2k
220  * k=1 2k x
221  *
222  * where the B2k are Bernoulli numbers.
223  *
224  * ACCURACY (float):
225  * Relative error (except absolute when |psi| < 1):
226  * arithmetic domain # trials peak rms
227  * IEEE 0,30 30000 1.3e-15 1.4e-16
228  * IEEE -30,0 40000 1.5e-15 2.2e-16
229  *
230  * ACCURACY (double):
231  * Absolute error, relative when |psi| > 1 :
232  * arithmetic domain # trials peak rms
233  * IEEE -33,0 30000 8.2e-7 1.2e-7
234  * IEEE 0,33 100000 7.3e-7 7.7e-8
235  *
236  * ERROR MESSAGES:
237  * message condition value returned
238  * psi singularity x integer <=0 INFINITY
239  */
240 
241  Scalar p, q, nz, s, w, y;
242  bool negative = false;
243 
244  const Scalar nan = NumTraits<Scalar>::quiet_NaN();
245  const Scalar m_pi = Scalar(EIGEN_PI);
246 
247  const Scalar zero = Scalar(0);
248  const Scalar one = Scalar(1);
249  const Scalar half = Scalar(0.5);
250  nz = zero;
251 
252  if (x <= zero) {
253  negative = true;
254  q = x;
255  p = numext::floor(q);
256  if (p == q) {
257  return nan;
258  }
259  /* Remove the zeros of tan(m_pi x)
260  * by subtracting the nearest integer from x
261  */
262  nz = q - p;
263  if (nz != half) {
264  if (nz > half) {
265  p += one;
266  nz = q - p;
267  }
268  nz = m_pi / numext::tan(m_pi * nz);
269  }
270  else {
271  nz = zero;
272  }
273  x = one - x;
274  }
275 
276  /* use the recurrence psi(x+1) = psi(x) + 1/x. */
277  s = x;
278  w = zero;
279  while (s < Scalar(10)) {
280  w += one / s;
281  s += one;
282  }
283 
285 
286  y = numext::log(s) - (half / s) - y - w;
287 
288  return (negative) ? y - nz : y;
289  }
290 };
291 
292 /****************************************************************************
293  * Implementation of erf, requires C++11/C99 *
294  ****************************************************************************/
295 
303 template <typename T>
305  // Clamp the inputs to the range [-4, 4] since anything outside
306  // this range is +/-1.0f in single-precision.
307  const T plus_4 = pset1<T>(4.f);
308  const T minus_4 = pset1<T>(-4.f);
309  const T x = pmax(pmin(a_x, plus_4), minus_4);
310  // The monomial coefficients of the numerator polynomial (odd).
311  const T alpha_1 = pset1<T>(-1.60960333262415e-02f);
312  const T alpha_3 = pset1<T>(-2.95459980854025e-03f);
313  const T alpha_5 = pset1<T>(-7.34990630326855e-04f);
314  const T alpha_7 = pset1<T>(-5.69250639462346e-05f);
315  const T alpha_9 = pset1<T>(-2.10102402082508e-06f);
316  const T alpha_11 = pset1<T>(2.77068142495902e-08f);
317  const T alpha_13 = pset1<T>(-2.72614225801306e-10f);
318 
319  // The monomial coefficients of the denominator polynomial (even).
320  const T beta_0 = pset1<T>(-1.42647390514189e-02f);
321  const T beta_2 = pset1<T>(-7.37332916720468e-03f);
322  const T beta_4 = pset1<T>(-1.68282697438203e-03f);
323  const T beta_6 = pset1<T>(-2.13374055278905e-04f);
324  const T beta_8 = pset1<T>(-1.45660718464996e-05f);
325 
326  // Since the polynomials are odd/even, we need x^2.
327  const T x2 = pmul(x, x);
328 
329  // Evaluate the numerator polynomial p.
330  T p = pmadd(x2, alpha_13, alpha_11);
331  p = pmadd(x2, p, alpha_9);
332  p = pmadd(x2, p, alpha_7);
333  p = pmadd(x2, p, alpha_5);
334  p = pmadd(x2, p, alpha_3);
335  p = pmadd(x2, p, alpha_1);
336  p = pmul(x, p);
337 
338  // Evaluate the denominator polynomial p.
339  T q = pmadd(x2, beta_8, beta_6);
340  q = pmadd(x2, q, beta_4);
341  q = pmadd(x2, q, beta_2);
342  q = pmadd(x2, q, beta_0);
343 
344  // Divide the numerator by the denominator.
345  return pdiv(p, q);
346 }
347 
348 template <typename T>
349 struct erf_impl {
351  static EIGEN_STRONG_INLINE T run(const T& x) {
352  return generic_fast_erf_float(x);
353  }
354 };
355 
356 template <typename Scalar>
357 struct erf_retval {
358  typedef Scalar type;
359 };
360 
361 #if EIGEN_HAS_C99_MATH
362 template <>
363 struct erf_impl<float> {
365  static EIGEN_STRONG_INLINE float run(float x) {
366 #if defined(SYCL_DEVICE_ONLY)
367  return cl::sycl::erf(x);
368 #else
369  return generic_fast_erf_float(x);
370 #endif
371  }
372 };
373 
374 template <>
375 struct erf_impl<double> {
377  static EIGEN_STRONG_INLINE double run(double x) {
378 #if defined(SYCL_DEVICE_ONLY)
379  return cl::sycl::erf(x);
380 #else
381  return ::erf(x);
382 #endif
383  }
384 };
385 #endif // EIGEN_HAS_C99_MATH
386 
387 /***************************************************************************
388 * Implementation of erfc, requires C++11/C99 *
389 ****************************************************************************/
390 
391 template <typename Scalar>
392 struct erfc_impl {
394  static EIGEN_STRONG_INLINE Scalar run(const Scalar) {
396  THIS_TYPE_IS_NOT_SUPPORTED);
397  return Scalar(0);
398  }
399 };
400 
401 template <typename Scalar>
402 struct erfc_retval {
403  typedef Scalar type;
404 };
405 
406 #if EIGEN_HAS_C99_MATH
407 template <>
408 struct erfc_impl<float> {
410  static EIGEN_STRONG_INLINE float run(const float x) {
411 #if defined(SYCL_DEVICE_ONLY)
412  return cl::sycl::erfc(x);
413 #else
414  return ::erfcf(x);
415 #endif
416  }
417 };
418 
419 template <>
420 struct erfc_impl<double> {
422  static EIGEN_STRONG_INLINE double run(const double x) {
423 #if defined(SYCL_DEVICE_ONLY)
424  return cl::sycl::erfc(x);
425 #else
426  return ::erfc(x);
427 #endif
428  }
429 };
430 #endif // EIGEN_HAS_C99_MATH
431 
432 
433 /***************************************************************************
434 * Implementation of ndtri. *
435 ****************************************************************************/
436 
437 /* Inverse of Normal distribution function (modified for Eigen).
438  *
439  *
440  * SYNOPSIS:
441  *
442  * double x, y, ndtri();
443  *
444  * x = ndtri( y );
445  *
446  *
447  *
448  * DESCRIPTION:
449  *
450  * Returns the argument, x, for which the area under the
451  * Gaussian probability density function (integrated from
452  * minus infinity to x) is equal to y.
453  *
454  *
455  * For small arguments 0 < y < exp(-2), the program computes
456  * z = sqrt( -2.0 * log(y) ); then the approximation is
457  * x = z - log(z)/z - (1/z) P(1/z) / Q(1/z).
458  * There are two rational functions P/Q, one for 0 < y < exp(-32)
459  * and the other for y up to exp(-2). For larger arguments,
460  * w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
461  *
462  *
463  * ACCURACY:
464  *
465  * Relative error:
466  * arithmetic domain # trials peak rms
467  * DEC 0.125, 1 5500 9.5e-17 2.1e-17
468  * DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17
469  * IEEE 0.125, 1 20000 7.2e-16 1.3e-16
470  * IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
471  *
472  *
473  * ERROR MESSAGES:
474  *
475  * message condition value returned
476  * ndtri domain x <= 0 -MAXNUM
477  * ndtri domain x >= 1 MAXNUM
478  *
479  */
480  /*
481  Cephes Math Library Release 2.2: June, 1992
482  Copyright 1985, 1987, 1992 by Stephen L. Moshier
483  Direct inquiries to 30 Frost Street, Cambridge, MA 02140
484  */
485 
486 
487 // TODO: Add a cheaper approximation for float.
488 
489 
490 template<typename T>
492  const T& should_flipsign, const T& x) {
493  typedef typename unpacket_traits<T>::type Scalar;
494  const T sign_mask = pset1<T>(Scalar(-0.0));
495  T sign_bit = pand<T>(should_flipsign, sign_mask);
496  return pxor<T>(sign_bit, x);
497 }
498 
499 template<>
501  const double& should_flipsign, const double& x) {
502  return should_flipsign == 0 ? x : -x;
503 }
504 
505 template<>
507  const float& should_flipsign, const float& x) {
508  return should_flipsign == 0 ? x : -x;
509 }
510 
511 // We split this computation in to two so that in the scalar path
512 // only one branch is evaluated (due to our template specialization of pselect
513 // being an if statement.)
514 
515 template <typename T, typename ScalarType>
517  const ScalarType p0[] = {
518  ScalarType(-5.99633501014107895267e1),
519  ScalarType(9.80010754185999661536e1),
520  ScalarType(-5.66762857469070293439e1),
521  ScalarType(1.39312609387279679503e1),
522  ScalarType(-1.23916583867381258016e0)
523  };
524  const ScalarType q0[] = {
525  ScalarType(1.0),
526  ScalarType(1.95448858338141759834e0),
527  ScalarType(4.67627912898881538453e0),
528  ScalarType(8.63602421390890590575e1),
529  ScalarType(-2.25462687854119370527e2),
530  ScalarType(2.00260212380060660359e2),
531  ScalarType(-8.20372256168333339912e1),
532  ScalarType(1.59056225126211695515e1),
533  ScalarType(-1.18331621121330003142e0)
534  };
535  const T sqrt2pi = pset1<T>(ScalarType(2.50662827463100050242e0));
536  const T half = pset1<T>(ScalarType(0.5));
537  T c, c2, ndtri_gt_exp_neg_two;
538 
539  c = psub(b, half);
540  c2 = pmul(c, c);
541  ndtri_gt_exp_neg_two = pmadd(c, pmul(
542  c2, pdiv(
545  return pmul(ndtri_gt_exp_neg_two, sqrt2pi);
546 }
547 
548 template <typename T, typename ScalarType>
550  const T& b, const T& should_flipsign) {
551  /* Approximation for interval z = sqrt(-2 log a ) between 2 and 8
552  * i.e., a between exp(-2) = .135 and exp(-32) = 1.27e-14.
553  */
554  const ScalarType p1[] = {
555  ScalarType(4.05544892305962419923e0),
556  ScalarType(3.15251094599893866154e1),
557  ScalarType(5.71628192246421288162e1),
558  ScalarType(4.40805073893200834700e1),
559  ScalarType(1.46849561928858024014e1),
560  ScalarType(2.18663306850790267539e0),
561  ScalarType(-1.40256079171354495875e-1),
562  ScalarType(-3.50424626827848203418e-2),
563  ScalarType(-8.57456785154685413611e-4)
564  };
565  const ScalarType q1[] = {
566  ScalarType(1.0),
567  ScalarType(1.57799883256466749731e1),
568  ScalarType(4.53907635128879210584e1),
569  ScalarType(4.13172038254672030440e1),
570  ScalarType(1.50425385692907503408e1),
571  ScalarType(2.50464946208309415979e0),
572  ScalarType(-1.42182922854787788574e-1),
573  ScalarType(-3.80806407691578277194e-2),
574  ScalarType(-9.33259480895457427372e-4)
575  };
576  /* Approximation for interval z = sqrt(-2 log a ) between 8 and 64
577  * i.e., a between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
578  */
579  const ScalarType p2[] = {
580  ScalarType(3.23774891776946035970e0),
581  ScalarType(6.91522889068984211695e0),
582  ScalarType(3.93881025292474443415e0),
583  ScalarType(1.33303460815807542389e0),
584  ScalarType(2.01485389549179081538e-1),
585  ScalarType(1.23716634817820021358e-2),
586  ScalarType(3.01581553508235416007e-4),
587  ScalarType(2.65806974686737550832e-6),
588  ScalarType(6.23974539184983293730e-9)
589  };
590  const ScalarType q2[] = {
591  ScalarType(1.0),
592  ScalarType(6.02427039364742014255e0),
593  ScalarType(3.67983563856160859403e0),
594  ScalarType(1.37702099489081330271e0),
595  ScalarType(2.16236993594496635890e-1),
596  ScalarType(1.34204006088543189037e-2),
597  ScalarType(3.28014464682127739104e-4),
598  ScalarType(2.89247864745380683936e-6),
599  ScalarType(6.79019408009981274425e-9)
600  };
601  const T eight = pset1<T>(ScalarType(8.0));
602  const T one = pset1<T>(ScalarType(1));
603  const T neg_two = pset1<T>(ScalarType(-2));
604  T x, x0, x1, z;
605 
606  x = psqrt(pmul(neg_two, plog(b)));
607  x0 = psub(x, pdiv(plog(x), x));
608  z = pdiv(one, x);
609  x1 = pmul(
610  z, pselect(
611  pcmp_lt(x, eight),
616  return flipsign(should_flipsign, psub(x0, x1));
617 }
618 
619 template <typename T, typename ScalarType>
621 T generic_ndtri(const T& a) {
622  const T maxnum = pset1<T>(NumTraits<ScalarType>::infinity());
623  const T neg_maxnum = pset1<T>(-NumTraits<ScalarType>::infinity());
624 
625  const T zero = pset1<T>(ScalarType(0));
626  const T one = pset1<T>(ScalarType(1));
627  // exp(-2)
628  const T exp_neg_two = pset1<T>(ScalarType(0.13533528323661269189));
629  T b, ndtri, should_flipsign;
630 
631  should_flipsign = pcmp_le(a, psub(one, exp_neg_two));
632  b = pselect(should_flipsign, a, psub(one, a));
633 
634  ndtri = pselect(
635  pcmp_lt(exp_neg_two, b),
636  generic_ndtri_gt_exp_neg_two<T, ScalarType>(b),
637  generic_ndtri_lt_exp_neg_two<T, ScalarType>(b, should_flipsign));
638 
639  return pselect(
640  pcmp_le(a, zero), neg_maxnum,
641  pselect(pcmp_le(one, a), maxnum, ndtri));
642 }
643 
644 template <typename Scalar>
645 struct ndtri_retval {
646  typedef Scalar type;
647 };
648 
649 #if !EIGEN_HAS_C99_MATH
650 
651 template <typename Scalar>
652 struct ndtri_impl {
656  THIS_TYPE_IS_NOT_SUPPORTED);
657  return Scalar(0);
658  }
659 };
660 
661 # else
662 
663 template <typename Scalar>
664 struct ndtri_impl {
666  static EIGEN_STRONG_INLINE Scalar run(const Scalar x) {
667  return generic_ndtri<Scalar, Scalar>(x);
668  }
669 };
670 
671 #endif // EIGEN_HAS_C99_MATH
672 
673 
674 /**************************************************************************************************************
675  * Implementation of igammac (complemented incomplete gamma integral), based on Cephes but requires C++11/C99 *
676  **************************************************************************************************************/
677 
678 template <typename Scalar>
680  typedef Scalar type;
681 };
682 
683 // NOTE: cephes_helper is also used to implement zeta
684 template <typename Scalar>
687  static EIGEN_STRONG_INLINE Scalar machep() { assert(false && "machep not supported for this type"); return 0.0; }
689  static EIGEN_STRONG_INLINE Scalar big() { assert(false && "big not supported for this type"); return 0.0; }
691  static EIGEN_STRONG_INLINE Scalar biginv() { assert(false && "biginv not supported for this type"); return 0.0; }
692 };
693 
694 template <>
697  static EIGEN_STRONG_INLINE float machep() {
698  return NumTraits<float>::epsilon() / 2; // 1.0 - machep == 1.0
699  }
701  static EIGEN_STRONG_INLINE float big() {
702  // use epsneg (1.0 - epsneg == 1.0)
703  return 1.0f / (NumTraits<float>::epsilon() / 2);
704  }
706  static EIGEN_STRONG_INLINE float biginv() {
707  // epsneg
708  return machep();
709  }
710 };
711 
712 template <>
713 struct cephes_helper<double> {
715  static EIGEN_STRONG_INLINE double machep() {
716  return NumTraits<double>::epsilon() / 2; // 1.0 - machep == 1.0
717  }
719  static EIGEN_STRONG_INLINE double big() {
720  return 1.0 / NumTraits<double>::epsilon();
721  }
723  static EIGEN_STRONG_INLINE double biginv() {
724  // inverse of eps
726  }
727 };
728 
730 
731 template <typename Scalar>
734  /* Compute x**a * exp(-x) / gamma(a) */
735  Scalar logax = a * numext::log(x) - x - lgamma_impl<Scalar>::run(a);
736  if (logax < -numext::log(NumTraits<Scalar>::highest()) ||
737  // Assuming x and a aren't Nan.
738  (numext::isnan)(logax)) {
739  return Scalar(0);
740  }
741  return numext::exp(logax);
742 }
743 
744 template <typename Scalar, IgammaComputationMode mode>
747  /* Returns the maximum number of internal iterations for igamma computation.
748  */
749  if (mode == VALUE) {
750  return 2000;
751  }
752 
754  return 200;
756  return 500;
757  } else {
758  return 2000;
759  }
760 }
761 
762 template <typename Scalar, IgammaComputationMode mode>
764  /* Computes igamc(a, x) or derivative (depending on the mode)
765  * using the continued fraction expansion of the complementary
766  * incomplete Gamma function.
767  *
768  * Preconditions:
769  * a > 0
770  * x >= 1
771  * x >= a
772  */
774  static Scalar run(Scalar a, Scalar x) {
775  const Scalar zero = 0;
776  const Scalar one = 1;
777  const Scalar two = 2;
778  const Scalar machep = cephes_helper<Scalar>::machep();
781 
782  if ((numext::isinf)(x)) {
783  return zero;
784  }
785 
786  Scalar ax = main_igamma_term<Scalar>(a, x);
787  // This is independent of mode. If this value is zero,
788  // then the function value is zero. If the function value is zero,
789  // then we are in a neighborhood where the function value evalutes to zero,
790  // so the derivative is zero.
791  if (ax == zero) {
792  return zero;
793  }
794 
795  // continued fraction
796  Scalar y = one - a;
797  Scalar z = x + y + one;
798  Scalar c = zero;
799  Scalar pkm2 = one;
800  Scalar qkm2 = x;
801  Scalar pkm1 = x + one;
802  Scalar qkm1 = z * x;
803  Scalar ans = pkm1 / qkm1;
804 
805  Scalar dpkm2_da = zero;
806  Scalar dqkm2_da = zero;
807  Scalar dpkm1_da = zero;
808  Scalar dqkm1_da = -x;
809  Scalar dans_da = (dpkm1_da - ans * dqkm1_da) / qkm1;
810 
811  for (int i = 0; i < igamma_num_iterations<Scalar, mode>(); i++) {
812  c += one;
813  y += one;
814  z += two;
815 
816  Scalar yc = y * c;
817  Scalar pk = pkm1 * z - pkm2 * yc;
818  Scalar qk = qkm1 * z - qkm2 * yc;
819 
820  Scalar dpk_da = dpkm1_da * z - pkm1 - dpkm2_da * yc + pkm2 * c;
821  Scalar dqk_da = dqkm1_da * z - qkm1 - dqkm2_da * yc + qkm2 * c;
822 
823  if (qk != zero) {
824  Scalar ans_prev = ans;
825  ans = pk / qk;
826 
827  Scalar dans_da_prev = dans_da;
828  dans_da = (dpk_da - ans * dqk_da) / qk;
829 
830  if (mode == VALUE) {
831  if (numext::abs(ans_prev - ans) <= machep * numext::abs(ans)) {
832  break;
833  }
834  } else {
835  if (numext::abs(dans_da - dans_da_prev) <= machep) {
836  break;
837  }
838  }
839  }
840 
841  pkm2 = pkm1;
842  pkm1 = pk;
843  qkm2 = qkm1;
844  qkm1 = qk;
845 
846  dpkm2_da = dpkm1_da;
847  dpkm1_da = dpk_da;
848  dqkm2_da = dqkm1_da;
849  dqkm1_da = dqk_da;
850 
851  if (numext::abs(pk) > big) {
852  pkm2 *= biginv;
853  pkm1 *= biginv;
854  qkm2 *= biginv;
855  qkm1 *= biginv;
856 
857  dpkm2_da *= biginv;
858  dpkm1_da *= biginv;
859  dqkm2_da *= biginv;
860  dqkm1_da *= biginv;
861  }
862  }
863 
864  /* Compute x**a * exp(-x) / gamma(a) */
866  Scalar dax_da = ax * dlogax_da;
867 
868  switch (mode) {
869  case VALUE:
870  return ans * ax;
871  case DERIVATIVE:
872  return ans * dax_da + dans_da * ax;
873  case SAMPLE_DERIVATIVE:
874  default: // this is needed to suppress clang warning
875  return -(dans_da + ans * dlogax_da) * x;
876  }
877  }
878 };
879 
880 template <typename Scalar, IgammaComputationMode mode>
882  /* Computes igam(a, x) or its derivative (depending on the mode)
883  * using the series expansion of the incomplete Gamma function.
884  *
885  * Preconditions:
886  * x > 0
887  * a > 0
888  * !(x > 1 && x > a)
889  */
891  static Scalar run(Scalar a, Scalar x) {
892  const Scalar zero = 0;
893  const Scalar one = 1;
894  const Scalar machep = cephes_helper<Scalar>::machep();
895 
896  Scalar ax = main_igamma_term<Scalar>(a, x);
897 
898  // This is independent of mode. If this value is zero,
899  // then the function value is zero. If the function value is zero,
900  // then we are in a neighborhood where the function value evalutes to zero,
901  // so the derivative is zero.
902  if (ax == zero) {
903  return zero;
904  }
905 
906  ax /= a;
907 
908  /* power series */
909  Scalar r = a;
910  Scalar c = one;
911  Scalar ans = one;
912 
913  Scalar dc_da = zero;
914  Scalar dans_da = zero;
915 
916  for (int i = 0; i < igamma_num_iterations<Scalar, mode>(); i++) {
917  r += one;
918  Scalar term = x / r;
919  Scalar dterm_da = -x / (r * r);
920  dc_da = term * dc_da + dterm_da * c;
921  dans_da += dc_da;
922  c *= term;
923  ans += c;
924 
925  if (mode == VALUE) {
926  if (c <= machep * ans) {
927  break;
928  }
929  } else {
930  if (numext::abs(dc_da) <= machep * numext::abs(dans_da)) {
931  break;
932  }
933  }
934  }
935 
936  Scalar dlogax_da = numext::log(x) - digamma_impl<Scalar>::run(a + one);
937  Scalar dax_da = ax * dlogax_da;
938 
939  switch (mode) {
940  case VALUE:
941  return ans * ax;
942  case DERIVATIVE:
943  return ans * dax_da + dans_da * ax;
944  case SAMPLE_DERIVATIVE:
945  default: // this is needed to suppress clang warning
946  return -(dans_da + ans * dlogax_da) * x / a;
947  }
948  }
949 };
950 
951 #if !EIGEN_HAS_C99_MATH
952 
953 template <typename Scalar>
954 struct igammac_impl {
956  static Scalar run(Scalar a, Scalar x) {
958  THIS_TYPE_IS_NOT_SUPPORTED);
959  return Scalar(0);
960  }
961 };
962 
963 #else
964 
965 template <typename Scalar>
966 struct igammac_impl {
968  static Scalar run(Scalar a, Scalar x) {
969  /* igamc()
970  *
971  * Incomplete gamma integral (modified for Eigen)
972  *
973  *
974  *
975  * SYNOPSIS:
976  *
977  * double a, x, y, igamc();
978  *
979  * y = igamc( a, x );
980  *
981  * DESCRIPTION:
982  *
983  * The function is defined by
984  *
985  *
986  * igamc(a,x) = 1 - igam(a,x)
987  *
988  * inf.
989  * -
990  * 1 | | -t a-1
991  * = ----- | e t dt.
992  * - | |
993  * | (a) -
994  * x
995  *
996  *
997  * In this implementation both arguments must be positive.
998  * The integral is evaluated by either a power series or
999  * continued fraction expansion, depending on the relative
1000  * values of a and x.
1001  *
1002  * ACCURACY (float):
1003  *
1004  * Relative error:
1005  * arithmetic domain # trials peak rms
1006  * IEEE 0,30 30000 7.8e-6 5.9e-7
1007  *
1008  *
1009  * ACCURACY (double):
1010  *
1011  * Tested at random a, x.
1012  * a x Relative error:
1013  * arithmetic domain domain # trials peak rms
1014  * IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
1015  * IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
1016  *
1017  */
1018  /*
1019  Cephes Math Library Release 2.2: June, 1992
1020  Copyright 1985, 1987, 1992 by Stephen L. Moshier
1021  Direct inquiries to 30 Frost Street, Cambridge, MA 02140
1022  */
1023  const Scalar zero = 0;
1024  const Scalar one = 1;
1025  const Scalar nan = NumTraits<Scalar>::quiet_NaN();
1026 
1027  if ((x < zero) || (a <= zero)) {
1028  // domain error
1029  return nan;
1030  }
1031 
1032  if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans
1033  return nan;
1034  }
1035 
1036  if ((x < one) || (x < a)) {
1037  return (one - igamma_series_impl<Scalar, VALUE>::run(a, x));
1038  }
1039 
1041  }
1042 };
1043 
1044 #endif // EIGEN_HAS_C99_MATH
1045 
1046 /************************************************************************************************
1047  * Implementation of igamma (incomplete gamma integral), based on Cephes but requires C++11/C99 *
1048  ************************************************************************************************/
1049 
1050 #if !EIGEN_HAS_C99_MATH
1051 
1052 template <typename Scalar, IgammaComputationMode mode>
1057  THIS_TYPE_IS_NOT_SUPPORTED);
1058  return Scalar(0);
1059  }
1060 };
1061 
1062 #else
1063 
1064 template <typename Scalar, IgammaComputationMode mode>
1065 struct igamma_generic_impl {
1067  static Scalar run(Scalar a, Scalar x) {
1068  /* Depending on the mode, returns
1069  * - VALUE: incomplete Gamma function igamma(a, x)
1070  * - DERIVATIVE: derivative of incomplete Gamma function d/da igamma(a, x)
1071  * - SAMPLE_DERIVATIVE: implicit derivative of a Gamma random variable
1072  * x ~ Gamma(x | a, 1), dx/da = -1 / Gamma(x | a, 1) * d igamma(a, x) / dx
1073  *
1074  * Derivatives are implemented by forward-mode differentiation.
1075  */
1076  const Scalar zero = 0;
1077  const Scalar one = 1;
1078  const Scalar nan = NumTraits<Scalar>::quiet_NaN();
1079 
1080  if (x == zero) return zero;
1081 
1082  if ((x < zero) || (a <= zero)) { // domain error
1083  return nan;
1084  }
1085 
1086  if ((numext::isnan)(a) || (numext::isnan)(x)) { // propagate nans
1087  return nan;
1088  }
1089 
1090  if ((x > one) && (x > a)) {
1092  if (mode == VALUE) {
1093  return one - ret;
1094  } else {
1095  return -ret;
1096  }
1097  }
1098 
1100  }
1101 };
1102 
1103 #endif // EIGEN_HAS_C99_MATH
1104 
1105 template <typename Scalar>
1107  typedef Scalar type;
1108 };
1109 
1110 template <typename Scalar>
1111 struct igamma_impl : igamma_generic_impl<Scalar, VALUE> {
1112  /* igam()
1113  * Incomplete gamma integral.
1114  *
1115  * The CDF of Gamma(a, 1) random variable at the point x.
1116  *
1117  * Accuracy estimation. For each a in [10^-2, 10^-1...10^3] we sample
1118  * 50 Gamma random variables x ~ Gamma(x | a, 1), a total of 300 points.
1119  * The ground truth is computed by mpmath. Mean absolute error:
1120  * float: 1.26713e-05
1121  * double: 2.33606e-12
1122  *
1123  * Cephes documentation below.
1124  *
1125  * SYNOPSIS:
1126  *
1127  * double a, x, y, igam();
1128  *
1129  * y = igam( a, x );
1130  *
1131  * DESCRIPTION:
1132  *
1133  * The function is defined by
1134  *
1135  * x
1136  * -
1137  * 1 | | -t a-1
1138  * igam(a,x) = ----- | e t dt.
1139  * - | |
1140  * | (a) -
1141  * 0
1142  *
1143  *
1144  * In this implementation both arguments must be positive.
1145  * The integral is evaluated by either a power series or
1146  * continued fraction expansion, depending on the relative
1147  * values of a and x.
1148  *
1149  * ACCURACY (double):
1150  *
1151  * Relative error:
1152  * arithmetic domain # trials peak rms
1153  * IEEE 0,30 200000 3.6e-14 2.9e-15
1154  * IEEE 0,100 300000 9.9e-14 1.5e-14
1155  *
1156  *
1157  * ACCURACY (float):
1158  *
1159  * Relative error:
1160  * arithmetic domain # trials peak rms
1161  * IEEE 0,30 20000 7.8e-6 5.9e-7
1162  *
1163  */
1164  /*
1165  Cephes Math Library Release 2.2: June, 1992
1166  Copyright 1985, 1987, 1992 by Stephen L. Moshier
1167  Direct inquiries to 30 Frost Street, Cambridge, MA 02140
1168  */
1169 
1170  /* left tail of incomplete gamma function:
1171  *
1172  * inf. k
1173  * a -x - x
1174  * x e > ----------
1175  * - -
1176  * k=0 | (a+k+1)
1177  *
1178  */
1179 };
1180 
1181 template <typename Scalar>
1183 
1184 template <typename Scalar>
1185 struct igamma_der_a_impl : igamma_generic_impl<Scalar, DERIVATIVE> {
1186  /* Derivative of the incomplete Gamma function with respect to a.
1187  *
1188  * Computes d/da igamma(a, x) by forward differentiation of the igamma code.
1189  *
1190  * Accuracy estimation. For each a in [10^-2, 10^-1...10^3] we sample
1191  * 50 Gamma random variables x ~ Gamma(x | a, 1), a total of 300 points.
1192  * The ground truth is computed by mpmath. Mean absolute error:
1193  * float: 6.17992e-07
1194  * double: 4.60453e-12
1195  *
1196  * Reference:
1197  * R. Moore. "Algorithm AS 187: Derivatives of the incomplete gamma
1198  * integral". Journal of the Royal Statistical Society. 1982
1199  */
1200 };
1201 
1202 template <typename Scalar>
1204 
1205 template <typename Scalar>
1207  : igamma_generic_impl<Scalar, SAMPLE_DERIVATIVE> {
1208  /* Derivative of a Gamma random variable sample with respect to alpha.
1209  *
1210  * Consider a sample of a Gamma random variable with the concentration
1211  * parameter alpha: sample ~ Gamma(alpha, 1). The reparameterization
1212  * derivative that we want to compute is dsample / dalpha =
1213  * d igammainv(alpha, u) / dalpha, where u = igamma(alpha, sample).
1214  * However, this formula is numerically unstable and expensive, so instead
1215  * we use implicit differentiation:
1216  *
1217  * igamma(alpha, sample) = u, where u ~ Uniform(0, 1).
1218  * Apply d / dalpha to both sides:
1219  * d igamma(alpha, sample) / dalpha
1220  * + d igamma(alpha, sample) / dsample * dsample/dalpha = 0
1221  * d igamma(alpha, sample) / dalpha
1222  * + Gamma(sample | alpha, 1) dsample / dalpha = 0
1223  * dsample/dalpha = - (d igamma(alpha, sample) / dalpha)
1224  * / Gamma(sample | alpha, 1)
1225  *
1226  * Here Gamma(sample | alpha, 1) is the PDF of the Gamma distribution
1227  * (note that the derivative of the CDF w.r.t. sample is the PDF).
1228  * See the reference below for more details.
1229  *
1230  * The derivative of igamma(alpha, sample) is computed by forward
1231  * differentiation of the igamma code. Division by the Gamma PDF is performed
1232  * in the same code, increasing the accuracy and speed due to cancellation
1233  * of some terms.
1234  *
1235  * Accuracy estimation. For each alpha in [10^-2, 10^-1...10^3] we sample
1236  * 50 Gamma random variables sample ~ Gamma(sample | alpha, 1), a total of 300
1237  * points. The ground truth is computed by mpmath. Mean absolute error:
1238  * float: 2.1686e-06
1239  * double: 1.4774e-12
1240  *
1241  * Reference:
1242  * M. Figurnov, S. Mohamed, A. Mnih "Implicit Reparameterization Gradients".
1243  * 2018
1244  */
1245 };
1246 
1247 /*****************************************************************************
1248  * Implementation of Riemann zeta function of two arguments, based on Cephes *
1249  *****************************************************************************/
1250 
1251 template <typename Scalar>
1252 struct zeta_retval {
1253  typedef Scalar type;
1254 };
1256 template <typename Scalar>
1261  THIS_TYPE_IS_NOT_SUPPORTED);
1262  return Scalar(0);
1263  }
1264 };
1265 
1266 template <>
1269  static EIGEN_STRONG_INLINE bool run(float& a, float& b, float& s, const float x, const float machep) {
1270  int i = 0;
1271  while(i < 9)
1272  {
1273  i += 1;
1274  a += 1.0f;
1275  b = numext::pow( a, -x );
1276  s += b;
1277  if( numext::abs(b/s) < machep )
1278  return true;
1279  }
1280 
1281  //Return whether we are done
1282  return false;
1283  }
1284 };
1285 
1286 template <>
1287 struct zeta_impl_series<double> {
1289  static EIGEN_STRONG_INLINE bool run(double& a, double& b, double& s, const double x, const double machep) {
1290  int i = 0;
1291  while( (i < 9) || (a <= 9.0) )
1292  {
1293  i += 1;
1294  a += 1.0;
1295  b = numext::pow( a, -x );
1296  s += b;
1297  if( numext::abs(b/s) < machep )
1298  return true;
1299  }
1300 
1301  //Return whether we are done
1302  return false;
1303  }
1304 };
1305 
1306 template <typename Scalar>
1307 struct zeta_impl {
1309  static Scalar run(Scalar x, Scalar q) {
1310  /* zeta.c
1311  *
1312  * Riemann zeta function of two arguments
1313  *
1314  *
1315  *
1316  * SYNOPSIS:
1317  *
1318  * double x, q, y, zeta();
1319  *
1320  * y = zeta( x, q );
1321  *
1322  *
1323  *
1324  * DESCRIPTION:
1325  *
1326  *
1327  *
1328  * inf.
1329  * - -x
1330  * zeta(x,q) = > (k+q)
1331  * -
1332  * k=0
1333  *
1334  * where x > 1 and q is not a negative integer or zero.
1335  * The Euler-Maclaurin summation formula is used to obtain
1336  * the expansion
1337  *
1338  * n
1339  * - -x
1340  * zeta(x,q) = > (k+q)
1341  * -
1342  * k=1
1343  *
1344  * 1-x inf. B x(x+1)...(x+2j)
1345  * (n+q) 1 - 2j
1346  * + --------- - ------- + > --------------------
1347  * x-1 x - x+2j+1
1348  * 2(n+q) j=1 (2j)! (n+q)
1349  *
1350  * where the B2j are Bernoulli numbers. Note that (see zetac.c)
1351  * zeta(x,1) = zetac(x) + 1.
1352  *
1353  *
1354  *
1355  * ACCURACY:
1356  *
1357  * Relative error for single precision:
1358  * arithmetic domain # trials peak rms
1359  * IEEE 0,25 10000 6.9e-7 1.0e-7
1360  *
1361  * Large arguments may produce underflow in powf(), in which
1362  * case the results are inaccurate.
1363  *
1364  * REFERENCE:
1365  *
1366  * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals,
1367  * Series, and Products, p. 1073; Academic Press, 1980.
1368  *
1369  */
1370 
1371  int i;
1372  Scalar p, r, a, b, k, s, t, w;
1373 
1374  const Scalar A[] = {
1375  Scalar(12.0),
1376  Scalar(-720.0),
1377  Scalar(30240.0),
1378  Scalar(-1209600.0),
1379  Scalar(47900160.0),
1380  Scalar(-1.8924375803183791606e9), /*1.307674368e12/691*/
1381  Scalar(7.47242496e10),
1382  Scalar(-2.950130727918164224e12), /*1.067062284288e16/3617*/
1383  Scalar(1.1646782814350067249e14), /*5.109094217170944e18/43867*/
1384  Scalar(-4.5979787224074726105e15), /*8.028576626982912e20/174611*/
1385  Scalar(1.8152105401943546773e17), /*1.5511210043330985984e23/854513*/
1386  Scalar(-7.1661652561756670113e18) /*1.6938241367317436694528e27/236364091*/
1387  };
1388 
1389  const Scalar maxnum = NumTraits<Scalar>::infinity();
1390  const Scalar zero = 0.0, half = 0.5, one = 1.0;
1391  const Scalar machep = cephes_helper<Scalar>::machep();
1392  const Scalar nan = NumTraits<Scalar>::quiet_NaN();
1393 
1394  if( x == one )
1395  return maxnum;
1396 
1397  if( x < one )
1398  {
1399  return nan;
1400  }
1401 
1402  if( q <= zero )
1403  {
1404  if(q == numext::floor(q))
1405  {
1406  if (x == numext::floor(x) && long(x) % 2 == 0) {
1407  return maxnum;
1408  }
1409  else {
1410  return nan;
1411  }
1412  }
1413  p = x;
1414  r = numext::floor(p);
1415  if (p != r)
1416  return nan;
1417  }
1418 
1419  /* Permit negative q but continue sum until n+q > +9 .
1420  * This case should be handled by a reflection formula.
1421  * If q<0 and x is an integer, there is a relation to
1422  * the polygamma function.
1423  */
1424  s = numext::pow( q, -x );
1425  a = q;
1426  b = zero;
1427  // Run the summation in a helper function that is specific to the floating precision
1428  if (zeta_impl_series<Scalar>::run(a, b, s, x, machep)) {
1429  return s;
1430  }
1431 
1432  w = a;
1433  s += b*w/(x-one);
1434  s -= half * b;
1435  a = one;
1436  k = zero;
1437  for( i=0; i<12; i++ )
1438  {
1439  a *= x + k;
1440  b /= w;
1441  t = a*b/A[i];
1442  s = s + t;
1443  t = numext::abs(t/s);
1444  if( t < machep ) {
1445  break;
1446  }
1447  k += one;
1448  a *= x + k;
1449  b /= w;
1450  k += one;
1451  }
1452  return s;
1453  }
1454 };
1455 
1456 /****************************************************************************
1457  * Implementation of polygamma function, requires C++11/C99 *
1458  ****************************************************************************/
1459 
1460 template <typename Scalar>
1462  typedef Scalar type;
1463 };
1465 #if !EIGEN_HAS_C99_MATH
1466 
1467 template <typename Scalar>
1472  THIS_TYPE_IS_NOT_SUPPORTED);
1473  return Scalar(0);
1474  }
1475 };
1476 
1477 #else
1478 
1479 template <typename Scalar>
1480 struct polygamma_impl {
1482  static Scalar run(Scalar n, Scalar x) {
1483  Scalar zero = 0.0, one = 1.0;
1484  Scalar nplus = n + one;
1485  const Scalar nan = NumTraits<Scalar>::quiet_NaN();
1486 
1487  // Check that n is a non-negative integer
1488  if (numext::floor(n) != n || n < zero) {
1489  return nan;
1490  }
1491  // Just return the digamma function for n = 0
1492  else if (n == zero) {
1493  return digamma_impl<Scalar>::run(x);
1494  }
1495  // Use the same implementation as scipy
1496  else {
1497  Scalar factorial = numext::exp(lgamma_impl<Scalar>::run(nplus));
1498  return numext::pow(-one, nplus) * factorial * zeta_impl<Scalar>::run(nplus, x);
1499  }
1500  }
1501 };
1502 
1503 #endif // EIGEN_HAS_C99_MATH
1504 
1505 /************************************************************************************************
1506  * Implementation of betainc (incomplete beta integral), based on Cephes but requires C++11/C99 *
1507  ************************************************************************************************/
1508 
1509 template <typename Scalar>
1511  typedef Scalar type;
1512 };
1514 #if !EIGEN_HAS_C99_MATH
1515 
1516 template <typename Scalar>
1521  THIS_TYPE_IS_NOT_SUPPORTED);
1522  return Scalar(0);
1523  }
1524 };
1525 
1526 #else
1527 
1528 template <typename Scalar>
1529 struct betainc_impl {
1532  /* betaincf.c
1533  *
1534  * Incomplete beta integral
1535  *
1536  *
1537  * SYNOPSIS:
1538  *
1539  * float a, b, x, y, betaincf();
1540  *
1541  * y = betaincf( a, b, x );
1542  *
1543  *
1544  * DESCRIPTION:
1545  *
1546  * Returns incomplete beta integral of the arguments, evaluated
1547  * from zero to x. The function is defined as
1548  *
1549  * x
1550  * - -
1551  * | (a+b) | | a-1 b-1
1552  * ----------- | t (1-t) dt.
1553  * - - | |
1554  * | (a) | (b) -
1555  * 0
1556  *
1557  * The domain of definition is 0 <= x <= 1. In this
1558  * implementation a and b are restricted to positive values.
1559  * The integral from x to 1 may be obtained by the symmetry
1560  * relation
1561  *
1562  * 1 - betainc( a, b, x ) = betainc( b, a, 1-x ).
1563  *
1564  * The integral is evaluated by a continued fraction expansion.
1565  * If a < 1, the function calls itself recursively after a
1566  * transformation to increase a to a+1.
1567  *
1568  * ACCURACY (float):
1569  *
1570  * Tested at random points (a,b,x) with a and b in the indicated
1571  * interval and x between 0 and 1.
1572  *
1573  * arithmetic domain # trials peak rms
1574  * Relative error:
1575  * IEEE 0,30 10000 3.7e-5 5.1e-6
1576  * IEEE 0,100 10000 1.7e-4 2.5e-5
1577  * The useful domain for relative error is limited by underflow
1578  * of the single precision exponential function.
1579  * Absolute error:
1580  * IEEE 0,30 100000 2.2e-5 9.6e-7
1581  * IEEE 0,100 10000 6.5e-5 3.7e-6
1582  *
1583  * Larger errors may occur for extreme ratios of a and b.
1584  *
1585  * ACCURACY (double):
1586  * arithmetic domain # trials peak rms
1587  * IEEE 0,5 10000 6.9e-15 4.5e-16
1588  * IEEE 0,85 250000 2.2e-13 1.7e-14
1589  * IEEE 0,1000 30000 5.3e-12 6.3e-13
1590  * IEEE 0,10000 250000 9.3e-11 7.1e-12
1591  * IEEE 0,100000 10000 8.7e-10 4.8e-11
1592  * Outputs smaller than the IEEE gradual underflow threshold
1593  * were excluded from these statistics.
1594  *
1595  * ERROR MESSAGES:
1596  * message condition value returned
1597  * incbet domain x<0, x>1 nan
1598  * incbet underflow nan
1599  */
1600 
1602  THIS_TYPE_IS_NOT_SUPPORTED);
1603  return Scalar(0);
1604  }
1605 };
1606 
1607 /* Continued fraction expansion #1 for incomplete beta integral (small_branch = True)
1608  * Continued fraction expansion #2 for incomplete beta integral (small_branch = False)
1609  */
1610 template <typename Scalar>
1611 struct incbeta_cfe {
1613  static EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x, bool small_branch) {
1616  THIS_TYPE_IS_NOT_SUPPORTED);
1618  const Scalar machep = cephes_helper<Scalar>::machep();
1620 
1621  const Scalar zero = 0;
1622  const Scalar one = 1;
1623  const Scalar two = 2;
1624 
1625  Scalar xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
1626  Scalar k1, k2, k3, k4, k5, k6, k7, k8, k26update;
1627  Scalar ans;
1628  int n;
1629 
1630  const int num_iters = (internal::is_same<Scalar, float>::value) ? 100 : 300;
1631  const Scalar thresh =
1632  (internal::is_same<Scalar, float>::value) ? machep : Scalar(3) * machep;
1634 
1635  if (small_branch) {
1636  k1 = a;
1637  k2 = a + b;
1638  k3 = a;
1639  k4 = a + one;
1640  k5 = one;
1641  k6 = b - one;
1642  k7 = k4;
1643  k8 = a + two;
1644  k26update = one;
1645  } else {
1646  k1 = a;
1647  k2 = b - one;
1648  k3 = a;
1649  k4 = a + one;
1650  k5 = one;
1651  k6 = a + b;
1652  k7 = a + one;
1653  k8 = a + two;
1654  k26update = -one;
1655  x = x / (one - x);
1656  }
1657 
1658  pkm2 = zero;
1659  qkm2 = one;
1660  pkm1 = one;
1661  qkm1 = one;
1662  ans = one;
1663  n = 0;
1664 
1665  do {
1666  xk = -(x * k1 * k2) / (k3 * k4);
1667  pk = pkm1 + pkm2 * xk;
1668  qk = qkm1 + qkm2 * xk;
1669  pkm2 = pkm1;
1670  pkm1 = pk;
1671  qkm2 = qkm1;
1672  qkm1 = qk;
1673 
1674  xk = (x * k5 * k6) / (k7 * k8);
1675  pk = pkm1 + pkm2 * xk;
1676  qk = qkm1 + qkm2 * xk;
1677  pkm2 = pkm1;
1678  pkm1 = pk;
1679  qkm2 = qkm1;
1680  qkm1 = qk;
1681 
1682  if (qk != zero) {
1683  r = pk / qk;
1684  if (numext::abs(ans - r) < numext::abs(r) * thresh) {
1685  return r;
1686  }
1687  ans = r;
1688  }
1689 
1690  k1 += one;
1691  k2 += k26update;
1692  k3 += two;
1693  k4 += two;
1694  k5 += one;
1695  k6 -= k26update;
1696  k7 += two;
1697  k8 += two;
1698 
1699  if ((numext::abs(qk) + numext::abs(pk)) > big) {
1700  pkm2 *= biginv;
1701  pkm1 *= biginv;
1702  qkm2 *= biginv;
1703  qkm1 *= biginv;
1704  }
1705  if ((numext::abs(qk) < biginv) || (numext::abs(pk) < biginv)) {
1706  pkm2 *= big;
1707  pkm1 *= big;
1708  qkm2 *= big;
1709  qkm1 *= big;
1710  }
1711  } while (++n < num_iters);
1712 
1713  return ans;
1714  }
1715 };
1716 
1717 /* Helper functions depending on the Scalar type */
1718 template <typename Scalar>
1719 struct betainc_helper {};
1720 
1721 template <>
1722 struct betainc_helper<float> {
1723  /* Core implementation, assumes a large (> 1.0) */
1724  EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE float incbsa(float aa, float bb,
1725  float xx) {
1726  float ans, a, b, t, x, onemx;
1727  bool reversed_a_b = false;
1728 
1729  onemx = 1.0f - xx;
1730 
1731  /* see if x is greater than the mean */
1732  if (xx > (aa / (aa + bb))) {
1733  reversed_a_b = true;
1734  a = bb;
1735  b = aa;
1736  t = xx;
1737  x = onemx;
1738  } else {
1739  a = aa;
1740  b = bb;
1741  t = onemx;
1742  x = xx;
1743  }
1744 
1745  /* Choose expansion for optimal convergence */
1746  if (b > 10.0f) {
1747  if (numext::abs(b * x / a) < 0.3f) {
1748  t = betainc_helper<float>::incbps(a, b, x);
1749  if (reversed_a_b) t = 1.0f - t;
1750  return t;
1751  }
1752  }
1753 
1754  ans = x * (a + b - 2.0f) / (a - 1.0f);
1755  if (ans < 1.0f) {
1756  ans = incbeta_cfe<float>::run(a, b, x, true /* small_branch */);
1757  t = b * numext::log(t);
1758  } else {
1759  ans = incbeta_cfe<float>::run(a, b, x, false /* small_branch */);
1760  t = (b - 1.0f) * numext::log(t);
1761  }
1762 
1763  t += a * numext::log(x) + lgamma_impl<float>::run(a + b) -
1765  t += numext::log(ans / a);
1766  t = numext::exp(t);
1767 
1768  if (reversed_a_b) t = 1.0f - t;
1769  return t;
1770  }
1771 
1773  static EIGEN_STRONG_INLINE float incbps(float a, float b, float x) {
1774  float t, u, y, s;
1775  const float machep = cephes_helper<float>::machep();
1776 
1777  y = a * numext::log(x) + (b - 1.0f) * numext::log1p(-x) - numext::log(a);
1779  y += lgamma_impl<float>::run(a + b);
1780 
1781  t = x / (1.0f - x);
1782  s = 0.0f;
1783  u = 1.0f;
1784  do {
1785  b -= 1.0f;
1786  if (b == 0.0f) {
1787  break;
1788  }
1789  a += 1.0f;
1790  u *= t * b / a;
1791  s += u;
1792  } while (numext::abs(u) > machep);
1793 
1794  return numext::exp(y) * (1.0f + s);
1795  }
1796 };
1797 
1798 template <>
1799 struct betainc_impl<float> {
1801  static float run(float a, float b, float x) {
1802  const float nan = NumTraits<float>::quiet_NaN();
1803  float ans, t;
1804 
1805  if (a <= 0.0f) return nan;
1806  if (b <= 0.0f) return nan;
1807  if ((x <= 0.0f) || (x >= 1.0f)) {
1808  if (x == 0.0f) return 0.0f;
1809  if (x == 1.0f) return 1.0f;
1810  // mtherr("betaincf", DOMAIN);
1811  return nan;
1812  }
1813 
1814  /* transformation for small aa */
1815  if (a <= 1.0f) {
1816  ans = betainc_helper<float>::incbsa(a + 1.0f, b, x);
1817  t = a * numext::log(x) + b * numext::log1p(-x) +
1820  return (ans + numext::exp(t));
1821  } else {
1822  return betainc_helper<float>::incbsa(a, b, x);
1823  }
1824  }
1825 };
1826 
1827 template <>
1828 struct betainc_helper<double> {
1830  static EIGEN_STRONG_INLINE double incbps(double a, double b, double x) {
1831  const double machep = cephes_helper<double>::machep();
1832 
1833  double s, t, u, v, n, t1, z, ai;
1834 
1835  ai = 1.0 / a;
1836  u = (1.0 - b) * x;
1837  v = u / (a + 1.0);
1838  t1 = v;
1839  t = u;
1840  n = 2.0;
1841  s = 0.0;
1842  z = machep * ai;
1843  while (numext::abs(v) > z) {
1844  u = (n - b) * x / n;
1845  t *= u;
1846  v = t / (a + n);
1847  s += v;
1848  n += 1.0;
1849  }
1850  s += t1;
1851  s += ai;
1852 
1853  u = a * numext::log(x);
1854  // TODO: gamma() is not directly implemented in Eigen.
1855  /*
1856  if ((a + b) < maxgam && numext::abs(u) < maxlog) {
1857  t = gamma(a + b) / (gamma(a) * gamma(b));
1858  s = s * t * pow(x, a);
1859  }
1860  */
1863  return s = numext::exp(t);
1864  }
1865 };
1866 
1867 template <>
1868 struct betainc_impl<double> {
1870  static double run(double aa, double bb, double xx) {
1871  const double nan = NumTraits<double>::quiet_NaN();
1872  const double machep = cephes_helper<double>::machep();
1873  // const double maxgam = 171.624376956302725;
1874 
1875  double a, b, t, x, xc, w, y;
1876  bool reversed_a_b = false;
1877 
1878  if (aa <= 0.0 || bb <= 0.0) {
1879  return nan; // goto domerr;
1880  }
1881 
1882  if ((xx <= 0.0) || (xx >= 1.0)) {
1883  if (xx == 0.0) return (0.0);
1884  if (xx == 1.0) return (1.0);
1885  // mtherr("incbet", DOMAIN);
1886  return nan;
1887  }
1888 
1889  if ((bb * xx) <= 1.0 && xx <= 0.95) {
1890  return betainc_helper<double>::incbps(aa, bb, xx);
1891  }
1892 
1893  w = 1.0 - xx;
1894 
1895  /* Reverse a and b if x is greater than the mean. */
1896  if (xx > (aa / (aa + bb))) {
1897  reversed_a_b = true;
1898  a = bb;
1899  b = aa;
1900  xc = xx;
1901  x = w;
1902  } else {
1903  a = aa;
1904  b = bb;
1905  xc = w;
1906  x = xx;
1907  }
1908 
1909  if (reversed_a_b && (b * x) <= 1.0 && x <= 0.95) {
1910  t = betainc_helper<double>::incbps(a, b, x);
1911  if (t <= machep) {
1912  t = 1.0 - machep;
1913  } else {
1914  t = 1.0 - t;
1915  }
1916  return t;
1917  }
1918 
1919  /* Choose expansion for better convergence. */
1920  y = x * (a + b - 2.0) - (a - 1.0);
1921  if (y < 0.0) {
1922  w = incbeta_cfe<double>::run(a, b, x, true /* small_branch */);
1923  } else {
1924  w = incbeta_cfe<double>::run(a, b, x, false /* small_branch */) / xc;
1925  }
1926 
1927  /* Multiply w by the factor
1928  a b _ _ _
1929  x (1-x) | (a+b) / ( a | (a) | (b) ) . */
1930 
1931  y = a * numext::log(x);
1932  t = b * numext::log(xc);
1933  // TODO: gamma is not directly implemented in Eigen.
1934  /*
1935  if ((a + b) < maxgam && numext::abs(y) < maxlog && numext::abs(t) < maxlog)
1936  {
1937  t = pow(xc, b);
1938  t *= pow(x, a);
1939  t /= a;
1940  t *= w;
1941  t *= gamma(a + b) / (gamma(a) * gamma(b));
1942  } else {
1943  */
1944  /* Resort to logarithms. */
1947  y += numext::log(w / a);
1948  t = numext::exp(y);
1949 
1950  /* } */
1951  // done:
1952 
1953  if (reversed_a_b) {
1954  if (t <= machep) {
1955  t = 1.0 - machep;
1956  } else {
1957  t = 1.0 - t;
1958  }
1959  }
1960  return t;
1961  }
1962 };
1963 
1964 #endif // EIGEN_HAS_C99_MATH
1965 
1966 } // end namespace internal
1967 
1968 namespace numext {
1969 
1970 template <typename Scalar>
1972  lgamma(const Scalar& x) {
1973  return EIGEN_MATHFUNC_IMPL(lgamma, Scalar)::run(x);
1974 }
1975 
1976 template <typename Scalar>
1978  digamma(const Scalar& x) {
1979  return EIGEN_MATHFUNC_IMPL(digamma, Scalar)::run(x);
1980 }
1981 
1982 template <typename Scalar>
1984 zeta(const Scalar& x, const Scalar& q) {
1985  return EIGEN_MATHFUNC_IMPL(zeta, Scalar)::run(x, q);
1986 }
1987 
1988 template <typename Scalar>
1990 polygamma(const Scalar& n, const Scalar& x) {
1991  return EIGEN_MATHFUNC_IMPL(polygamma, Scalar)::run(n, x);
1992 }
1993 
1994 template <typename Scalar>
1996  erf(const Scalar& x) {
1997  return EIGEN_MATHFUNC_IMPL(erf, Scalar)::run(x);
1998 }
1999 
2000 template <typename Scalar>
2002  erfc(const Scalar& x) {
2003  return EIGEN_MATHFUNC_IMPL(erfc, Scalar)::run(x);
2004 }
2005 
2006 template <typename Scalar>
2008  ndtri(const Scalar& x) {
2009  return EIGEN_MATHFUNC_IMPL(ndtri, Scalar)::run(x);
2010 }
2011 
2012 template <typename Scalar>
2014  igamma(const Scalar& a, const Scalar& x) {
2015  return EIGEN_MATHFUNC_IMPL(igamma, Scalar)::run(a, x);
2016 }
2017 
2018 template <typename Scalar>
2020  igamma_der_a(const Scalar& a, const Scalar& x) {
2022 }
2023 
2024 template <typename Scalar>
2026  gamma_sample_der_alpha(const Scalar& a, const Scalar& x) {
2028 }
2029 
2030 template <typename Scalar>
2032  igammac(const Scalar& a, const Scalar& x) {
2033  return EIGEN_MATHFUNC_IMPL(igammac, Scalar)::run(a, x);
2034 }
2035 
2036 template <typename Scalar>
2038  betainc(const Scalar& a, const Scalar& b, const Scalar& x) {
2039  return EIGEN_MATHFUNC_IMPL(betainc, Scalar)::run(a, b, x);
2040 }
2041 
2042 } // end namespace numext
2043 } // end namespace Eigen
2044 
2045 #endif // EIGEN_SPECIAL_FUNCTIONS_H
Eigen::numext::x
EIGEN_DEVICE_FUNC const Scalar & x
Definition: SpecialFunctionsImpl.h:1990
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
Eigen::internal::igamma_num_iterations
EIGEN_DEVICE_FUNC int igamma_num_iterations()
Definition: SpecialFunctionsImpl.h:746
Eigen::internal::cephes_helper< double >::biginv
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double biginv()
Definition: SpecialFunctionsImpl.h:723
Eigen::internal::gamma_sample_der_alpha_impl
Definition: SpecialFunctionsImpl.h:1206
EIGEN_DEVICE_FUNC
#define EIGEN_DEVICE_FUNC
Definition: Macros.h:976
Eigen::internal::gamma_sample_der_alpha_retval
Definition: SpecialFunctionsImpl.h:1203
Eigen::internal::ppolevl::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet run(const Packet &x, const typename unpacket_traits< Packet >::type coeff[])
Definition: GenericPacketMathFunctions.h:1561
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::internal::psqrt
EIGEN_STRONG_INLINE Packet4f psqrt(const Packet4f &a)
Definition: MSA/PacketMath.h:723
Eigen::igamma_der_a
const EIGEN_STRONG_INLINE Eigen::CwiseBinaryOp< Eigen::internal::scalar_igamma_der_a_op< typename Derived::Scalar >, const Derived, const ExponentDerived > igamma_der_a(const Eigen::ArrayBase< Derived > &a, const Eigen::ArrayBase< ExponentDerived > &x)
Definition: SpecialFunctionsArrayAPI.h:51
Eigen::internal::igammac_impl::run
static EIGEN_DEVICE_FUNC Scalar run(Scalar a, Scalar x)
Definition: SpecialFunctionsImpl.h:956
Eigen::betainc
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE 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
EIGEN_PI
#define EIGEN_PI
Definition: Eigen/src/Core/MathFunctions.h:16
Eigen::internal::cephes_helper< double >::big
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double big()
Definition: SpecialFunctionsImpl.h:719
Eigen::internal::is_same::value
@ value
Definition: Meta.h:148
s
RealScalar s
Definition: level1_cplx_impl.h:126
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Eigen::internal::digamma_retval
Definition: SpecialFunctionsImpl.h:110
Eigen::internal::igamma_generic_impl
Definition: SpecialFunctionsImpl.h:1053
Eigen::internal::digamma_impl_maybe_poly::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(const Scalar)
Definition: SpecialFunctionsImpl.h:130
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Eigen::internal::pcmp_lt
EIGEN_STRONG_INLINE Packet4f pcmp_lt(const Packet4f &a, const Packet4f &b)
Definition: AltiVec/PacketMath.h:868
Eigen::internal::erf_retval
Definition: SpecialFunctionsImpl.h:357
Eigen::numext::isinf
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool() isinf(const Eigen::bfloat16 &h)
Definition: BFloat16.h:665
Eigen::gamma_sample_der_alpha
const EIGEN_STRONG_INLINE Eigen::CwiseBinaryOp< Eigen::internal::scalar_gamma_sample_der_alpha_op< typename AlphaDerived::Scalar >, const AlphaDerived, const SampleDerived > gamma_sample_der_alpha(const Eigen::ArrayBase< AlphaDerived > &alpha, const Eigen::ArrayBase< SampleDerived > &sample)
Definition: SpecialFunctionsArrayAPI.h:72
ret
DenseIndex ret
Definition: level1_cplx_impl.h:44
Eigen::internal::betainc_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:1513
Eigen::internal::cephes_helper< float >::big
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float big()
Definition: SpecialFunctionsImpl.h:701
Eigen::internal::cephes_helper< float >::machep
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float machep()
Definition: SpecialFunctionsImpl.h:697
Eigen::internal::generic_ndtri_gt_exp_neg_two
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_ndtri_gt_exp_neg_two(const T &b)
Definition: SpecialFunctionsImpl.h:516
Eigen::internal::digamma_impl::run
static EIGEN_DEVICE_FUNC Scalar run(Scalar x)
Definition: SpecialFunctionsImpl.h:183
Eigen::internal::lgamma_impl
Definition: SpecialFunctionsImpl.h:45
Eigen::internal::pdiv
EIGEN_DEVICE_FUNC Packet pdiv(const Packet &a, const Packet &b)
Definition: GenericPacketMath.h:244
Eigen::internal::IgammaComputationMode
IgammaComputationMode
Definition: SpecialFunctionsImpl.h:729
Eigen::internal::zeta_impl
Definition: SpecialFunctionsImpl.h:1307
Eigen::internal::generic_ndtri
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T generic_ndtri(const T &a)
Definition: SpecialFunctionsImpl.h:621
Eigen::internal::erfc_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:403
pybind_wrapper_test_script.dummy
dummy
Definition: pybind_wrapper_test_script.py:42
Eigen::internal::igamma_impl
Definition: SpecialFunctionsImpl.h:1111
Eigen::internal::digamma_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:113
Eigen::numext::exp
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T exp(const T &x)
Definition: Eigen/src/Core/MathFunctions.h:1548
digamma
const EIGEN_DEVICE_FUNC DigammaReturnType digamma() const
Definition: ArrayCwiseUnaryOps.h:637
Eigen::internal::igamma_der_a_retval
Definition: SpecialFunctionsImpl.h:1182
n
int n
Definition: BiCGSTAB_simple.cpp:1
epsilon
static double epsilon
Definition: testRot3.cpp:37
Eigen::internal::digamma_impl_maybe_poly< float >::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float run(const float s)
Definition: SpecialFunctionsImpl.h:141
Eigen::internal::zeta_impl_series
Definition: SpecialFunctionsImpl.h:1257
simple::p2
static Point3 p2
Definition: testInitializePose3.cpp:51
Eigen::internal::flipsign< float >
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float flipsign< float >(const float &should_flipsign, const float &x)
Definition: SpecialFunctionsImpl.h:506
Eigen::internal::VALUE
@ VALUE
Definition: SpecialFunctionsImpl.h:729
Eigen::internal::main_igamma_term
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar main_igamma_term(Scalar a, Scalar x)
Definition: SpecialFunctionsImpl.h:733
Eigen::internal::zeta_retval
Definition: SpecialFunctionsImpl.h:1252
A
Definition: test_numpy_dtypes.cpp:298
Eigen::internal::zeta_impl_series< double >::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool run(double &a, double &b, double &s, const double x, const double machep)
Definition: SpecialFunctionsImpl.h:1289
Eigen::internal::ndtri_impl::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(const Scalar)
Definition: SpecialFunctionsImpl.h:654
Eigen::polygamma
const EIGEN_STRONG_INLINE 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:112
lgamma
const EIGEN_DEVICE_FUNC LgammaReturnType lgamma() const
Definition: ArrayCwiseUnaryOps.h:620
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
x1
Pose3 x1
Definition: testPose3.cpp:663
Eigen::internal::betainc_impl::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar b, Scalar x)
Definition: SpecialFunctionsImpl.h:1519
Eigen::internal::lgamma_retval
Definition: SpecialFunctionsImpl.h:55
k1
double k1(double x)
Definition: k1.c:133
Eigen::internal::igamma_der_a_impl
Definition: SpecialFunctionsImpl.h:1185
Eigen::internal::generic_fast_erf_float
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_fast_erf_float(const T &a_x)
Definition: SpecialFunctionsImpl.h:304
Eigen::internal::pselect
EIGEN_STRONG_INLINE Packet4f pselect(const Packet4f &mask, const Packet4f &a, const Packet4f &b)
Definition: AltiVec/PacketMath.h:917
zero
EIGEN_DONT_INLINE Scalar zero()
Definition: svd_common.h:296
Eigen::internal::ndtri_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:646
EIGEN_STRONG_INLINE
#define EIGEN_STRONG_INLINE
Definition: Macros.h:917
log1p
double log1p(double x)
Definition: unity.c:49
x0
static Symbol x0('x', 0)
Eigen::internal::pmax
EIGEN_DEVICE_FUNC Packet pmax(const Packet &a, const Packet &b)
Definition: GenericPacketMath.h:524
Eigen::numext::isnan
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool() isnan(const Eigen::bfloat16 &h)
Definition: BFloat16.h:659
biginv
static double biginv
Definition: igam.c:120
gtsam.examples.DogLegOptimizerExample.run
def run(args)
Definition: DogLegOptimizerExample.py:21
Eigen::internal::cephes_helper::big
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar big()
Definition: SpecialFunctionsImpl.h:689
Eigen::internal::Packet
Definition: ZVector/PacketMath.h:47
Eigen::internal::psub
EIGEN_DEVICE_FUNC Packet psub(const Packet &a, const Packet &b)
Definition: GenericPacketMath.h:222
Eigen::internal::generic_ndtri_lt_exp_neg_two
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T generic_ndtri_lt_exp_neg_two(const T &b, const T &should_flipsign)
Definition: SpecialFunctionsImpl.h:549
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
Eigen::internal::lgamma_impl::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(const Scalar)
Definition: SpecialFunctionsImpl.h:49
EIGEN_ALWAYS_INLINE
#define EIGEN_ALWAYS_INLINE
Definition: Macros.h:932
Eigen::Triplet< double >
Eigen::internal::digamma_impl_maybe_poly
Definition: SpecialFunctionsImpl.h:128
Eigen::internal::zeta_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:1255
Eigen::internal::betainc_impl
Definition: SpecialFunctionsImpl.h:1517
Eigen::internal::erf_impl
Definition: SpecialFunctionsImpl.h:349
Eigen::internal::erf_impl::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T run(const T &x)
Definition: SpecialFunctionsImpl.h:351
Eigen::internal::igamma_generic_impl::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(Scalar a, Scalar x)
Definition: SpecialFunctionsImpl.h:1057
p1
Vector3f p1
Definition: MatrixBase_all.cpp:2
Eigen::internal::zeta_impl_series::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(const Scalar)
Definition: SpecialFunctionsImpl.h:1259
Eigen::numext::pow
EIGEN_DEVICE_FUNC internal::pow_impl< ScalarX, ScalarY >::result_type pow(const ScalarX &x, const ScalarY &y)
Definition: Eigen/src/Core/MathFunctions.h:1363
Eigen::internal::zeta_impl::run
static EIGEN_DEVICE_FUNC Scalar run(Scalar x, Scalar q)
Definition: SpecialFunctionsImpl.h:1309
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
Eigen::internal::digamma_impl
Definition: SpecialFunctionsImpl.h:181
Eigen::internal::flipsign< double >
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double flipsign< double >(const double &should_flipsign, const double &x)
Definition: SpecialFunctionsImpl.h:500
Eigen::internal::flipsign
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T flipsign(const T &should_flipsign, const T &x)
Definition: SpecialFunctionsImpl.h:491
Eigen::internal::cephes_helper< double >::machep
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double machep()
Definition: SpecialFunctionsImpl.h:715
Eigen::igamma
const EIGEN_STRONG_INLINE 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
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
Eigen::internal::pmul
EIGEN_DEVICE_FUNC Packet pmul(const Packet &a, const Packet &b)
Definition: GenericPacketMath.h:237
Eigen::internal::igamma_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:1107
Eigen::internal::polygamma_retval
Definition: SpecialFunctionsImpl.h:1461
Eigen::internal::ndtri_retval
Definition: SpecialFunctionsImpl.h:645
Eigen::internal::erfc_impl::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(const Scalar)
Definition: SpecialFunctionsImpl.h:396
Eigen::internal::DERIVATIVE
@ DERIVATIVE
Definition: SpecialFunctionsImpl.h:729
EIGEN_STATIC_ASSERT
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:127
Eigen::numext::EIGEN_MATHFUNC_RETVAL
EIGEN_DEVICE_FUNC EIGEN_MATHFUNC_RETVAL(real, Scalar) real(const Scalar &x)
Definition: Eigen/src/Core/MathFunctions.h:1232
Eigen::internal::igamma_retval
Definition: SpecialFunctionsImpl.h:1106
p0
Vector3f p0
Definition: MatrixBase_all.cpp:2
Eigen::internal::y
const Scalar & y
Definition: Eigen/src/Core/MathFunctions.h:821
Eigen::numext::floor
EIGEN_DEVICE_FUNC T() floor(const T &x)
Definition: Eigen/src/Core/MathFunctions.h:1402
Eigen::internal::igamma_series_impl::run
static EIGEN_DEVICE_FUNC Scalar run(Scalar a, Scalar x)
Definition: SpecialFunctionsImpl.h:891
big
static double big
Definition: igam.c:119
Eigen::internal::pmadd
EIGEN_STRONG_INLINE Packet4f pmadd(const Packet4f &a, const Packet4f &b, const Packet4f &c)
Definition: AltiVec/PacketMath.h:827
Eigen::internal::cephes_helper::machep
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar machep()
Definition: SpecialFunctionsImpl.h:687
Eigen::internal::plog
EIGEN_DECLARE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS Packet plog(const Packet &a)
Definition: GenericPacketMath.h:796
Eigen::internal::polygamma_impl::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar run(Scalar n, Scalar x)
Definition: SpecialFunctionsImpl.h:1470
p
float * p
Definition: Tutorial_Map_using.cpp:9
c2
static double c2
Definition: airy.c:55
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
Eigen::zeta
const EIGEN_STRONG_INLINE 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:156
Eigen::internal::pmin
EIGEN_DEVICE_FUNC Packet pmin(const Packet &a, const Packet &b)
Definition: GenericPacketMath.h:512
erf
double erf(double x)
Definition: ndtr.c:285
Eigen::internal::cephes_helper::biginv
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar biginv()
Definition: SpecialFunctionsImpl.h:691
gtsam.examples.DogLegOptimizerExample.float
float
Definition: DogLegOptimizerExample.py:113
mode
static const DiscreteKey mode(modeKey, 2)
Eigen::internal::is_same
Definition: Meta.h:148
ndtri
double ndtri(double y0)
Definition: ndtri.c:134
Eigen::numext::abs
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE internal::enable_if< NumTraits< T >::IsSigned||NumTraits< T >::IsComplex, typename NumTraits< T >::Real >::type abs(const T &x)
Definition: Eigen/src/Core/MathFunctions.h:1511
Eigen::internal::polygamma_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:1464
Eigen::internal::ppolevl
Definition: GenericPacketMathFunctions.h:1560
Eigen::internal::igammac_impl
Definition: SpecialFunctionsImpl.h:954
internal
Definition: BandTriangularSolver.h:13
Eigen::internal::betainc_retval
Definition: SpecialFunctionsImpl.h:1510
Eigen::internal::zeta_impl_series< float >::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool run(float &a, float &b, float &s, const float x, const float machep)
Definition: SpecialFunctionsImpl.h:1269
Eigen::internal::cephes_helper< float >::biginv
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float biginv()
Definition: SpecialFunctionsImpl.h:706
erfc
double erfc(double a)
Definition: ndtr.c:227
Eigen::internal::igamma_series_impl
Definition: SpecialFunctionsImpl.h:881
Eigen::internal::igammac_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:682
Eigen::internal::erf_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:358
Eigen::half
Definition: Half.h:142
align_3::t
Point2 t(10, 10)
Eigen::internal::erfc_impl
Definition: SpecialFunctionsImpl.h:392
Eigen::internal::ndtri_impl
Definition: SpecialFunctionsImpl.h:652
Eigen::internal::lgamma_retval::type
Scalar type
Definition: SpecialFunctionsImpl.h:56
Eigen::internal::polygamma_impl
Definition: SpecialFunctionsImpl.h:1468
Eigen::numext::log
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T log(const T &x)
Definition: Eigen/src/Core/MathFunctions.h:1491
Eigen::internal::SAMPLE_DERIVATIVE
@ SAMPLE_DERIVATIVE
Definition: SpecialFunctionsImpl.h:729
Eigen::internal::erfc_retval
Definition: SpecialFunctionsImpl.h:402
Eigen::internal::pcmp_le
EIGEN_STRONG_INLINE Packet4f pcmp_le(const Packet4f &a, const Packet4f &b)
Definition: AltiVec/PacketMath.h:867
Eigen::numext::tan
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T tan(const T &x)
Definition: Eigen/src/Core/MathFunctions.h:1640
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
x2
Pose3 x2(Rot3::Ypr(0.0, 0.0, 0.0), l2)
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::internal::igammac_cf_impl
Definition: SpecialFunctionsImpl.h:763
Eigen::internal::igammac_retval
Definition: SpecialFunctionsImpl.h:679
Eigen::numext::b
EIGEN_DEVICE_FUNC const Scalar & b
Definition: SpecialFunctionsImpl.h:2038
Eigen::internal::cephes_helper
Definition: SpecialFunctionsImpl.h:685
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
EIGEN_MATHFUNC_IMPL
#define EIGEN_MATHFUNC_IMPL(func, scalar)
Definition: Eigen/src/Core/MathFunctions.h:70
Eigen::internal::digamma_impl_maybe_poly< double >::run
static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE double run(const double s)
Definition: SpecialFunctionsImpl.h:160
Eigen::igammac
const EIGEN_STRONG_INLINE 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:90
Eigen::internal::igammac_cf_impl::run
static EIGEN_DEVICE_FUNC Scalar run(Scalar a, Scalar x)
Definition: SpecialFunctionsImpl.h:774


gtsam
Author(s):
autogenerated on Thu Jun 13 2024 03:06:31