Full_Name: M Welinder Version: 1.4 OS: (src) Submission from: (NULL) (192.5.35.38) It seems to me that pexp can be improved in the lower_tail=TRUE and log_p=FALSE case by using expm1. Something like -expm1 (-x / scale); I think. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>>>>> "terra" == terra <terra@diku.dk> writes:terra> Full_Name: M Welinder Version: 1.4 OS: (src) terra> Submission from: (NULL) (192.5.35.38) terra> It seems to me that pexp can be improved in the terra> lower_tail=TRUE and log_p=FALSE case by using expm1. terra> Something like terra> -expm1 (-x / scale); terra> I think. Thank you, "terra". Similar thoughts had crossed my mind quite a few months ago, e.g. you find in src/nmath/qt.c : /* FIXME: Following cutoff is machine-precision dependent ----- Really, use stable impl. of expm1(y) == exp(y) - 1, as it is in GNU's mathlib ..*/ if (y > 1e-6) /* was (y > 0.002) */ y = exp(y) - 1; else { /* Taylor of e^y -1 : */ y = (0.5 * y + 1) * y; } Unfortunately, expm1() is NOT ANSI C / ISO C as far as I know. For the `complementary' function, log1p(), we have been having our own builtin (src/nmath/log1p.c) to use (unfortunately, we've used it unconditionally; from 1.5.0 on, we will only use it when there's no system log1p() available). For expm1() it's the same: We (i.e. the "configure" script) should check for availability of a system "exp1m()" and when it is not available, for portability, we need to provide our own version. Since GNU's mathlib does provide it, the source code is available under GPL, i.e. directly usable for R. Would you provide a src/nmath/expm1.c file, modeled after R-devel's (not R 1.4.x!) src/nmath/log1p.c one? Feel free to ask me for more advice, thank you in advance! Martin Maechler <maechler@stat.math.ethz.ch> http://stat.ethz.ch/~maechler/ Seminar fuer Statistik, ETH-Zentrum LEO C16 Leonhardstr. 27 ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND phone: x-41-1-632-3408 fax: ...-1228 <>< -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>>>>> "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...) Morten> Morten looks very reasonable; this might make it (no promises yet). Martin Morten> double Morten> myexpm1 (double x) Morten> { Morten> double y; Morten> if (fabs (x) > 1e-6) Morten> { Morten> y = exp (x) - 1; Morten> if (y > 10000) Morten> return y; Morten> } Morten> else Morten> y = (x / 2 + 1) * x; /* Taylor expansion */ Morten> /* Newton step. */ Morten> y -= (1 + y) * (log1p (y) - x); Morten> return y; Morten> } -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>>>>> "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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._