GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
snorm.f
Go to the documentation of this file.
1  REAL FUNCTION snorm()
2 C**********************************************************************C
3 C C
4 C C
5 C (STANDARD-) N O R M A L DISTRIBUTION C
6 C C
7 C C
8 C**********************************************************************C
9 C**********************************************************************C
10 C C
11 C FOR DETAILS SEE: C
12 C C
13 C AHRENS, J.H. AND DIETER, U. C
14 C EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM C
15 C SAMPLING FROM THE NORMAL DISTRIBUTION. C
16 C MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937. C
17 C C
18 C ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL' C
19 C (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) C
20 C C
21 C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C
22 C SUNIF. The argument IR thus goes away. C
23 C C
24 C**********************************************************************C
25 C
26 C
27 C THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
28 C H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
29 C
30 C .. Local Scalars ..
31  REAL aa,s,tt,u,ustar,w,y
32  INTEGER i
33 C ..
34 C .. Local Arrays ..
35  REAL a(32),d(31),h(31),t(31)
36 C ..
37 C .. External Functions ..
38  REAL ranf
39  EXTERNAL ranf
40 C ..
41 C .. Intrinsic Functions ..
42  INTRINSIC float,int
43 C ..
44 C .. Save statement ..
45 C JJV added a Save statement for arrays initialized in Data statmts
46  SAVE a,d,t,h
47 C ..
48 C .. Data statements ..
49  DATA a/0.0,.3917609e-1,.7841241e-1,.1177699,.1573107,.1970991,
50  + .2372021,.2776904,.3186394,.3601299,.4022501,.4450965,
51  + .4887764,.5334097,.5791322,.6260990,.6744898,.7245144,
52  + .7764218,.8305109,.8871466,.9467818,1.009990,1.077516,
53  + 1.150349,1.229859,1.318011,1.417797,1.534121,1.675940,
54  + 1.862732,2.153875/
55  DATA d/5*0.0,.2636843,.2425085,.2255674,.2116342,.1999243,
56  + .1899108,.1812252,.1736014,.1668419,.1607967,.1553497,
57  + .1504094,.1459026,.1417700,.1379632,.1344418,.1311722,
58  + .1281260,.1252791,.1226109,.1201036,.1177417,.1155119,
59  + .1134023,.1114027,.1095039/
60  DATA t/.7673828e-3,.2306870e-2,.3860618e-2,.5438454e-2,
61  + .7050699e-2,.8708396e-2,.1042357e-1,.1220953e-1,.1408125e-1,
62  + .1605579e-1,.1815290e-1,.2039573e-1,.2281177e-1,.2543407e-1,
63  + .2830296e-1,.3146822e-1,.3499233e-1,.3895483e-1,.4345878e-1,
64  + .4864035e-1,.5468334e-1,.6184222e-1,.7047983e-1,.8113195e-1,
65  + .9462444e-1,.1123001,.1364980,.1716886,.2276241,.3304980,
66  + .5847031/
67  DATA h/.3920617e-1,.3932705e-1,.3950999e-1,.3975703e-1,
68  + .4007093e-1,.4045533e-1,.4091481e-1,.4145507e-1,.4208311e-1,
69  + .4280748e-1,.4363863e-1,.4458932e-1,.4567523e-1,.4691571e-1,
70  + .4833487e-1,.4996298e-1,.5183859e-1,.5401138e-1,.5654656e-1,
71  + .5953130e-1,.6308489e-1,.6737503e-1,.7264544e-1,.7926471e-1,
72  + .8781922e-1,.9930398e-1,.1155599,.1404344,.1836142,.2790016,
73  + .7010474/
74 C ..
75 C .. Executable Statements ..
76 C
77  10 u = ranf()
78  s = 0.0
79  IF (u.GT.0.5) s = 1.0
80  u = u + u - s
81  20 u = 32.0*u
82  i = int(u)
83  IF (i.EQ.32) i = 31
84  IF (i.EQ.0) go to 100
85 C
86 C START CENTER
87 C
88  30 ustar = u - float(i)
89  aa = a(i)
90  40 IF (ustar.LE.t(i)) go to 60
91  w = (ustar-t(i))*h(i)
92 C
93 C EXIT (BOTH CASES)
94 C
95  50 y = aa + w
96  snorm = y
97  IF (s.EQ.1.0) snorm = -y
98  RETURN
99 C
100 C CENTER CONTINUED
101 C
102  60 u = ranf()
103  w = u* (a(i+1)-aa)
104  tt = (0.5*w+aa)*w
105  go to 80
106 
107  70 tt = u
108  ustar = ranf()
109  80 IF (ustar.GT.tt) go to 50
110  90 u = ranf()
111  IF (ustar.GE.u) go to 70
112  ustar = ranf()
113  go to 40
114 C
115 C START TAIL
116 C
117  100 i = 6
118  aa = a(32)
119  go to 120
120 
121  110 aa = aa + d(i)
122  i = i + 1
123  120 u = u + u
124  IF (u.LT.1.0) go to 110
125  130 u = u - 1.0
126  140 w = u*d(i)
127  tt = (0.5*w+aa)*w
128  go to 160
129 
130  150 tt = u
131  160 ustar = ranf()
132  IF (ustar.GT.tt) go to 50
133  170 u = ranf()
134  IF (ustar.GE.u) go to 150
135  u = ranf()
136  go to 140
137 
138  END