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 / sexpo.f < prev    next >
Text File  |  1996-09-28  |  2KB  |  58 lines
  1.       REAL FUNCTION sexpo()
  2. C**********************************************************************C
  3. C                                                                      C
  4. C                                                                      C
  5. C     (STANDARD-)  E X P O N E N T I 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               COMPUTER METHODS FOR SAMPLING FROM THE                 C
  15. C               EXPONENTIAL AND NORMAL DISTRIBUTIONS.                  C
  16. C               COMM. ACM, 15,10 (OCT. 1972), 873 - 882.               C
  17. C                                                                      C
  18. C     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM       C
  19. C     'SA' 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 q(8)
  27.       EQUIVALENCE (q(1),q1)
  28. C
  29. C     Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
  30. C     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
  31. C
  32.       DATA q/.6931472,.9333737,.9888778,.9984959,.9998293,.9999833,
  33.      +     .9999986,.9999999/
  34. C
  35.    10 a = 0.0
  36.       u = ranf()
  37.       GO TO 30
  38.  
  39.    20 a = a + q1
  40.    30 u = u + u
  41.       IF (u.LE.1.0) GO TO 20
  42.    40 u = u - 1.0
  43.       IF (u.GT.q1) GO TO 60
  44.    50 sexpo = a + u
  45.       RETURN
  46.  
  47.    60 i = 1
  48.       ustar = ranf()
  49.       umin = ustar
  50.    70 ustar = ranf()
  51.       IF (ustar.LT.umin) umin = ustar
  52.    80 i = i + 1
  53.       IF (u.GT.q(i)) GO TO 70
  54.    90 sexpo = a + umin*q1
  55.       RETURN
  56.  
  57.       END
  58.