| /* |
| * 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*/ |
| } |