The R Project SVN R

Rev

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

Rev 85634 Rev 87047
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 2000-2023 The R Core Team
 3 
 *  Copyright (C) 2000-2024 The R Core Team
4 
 *  Copyright (C) 2003	    The R Foundation
 4 
 *  Copyright (C) 2003	    The R Foundation
5 
 *  Copyright (C) 1998	    Ross Ihaka
 5 
 *  Copyright (C) 1998	    Ross Ihaka
6 
 *
 6 
 *
7 
 *  This program is free software; you can redistribute it and/or modify
 7 
 *  This program is free software; you can redistribute it and/or modify
8 
 *  it under the terms of the GNU General Public License as published by
 8 
 *  it under the terms of the GNU General Public License as published by
Line 221... Line 221...
221 
 * Cody had (-37.5193 < x  &&  x < 8.2924) ; R originally had y < 50
 221 
 * Cody had (-37.5193 < x  &&  x < 8.2924) ; R originally had y < 50
222 
 *
 222 
 *
223 
 * Note that we do want symmetry(0), lower/upper -> hence use y
 223 
 * Note that we do want symmetry(0), lower/upper -> hence use y
224 
 */
 224 
 */
225 
    else if((log_p && y < 1e170) /* avoid underflow below */
 225 
    else if((log_p && y < 1e170) /* avoid underflow below */
226 
	/*  ^^^^^ MM FIXME: could speed up for log_p and  |x| >> 5.657 !
 226 
	/*  ^^^^^ MM FIXME: could speed up for log_p and  y := |x| >> 5.657 !
227 
	 * Then, make use of  Abramowitz & Stegun, 26.2.13, p.932,  something like
 227 
	 * Then, make use of  Abramowitz & Stegun, 26.2.13, p.932,  something like
228 
 
 228 
 
229 
	 * Even smarter: work with   example(pnormAsymp, package="DPQ")
 229 
	 * Even smarter: work with   example(pnormAsymp, package="DPQ")
230 
 
 230 
 
231 
	 xsq = x*x;
 231 
	 xsq = x*x;
Line 241... Line 241...
241 
 
 241 
 
242 
	 Yes, but xsq might be infinite;
 242 
	 Yes, but xsq might be infinite;
243 
 	 well, actually  x = -1.34..e154 = -sqrt(DBL_MAX) already overflows x^2
 243 
 	 well, actually  x = -1.34..e154 = -sqrt(DBL_MAX) already overflows x^2
244 
	 The largest x for which  x/2*x is finite is
 244 
	 The largest x for which  x/2*x is finite is
245 
	 x = +/- 1.89615038e154 ~= sqrt(2) * sqrt(.Machine$double.xmax)
 245 
	 x = +/- 1.89615038e154 ~= sqrt(2) * sqrt(.Machine$double.xmax)
- 
 
 246 
 
- 
 
 247 
	 NB: allowing "DENORMS" ==> boundaries at +/- 38.4674  <--> qnorm(log(2^-1074), log.p=TRUE)
- 
 
 248 
	 --                               rather than 37.5193 (up to R 4.4.x)
246 
	*/
 249 
	*/
247 
	    || (lower && -37.5193 < x  &&  x < 8.2924)
 250 
	    || (lower && -38.4674 < x  &&  x < 8.2924)
248 
	    || (upper && -8.2924  < x  &&  x < 37.5193)
 251 
	    || (upper && -8.2924  < x  &&  x < 38.4674)
249 
	) {
 252 
	) {
250 
 
 253 
 
251 
	/* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
 254 
	/* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
252 
	xsq = 1.0 / (x * x); /* (1./x)*(1./x) might be better */
 255 
	xsq = 1.0 / (x * x); /* (1./x)*(1./x) might be better */
253 
	xnum = p[5] * xsq;
 256 
	xnum = p[5] * xsq;
Line 259... Line 262...
259 
	temp = xsq * (xnum + p[4]) / (xden + q[4]);
 262 
	temp = xsq * (xnum + p[4]) / (xden + q[4]);
260 
	temp = (M_1_SQRT_2PI - temp) / y;
 263 
	temp = (M_1_SQRT_2PI - temp) / y;
261 
 
 264 
 
262 
	do_del(x);
 265 
	do_del(x);
263 
	swap_tail;
 266 
	swap_tail;
264 
    } else { /* large x such that probs are 0 or 1 */
 267 
    } else { /* large |x| such that probs are 0 or 1 */
265 
	if(x > 0) {	*cum = R_D__1; *ccum = R_D__0;	}
 268 
	if(x > 0) {	*cum = R_D__1; *ccum = R_D__0;	}
266 
	else {	        *cum = R_D__0; *ccum = R_D__1;	}
 269 
	else {	        *cum = R_D__0; *ccum = R_D__1;	}
267 
    }
 270 
    }
268 
 
 271 
 
269 
 
 272