GNU Octave: libcruft/ranlib/snorm.f Source File

Go to the documentation of this file.
00001       REAL FUNCTION snorm()
00002 C**********************************************************************C
00003 C                                                                      C
00004 C                                                                      C
00005 C     (STANDARD-)  N O R M A L  DISTRIBUTION                           C
00006 C                                                                      C
00007 C                                                                      C
00008 C**********************************************************************C
00009 C**********************************************************************C
00010 C                                                                      C
00011 C     FOR DETAILS SEE:                                                 C
00012 C                                                                      C
00013 C               AHRENS, J.H. AND DIETER, U.                            C
00014 C               EXTENSIONS OF FORSYTHE'S METHOD FOR RANDOM             C
00015 C               SAMPLING FROM THE NORMAL DISTRIBUTION.                 C
00016 C               MATH. COMPUT., 27,124 (OCT. 1973), 927 - 937.          C
00017 C                                                                      C
00018 C     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM 'FL'  C
00019 C     (M=5) IN THE ABOVE PAPER     (SLIGHTLY MODIFIED IMPLEMENTATION)  C
00020 C                                                                      C
00021 C     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   C
00022 C     SUNIF.  The argument IR thus goes away.                          C
00023 C                                                                      C
00024 C**********************************************************************C
00025 C
00026 C
00027 C     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
00028 C     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
00029 C
00030 C     .. Local Scalars ..
00031       REAL aa,s,tt,u,ustar,w,y
00032       INTEGER i
00033 C     ..
00034 C     .. Local Arrays ..
00035       REAL a(32),d(31),h(31),t(31)
00036 C     ..
00037 C     .. External Functions ..
00038       REAL ranf
00039       EXTERNAL ranf
00040 C     ..
00041 C     .. Intrinsic Functions ..
00042       INTRINSIC float,int
00043 C     ..
00044 C     .. Save statement ..
00045 C     JJV added a Save statement for arrays initialized in Data statmts
00046       SAVE a,d,t,h
00047 C     ..
00048 C     .. Data statements ..
00049       DATA a/0.0,.3917609E-1,.7841241E-1,.1177699,.1573107,.1970991,
00050      +     .2372021,.2776904,.3186394,.3601299,.4022501,.4450965,
00051      +     .4887764,.5334097,.5791322,.6260990,.6744898,.7245144,
00052      +     .7764218,.8305109,.8871466,.9467818,1.009990,1.077516,
00053      +     1.150349,1.229859,1.318011,1.417797,1.534121,1.675940,
00054      +     1.862732,2.153875/
00055       DATA d/5*0.0,.2636843,.2425085,.2255674,.2116342,.1999243,
00056      +     .1899108,.1812252,.1736014,.1668419,.1607967,.1553497,
00057      +     .1504094,.1459026,.1417700,.1379632,.1344418,.1311722,
00058      +     .1281260,.1252791,.1226109,.1201036,.1177417,.1155119,
00059      +     .1134023,.1114027,.1095039/
00060       DATA t/.7673828E-3,.2306870E-2,.3860618E-2,.5438454E-2,
00061      +     .7050699E-2,.8708396E-2,.1042357E-1,.1220953E-1,.1408125E-1,
00062      +     .1605579E-1,.1815290E-1,.2039573E-1,.2281177E-1,.2543407E-1,
00063      +     .2830296E-1,.3146822E-1,.3499233E-1,.3895483E-1,.4345878E-1,
00064      +     .4864035E-1,.5468334E-1,.6184222E-1,.7047983E-1,.8113195E-1,
00065      +     .9462444E-1,.1123001,.1364980,.1716886,.2276241,.3304980,
00066      +     .5847031/
00067       DATA h/.3920617E-1,.3932705E-1,.3950999E-1,.3975703E-1,
00068      +     .4007093E-1,.4045533E-1,.4091481E-1,.4145507E-1,.4208311E-1,
00069      +     .4280748E-1,.4363863E-1,.4458932E-1,.4567523E-1,.4691571E-1,
00070      +     .4833487E-1,.4996298E-1,.5183859E-1,.5401138E-1,.5654656E-1,
00071      +     .5953130E-1,.6308489E-1,.6737503E-1,.7264544E-1,.7926471E-1,
00072      +     .8781922E-1,.9930398E-1,.1155599,.1404344,.1836142,.2790016,
00073      +     .7010474/
00074 C     ..
00075 C     .. Executable Statements ..
00076 C
00077    10 u = ranf()
00078       s = 0.0
00079       IF (u.GT.0.5) s = 1.0
00080       u = u + u - s
00081    20 u = 32.0*u
00082       i = int(u)
00083       IF (i.EQ.32) i = 31
00084       IF (i.EQ.0) GO TO 100
00085 C
00086 C                                START CENTER
00087 C
00088    30 ustar = u - float(i)
00089       aa = a(i)
00090    40 IF (ustar.LE.t(i)) GO TO 60
00091       w = (ustar-t(i))*h(i)
00092 C
00093 C                                EXIT   (BOTH CASES)
00094 C
00095    50 y = aa + w
00096       snorm = y
00097       IF (s.EQ.1.0) snorm = -y
00098       RETURN
00099 C
00100 C                                CENTER CONTINUED
00101 C
00102    60 u = ranf()
00103       w = u* (a(i+1)-aa)
00104       tt = (0.5*w+aa)*w
00105       GO TO 80
00106 
00107    70 tt = u
00108       ustar = ranf()
00109    80 IF (ustar.GT.tt) GO TO 50
00110    90 u = ranf()
00111       IF (ustar.GE.u) GO TO 70
00112       ustar = ranf()
00113       GO TO 40
00114 C
00115 C                                START TAIL
00116 C
00117   100 i = 6
00118       aa = a(32)
00119       GO TO 120
00120 
00121   110 aa = aa + d(i)
00122       i = i + 1
00123   120 u = u + u
00124       IF (u.LT.1.0) GO TO 110
00125   130 u = u - 1.0
00126   140 w = u*d(i)
00127       tt = (0.5*w+aa)*w
00128       GO TO 160
00129 
00130   150 tt = u
00131   160 ustar = ranf()
00132       IF (ustar.GT.tt) GO TO 50
00133   170 u = ranf()
00134       IF (ustar.GE.u) GO TO 150
00135       u = ranf()
00136       GO TO 140
00137 
00138       END