Eigen::internal::zeta_impl< Scalar > Struct Template Reference

#include <SpecialFunctionsImpl.h>

Static Public Member Functions
static EIGEN_DEVICE_FUNC Scalar	run (Scalar x, Scalar q)

Member Function Documentation

run()

template<typename Scalar >

static EIGEN_DEVICE_FUNC Scalar Eigen::internal::zeta_impl< Scalar >::run ( Scalar  x, Scalar  q  )

inlinestatic

                                                          {
    /*                          zeta.c
     *
     *  Riemann zeta function of two arguments
     *
     *
     *
     * SYNOPSIS:
     *
     * double x, q, y, zeta();
     *
     * y = zeta( x, q );
     *
     *
     *
     * DESCRIPTION:
     *
     *
     *
     *                 inf.
     *                  -        -x
     *   zeta(x,q)  =   >   (k+q)
     *                  -
     *                 k=0
     *
     * where x > 1 and q is not a negative integer or zero.
     * The Euler-Maclaurin summation formula is used to obtain
     * the expansion
     *
     *                n
     *                -       -x
     * zeta(x,q)  =   >  (k+q)
     *                -
     *               k=1
     *
     *           1-x                 inf.  B   x(x+1)...(x+2j)
     *      (n+q)           1         -     2j
     *  +  ---------  -  -------  +   >    --------------------
     *        x-1              x      -                   x+2j+1
     *                   2(n+q)      j=1       (2j)! (n+q)
     *
     * where the B2j are Bernoulli numbers.  Note that (see zetac.c)
     * zeta(x,1) = zetac(x) + 1.
     *
     *
     *
     * ACCURACY:
     *
     * Relative error for single precision:
     * arithmetic   domain     # trials      peak         rms
     *    IEEE      0,25        10000       6.9e-7      1.0e-7
     *
     * Large arguments may produce underflow in powf(), in which
     * case the results are inaccurate.
     *
     * REFERENCE:
     *
     * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals,
     * Series, and Products, p. 1073; Academic Press, 1980.
     *
     */
 
    int i;
    Scalar p, r, a, b, k, s, t, w;
 
    const Scalar A[] = {
        Scalar(12.0),
        Scalar(-720.0),
        Scalar(30240.0),
        Scalar(-1209600.0),
        Scalar(47900160.0),
        Scalar(-1.8924375803183791606e9), /*1.307674368e12/691*/
        Scalar(7.47242496e10),
        Scalar(-2.950130727918164224e12),  /*1.067062284288e16/3617*/
        Scalar(1.1646782814350067249e14),  /*5.109094217170944e18/43867*/
        Scalar(-4.5979787224074726105e15), /*8.028576626982912e20/174611*/
        Scalar(1.8152105401943546773e17),  /*1.5511210043330985984e23/854513*/
        Scalar(-7.1661652561756670113e18)  /*1.6938241367317436694528e27/236364091*/
    };
 
    const Scalar maxnum = NumTraits<Scalar>::infinity();
    const Scalar zero = Scalar(0.0), half = Scalar(0.5), one = Scalar(1.0);
    const Scalar machep = cephes_helper<Scalar>::machep();
    const Scalar nan = NumTraits<Scalar>::quiet_NaN();
 
    if (x == one) return maxnum;
 
    if (x < one) {
      return nan;
    }
 
    if (q <= zero) {
      if (q == numext::floor(q)) {
        if (x == numext::floor(x) && long(x) % 2 == 0) {
          return maxnum;
        } else {
          return nan;
        }
      }
      p = x;
      r = numext::floor(p);
      if (p != r) return nan;
    }
 
    /* Permit negative q but continue sum until n+q > +9 .
     * This case should be handled by a reflection formula.
     * If q<0 and x is an integer, there is a relation to
     * the polygamma function.
     */
    s = numext::pow(q, -x);
    a = q;
    b = zero;
    // Run the summation in a helper function that is specific to the floating precision
    if (zeta_impl_series<Scalar>::run(a, b, s, x, machep)) {
      return s;
    }
 
    // If b is zero, then the tail sum will also end up being zero.
    // Exiting early here can prevent NaNs for some large inputs, where
    // the tail sum computed below has term `a` which can overflow to `inf`.
    if (numext::equal_strict(b, zero)) {
      return s;
    }
 
    w = a;
    s += b * w / (x - one);
    s -= half * b;
    a = one;
    k = zero;
 
    for (i = 0; i < 12; i++) {
      a *= x + k;
      b /= w;
      t = a * b / A[i];
      s = s + t;
      t = numext::abs(t / s);
      if (t < machep) {
        break;
      }
      k += one;
      a *= x + k;
      b /= w;
      k += one;
    }
    return s;
  }

References a, Eigen::numext::abs(), Eigen::numext::b, Eigen::numext::equal_strict(), Eigen::numext::floor(), i, k, Eigen::internal::cephes_helper< Scalar >::machep(), p, Eigen::numext::pow(), Eigen::numext::q, UniformPSDSelfTest::r, s, plotPSD::t, w, Eigen::numext::x, and zero().

---

The documentation for this struct was generated from the following file:

• SpecialFunctionsImpl.h