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