[Rd] pexp.c (PR#1335)

[maechler@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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._