subroutine dranini(iseed0) implicit double precision(a-h,o-z) 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(randxx(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) do 856 i=0,m 856 g(i)=g(i)/u2th return end ! subroutine dranread(iunit) ! implicit double precision(a-h,o-z) ! 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 ! open(UNIT=iunit) ! read (iunit,*) ic ! read (iunit,*) (ix(i),i=1,ip) ! read (iunit,*) (g(i),i=0,m) ! close(iunit) ! return ! end ! ! subroutine dranwrite(iunit) ! implicit double precision(a-h,o-z) ! 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 ! open(UNIT=iunit) ! write (iunit,*) ic ! write (iunit,*) (ix(i),i=1,ip) ! write (iunit,*) (g(i),i=0,m) ! close(iunit) ! return ! end ! ! subroutine dranread2(iunit,cadena) ! implicit double precision(a-h,o-z) ! character*100 cadena ! 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 ! open(UNIT=iunit,FILE=cadena) ! read (iunit,*) ic ! read (iunit,*) (ix(i),i=1,ip) ! read (iunit,*) (g(i),i=0,m) ! close(iunit) ! return ! end ! ! subroutine dranwrite2(iunit,cadena) ! implicit double precision(a-h,o-z) ! character*100 cadena ! 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 ! open(UNIT=iunit,FILE=cadena,STATUS="UNKNOWN") ! write (iunit,*) ic ! write (iunit,*) (ix(i),i=1,ip) ! write (iunit,*) (g(i),i=0,m) ! close(iunit) ! return ! end function idran(n) implicit double precision(a-h,o-z) 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 idran=ix(ic) if (n.gt.0) idran=mod(idran,n)+1 return end function dranu() implicit double precision(a-h,o-z) 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 dranu=real(ix(ic))/rmax return end function drang() implicit double precision(a-h,o-z) 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) drang=i2*g(i+1)+(nn1-i2)*g(i) return end subroutine dranbm(x1,x2) implicit double precision(a-h,o-z) 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 drangv(u,n) implicit double precision(a-h,o-z) 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 drangvv(iv,u,n) implicit double precision(a-h,o-z) 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 randxx(dseed) double precision a,c,xm,rm,dseed,randxx parameter (xm=2.d0**32,rm=1.d0/xm,a=69069.d0,c=1.d0) dseed=mod(dseed*a+c,xm) randxx=dseed*rm return end