GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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*4 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
ranf
real function ranf()
Definition: ranf.f:2
snorm
real function snorm()
Definition: snorm.f:2