The R Project SVN R

Rev

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

Rev 23377 Rev 24322
Line 117... Line 117...
117 
     *    JASA 71, 893-896.
 117 
     *    JASA 71, 893-896.
118 
     */
 118 
     */
119 
 
 119 
 
120 
#define C1		0.398942280401433
 120 
#define C1		0.398942280401433
121 
#define C2		0.180025191068563
 121 
#define C2		0.180025191068563
122 
#define g(x)		(C1*exp(-x*x/2.0)-C2*(A-fabs(x)))
 122 
#define g(x)		(C1*exp(-x*x/2.0)-C2*(A-x))
123 
 
 123 
 
124 
    const double A =  2.216035867166471;
 124 
    const double A =  2.216035867166471;
125 
 
 125 
 
126 
    double s, u1, w, y, u2, u3, aa, tt, theta, R;
 126 
    double s, u1, w, y, u2, u3, aa, tt, theta, R;
127 
    int i;
 127 
    int i;
Line 191... Line 191...
191 
	y = aa + w;
 191 
	y = aa + w;
192 
	return (s == 1.0) ? -y : y;
 192 
	return (s == 1.0) ? -y : y;
193 
 
 193 
 
194 
	/*-----------------------------------------------------------*/
 194 
	/*-----------------------------------------------------------*/
195
195
196 
    case KINDERMAN_RAMAGE: /* see Reference above */
 196 
    case BUGGY_KINDERMAN_RAMAGE: /* see Reference above */
- 
 
 197 
	/* note: this has problems, but is retained for
- 
 
 198 
	 * reproducibility of older codes, with the same
- 
 
 199 
	 * numeric code */
197 
	u1 = unif_rand();
 200 
	u1 = unif_rand();
198 
	if(u1 < 0.884070402298758) {
 201 
	if(u1 < 0.884070402298758) {
199 
	    u2 = unif_rand();
 202 
	    u2 = unif_rand();
200 
	    return A*(1.13113163544180*u1+u2-1);
 203 
	    return A*(1.13113163544180*u1+u2-1);
201 
	}
 204 
	}
Line 261... Line 264...
261 
#define BIG 134217728 /* 2^27 */
 264 
#define BIG 134217728 /* 2^27 */
262 
	/* unif_rand() alone is not of high enough precision */
 265 
	/* unif_rand() alone is not of high enough precision */
263 
	u1 = unif_rand();
 266 
	u1 = unif_rand();
264 
	u1 = (int)(BIG*u1) + unif_rand();
 267 
	u1 = (int)(BIG*u1) + unif_rand();
265 
	return qnorm5(u1/BIG, 0.0, 1.0, 1, 0);
 268 
	return qnorm5(u1/BIG, 0.0, 1.0, 1, 0);
- 
 
 269 
    case KINDERMAN_RAMAGE: /* see Reference above */
- 
 
 270 
	/* corrected version from Josef Leydold
- 
 
 271 
	 * */
- 
 
 272 
	u1 = unif_rand();
- 
 
 273 
	if(u1 < 0.884070402298758) {
- 
 
 274 
	    u2 = unif_rand();
- 
 
 275 
	    return A*(1.131131635444180*u1+u2-1);
- 
 
 276 
	}
- 
 
 277
- 
 
 278 
	if(u1 >= 0.973310954173898) { /* tail: */
- 
 
 279 
	    repeat {
- 
 
 280 
		u2 = unif_rand();
- 
 
 281 
		u3 = unif_rand();
- 
 
 282 
		tt = (A*A-2*log(u3));
- 
 
 283 
		if( u2*u2<(A*A)/tt )
- 
 
 284 
		    return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
- 
 
 285 
	    }
- 
 
 286 
	}
- 
 
 287
- 
 
 288 
	if(u1 >= 0.958720824790463) { /* region3: */
- 
 
 289 
	    repeat {
- 
 
 290 
		u2 = unif_rand();
- 
 
 291 
		u3 = unif_rand();
- 
 
 292 
		tt = A - 0.630834801921960* fmin2(u2,u3);
- 
 
 293 
		if(fmax2(u2,u3) <= 0.755591531667601)
- 
 
 294 
		    return (u2<u3) ? tt : -tt;
- 
 
 295 
		if(0.034240503750111*fabs(u2-u3) <= g(tt))
- 
 
 296 
		    return (u2<u3) ? tt : -tt;
- 
 
 297 
	    }
- 
 
 298 
	}
- 
 
 299
- 
 
 300 
	if(u1 >= 0.911312780288703) { /* region2: */
- 
 
 301 
	    repeat {
- 
 
 302 
		u2 = unif_rand();
- 
 
 303 
		u3 = unif_rand();
- 
 
 304 
		tt = 0.479727404222441+1.105473661022070*fmin2(u2,u3);
- 
 
 305 
		if( fmax2(u2,u3)<=0.872834976671790 )
- 
 
 306 
		    return (u2<u3) ? tt : -tt;
- 
 
 307 
		if( 0.049264496373128*fabs(u2-u3)<=g(tt) )
- 
 
 308 
		    return (u2<u3) ? tt : -tt;
- 
 
 309 
	    }
- 
 
 310 
	}
- 
 
 311 
 
- 
 
 312 
	/* ELSE	 region1: */
- 
 
 313 
	repeat {
- 
 
 314 
	    u2 = unif_rand();
- 
 
 315 
	    u3 = unif_rand();
- 
 
 316 
	    tt = 0.479727404222441-0.595507138015940*fmin2(u2,u3);
- 
 
 317 
	    if (tt < 0.) continue;
- 
 
 318 
	    if(fmax2(u2,u3) <= 0.805577924423817)
- 
 
 319 
		return (u2<u3) ? tt : -tt;
- 
 
 320 
     	    if(0.053377549506886*fabs(u2-u3) <= g(tt))
- 
 
 321 
		return (u2<u3) ? tt : -tt;
- 
 
 322 
	}
266 
    default:
 323 
    default:
267 
	MATHLIB_ERROR("norm_rand(): invalid N01_kind: %d\n", N01_kind)
 324 
	MATHLIB_ERROR("norm_rand(): invalid N01_kind: %d\n", N01_kind)
268 
	    return 0.0;/*- -Wall */
 325 
	    return 0.0;/*- -Wall */
269 
    }/*switch*/
 326 
    }/*switch*/
270 
}
 327 
}