[Rd] pexp.c (PR#1335)

terra@diku.dk terra@diku.dk
Fri, 1 Mar 2002 15:58:08 +0100 (MET)



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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._