_igamc_helper_continued_fraction Class — pytorch Architecture

Source Code

aten/src/ATen/native/Math.h lines 1002–1061

template <typename scalar_t>
static scalar_t _igamc_helper_continued_fraction(scalar_t a, scalar_t x) {
  // Compute igamc using DLMF 8.9.2. [igam1]
  int i;
  scalar_t ans, ax, c, yc, r, t, y, z;
  scalar_t pk, pkm1, pkm2, qk, qkm1, qkm2;
  int MAXITER = 2000;
  static scalar_t MACHEP = std::is_same_v<scalar_t, double> ?
    1.11022302462515654042E-16 : 5.9604644775390625E-8;
  static scalar_t BIG = std::is_same_v<scalar_t,double> ?
    4.503599627370496e15 : 16777216.;
  static scalar_t BIGINV = std::is_same_v<scalar_t,double> ?
    2.22044604925031308085e-16 : 5.9604644775390625E-8;

  ax = _igam_helper_fac(a, x);
  if (ax == 0.0) {
    return 0.0;
  }

  /* continued fraction */
  y = 1.0 - a;
  z = x + y + 1.0;
  c = 0.0;
  pkm2 = 1.0;
  qkm2 = x;
  pkm1 = x + 1.0;
  qkm1 = z * x;
  ans = pkm1 / qkm1;

  for (i = 0; i < MAXITER; i++) {
    c += 1.0;
    y += 1.0;
    z += 2.0;
    yc = y * c;
    pk = pkm1 * z - pkm2 * yc;
    qk = qkm1 * z - qkm2 * yc;
    if (qk != 0) {
      r = pk / qk;
      t = std::fabs((ans - r) / r);
      ans = r;
    }
    else {
      t = 1.0;
    }
    pkm2 = pkm1;
    pkm1 = pk;
    qkm2 = qkm1;
    qkm1 = qk;
    if (std::fabs(pk) > BIG) {
      pkm2 *= BIGINV;
      pkm1 *= BIGINV;
      qkm2 *= BIGINV;
      qkm1 *= BIGINV;
    }
    if (t <= MACHEP) {
      break;
    }
  }
  return ans * ax;
}