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


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:44:49