home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / ranlib / snorm.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  119 lines
  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.       DIMENSION a(32),d(31),t(31),h(31)
  27. C
  28. C     THE DEFINITIONS OF THE CONSTANTS A(K), D(K), T(K) AND
  29. C     H(K) ARE ACCORDING TO THE ABOVEMENTIONED ARTICLE
  30. C
  31.       DATA a/0.0,.3917609E-1,.7841241E-1,.1177699,.1573107,.1970991,
  32.      +     .2372021,.2776904,.3186394,.3601299,.4022501,.4450965,
  33.      +     .4887764,.5334097,.5791322,.6260990,.6744898,.7245144,
  34.      +     .7764218,.8305109,.8871466,.9467818,1.009990,1.077516,
  35.      +     1.150349,1.229859,1.318011,1.417797,1.534121,1.675940,
  36.      +     1.862732,2.153875/
  37.       DATA d/5*0.0,.2636843,.2425085,.2255674,.2116342,.1999243,
  38.      +     .1899108,.1812252,.1736014,.1668419,.1607967,.1553497,
  39.      +     .1504094,.1459026,.1417700,.1379632,.1344418,.1311722,
  40.      +     .1281260,.1252791,.1226109,.1201036,.1177417,.1155119,
  41.      +     .1134023,.1114027,.1095039/
  42.       DATA t/.7673828E-3,.2306870E-2,.3860618E-2,.5438454E-2,
  43.      +     .7050699E-2,.8708396E-2,.1042357E-1,.1220953E-1,.1408125E-1,
  44.      +     .1605579E-1,.1815290E-1,.2039573E-1,.2281177E-1,.2543407E-1,
  45.      +     .2830296E-1,.3146822E-1,.3499233E-1,.3895483E-1,.4345878E-1,
  46.      +     .4864035E-1,.5468334E-1,.6184222E-1,.7047983E-1,.8113195E-1,
  47.      +     .9462444E-1,.1123001,.1364980,.1716886,.2276241,.3304980,
  48.      +     .5847031/
  49.       DATA h/.3920617E-1,.3932705E-1,.3950999E-1,.3975703E-1,
  50.      +     .4007093E-1,.4045533E-1,.4091481E-1,.4145507E-1,.4208311E-1,
  51.      +     .4280748E-1,.4363863E-1,.4458932E-1,.4567523E-1,.4691571E-1,
  52.      +     .4833487E-1,.4996298E-1,.5183859E-1,.5401138E-1,.5654656E-1,
  53.      +     .5953130E-1,.6308489E-1,.6737503E-1,.7264544E-1,.7926471E-1,
  54.      +     .8781922E-1,.9930398E-1,.1155599,.1404344,.1836142,.2790016,
  55.      +     .7010474/
  56. C
  57.    10 u = ranf()
  58.       s = 0.0
  59.       IF (u.GT.0.5) s = 1.0
  60.       u = u + u - s
  61.    20 u = 32.0*u
  62.       i = int(u)
  63.       IF (i.EQ.32) i = 31
  64.       IF (i.EQ.0) GO TO 100
  65. C
  66. C                                START CENTER
  67. C
  68.    30 ustar = u - float(i)
  69.       aa = a(i)
  70.    40 IF (ustar.LE.t(i)) GO TO 60
  71.       w = (ustar-t(i))*h(i)
  72. C
  73. C                                EXIT   (BOTH CASES)
  74. C
  75.    50 y = aa + w
  76.       snorm = y
  77.       IF (s.EQ.1.0) snorm = -y
  78.       RETURN
  79. C
  80. C                                CENTER CONTINUED
  81. C
  82.    60 u = ranf()
  83.       w = u* (a(i+1)-aa)
  84.       tt = (0.5*w+aa)*w
  85.       GO TO 80
  86.  
  87.    70 tt = u
  88.       ustar = ranf()
  89.    80 IF (ustar.GT.tt) GO TO 50
  90.    90 u = ranf()
  91.       IF (ustar.GE.u) GO TO 70
  92.       ustar = ranf()
  93.       GO TO 40
  94. C
  95. C                                START TAIL
  96. C
  97.   100 i = 6
  98.       aa = a(32)
  99.       GO TO 120
  100.  
  101.   110 aa = aa + d(i)
  102.       i = i + 1
  103.   120 u = u + u
  104.       IF (u.LT.1.0) GO TO 110
  105.   130 u = u - 1.0
  106.   140 w = u*d(i)
  107.       tt = (0.5*w+aa)*w
  108.       GO TO 160
  109.  
  110.   150 tt = u
  111.   160 ustar = ranf()
  112.       IF (ustar.GT.tt) GO TO 50
  113.   170 u = ranf()
  114.       IF (ustar.GE.u) GO TO 150
  115.       u = ranf()
  116.       GO TO 140
  117.  
  118.       END
  119.