The R Project SVN – R – Diff – /trunk/src/nmath/pnchisq.c



 
 * -- Feb.6, 2000 (R pre0.99); M.Maechler:  still have
 * bad precision & non-convergence in some cases (x ~= f, both LARGE)
 */
 
#ifdef HAVE_LONG_DOUBLE
# define EXP expl
# define FABS fabsl
#else
# define EXP exp
# define FABS fabs
#endif
 
double pnchisq(double x, double df, double ncp, int lower_tail, int log_p)
{
    double ans;
#ifdef IEEE_754

{
    double lam, x2, f2, term, bound, f_x_2n, f_2n;
    double l_lam = -1., l_x = -1.; /* initialized for -Wall */
    int n;
    Rboolean lamSml, tSml, is_r, is_b, is_it;
    LDOUBLE ans, u, v, t, lt, lu =-1;
 
    static const double _dbl_min_exp = M_LN2 * DBL_MIN_EXP;
    /*= -708.3964 for IEEE double precision */
 
    if (x <= 0.) {

#ifndef MATHLIB_STANDALONE
    R_CheckUserInterrupt();
#endif
 
    if(theta < 80) { /* use 110 for Inf, as ppois(110, 80/2, lower.tail=FALSE) is 2e-20 */
	LDOUBLE sum = 0, sum2 = 0, lambda = 0.5*theta, 
	    pr = EXP(-lambda); // does this need a feature test?
	double ans;
	int i;
	/* we need to renormalize here: the result could be very close to 1 */
	for(i = 0; i < 110;  pr *= lambda/++i) {
	    sum2 += pr;

#ifdef DEBUG_pnch
    REprintf("-- v=exp(-th/2)=%g, x/2= %g, f/2= %g\n",v,x2,f2);
#endif
 
    if(f2 * DBL_EPSILON > 0.125 && /* very large f and x ~= f: probably needs */
       FABS(t = x2 - f2) <         /* another algorithm anyway */
       sqrt(DBL_EPSILON) * f2) {
	/* evade cancellation error */
	/* t = exp((1 - t)*(2 - t/(f2 + 1))) / sqrt(2*M_PI*(f2 + 1));*/
        lt = (1 - t)*(2 - t/(f2 + 1)) - 0.5 * log(2*M_PI*(f2 + 1));
#ifdef DEBUG_pnch

	} /* else */
	l_x = log(x);
	ans = term = 0.; t = 0;
    }
    else {
	t = EXP(lt);
#ifdef DEBUG_pnch
 	REprintf(", t=exp(lt)= %g\n", t);
#endif
	ans = term = (double) (v * t);
    }

		/* no underflow anymore ==> change regime */
#ifdef DEBUG_pnch
                REprintf(" n=%d; nomore underflow in u = exp(lu) ==> change\n",
			 n);
#endif
                v = u = EXP(lu); /* the first non-0 'u' */
                lamSml = FALSE;
            }
        } else {
	    u *= lam / n;
	    v += u;

		/* no underflow anymore ==> change regime */
#ifdef DEBUG_pnch
                REprintf("  n=%d; nomore underflow in t = exp(lt) ==> change\n",
			 n);
#endif
                t = EXP(lt); /* the first non-0 't' */
                tSml = FALSE;
            }
        } else {
	    t *= x / f_2n;
	}

Rev 59171	Rev 60268
Line 24...	Line 24...
24		24
25	`* -- Feb.6, 2000 (R pre0.99); M.Maechler: still have`	25	`* -- Feb.6, 2000 (R pre0.99); M.Maechler: still have`
26	`* bad precision & non-convergence in some cases (x ~= f, both LARGE)`	26	`* bad precision & non-convergence in some cases (x ~= f, both LARGE)`
27	`*/`	27	`*/`
28		28
-		29	`#ifdef HAVE_LONG_DOUBLE`
-		30	`# define EXP expl`
-		31	`# define FABS fabsl`
-		32	`#else`
-		33	`# define EXP exp`
-		34	`# define FABS fabs`
-		35	`#endif`
29		36
30	`double pnchisq(double x, double df, double ncp, int lower_tail, int log_p)`	37	`double pnchisq(double x, double df, double ncp, int lower_tail, int log_p)`
31	`{`	38	`{`
32	`double ans;`	39	`double ans;`
33	`#ifdef IEEE_754`	40	`#ifdef IEEE_754`
Line 62...	Line 69...
62	`{`	69	`{`
63	`double lam, x2, f2, term, bound, f_x_2n, f_2n;`	70	`double lam, x2, f2, term, bound, f_x_2n, f_2n;`
64	`double l_lam = -1., l_x = -1.; /* initialized for -Wall */`	71	`double l_lam = -1., l_x = -1.; /* initialized for -Wall */`
65	`int n;`	72	`int n;`
66	`Rboolean lamSml, tSml, is_r, is_b, is_it;`	73	`Rboolean lamSml, tSml, is_r, is_b, is_it;`
67	`long double ans, u, v, t, lt, lu =-1;`	74	`LDOUBLE ans, u, v, t, lt, lu =-1;`
68		75
69	`static const double _dbl_min_exp = M_LN2 * DBL_MIN_EXP;`	76	`static const double _dbl_min_exp = M_LN2 * DBL_MIN_EXP;`
70	`/= -708.3964 for IEEE double precision /`	77	`/= -708.3964 for IEEE double precision /`
71		78
72	`if (x <= 0.) {`	79	`if (x <= 0.) {`
Line 81...	Line 88...
81	`#ifndef MATHLIB_STANDALONE`	88	`#ifndef MATHLIB_STANDALONE`
82	`R_CheckUserInterrupt();`	89	`R_CheckUserInterrupt();`
83	`#endif`	90	`#endif`
84		91
85	`if(theta < 80) { /* use 110 for Inf, as ppois(110, 80/2, lower.tail=FALSE) is 2e-20 */`	92	`if(theta < 80) { /* use 110 for Inf, as ppois(110, 80/2, lower.tail=FALSE) is 2e-20 */`
86	`long double sum = 0, sum2 = 0, lambda = 0.5*theta,`	93	`LDOUBLE sum = 0, sum2 = 0, lambda = 0.5*theta,`
87	`pr = expl(-lambda); // does this need a feature test?`	94	`pr = EXP(-lambda); // does this need a feature test?`
88	`double ans;`	95	`double ans;`
89	`int i;`	96	`int i;`
90	`/* we need to renormalize here: the result could be very close to 1 */`	97	`/* we need to renormalize here: the result could be very close to 1 */`
91	`for(i = 0; i < 110; pr *= lambda/++i) {`	98	`for(i = 0; i < 110; pr *= lambda/++i) {`
92	`sum2 += pr;`	99	`sum2 += pr;`
Line 123...	Line 130...
123	`#ifdef DEBUG_pnch`	130	`#ifdef DEBUG_pnch`
124	`REprintf("-- v=exp(-th/2)=%g, x/2= %g, f/2= %g\n",v,x2,f2);`	131	`REprintf("-- v=exp(-th/2)=%g, x/2= %g, f/2= %g\n",v,x2,f2);`
125	`#endif`	132	`#endif`
126		133
127	`if(f2 * DBL_EPSILON > 0.125 && /* very large f and x ~= f: probably needs */`	134	`if(f2 * DBL_EPSILON > 0.125 && /* very large f and x ~= f: probably needs */`
128	`fabsl(t = x2 - f2) < /* another algorithm anyway */`	135	`FABS(t = x2 - f2) < /* another algorithm anyway */`
129	`sqrt(DBL_EPSILON) * f2) {`	136	`sqrt(DBL_EPSILON) * f2) {`
130	`/* evade cancellation error */`	137	`/* evade cancellation error */`
131	`/* t = exp((1 - t)(2 - t/(f2 + 1))) / sqrt(2M_PI(f2 + 1));/`	138	`/* t = exp((1 - t)(2 - t/(f2 + 1))) / sqrt(2M_PI(f2 + 1));/`
132	`lt = (1 - t)(2 - t/(f2 + 1)) - 0.5 log(2M_PI(f2 + 1));`	139	`lt = (1 - t)(2 - t/(f2 + 1)) - 0.5 log(2M_PI(f2 + 1));`
133	`#ifdef DEBUG_pnch`	140	`#ifdef DEBUG_pnch`
Line 150...	Line 157...
150	`} /* else */`	157	`} /* else */`
151	`l_x = log(x);`	158	`l_x = log(x);`
152	`ans = term = 0.; t = 0;`	159	`ans = term = 0.; t = 0;`
153	`}`	160	`}`
154	`else {`	161	`else {`
155	`t = expl(lt);`	162	`t = EXP(lt);`
156	`#ifdef DEBUG_pnch`	163	`#ifdef DEBUG_pnch`
157	`REprintf(", t=exp(lt)= %g\n", t);`	164	`REprintf(", t=exp(lt)= %g\n", t);`
158	`#endif`	165	`#endif`
159	`ans = term = (double) (v * t);`	166	`ans = term = (double) (v * t);`
160	`}`	167	`}`
Line 200...	Line 207...
200	`/* no underflow anymore ==> change regime */`	207	`/* no underflow anymore ==> change regime */`
201	`#ifdef DEBUG_pnch`	208	`#ifdef DEBUG_pnch`
202	`REprintf(" n=%d; nomore underflow in u = exp(lu) ==> change\n",`	209	`REprintf(" n=%d; nomore underflow in u = exp(lu) ==> change\n",`
203	`n);`	210	`n);`
204	`#endif`	211	`#endif`
205	`v = u = expl(lu); /* the first non-0 'u' */`	212	`v = u = EXP(lu); /* the first non-0 'u' */`
206	`lamSml = FALSE;`	213	`lamSml = FALSE;`
207	`}`	214	`}`
208	`} else {`	215	`} else {`
209	`u *= lam / n;`	216	`u *= lam / n;`
210	`v += u;`	217	`v += u;`
Line 215...	Line 222...
215	`/* no underflow anymore ==> change regime */`	222	`/* no underflow anymore ==> change regime */`
216	`#ifdef DEBUG_pnch`	223	`#ifdef DEBUG_pnch`
217	`REprintf(" n=%d; nomore underflow in t = exp(lt) ==> change\n",`	224	`REprintf(" n=%d; nomore underflow in t = exp(lt) ==> change\n",`
218	`n);`	225	`n);`
219	`#endif`	226	`#endif`
220	`t = expl(lt); /* the first non-0 't' */`	227	`t = EXP(lt); /* the first non-0 't' */`
221	`tSml = FALSE;`	228	`tSml = FALSE;`
222	`}`	229	`}`
223	`} else {`	230	`} else {`
224	`t *= x / f_2n;`	231	`t *= x / f_2n;`
225	`}`	232	`}`

The R Project SVN R

(root)/trunk/src/nmath/pnchisq.c – Rev 59171 → 60268