subroutine ran_ini(iseed0) double precision dseed parameter(ip=1279) parameter(np=14) parameter(nbit=31) parameter(m=2**np,np1=nbit-np,nn=2**np1-1,nn1=nn+1) integer ix(ip) dimension g(0:m) data c0,c1,c2/2.515517,0.802853,0.010328/ data d1,d2,d3/1.432788,0.189269,0.001308/ c common /ixx/ ix common /icc/ ic common /gg/ g c dseed=iseed0 do 200 i=1,ip ix(i)=0 do 200 j=0,nbit-1 if(rand_xx(dseed).lt.0.5) ix(i)=ibset(ix(i),j) 200 continue ic=0 c pi=4.0d0*datan(1.0d0) do 1 i=m/2,m p=1.0-real(i+1)/(m+2) t=sqrt(-2.0*log(p)) x=t-(c0+t*(c1+c2*t))/(1.0+t*(d1+t*(d2+t*d3))) g(i)=x g(m-i)=-x 1 continue u2th=1.0-real(m+2)/m*sqrt(2.0/pi)*g(m)*exp(-g(m)*g(m)/2) u2th=nn1*sqrt(u2th) write(6,*) u2th do 856 i=0,m 856 g(i)=g(i)/u2th return end subroutine ran_read(iunit) parameter(ip=1279) parameter(np=14) parameter(m=2**np) integer ix(ip) dimension g(0:m) common /ixx/ ix common /icc/ ic common /gg/ g read (iunit,*) ic read (iunit,*) (ix(i),i=1,ip) read (iunit,*) (g(i),i=0,m) return end subroutine ran_write(iunit) parameter(ip=1279) parameter(np=14) parameter(m=2**np) integer ix(ip) dimension g(0:m) common /ixx/ ix common /icc/ ic common /gg/ g write (iunit,*) ic write (iunit,*) (ix(i),i=1,ip) write (iunit,*) (g(i),i=0,m) return end function i_ran(n) parameter(ip=1279) parameter(iq=418) parameter(is=ip-iq) integer ix(ip) common /ixx/ ix common /icc/ ic ic=ic+1 if(ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif i_ran=ix(ic) if (n.gt.0) i_ran=mod(i_ran,n)+1 return end function ran_u() parameter(ip=1279) parameter(iq=418) parameter(is=ip-iq) parameter (rmax=2147483647.0) integer ix(ip) common /ixx/ ix common /icc/ ic ic=ic+1 if(ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif ran_u=real(ix(ic))/rmax return end function ran_g() parameter(ip=1279) parameter(iq=418) parameter(np=14) parameter(nbit=31) parameter(m=2**np,np1=nbit-np,nn=2**np1-1,nn1=nn+1) parameter(is=ip-iq) integer ix(ip) dimension g(0:m) common /ixx/ ix common /icc/ ic common /gg/ g ic=ic+1 if(ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif i=ishft(ix(ic),-np1) i2=iand(ix(ic),nn) ran_g=i2*g(i+1)+(nn1-i2)*g(i) return end subroutine ran_bm(x1,x2) parameter(ip=1279) parameter(iq=418) parameter(is=ip-iq) parameter (rmax=2147483647.0) integer ix(ip) common /ixx/ ix common /icc/ ic data pi2 /6.2831853072/ ic=ic+1 if (ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif u=pi2*real(ix(ic))/rmax ic=ic+1 if(ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif v=(real(ix(ic))+0.5)/rmax v=sqrt(-2.0*log(v)) x1=cos(u)*v x2=sin(u)*v return end subroutine ran_gv(u,n) parameter(ip=1279) parameter(iq=418) parameter(np=14) parameter(nbit=31) parameter(m=2**np,np1=nbit-np,nn=2**np1-1,nn1=nn+1) parameter(is=ip-iq) dimension g(0:m) dimension u(n) dimension ix(ip) common /gg/ g common /ixx/ ix common /icc/ic n1=0 10 if (ic.lt.iq) then kmax=min(n-n1,iq-ic) do 99 k=1,kmax ic=ic+1 ix(ic)=ieor(ix(ic),ix(ic+is)) i=ishft(ix(ic),-np1) i2=iand(ix(ic),nn) u(n1+k)=i2*g(i+1)+(nn1-i2)*g(i) 99 continue else kmax=min(n-n1,ip-ic) do 98 k=1,kmax ic=ic+1 ix(ic)=ieor(ix(ic),ix(ic-iq)) i=ishft(ix(ic),-np1) i2=iand(ix(ic),nn) u(n1+k)=i2*g(i+1)+(nn1-i2)*g(i) 98 enddo endif if(ic.ge.ip) ic=0 n1=n1+kmax if (n1.lt.n) goto 10 return end subroutine ran_gvv(iv,u,n) parameter(ip=1279) parameter(iq=418) parameter(np=14) parameter(nbit=31) parameter(m=2**np,np1=nbit-np,nn=2**np1-1,nn1=nn+1) parameter(is=ip-iq) dimension iv(n+ip) dimension g(0:m) dimension u(n) dimension ix(ip) common /gg/ g common /ixx/ ix common /icc/ic do 10 i=ic+1,ip 10 iv(i-ic)=ix(i) do 20 i=1,ic 20 iv(ip-ic+i)=ix(i) cCDEC$ INIT_DEP_FWD do 99 k=1,n ir=ieor(iv(k+is),iv(k)) i=ishft(ir,-np1) i2=iand(ir,nn) u(k)=i2*g(i+1)+(nn1-i2)*g(i) iv(k+ip)=ir 99 continue do 11 i=ic+1,ip 11 ix(i)=iv(n+i-ic) do 21 i=1,ic 21 ix(i)=iv(n+ip-ic+i) return end function rand_xx(dseed) double precision a,c,xm,rm,dseed parameter (xm=2.d0**32,rm=1.d0/xm,a=69069.d0,c=1.d0) dseed=mod(dseed*a+c,xm) rand_xx=dseed*rm return end