[Rd] pexp.c (PR#1335)

[terra@diku.dk]

If you want an expm1, this will do just fine, given the availability
of log1p.  It compares well with Solaris' expm1.  (And, luckily, is
much better than just exp(x)-1...)

Morten



double
myexpm1 (double x)
{
  double y;

  if (fabs (x) > 1e-6)
    {
      y = exp (x) - 1;
      if (y > 10000)
      return y;
    }
  else
    y = (x / 2 + 1) * x; /* Taylor expansion */

  /* Newton step.  */
  y -= (1 + y) * (log1p (y) - x);
  return y;
}

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._