Mercurial > hg > octave-lyh

annotate liboctave/cruft/ranlib/snorm.f @ 17428:577a19afdaf5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Backed out changeset f81401b6b1f7. * test/parser.tst: bug #33304 should remain an error until it is fixed.
author Rik <rik@octave.org>
date Sat, 14 Sep 2013 17:07:05 -0700
parents 648dabbb4c6b
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
1 REAL FUNCTION snorm()
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
2 C**********************************************************************C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
3 C C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
4 C C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
5 C (STANDARD-) N O R M A L DISTRIBUTION C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
6 C C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
7 C C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
8 C**********************************************************************C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
9 C**********************************************************************C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
10 C C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
11 C FOR DETAILS SEE: C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
12 C C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
13 C AHRENS, J.H. AND DIETER, U. C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
14 C EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
15 C SAMPLING FROM THE NORMAL DISTRIBUTION. C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
16 C MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937. C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
17 C C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
18 C ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL' C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
19 C (M=5) IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION) C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
20 C C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
21 C Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
22 C SUNIF. The argument IR thus goes away. C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
23 C C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
24 C**********************************************************************C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
25 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
26 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
27 C THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
28 C H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
29 C
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
30 C .. Local Scalars ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
31 REAL aa,s,tt,u,ustar,w,y
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
32 INTEGER i
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
33 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
34 C .. Local Arrays ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
35 REAL a(32),d(31),h(31),t(31)
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
36 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
37 C .. External Functions ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
38 REAL ranf
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
39 EXTERNAL ranf
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
40 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
41 C .. Intrinsic Functions ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
42 INTRINSIC float,int
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
43 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
44 C .. Save statement ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
45 C JJV added a Save statement for arrays initialized in Data statmts
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
46 SAVE a,d,t,h
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
47 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
48 C .. Data statements ..
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
49 DATA a/0.0,.3917609E-1,.7841241E-1,.1177699,.1573107,.1970991,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
50 + .2372021,.2776904,.3186394,.3601299,.4022501,.4450965,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
51 + .4887764,.5334097,.5791322,.6260990,.6744898,.7245144,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
52 + .7764218,.8305109,.8871466,.9467818,1.009990,1.077516,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
53 + 1.150349,1.229859,1.318011,1.417797,1.534121,1.675940,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
54 + 1.862732,2.153875/
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
55 DATA d/5*0.0,.2636843,.2425085,.2255674,.2116342,.1999243,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
56 + .1899108,.1812252,.1736014,.1668419,.1607967,.1553497,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
57 + .1504094,.1459026,.1417700,.1379632,.1344418,.1311722,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
58 + .1281260,.1252791,.1226109,.1201036,.1177417,.1155119,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
59 + .1134023,.1114027,.1095039/
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
60 DATA t/.7673828E-3,.2306870E-2,.3860618E-2,.5438454E-2,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
61 + .7050699E-2,.8708396E-2,.1042357E-1,.1220953E-1,.1408125E-1,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
62 + .1605579E-1,.1815290E-1,.2039573E-1,.2281177E-1,.2543407E-1,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
63 + .2830296E-1,.3146822E-1,.3499233E-1,.3895483E-1,.4345878E-1,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
64 + .4864035E-1,.5468334E-1,.6184222E-1,.7047983E-1,.8113195E-1,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
65 + .9462444E-1,.1123001,.1364980,.1716886,.2276241,.3304980,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
66 + .5847031/
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
67 DATA h/.3920617E-1,.3932705E-1,.3950999E-1,.3975703E-1,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
68 + .4007093E-1,.4045533E-1,.4091481E-1,.4145507E-1,.4208311E-1,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
69 + .4280748E-1,.4363863E-1,.4458932E-1,.4567523E-1,.4691571E-1,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
70 + .4833487E-1,.4996298E-1,.5183859E-1,.5401138E-1,.5654656E-1,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
71 + .5953130E-1,.6308489E-1,.6737503E-1,.7264544E-1,.7926471E-1,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
72 + .8781922E-1,.9930398E-1,.1155599,.1404344,.1836142,.2790016,
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
73 + .7010474/
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
74 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents: 2329
diff changeset
75 C .. Executable Statements ..
2329
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
76 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
77 10 u = ranf()
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
78 s = 0.0
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
79 IF (u.GT.0.5) s = 1.0
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
80 u = u + u - s
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
81 20 u = 32.0*u
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
82 i = int(u)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
83 IF (i.EQ.32) i = 31
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
84 IF (i.EQ.0) GO TO 100
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
85 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
86 C START CENTER
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
87 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
88 30 ustar = u - float(i)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
89 aa = a(i)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
90 40 IF (ustar.LE.t(i)) GO TO 60
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
91 w = (ustar-t(i))*h(i)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
92 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
93 C EXIT (BOTH CASES)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
94 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
95 50 y = aa + w
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
96 snorm = y
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
97 IF (s.EQ.1.0) snorm = -y
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
98 RETURN
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
99 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
100 C CENTER CONTINUED
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
101 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
102 60 u = ranf()
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
103 w = u* (a(i+1)-aa)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
104 tt = (0.5*w+aa)*w
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
105 GO TO 80
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
106
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
107 70 tt = u
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
108 ustar = ranf()
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
109 80 IF (ustar.GT.tt) GO TO 50
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
110 90 u = ranf()
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
111 IF (ustar.GE.u) GO TO 70
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
112 ustar = ranf()
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
113 GO TO 40
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
114 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
115 C START TAIL
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
116 C
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
117 100 i = 6
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
118 aa = a(32)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
119 GO TO 120
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
120
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
121 110 aa = aa + d(i)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
122 i = i + 1
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
123 120 u = u + u
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
124 IF (u.LT.1.0) GO TO 110
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
125 130 u = u - 1.0
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
126 140 w = u*d(i)
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
127 tt = (0.5*w+aa)*w
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
128 GO TO 160
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
129
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
130 150 tt = u
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
131 160 ustar = ranf()
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
132 IF (ustar.GT.tt) GO TO 50
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
133 170 u = ranf()
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
134 IF (ustar.GE.u) GO TO 150
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
135 u = ranf()
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
136 GO TO 140
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
137
30c606bec7a8 [project @ 1996-07-19 01:29:05 by jwe]
jwe
parents:
diff changeset
138 END