src/nmath/snorm.c - R - Git at Google

 /*
  *  Mathlib : A C Library of Special Functions
  *  Copyright (C) 1998   Ross Ihaka
  *  Copyright (C) 2000-9 The R Core Team
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation; either version 2 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, a copy is available at
  *  https://www.R-project.org/Licenses/
  *
  *  SYNOPSIS
  *
  *    #include <Rmath.h>
  *    double norm_rand(void);
  *
  *  DESCRIPTION
  *
  *    Random variates from the STANDARD normal distribution  N(0,1).
  *
  * Is called from  rnorm(..), but also rt(), rf(), rgamma(), ...
  */

 #include <R_ext/Random.h>
 #include "nmath.h"

 #define repeat for(;;)

 #ifdef MATHLIB_STANDALONE
 static
 #else
 attribute_hidden
 #endif
 double BM_norm_keep = 0.0;

 N01type N01_kind = INVERSION;

 #ifndef MATHLIB_STANDALONE
 typedef void * (*DL_FUNC)();
 extern DL_FUNC  User_norm_fun; /* declared and set in ../main/RNG.c */
 #endif

 /*
  *  REFERENCE
  *
  *    Ahrens, J.H. and Dieter, U.
  *    Extensions of Forsythe's method for random sampling from
  *    the normal distribution.
  *    Math. Comput. 27, 927-937.
  *
  *    The definitions of the constants a[k], d[k], t[k] and
  *    h[k] are according to the abovementioned article
  */
 double norm_rand(void)
 {

     const static double a[32] =
     {
 	0.0000000, 0.03917609, 0.07841241, 0.1177699,
 	0.1573107, 0.19709910, 0.23720210, 0.2776904,
 	0.3186394, 0.36012990, 0.40225010, 0.4450965,
 	0.4887764, 0.53340970, 0.57913220, 0.6260990,
 	0.6744898, 0.72451440, 0.77642180, 0.8305109,
 	0.8871466, 0.94678180, 1.00999000, 1.0775160,
 	1.1503490, 1.22985900, 1.31801100, 1.4177970,
 	1.5341210, 1.67594000, 1.86273200, 2.1538750
     };

     const static double d[31] =
     {
 	0.0000000, 0.0000000, 0.0000000, 0.0000000,
 	0.0000000, 0.2636843, 0.2425085, 0.2255674,
 	0.2116342, 0.1999243, 0.1899108, 0.1812252,
 	0.1736014, 0.1668419, 0.1607967, 0.1553497,
 	0.1504094, 0.1459026, 0.1417700, 0.1379632,
 	0.1344418, 0.1311722, 0.1281260, 0.1252791,
 	0.1226109, 0.1201036, 0.1177417, 0.1155119,
 	0.1134023, 0.1114027, 0.1095039
     };

     const static double t[31] =
     {
 	7.673828e-4, 0.002306870, 0.003860618, 0.005438454,
 	0.007050699, 0.008708396, 0.010423570, 0.012209530,
 	0.014081250, 0.016055790, 0.018152900, 0.020395730,
 	0.022811770, 0.025434070, 0.028302960, 0.031468220,
 	0.034992330, 0.038954830, 0.043458780, 0.048640350,
 	0.054683340, 0.061842220, 0.070479830, 0.081131950,
 	0.094624440, 0.112300100, 0.136498000, 0.171688600,
 	0.227624100, 0.330498000, 0.584703100
     };

     const static double h[31] =
     {
 	0.03920617, 0.03932705, 0.03950999, 0.03975703,
 	0.04007093, 0.04045533, 0.04091481, 0.04145507,
 	0.04208311, 0.04280748, 0.04363863, 0.04458932,
 	0.04567523, 0.04691571, 0.04833487, 0.04996298,
 	0.05183859, 0.05401138, 0.05654656, 0.05953130,
 	0.06308489, 0.06737503, 0.07264544, 0.07926471,
 	0.08781922, 0.09930398, 0.11555990, 0.14043440,
 	0.18361420, 0.27900160, 0.70104740
     };

     /*----------- Constants and definitions for  Kinderman - Ramage --- */
     /*
      *  REFERENCE
      *
      *    Kinderman A. J. and Ramage J. G. (1976).
      *    Computer generation of normal random variables.
      *    JASA 71, 893-896.
      */

 #define C1		0.398942280401433
 #define C2		0.180025191068563
 #define g(x)		(C1*exp(-x*x/2.0)-C2*(A-x))

     const static double A =  2.216035867166471;

     double s, u1, w, y, u2, u3, aa, tt, theta, R;
     int i;

     switch(N01_kind) {

     case  AHRENS_DIETER: /* see Reference above */

 	u1 = unif_rand();
 	s = 0.0;
 	if (u1 > 0.5)
 	    s = 1.0;
 	u1 = u1 + u1 - s;
 	u1 *= 32.0;
 	i = (int) u1;
 	if (i == 32)
 	    i = 31;
 	if (i != 0) {
 	    u2 = u1 - i;
 	    aa = a[i - 1];
 	    while (u2 <= t[i - 1]) {
 		u1 = unif_rand();
 		w = u1 * (a[i] - aa);
 		tt = (w * 0.5 + aa) * w;
 		repeat {
 		    if (u2 > tt)
 			goto deliver;
 		    u1 = unif_rand();
 		    if (u2 < u1)
 			break;
 		    tt = u1;
 		    u2 = unif_rand();
 		}
 		u2 = unif_rand();
 	    }
 	    w = (u2 - t[i - 1]) * h[i - 1];
 	}
 	else {
 	    i = 6;
 	    aa = a[31];
 	    repeat {
 		u1 = u1 + u1;
 		if (u1 >= 1.0)
 		    break;
 		aa = aa + d[i - 1];
 		i = i + 1;
 	    }
 	    u1 = u1 - 1.0;
 	    repeat {
 		w = u1 * d[i - 1];
 		tt = (w * 0.5 + aa) * w;
 		repeat {
 		    u2 = unif_rand();
 		    if (u2 > tt)
 			goto jump;
 		    u1 = unif_rand();
 		    if (u2 < u1)
 			break;
 		    tt = u1;
 		}
 		u1 = unif_rand();
 	    }
 	  jump:;
 	}

       deliver:
 	y = aa + w;
 	return (s == 1.0) ? -y : y;

 	/*-----------------------------------------------------------*/

     case BUGGY_KINDERMAN_RAMAGE: /* see Reference above */
 	/* note: this has problems, but is retained for
 	 * reproducibility of older codes, with the same
 	 * numeric code */
 	u1 = unif_rand();
 	if(u1 < 0.884070402298758) {
 	    u2 = unif_rand();
 	    return A*(1.13113163544180*u1+u2-1);
 	}

 	if(u1 >= 0.973310954173898) { /* tail: */
 	    repeat {
 		u2 = unif_rand();
 		u3 = unif_rand();
 		tt = (A*A-2*log(u3));
 		if( u2*u2<(A*A)/tt )
 		    return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
 	    }
 	}

 	if(u1 >= 0.958720824790463) { /* region3: */
 	    repeat {
 		u2 = unif_rand();
 		u3 = unif_rand();
 		tt = A - 0.630834801921960* fmin2(u2,u3);
 		if(fmax2(u2,u3) <= 0.755591531667601)
 		    return (u2<u3) ? tt : -tt;
 		if(0.034240503750111*fabs(u2-u3) <= g(tt))
 		    return (u2<u3) ? tt : -tt;
 	    }
 	}

 	if(u1 >= 0.911312780288703) { /* region2: */
 	    repeat {
 		u2 = unif_rand();
 		u3 = unif_rand();
 		tt = 0.479727404222441+1.105473661022070*fmin2(u2,u3);
 		if( fmax2(u2,u3)<=0.872834976671790 )
 		    return (u2<u3) ? tt : -tt;
 		if( 0.049264496373128*fabs(u2-u3)<=g(tt) )
 		    return (u2<u3) ? tt : -tt;
 	    }
 	}

 	/* ELSE	 region1: */
 	repeat {
 	    u2 = unif_rand();
 	    u3 = unif_rand();
 	    tt = 0.479727404222441-0.595507138015940*fmin2(u2,u3);
 	    if(fmax2(u2,u3) <= 0.805577924423817)
 		return (u2<u3) ? tt : -tt;
 	}
     case BOX_MULLER:
 	if(BM_norm_keep != 0.0) { /* An exact test is intentional */
 	    s = BM_norm_keep;
 	    BM_norm_keep = 0.0;
 	    return s;
 	} else {
 	    theta = 2 * M_PI * unif_rand();
 	    R = sqrt(-2 * log(unif_rand())) + 10*DBL_MIN; /* ensure non-zero */
 	    BM_norm_keep = R * sin(theta);
 	    return R * cos(theta);
 	}
 #ifndef MATHLIB_STANDALONE
     case USER_NORM:
 	return *((double *) User_norm_fun());
 #endif
     case INVERSION:
 #define BIG 134217728 /* 2^27 */
 	/* unif_rand() alone is not of high enough precision */
 	u1 = unif_rand();
 	u1 = (int)(BIG*u1) + unif_rand();
 	return qnorm5(u1/BIG, 0.0, 1.0, 1, 0);
     case KINDERMAN_RAMAGE: /* see Reference above */
 	/* corrected version from Josef Leydold
 	 * */
 	u1 = unif_rand();
 	if(u1 < 0.884070402298758) {
 	    u2 = unif_rand();
 	    return A*(1.131131635444180*u1+u2-1);
 	}

 	if(u1 >= 0.973310954173898) { /* tail: */
 	    repeat {
 		u2 = unif_rand();
 		u3 = unif_rand();
 		tt = (A*A-2*log(u3));
 		if( u2*u2<(A*A)/tt )
 		    return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
 	    }
 	}

 	if(u1 >= 0.958720824790463) { /* region3: */
 	    repeat {
 		u2 = unif_rand();
 		u3 = unif_rand();
 		tt = A - 0.630834801921960* fmin2(u2,u3);
 		if(fmax2(u2,u3) <= 0.755591531667601)
 		    return (u2<u3) ? tt : -tt;
 		if(0.034240503750111*fabs(u2-u3) <= g(tt))
 		    return (u2<u3) ? tt : -tt;
 	    }
 	}

 	if(u1 >= 0.911312780288703) { /* region2: */
 	    repeat {
 		u2 = unif_rand();
 		u3 = unif_rand();
 		tt = 0.479727404222441+1.105473661022070*fmin2(u2,u3);
 		if( fmax2(u2,u3)<=0.872834976671790 )
 		    return (u2<u3) ? tt : -tt;
 		if( 0.049264496373128*fabs(u2-u3)<=g(tt) )
 		    return (u2<u3) ? tt : -tt;
 	    }
 	}

 	/* ELSE	 region1: */
 	repeat {
 	    u2 = unif_rand();
 	    u3 = unif_rand();
 	    tt = 0.479727404222441-0.595507138015940*fmin2(u2,u3);
 	    if (tt < 0.) continue;
 	    if(fmax2(u2,u3) <= 0.805577924423817)
 		return (u2<u3) ? tt : -tt;
      	    if(0.053377549506886*fabs(u2-u3) <= g(tt))
 		return (u2<u3) ? tt : -tt;
 	}
     default:
 	MATHLIB_ERROR(_("norm_rand(): invalid N01_kind: %d\n"), N01_kind)
 	    return 0.0;/*- -Wall */
     }/*switch*/
 }
	/*
	* Mathlib : A C Library of Special Functions
	* Copyright (C) 1998 Ross Ihaka
	* Copyright (C) 2000-9 The R Core Team
	*
	* This program is free software; you can redistribute it and/or modify
	* it under the terms of the GNU General Public License as published by
	* the Free Software Foundation; either version 2 of the License, or
	* (at your option) any later version.
	*
	* This program is distributed in the hope that it will be useful,
	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	* GNU General Public License for more details.
	*
	* You should have received a copy of the GNU General Public License
	* along with this program; if not, a copy is available at
	* https://www.R-project.org/Licenses/
	*
	* SYNOPSIS
	*
	* #include <Rmath.h>
	* double norm_rand(void);
	*
	* DESCRIPTION
	*
	* Random variates from the STANDARD normal distribution N(0,1).
	*
	* Is called from rnorm(..), but also rt(), rf(), rgamma(), ...
	*/

	#include <R_ext/Random.h>
	#include "nmath.h"

	#define repeat for(;;)

	#ifdef MATHLIB_STANDALONE
	static
	#else
	attribute_hidden
	#endif
	double BM_norm_keep = 0.0;

	N01type N01_kind = INVERSION;

	#ifndef MATHLIB_STANDALONE
	typedef void * (*DL_FUNC)();
	extern DL_FUNC User_norm_fun; /* declared and set in ../main/RNG.c */
	#endif

	/*
	* REFERENCE
	*
	* Ahrens, J.H. and Dieter, U.
	* Extensions of Forsythe's method for random sampling from
	* the normal distribution.
	* Math. Comput. 27, 927-937.
	*
	* The definitions of the constants a[k], d[k], t[k] and
	* h[k] are according to the abovementioned article
	*/
	double norm_rand(void)
	{

	const static double a[32] =
	{
	0.0000000, 0.03917609, 0.07841241, 0.1177699,
	0.1573107, 0.19709910, 0.23720210, 0.2776904,
	0.3186394, 0.36012990, 0.40225010, 0.4450965,
	0.4887764, 0.53340970, 0.57913220, 0.6260990,
	0.6744898, 0.72451440, 0.77642180, 0.8305109,
	0.8871466, 0.94678180, 1.00999000, 1.0775160,
	1.1503490, 1.22985900, 1.31801100, 1.4177970,
	1.5341210, 1.67594000, 1.86273200, 2.1538750
	};

	const static double d[31] =
	{
	0.0000000, 0.0000000, 0.0000000, 0.0000000,
	0.0000000, 0.2636843, 0.2425085, 0.2255674,
	0.2116342, 0.1999243, 0.1899108, 0.1812252,
	0.1736014, 0.1668419, 0.1607967, 0.1553497,
	0.1504094, 0.1459026, 0.1417700, 0.1379632,
	0.1344418, 0.1311722, 0.1281260, 0.1252791,
	0.1226109, 0.1201036, 0.1177417, 0.1155119,
	0.1134023, 0.1114027, 0.1095039
	};

	const static double t[31] =
	{
	7.673828e-4, 0.002306870, 0.003860618, 0.005438454,
	0.007050699, 0.008708396, 0.010423570, 0.012209530,
	0.014081250, 0.016055790, 0.018152900, 0.020395730,
	0.022811770, 0.025434070, 0.028302960, 0.031468220,
	0.034992330, 0.038954830, 0.043458780, 0.048640350,
	0.054683340, 0.061842220, 0.070479830, 0.081131950,
	0.094624440, 0.112300100, 0.136498000, 0.171688600,
	0.227624100, 0.330498000, 0.584703100
	};

	const static double h[31] =
	{
	0.03920617, 0.03932705, 0.03950999, 0.03975703,
	0.04007093, 0.04045533, 0.04091481, 0.04145507,
	0.04208311, 0.04280748, 0.04363863, 0.04458932,
	0.04567523, 0.04691571, 0.04833487, 0.04996298,
	0.05183859, 0.05401138, 0.05654656, 0.05953130,
	0.06308489, 0.06737503, 0.07264544, 0.07926471,
	0.08781922, 0.09930398, 0.11555990, 0.14043440,
	0.18361420, 0.27900160, 0.70104740
	};

	/----------- Constants and definitions for Kinderman - Ramage --- /
	/*
	* REFERENCE
	*
	* Kinderman A. J. and Ramage J. G. (1976).
	* Computer generation of normal random variables.
	* JASA 71, 893-896.
	*/

	#define C1 0.398942280401433
	#define C2 0.180025191068563
	#define g(x) (C1exp(-xx/2.0)-C2*(A-x))

	const static double A = 2.216035867166471;

	double s, u1, w, y, u2, u3, aa, tt, theta, R;
	int i;

	switch(N01_kind) {

	case AHRENS_DIETER: /* see Reference above */

	u1 = unif_rand();
	s = 0.0;
	if (u1 > 0.5)
	s = 1.0;
	u1 = u1 + u1 - s;
	u1 *= 32.0;
	i = (int) u1;
	if (i == 32)
	i = 31;
	if (i != 0) {
	u2 = u1 - i;
	aa = a[i - 1];
	while (u2 <= t[i - 1]) {
	u1 = unif_rand();
	w = u1 * (a[i] - aa);
	tt = (w * 0.5 + aa) * w;
	repeat {
	if (u2 > tt)
	goto deliver;
	u1 = unif_rand();
	if (u2 < u1)
	break;
	tt = u1;
	u2 = unif_rand();
	}
	u2 = unif_rand();
	}
	w = (u2 - t[i - 1]) * h[i - 1];
	}
	else {
	i = 6;
	aa = a[31];
	repeat {
	u1 = u1 + u1;
	if (u1 >= 1.0)
	break;
	aa = aa + d[i - 1];
	i = i + 1;
	}
	u1 = u1 - 1.0;
	repeat {
	w = u1 * d[i - 1];
	tt = (w * 0.5 + aa) * w;
	repeat {
	u2 = unif_rand();
	if (u2 > tt)
	goto jump;
	u1 = unif_rand();
	if (u2 < u1)
	break;
	tt = u1;
	}
	u1 = unif_rand();
	}
	jump:;
	}

	deliver:
	y = aa + w;
	return (s == 1.0) ? -y : y;

	/-----------------------------------------------------------/

	case BUGGY_KINDERMAN_RAMAGE: /* see Reference above */
	/* note: this has problems, but is retained for
	* reproducibility of older codes, with the same
	* numeric code */
	u1 = unif_rand();
	if(u1 < 0.884070402298758) {
	u2 = unif_rand();
	return A(1.13113163544180u1+u2-1);
	}

	if(u1 >= 0.973310954173898) { /* tail: */
	repeat {
	u2 = unif_rand();
	u3 = unif_rand();
	tt = (AA-2log(u3));
	if( u2u2<(AA)/tt )
	return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
	}
	}

	if(u1 >= 0.958720824790463) { /* region3: */
	repeat {
	u2 = unif_rand();
	u3 = unif_rand();
	tt = A - 0.630834801921960* fmin2(u2,u3);
	if(fmax2(u2,u3) <= 0.755591531667601)
	return (u2<u3) ? tt : -tt;
	if(0.034240503750111*fabs(u2-u3) <= g(tt))
	return (u2<u3) ? tt : -tt;
	}
	}

	if(u1 >= 0.911312780288703) { /* region2: */
	repeat {
	u2 = unif_rand();
	u3 = unif_rand();
	tt = 0.479727404222441+1.105473661022070*fmin2(u2,u3);
	if( fmax2(u2,u3)<=0.872834976671790 )
	return (u2<u3) ? tt : -tt;
	if( 0.049264496373128*fabs(u2-u3)<=g(tt) )
	return (u2<u3) ? tt : -tt;
	}
	}

	/* ELSE region1: */
	repeat {
	u2 = unif_rand();
	u3 = unif_rand();
	tt = 0.479727404222441-0.595507138015940*fmin2(u2,u3);
	if(fmax2(u2,u3) <= 0.805577924423817)
	return (u2<u3) ? tt : -tt;
	}
	case BOX_MULLER:
	if(BM_norm_keep != 0.0) { /* An exact test is intentional */
	s = BM_norm_keep;
	BM_norm_keep = 0.0;
	return s;
	} else {
	theta = 2 * M_PI * unif_rand();
	R = sqrt(-2 * log(unif_rand())) + 10DBL_MIN; / ensure non-zero */
	BM_norm_keep = R * sin(theta);
	return R * cos(theta);
	}
	#ifndef MATHLIB_STANDALONE
	case USER_NORM:
	return ((double ) User_norm_fun());
	#endif
	case INVERSION:
	#define BIG 134217728 /* 2^27 */
	/* unif_rand() alone is not of high enough precision */
	u1 = unif_rand();
	u1 = (int)(BIG*u1) + unif_rand();
	return qnorm5(u1/BIG, 0.0, 1.0, 1, 0);
	case KINDERMAN_RAMAGE: /* see Reference above */
	/* corrected version from Josef Leydold
	* */
	u1 = unif_rand();
	if(u1 < 0.884070402298758) {
	u2 = unif_rand();
	return A(1.131131635444180u1+u2-1);
	}

	if(u1 >= 0.973310954173898) { /* tail: */
	repeat {
	u2 = unif_rand();
	u3 = unif_rand();
	tt = (AA-2log(u3));
	if( u2u2<(AA)/tt )
	return (u1 < 0.986655477086949) ? sqrt(tt) : -sqrt(tt);
	}
	}

	if(u1 >= 0.958720824790463) { /* region3: */
	repeat {
	u2 = unif_rand();
	u3 = unif_rand();
	tt = A - 0.630834801921960* fmin2(u2,u3);
	if(fmax2(u2,u3) <= 0.755591531667601)
	return (u2<u3) ? tt : -tt;
	if(0.034240503750111*fabs(u2-u3) <= g(tt))
	return (u2<u3) ? tt : -tt;
	}
	}

	if(u1 >= 0.911312780288703) { /* region2: */
	repeat {
	u2 = unif_rand();
	u3 = unif_rand();
	tt = 0.479727404222441+1.105473661022070*fmin2(u2,u3);
	if( fmax2(u2,u3)<=0.872834976671790 )
	return (u2<u3) ? tt : -tt;
	if( 0.049264496373128*fabs(u2-u3)<=g(tt) )
	return (u2<u3) ? tt : -tt;
	}
	}

	/* ELSE region1: */
	repeat {
	u2 = unif_rand();
	u3 = unif_rand();
	tt = 0.479727404222441-0.595507138015940*fmin2(u2,u3);
	if (tt < 0.) continue;
	if(fmax2(u2,u3) <= 0.805577924423817)
	return (u2<u3) ? tt : -tt;
	if(0.053377549506886*fabs(u2-u3) <= g(tt))
	return (u2<u3) ? tt : -tt;
	}
	default:
	MATHLIB_ERROR(_("norm_rand(): invalid N01_kind: %d\n"), N01_kind)
	return 0.0;/- -Wall /
	}/switch/
	}