home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / mt19937.zip / mt19937.f < prev    next >
Text File  |  2003-01-13  |  6KB  |  163 lines
  1. * A C-program for MT19937: Real number version
  2. *   genrand() generates one pseudorandom real number (double)
  3. * which is uniformly distributed on [0,1]-interval, for each
  4. * call. sgenrand(seed) set initial values to the working area
  5. * of 624 words. Before genrand(), sgenrand(seed) must be
  6. * called once. (seed is any 32-bit integer except for 0).
  7. * Integer generator is obtained by modifying two lines.
  8. *   Coded by Takuji Nishimura, considering the suggestions by
  9. * Topher Cooper and Marc Rieffel in July-Aug. 1997.
  10. *
  11. * This library is free software; you can redistribute it and/or
  12. * modify it under the terms of the GNU Library General Public
  13. * License as published by the Free Software Foundation; either
  14. * version 2 of the License, or (at your option) any later
  15. * version.
  16. * This library is distributed in the hope that it will be useful,
  17. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  19. * See the GNU Library General Public License for more details.
  20. * You should have received a copy of the GNU Library General
  21. * Public License along with this library; if not, write to the
  22. * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  23. * 02111-1307  USA
  24. *
  25. * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
  26. * When you use this, send an email to: matumoto@math.keio.ac.jp
  27. * with an appropriate reference to your work.
  28. *
  29. ************************************************************************
  30. * Fortran translation by Hiroshi Takano.  Jan. 13, 1999.
  31. *
  32. *   genrand()      -> double precision function grnd()
  33. *   sgenrand(seed) -> subroutine sgrnd(seed)
  34. *                     integer seed
  35. *
  36. * This program uses the following non-standard intrinsics.
  37. *   ishft(i,n): If n>0, shifts bits in i by n positions to left.
  38. *               If n<0, shifts bits in i by n positions to right.
  39. *   iand (i,j): Performs logical AND on corresponding bits of i and j.
  40. *   ior  (i,j): Performs inclusive OR on corresponding bits of i and j.
  41. *   ieor (i,j): Performs exclusive OR on corresponding bits of i and j.
  42. *
  43. ************************************************************************
  44. * this main() outputs first 1000 generated numbers
  45.       program main
  46. *
  47.       implicit integer(i-n)
  48.       implicit double precision(a-h,o-z)
  49. *
  50.       parameter(no=1000)
  51.       dimension r(0:7)
  52. *
  53. *      call sgrnd(4357)
  54. *                         any nonzero integer can be used as a seed
  55.       do 1000 j=0,no-1
  56.         r(mod(j,8))=grnd()
  57.         if(mod(j,8).eq.7) then
  58.           write(*,'(8(f8.6,'' ''))') (r(k),k=0,7)
  59.         else if(j.eq.no-1) then
  60.           write(*,'(8(f8.6,'' ''))') (r(k),k=0,mod(no-1,8))
  61.         endif
  62.  1000 continue
  63. *
  64.       stop
  65.       end
  66. ************************************************************************
  67.       subroutine sgrnd(seed)
  68. *
  69.       implicit integer(a-z)
  70. *
  71. * Period parameters
  72.       parameter(N     =  624)
  73. *
  74.       dimension mt(0:N-1)
  75. *                     the array for the state vector
  76.       common /block/mti,mt
  77.       save   /block/
  78. *
  79. *      setting initial seeds to mt[N] using
  80. *      the generator Line 25 of Table 1 in
  81. *      [KNUTH 1981, The Art of Computer Programming
  82. *         Vol. 2 (2nd Ed.), pp102]
  83. *
  84.       mt(0)= iand(seed,-1)
  85.       do 1000 mti=1,N-1
  86.         mt(mti) = iand(69069 * mt(mti-1),-1)
  87.  1000 continue
  88. *
  89.       return
  90.       end
  91. ************************************************************************
  92.       double precision function grnd()
  93. *
  94.       implicit integer(a-z)
  95. *
  96. * Period parameters
  97.       parameter(N     =  624)
  98.       parameter(N1    =  N+1)
  99.       parameter(M     =  397)
  100.       parameter(MATA  = -1727483681)
  101. *                                    constant vector a
  102.       parameter(UMASK = -2147483648)
  103. *                                    most significant w-r bits
  104.       parameter(LMASK =  2147483647)
  105. *                                    least significant r bits
  106. * Tempering parameters
  107.       parameter(TMASKB= -1658038656)
  108.       parameter(TMASKC= -272236544)
  109. *
  110.       dimension mt(0:N-1)
  111. *                     the array for the state vector
  112.       common /block/mti,mt
  113.       save   /block/
  114.       data   mti/N1/
  115. *                     mti==N+1 means mt[N] is not initialized
  116. *
  117.       dimension mag01(0:1)
  118.       data mag01/0, MATA/
  119.       save mag01
  120. *                        mag01(x) = x * MATA for x=0,1
  121. *
  122.       TSHFTU(y)=ishft(y,-11)
  123.       TSHFTS(y)=ishft(y,7)
  124.       TSHFTT(y)=ishft(y,15)
  125.       TSHFTL(y)=ishft(y,-18)
  126. *
  127.       if(mti.ge.N) then
  128. *                       generate N words at one time
  129.         if(mti.eq.N+1) then
  130. *                            if sgrnd() has not been called,
  131.           call sgrnd(4357)
  132. *                              a default initial seed is used
  133.         endif
  134. *
  135.         do 1000 kk=0,N-M-1
  136.             y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
  137.             mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1)))
  138.  1000   continue
  139.         do 1100 kk=N-M,N-2
  140.             y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK))
  141.             mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1)))
  142.  1100   continue
  143.         y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK))
  144.         mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1)))
  145.         mti = 0
  146.       endif
  147. *
  148.       y=mt(mti)
  149.       mti=mti+1
  150.       y=ieor(y,TSHFTU(y))
  151.       y=ieor(y,iand(TSHFTS(y),TMASKB))
  152.       y=ieor(y,iand(TSHFTT(y),TMASKC))
  153.       y=ieor(y,TSHFTL(y))
  154. *
  155.       if(y.lt.0) then
  156.         grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0)
  157.       else
  158.         grnd=dble(y)/(2.0d0**32-1.0d0)
  159.       endif
  160. *
  161.       return
  162.       end
  163.