Eigen::internal::igammac_cf_impl< Scalar, mode > Struct Template Reference

#include <SpecialFunctionsImpl.h>

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

Member Function Documentation

run()

template<typename Scalar , IgammaComputationMode mode>

static EIGEN_DEVICE_FUNC Scalar Eigen::internal::igammac_cf_impl< Scalar, mode >::run ( Scalar  a, Scalar  x  )

inlinestatic

                                                          {
    const Scalar zero = 0;
    const Scalar one = 1;
    const Scalar two = 2;
    const Scalar machep = cephes_helper<Scalar>::machep();
    const Scalar big = cephes_helper<Scalar>::big();
    const Scalar biginv = cephes_helper<Scalar>::biginv();
 
    if ((numext::isinf)(x)) {
      return zero;
    }
 
    Scalar ax = main_igamma_term<Scalar>(a, x);
    // This is independent of mode. If this value is zero,
    // then the function value is zero. If the function value is zero,
    // then we are in a neighborhood where the function value evaluates to zero,
    // so the derivative is zero.
    if (ax == zero) {
      return zero;
    }
 
    // continued fraction
    Scalar y = one - a;
    Scalar z = x + y + one;
    Scalar c = zero;
    Scalar pkm2 = one;
    Scalar qkm2 = x;
    Scalar pkm1 = x + one;
    Scalar qkm1 = z * x;
    Scalar ans = pkm1 / qkm1;
 
    Scalar dpkm2_da = zero;
    Scalar dqkm2_da = zero;
    Scalar dpkm1_da = zero;
    Scalar dqkm1_da = -x;
    Scalar dans_da = (dpkm1_da - ans * dqkm1_da) / qkm1;
 
    for (int i = 0; i < igamma_num_iterations<Scalar, mode>(); i++) {
      c += one;
      y += one;
      z += two;
 
      Scalar yc = y * c;
      Scalar pk = pkm1 * z - pkm2 * yc;
      Scalar qk = qkm1 * z - qkm2 * yc;
 
      Scalar dpk_da = dpkm1_da * z - pkm1 - dpkm2_da * yc + pkm2 * c;
      Scalar dqk_da = dqkm1_da * z - qkm1 - dqkm2_da * yc + qkm2 * c;
 
      if (qk != zero) {
        Scalar ans_prev = ans;
        ans = pk / qk;
 
        Scalar dans_da_prev = dans_da;
        dans_da = (dpk_da - ans * dqk_da) / qk;
 
        if (mode == VALUE) {
          if (numext::abs(ans_prev - ans) <= machep * numext::abs(ans)) {
            break;
          }
        } else {
          if (numext::abs(dans_da - dans_da_prev) <= machep) {
            break;
          }
        }
      }
 
      pkm2 = pkm1;
      pkm1 = pk;
      qkm2 = qkm1;
      qkm1 = qk;
 
      dpkm2_da = dpkm1_da;
      dpkm1_da = dpk_da;
      dqkm2_da = dqkm1_da;
      dqkm1_da = dqk_da;
 
      if (numext::abs(pk) > big) {
        pkm2 *= biginv;
        pkm1 *= biginv;
        qkm2 *= biginv;
        qkm1 *= biginv;
 
        dpkm2_da *= biginv;
        dpkm1_da *= biginv;
        dqkm2_da *= biginv;
        dqkm1_da *= biginv;
      }
    }
 
    /* Compute  x**a * exp(-x) / gamma(a)  */
    Scalar dlogax_da = numext::log(x) - digamma_impl<Scalar>::run(a);
    Scalar dax_da = ax * dlogax_da;
 
    switch (mode) {
      case VALUE:
        return ans * ax;
      case DERIVATIVE:
        return ans * dax_da + dans_da * ax;
      case SAMPLE_DERIVATIVE:
      default:  // this is needed to suppress clang warning
        return -(dans_da + ans * dlogax_da) * x;
    }
  }

References a, Eigen::numext::abs(), plotDoE::ax, Eigen::internal::cephes_helper< Scalar >::big(), Eigen::internal::cephes_helper< Scalar >::biginv(), calibrate::c, Eigen::internal::DERIVATIVE, i, Eigen::numext::isinf(), Eigen::numext::log(), Eigen::internal::cephes_helper< Scalar >::machep(), Eigen::internal::digamma_impl< Scalar >::run(), Eigen::internal::SAMPLE_DERIVATIVE, Eigen::internal::VALUE, Eigen::numext::x, Eigen::internal::y, and zero().

---

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

• SpecialFunctionsImpl.h