From maechler@stat.math.ethz.ch  Sat Mar  2 16:02:27 2002
From: maechler@stat.math.ethz.ch (maechler@stat.math.ethz.ch)
Date: Sat, 2 Mar 2002 17:02:27 +0100 (MET)
Subject: [Rd] pexp.c (PR#1335)
Message-ID: <200203021602.RAA19606@pubhealth.ku.dk>

>>>>> "Morten" == Morten Welinder <terra@diku.dk> writes:

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

yes, it looks fine; I've compared it with several other versions
(in R, of course).

I've added it for R-devel -- only when there's no system expm1().

    Morten> double myexpm1 (double x) { double y;

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

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

I've also improved pexp() following your suggestion -- but note my
answer to #1334 (pweibull) which applies here as well!

Martin

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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)
maechler@stat.math.ethz.ch