view release on metacpan or  search on metacpan

real.c  view on Meta::CPAN

        * x) * x) * x) * x) * x) * x) * x) * x) * x;

  } else if (x < 1) {
    /* This and the rest from Vazquez-Leal et al. (2019). */
    k1 = sqrtl(1.0L + 2.7182818284590452353603L * x);
    k2 = 0.33333333333333333333333L + 0.70710678118654752440084L / k1 - 0.058925565098878960366737L * k1 +
         (x + 0.36787944117144L) * (0.050248489761611L + (0.11138904851051 + 0.040744556245195L * x) * x)
         /
         (1.0L + (2.7090878606183L + (1.5510922597820L + 0.095477712183841L * x) * x) * x);
    w = -(k2-1)/k2;

  } else if (x < 40) {
    k1 = 1.0L + (5.950065500550155L + (13.96586471370701L + (10.52192021050505L + (3.065294254265870L + 0.1204576876518760L * x) * x) * x) * x) * x;
    w = 0.1600049638651493L * logl(k1);

  } else if (x < 20000) {
    k1 = 1.0L + (-3.16866642511229e11L + (3.420439800038598e10L +
         (-1.501433652432257e9L + (3.44887729947585e7L + (-4.453783741137856e5L +
         (3257.926478908996L + (-10.82545259305382L + (0.6898058947898353e-1L +
         0.4703653406071575e-4L * x) * x) * x) * x) * x) * x) * x) * x) * x;
    w = 0.9898045358731312e-1L * logl(k1);

  } else {
    k1 = 1.0L / (1.0L + logl(1.0L + x));
    k2 = 1.0L / k1;
    k3 = logl(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;
  }
  return w;
}

NV lambertw(NV x) {
  long double w;
  int i;

  if (x < -0.36787944117145L)
    croak("Invalid input to LambertW:  x must be >= -1/e");
  if (x == 0.0L) return 0.0L;

  /* Estimate initial value */
  w = _lambertw_approx(x);

  /* TODO: this section might not be best for quad precision */
  /* If input is too small, return .99999.... */
  /* if (w <= -1.0L) return -1.0L + LDBL_EPSILON; */
  /* For very small inputs, don't iterate, return approx directly. */
  if (x < -0.36768) return w;

#if 0  /* Halley */
  long double lastw = w;
  for (i = 0; i < 100; i++) {
    long double ew = expl(w);
    long double wew = w * ew;
    long double wewx = wew - x;
    long double w1 = w + 1;
    w = w - wewx / (ew * w1 - (w+2) * wewx/(2*w1));
    if (w != 0.0L && fabsl((w-lastw)/w) <= 8*LDBL_EPSILON) break;
    lastw = w;
  }
#else  /* Fritsch, see Veberic 2009.  1-2 iterations are enough. */
  for (i = 0; i < 6 && w != 0.0L; i++) {
    long double w1 = 1 + w;
    long double zn = logl((long double)x/w) - w;
    long double qn = 2 * w1 * (w1+(2.0L/3.0L)*zn);
    long double en = (zn/w1) * (qn-zn)/(qn-2.0L*zn);
    /* w *= 1.0L + en;  if (fabsl(en) <= 16*LDBL_EPSILON) break; */
    long double wen = w * en;
    if (isnan(wen)) return 0;
    w += wen;
    if (fabsl(wen) <= 64*LDBL_EPSILON) break;
  }
#endif

#if LNV_IS_QUAD /* For quadmath, one high precision correction */
  if (w != LNV_ZERO) {
    LNV lw = w;
    LNV w1 = LNV_ONE + lw;
    LNV zn = loglnv((LNV)x/lw) - lw;
    LNV qn = LNVCONST(2.0) * w1 * (w1+(LNVCONST(2.0)/LNVCONST(3.0))*zn);
    LNV en = (zn/w1) * (qn-zn)/(qn-LNVCONST(2.0)*zn);
    return lw + lw * en;
  }
#endif

  /* With long double = 64-bit double, we have 15 digits precision
   * near the branch point, and 16 over the rest of the range.
   * With long double = x86 extended precision, we have over 17 digits
   * over the entire range.
   * Correcting to the exact LDBL_EPSILON does not improve this. */

  return w;
}


/******************************************************************************/
/*                        Chebyshev PSI / THETA                               */
/******************************************************************************/

NV chebyshev_psi(UV n)
{
  UV k;
  SUM_INIT(sum);

  for (k = log2floor(n); k > 0; k--) {
    SUM_ADD(sum, chebyshev_theta(rootint(n,k)));
  }
  return SUM_FINAL(sum);
}

#if BITS_PER_WORD == 64
typedef struct {
  UV n;
  LNV theta;
} cheby_theta_t;
static const cheby_theta_t _cheby_theta[] = { /* >= quad math precision */
  { UVCONST(      67108864),LNVCONST(    67100507.6357700963903836828562472350035880) },
  { UVCONST(     100000000),LNVCONST(    99987730.0180220043832124342600487053812729) },
  { UVCONST(     134217728),LNVCONST(   134204014.5735572091791081610859055728165544) },
  { UVCONST(     268435456),LNVCONST(   268419741.6134308193112682817754501071404173) },
  { UVCONST(     536870912),LNVCONST(   536842885.8045763840625719515011160692495056) },