view release on metacpan or  search on metacpan

real.c  view on Meta::CPAN

  if (n < 46) {
    _bernfrac_comb(num, den, n, t);
  } else {
    _bernfrac_zeta(num, den, n, t);
  }
  mpz_gcd(t, num, den);
  mpz_divexact(num, num, t);
  mpz_divexact(den, den, t);
  mpz_clear(t);
}

/***************************       Lambert W       ***************************/

static double _lambertw_approx(double x) {
  double w, k1, k2, k3;

  if (isinf(x) || isnan(x))
    return 1000.0;

  if (x < -0.312) {
    /* Near the branch point.  See Fukushima (2013) section 2.5. */
    k2 = 2.0 * (1.0 + 2.71828182845904523536 * x);
    if (k2 <= 0) return -1.0 + 1*DBL_EPSILON;
    k1 = sqrt(k2);
    /* Use Puiseux series, e.g. Verberic 2009, Boost, Johannson (2020). */
    w = -1.0 + (1.0 + (-1.0/3.0 + (11.0/72.0 + (-43.0/540.0 + (769.0/17280.0 + (-221.0/8505.0 + (680863.0/43545600.0 + (-1963.0/204120.0 + 226287557.0/37623398400.0
    * k1) * k1) * k1) * k1) * k1) * k1) * k1) * k1) * k1;

  } else if (x > -0.14 && x < 0.085) {
    /* Around zero.  See Fukushima (2013) section 2.6. */
    w = (1.0 + (-1.0 + (3.0/2.0 + (-8.0/3.0 + (125.0/24.0 + (-54.0/5.0 + (16807.0/720.0 + (-16384.0/315.0 + 531441.0/4480.0
        * x) * x) * x) * x) * x) * x) * x) * x) * x;

  } else if (x < 1) {
    /* This and the rest from Vazquez-Leal et al. (2019). */
    k1 = sqrt(1.0 + 2.71828182845904523536 * x);
    k2 = 0.33333333333333333333 + 0.7071067811865475244 / k1 - 0.058925565098880 * k1 +
         (x + 0.36787944117144) * (0.050248489761611 + (0.11138904851051 + 0.040744556245195 * x) * x)
         /
         (1.0 + (2.7090878606183 + (1.5510922597820 + 0.095477712183841 * x) * x) * x);
    w = -(k2-1)/k2;

  } else if (x < 40) {
    k1 = 1.0 + (5.950065500550155 + (13.96586471370701 + (10.52192021050505 + (3.065294254265870 + 0.1204576876518760 * x) * x) * x) * x) * x;
    w = 0.1600049638651493 * log(k1);
  } else if (x < 20000) {
    k1 = 1.0 + (-3.16866642511229e11 + (3.420439800038598e10 +
         (-1.501433652432257e9 + (3.44887729947585e7 + (-4.453783741137856e5 +
         (3257.926478908996 + (-10.82545259305382 + (0.6898058947898353e-1 +
         0.4703653406071575e-4 * x) * x) * x) * x) * x) * x) * x) * x) * x;
    w = 0.9898045358731312e-1 * log(k1);

  } else {
    k1 = 1.0 / (1.0 + log(1.0 + x));
    k2 = 1.0 / k1;
    k3 = log(k2);
    w = k2-1-k3+(1+k3+(-1/2+(1/2)*k3*k3 +(-1/6+(-1+(-1/2+
        (1/3) * k3) * k3) * k3) * k1) * k1) * k1;
  }

  /* Improve the FP estimate using two simple Halley iterations. */
  if (x >= -0.36728) {
    if (w != 0) w = (w/(1.0+w)) * (1.0+log(x/w));
    if (w != 0) w = (w/(1.0+w)) * (1.0+log(x/w));
    if (isnan(w)) w = DBL_EPSILON;
  }

  /* Result is over 15 digits, 16 when x > 0 */
  return w;
}

/* Iteration (5) from https://arxiv.org/pdf/2008.06122. */
static void _mpf_lambertw_approx(mpf_t w, mpf_t x) {
  double dapprox;
  mpf_t T1, T2;
  const int apbits = 117;   /* 35 digits */
  int i, loops;

  mpf_init2(T1, apbits);
  mpf_init2(T2, apbits);

  loops = 2;
  dapprox = _lambertw_approx(mpf_get_d(x));
  if (dapprox < 800.0) {
    mpf_set_d(T1, dapprox);
  } else {
    /* Init using equation 21 */
    mpf_log(T1, x);
    mpf_log(T2, T1);
    mpf_sub(T1, T1, T2);    /* T1 = Beta0 = log(x)-log(log(x)) */
    loops = 3;
  }

  for (i = 0; i < loops; i++) {
    mpf_add_ui(T2, T1, 1);
    mpf_div(T2, T1, T2);    /* T2 = Bn/(1+Bn) */
    mpf_div(T1, x, T1);
    mpf_log(T1, T1);
    mpf_add_ui(T1,T1,1);    /* T1 = 1 + log(x/Bn) */
    mpf_mul(T1, T1, T2);    /* T1 = Beta_{n+1} */
  }
  mpf_set(w, T1);
  mpf_clear(T2); mpf_clear(T1);
}

static void _lambertw(mpf_t r, mpf_t x, unsigned long prec)
{
  int i;
  unsigned long bits = 96 + DIGS2BITS(prec);
  mpf_t w, t, tol, w1, zn, qn, en;

  if (mpf_cmp_d(x, -0.36787944117145) < 0)
    croak("Invalid input to LambertW:  x must be >= -1/e");
  if (mpf_sgn(x) == 0)
    { mpf_set(r, x); return; }

  /* Use Fritsch rather than Halley. */
  mpf_init2(w,   bits);
  mpf_init2(t,   bits);
  mpf_init2(tol, bits);
  mpf_init2(w1,  bits);