The R Project SVN R

Rev

Rev 19500 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 19500 Rev 19692
Line 47... Line 47...
47 
double expm1(double x)
 47 
double expm1(double x)
48 
{
 48 
{
49 
    double y, a = fabs(x);
 49 
    double y, a = fabs(x);
50 
 
 50 
 
51 
    if (a < DBL_EPSILON) return x;
 51 
    if (a < DBL_EPSILON) return x;
- 
 
 52 
    if (a > 0.697) return exp(x) - 1;  /* negligible cancellation */
52 
 
 53 
 
53 
    if (a > 1e-6) {
 54 
    if (a > 1e-8)
54 
	y = exp(x) - 1;
 55 
	y = exp(x) - 1;
55 
	if (y > 1.) /* no cancellation */
 - 
 
56 
	    return y;
 - 
 
57 
    }
 - 
 
58 
    else /* Taylor expansion */
 56 
    else /* Taylor expansion, more accurate in this range */
59 
	y = (x / 2 + 1) * x;
 57 
	y = (x / 2 + 1) * x;
60 
 
 58 
 
61 
    /* Newton step for solving   log(1 + y) = x   for y : */
 59 
    /* Newton step for solving   log(1 + y) = x   for y : */
- 
 
 60 
    /* WARNING: does not work for y ~ -1: bug in 1.5.0 */
62 
    y -= (1 + y) * (log1p (y) - x);
 61 
    y -= (1 + y) * (log1p (y) - x);
63 
    return y;
 62 
    return y;
64 
}
 63 
}
65 
#endif
 64 
#endif