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