The R Project SVN R

Rev

Rev 88193 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 88193 Rev 88208
Line 57... Line 57...
57 
		x_ = ldexp(x, -2),
 57 
		x_ = ldexp(x, -2),
58 
		n_ = ldexp(np,-2);
 58 
		n_ = ldexp(np,-2);
59 
	    v = (x_ - n_)/(x_ + n_);
 59 
	    v = (x_ - n_)/(x_ + n_);
60 
	}
 60 
	}
61 
	double s = ldexp(d, -1) * v; // was d * v
 61 
	double s = ldexp(d, -1) * v; // was d * v
62 
	if(fabs(s) < DBL_MIN*0.5) return ldexp(s, 1);// CARE: assumes subnormal numbers
 62 
	if(fabs(ldexp(s, 1)) < DBL_MIN) return ldexp(s, 1);
63 
	double ej = x * v; // as 2*x*v could overflow:  v > 1/2  <==> ej = 2xv > x
 63 
	double ej = x * v; // as 2*x*v could overflow:  v > 1/2  <==> ej = 2xv > x
64 
	v *= v; // "v = v^2"
 64 
	v *= v; // "v = v^2"
65 
	for (int j = 1; j < 1000; j++) { /* Taylor series; 1000: no infinite loop
 65 
	for (int j = 1; j < 1000; j++) { /* Taylor series; 1000: no infinite loop
66 
					    as |v| < .1,  v^2000 is "zero" */
 66 
					    as |v| < .1,  v^2000 is "zero" */
67 
	    ej *= v;// = x v^(2j+1)
 67 
	    ej *= v;// = x v^(2j+1)
Line 77... Line 77...
77 
	/* ---- the following should _never_ happen ------------ */
 77 
	/* ---- the following should _never_ happen ------------ */
78 
	MATHLIB_WARNING4("bd0(%g, %g): T.series failed to converge in 1000 it.; s=%g, ej/(2j+1)=%g\n",
 78 
	MATHLIB_WARNING4("bd0(%g, %g): T.series failed to converge in 1000 it.; s=%g, ej/(2j+1)=%g\n",
79 
			 x, np, s, ej/((1000<<1)+1));
 79 
			 x, np, s, ej/((1000<<1)+1));
80 
    }
 80 
    }
81 
    /* else:  | x - np |  is not too small */
 81 
    /* else:  | x - np |  is not too small */
- 
 
 82 
    /* NB: x/np |--> inf (overflow)  doesn't happen when called from rpois_raw() */
- 
 
 83 
#define lg_x_n (R_FINITE(x/np) ? log(x/np) : (log(x) - log(np)))
82 
    return (x > np) ? x*(log(x/np) -1.) + np
 84 
    return (x > np) ? x*(lg_x_n -1.) + np
83 
	            : x* log(x/np) + np -x;
 85 
	            : x* lg_x_n + np -x;
- 
 
 86 
#undef lg_x_n
84 
}
 87 
}
85 
 
 88 
 
86 
 
 89 
 
87 
// ebd0(): R Bugzilla PR#15628 -- proposed accuracy improvement by Morten Welinder
 90 
// ebd0(): R Bugzilla PR#15628 -- proposed accuracy improvement by Morten Welinder
88 
 
 91