home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 December / simtel1292_SIMTEL_1292_Walnut_Creek.iso / msdos / fortran / random.arc / RANDOM.FOR < prev   
Text File  |  1989-03-12  |  5KB  |  162 lines
  1. C This random number generator originally appeared in "Toward a Universal 
  2. C Random Number Generator" by George Marsaglia and Arif Zaman. 
  3. C Florida State University Report: FSU-SCRI-87-50 (1987)
  4. C 
  5. C It was later modified by F. James and published in "A Review of Pseudo-
  6. C random Number Generators" 
  7. C 
  8. C THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.
  9. C       (However, a newly discovered technique can yield 
  10. C         a period of 10^600. But that is still in the development stage.)
  11. C
  12. C It passes ALL of the tests for random number generators and has a period 
  13. C   of 2^144, is completely portable (gives bit identical results on all 
  14. C   machines with at least 24-bit mantissas in the floating point 
  15. C   representation). 
  16. C 
  17. C The algorithm is a combination of a Fibonacci sequence (with lags of 97
  18. C   and 33, and operation "subtraction plus one, modulo one") and an 
  19. C   "arithmetic sequence" (using subtraction).
  20. C
  21. C On a Vax 11/780, this random number generator can produce a number in 
  22. C    13 microseconds. 
  23. C======================================================================== 
  24.  
  25.       PROGRAM TstRAN
  26.       REAL temp(100)
  27.       INTEGER IJ, KL, len
  28. C These are the seeds needed to produce the test case results
  29.       IJ = 1802
  30.       KL = 9373
  31.  
  32.  
  33. C Do the initialization
  34.       call rmarin(ij,kl)
  35.  
  36. C Generate 20000 random numbers
  37.       len = 100
  38.       do 10 i = 1, 200
  39.          call RANMAR(temp, len)
  40. 10    continue
  41.  
  42. C If the random number generator is working properly, the next six random
  43. C numbers should be:
  44. C           6533892.0  14220222.0  7275067.0
  45. C           6172232.0  8354498.0   10633180.0
  46.  
  47.       len = 6
  48.       call RANMAR(temp, len)
  49.  
  50.         
  51.       write(6,20) (4096.0*4096.0*temp(I), I=1,6)
  52. 20    format (3f12.1)
  53.       end
  54.  
  55.       subroutine RMARIN(IJ,KL)
  56. C This is the initialization routine for the random number generator RANMAR()
  57. C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
  58. C                                                      0 <= KL <= 30081
  59. C The random number sequences created by these two seeds are of sufficient 
  60. C length to complete an entire calculation with. For example, if sveral 
  61. C different groups are working on different parts of the same calculation,
  62. C each group could be assigned its own IJ seed. This would leave each group
  63. C with 30000 choices for the second seed. That is to say, this random 
  64. C number generator can create 900 million different subsequences -- with 
  65. C each subsequence having a length of approximately 10^30.
  66. C 
  67. C Use IJ = 1802 & KL = 9373 to test the random number generator. The
  68. C subroutine RANMAR should be used to generate 20000 random numbers.
  69. C Then display the next six random numbers generated multiplied by 4096*4096
  70. C If the random number generator is working properly, the random numbers
  71. C should be:
  72. C           6533892.0  14220222.0  7275067.0
  73. C           6172232.0  8354498.0   10633180.0
  74.  
  75.  
  76.       real U(97), C, CD, CM
  77.       integer I97, J97
  78.       logical TEST
  79.       data TEST /.FALSE./
  80.       common /raset1/ U, C, CD, CM, I97, J97, TEST
  81.  
  82.       if( IJ .lt. 0  .or.  IJ .gt. 31328  .or.
  83.      *    KL .lt. 0  .or.  KL .gt. 30081 ) then
  84.           print '(A)', ' The first random number seed must have a value 
  85.      *between 0 and 31328'
  86.           print '(A)',' The second seed must have a value between 0 and         
  87.      *30081'
  88.             stop
  89.       endif
  90.  
  91.       i = mod(IJ/177, 177) + 2
  92.       j = mod(IJ    , 177) + 2
  93.       k = mod(KL/169, 178) + 1
  94.       l = mod(KL,     169) 
  95.  
  96.       do 2 ii = 1, 97
  97.          s = 0.0
  98.          t = 0.5
  99.          do 3 jj = 1, 24
  100.             m = mod(mod(i*j, 179)*k, 179)
  101.             i = j
  102.             j = k
  103.             k = m
  104.             l = mod(53*l+1, 169)
  105.             if (mod(l*m, 64) .ge. 32) then
  106.                s = s + t
  107.             endif
  108.             t = 0.5 * t
  109. 3        continue
  110.          U(ii) = s
  111. 2     continue
  112.  
  113.       C = 362436.0 / 16777216.0
  114.       CD = 7654321.0 / 16777216.0
  115.       CM = 16777213.0 /16777216.0
  116.  
  117.       I97 = 97
  118.       J97 = 33
  119.  
  120.       TEST = .TRUE.
  121.       return
  122.       end
  123.  
  124.       subroutine ranmar(RVEC, LEN)
  125. C This is the random number generator proposed by George Marsaglia in 
  126. C Florida State University Report: FSU-SCRI-87-50
  127. C It was slightly modified by F. James to produce an array of pseudorandom
  128. C numbers.
  129.       REAL RVEC(*)
  130.       real U(97), C, CD, CM
  131.       integer I97, J97
  132.       logical TEST
  133.       common /raset1/ U, C, CD, CM, I97, J97, TEST
  134.  
  135.       integer ivec
  136.  
  137.       if( .NOT. TEST ) then
  138.          print '(A)',' Call the init routine (RMARIN) before calling RAN        
  139.      *MAR'  
  140.          stop
  141.       endif
  142.  
  143.       do 100 ivec = 1, LEN
  144.          uni = U(I97) - U(J97)
  145.          if( uni .lt. 0.0 ) uni = uni + 1.0
  146.          U(I97) = uni
  147.          I97 = I97 - 1
  148.          if(I97 .eq. 0) I97 = 97
  149.          J97 = J97 - 1
  150.          if(J97 .eq. 0) J97 = 97
  151.          C = C - CD
  152.          if( C .lt. 0.0 ) C = C + CM
  153.          uni = uni - C
  154.          if( uni .lt. 0.0 ) uni = uni + 1.0
  155.          RVEC(ivec) = uni
  156. 100   continue
  157.       return
  158.       end
  159.  
  160. ------------------------------
  161.  
  162.