From terra@diku.dk  Fri Mar  1 14:58:08 2002
From: terra@diku.dk (terra@diku.dk)
Date: Fri, 1 Mar 2002 15:58:08 +0100 (MET)
Subject: [Rd] pexp.c (PR#1335)
Message-ID: <200203011458.PAA16714@pubhealth.ku.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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._


[Rd] pexp.c (PR#1335)
terra@diku.dk