* A C-program for MT19937: Integer version * genrand() generates one pseudorandom unsigned integer (32bit) * which is uniformly distributed among 0 to 2^32-1 for each * call. sgenrand(seed) set initial values to the working area * of 624 words. Before genrand(), sgenrand(seed) must be * called once. (seed is any 32-bit integer except for 0). * Coded by Takuji Nishimura, considering the suggestions by * Topher Cooper and Marc Rieffel in July-Aug. 1997. * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later * version. * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. * See the GNU Library General Public License for more details. * You should have received a copy of the GNU Library General * Public License along with this library; if not, write to the * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA * 02111-1307 USA * * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. * When you use this, send an email to: matumoto@math.keio.ac.jp * with an appropriate reference to your work. * ************************************************************************ * Fortran translation by Hiroshi Takano. Jan. 13, 1999. * * genrand() -> integer function grnd() * sgenrand(seed) -> subroutine sgrnd(seed) * integer seed * * This program uses the following non-standard intrinsics. * ishft(i,n): If n>0, shifts bits in i by n positions to left. * If n<0, shifts bits in i by n positions to right. * iand (i,j): Performs logical AND on corresponding bits of i and j. * ior (i,j): Performs inclusive OR on corresponding bits of i and j. * ieor (i,j): Performs exclusive OR on corresponding bits of i and j. * ************************************************************************ * Reverse Generator by Katsumi Hagita. Jun. 10, 2000. * * revrand() means integer functions of grnd() in reverse * ************************************************************************ * this main() outputs first 1000 generated numbers * and the 1000 numbers in reverse of the order program main * implicit integer(a-z) * * Period parameters parameter(N = 624) parameter(NO = 100000) * dimension mtn(0:NO) dimension mtr(0:NO) dimension mt(0:N-1) * the array for the state vector dimension i(0:7) * common /block/mti,mt save /block/ * call sgrnd(4357) * any nonzero integer can be used as a seed do 1000 j=0,no CC i(mod(j,8))=grnd() mtn(j)=grnd() i(mod(j,8))=mtn(j) if(mod(j,8).eq.7) then CCCC write(*,'(8(i12,'' ''))') (i(k),k=0,7) else if(j.eq.no) then CCCC write(*,'(8(i12,'' ''))') (i(k),k=0,mod(no-1,8)) endif 1000 continue * do 1100 j=0,no CC i(mod(j,8))=revgrnd() mtr(no-j)=revgrnd() i(mod(j,8))=mtr(no-j) if(mod(j,8).eq.7) then CCCC write(*,'(8(i12,'' ''))') (i(k),k=0,7) else if(j.eq.no) then CCCC write(*,'(8(i12,'' ''))') (i(k),k=0,mod(no-1,8)) endif 1100 continue do 1200 j=0,no if( mtn(j)-mtr(j).ne.0) then write(6,*) j,mtn(j),mtr(j), mtn(j)-mtr(j) endif 1200 continue stop end ************************************************************************ subroutine sgrnd(seed) * implicit integer(a-z) * * Period parameters parameter(N = 624) * dimension mt(0:N-1) * the array for the state vector common /block/mti,mt save /block/ * * setting initial seeds to mt[N] using * the generator Line 25 of Table 1 in * [KNUTH 1981, The Art of Computer Programming * Vol. 2 (2nd Ed.), pp102] * mt(0)= iand(seed,-1) do 1000 mti=1,N-1 mt(mti) = iand(69069 * mt(mti-1),-1) 1000 continue * return end ************************************************************************ integer function grnd() * implicit integer(a-z) * * Period parameters parameter(N = 624) parameter(N1 = N+1) parameter(M = 397) parameter(MATA = -1727483681) * constant vector a parameter(UMASK = -2147483648) * most significant w-r bits parameter(LMASK = 2147483647) * least significant r bits * Tempering parameters parameter(TMASKB= -1658038656) parameter(TMASKC= -272236544) * dimension mt(0:N-1) * the array for the state vector common /block/mti,mt save /block/ data mti/N1/ * mti==N+1 means mt[N] is not initialized * dimension mag01(0:1) data mag01/0, MATA/ save mag01 * mag01(x) = x * MATA for x=0,1 * TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) * if(mti.ge.N) then * generate N words at one time if(mti.eq.N+1) then * if sgrnd() has not been called, call sgrnd(4357) * a default initial seed is used endif * do 1000 kk=0,N-M-1 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) 1000 continue do 1100 kk=N-M,N-2 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) 1100 continue y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK)) mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti = 0 endif * y=mt(mti) mti=mti+1 y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) * grnd=y * return end ********************************** ********************************** integer function revgrnd() * implicit integer(a-z) * * Period parameters parameter(N = 624) parameter(N1 = N+1) parameter(M = 397) parameter(MATA = -1727483681) * constant vector a parameter(UMASK = -2147483648) * most significant w-r bits parameter(LMASK = 2147483647) * least significant r bits parameter(YMASK = -2) * Tempering parameters parameter(TMASKB= -1658038656) parameter(TMASKC= -272236544) * dimension mt(0:N-1) * the array for the state vector common /block/mti,mt save /block/ *** data mti/N1/ * mti==N+1 means mt[N] is not initialized * dimension mag01(0:1) data mag01/0, MATA/ save mag01 * mag01(x) = x * MATA for x=0,1 * TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) * mti=mti-1 y=mt(mti) if(mti.eq.0) then z=ieor(mt(N-1),mt(M-1)) x=ishft(z,-31) y=ieor(z,mag01(x)) p=ior(iand(ishft(y,1),YMASK),x) mt0L=iand(LMASK,p) mt(0)=ior(iand(mt(0),UMASK),mt0L) mt0=mt(0) do 2000 kk=N-1,N-M+1,-1 z=ieor(mt(kk-1),mt(kk-1+M-N)) x=ishft(z,-31) y=ieor(z,mag01(x)) q=ior(iand(ishft(y,1),YMASK),x) mt(kk)=ior(iand(UMASK,p),iand(LMASK,q)) p=q 2000 continue do 2010 kk=N-M,1,-1 z=ieor(mt(kk-1),mt(kk-1+M)) x=ishft(z,-31) y=ieor(z,mag01(x)) q=ior(iand(ishft(y,1),YMASK),x) mt(kk)=ior(iand(UMASK,p),iand(LMASK,q)) p=q 2010 continue mt(0)=ior(iand(UMASK,p),mt0L) y=mt0 mti = N endif * y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) * revgrnd=y return end