[Rd] pexp.c (PR#1335)
maechler@stat.math.ethz.ch
maechler@stat.math.ethz.ch
Sat, 2 Mar 2002 17:02:27 +0100 (MET)
>>>>> "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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._