Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (breakdown) Browse

fennecs/src/random.f
3eaf507ffec2breakdown

Last Change
Raw File

random.f
View Options

* Version 1.0 of random number routines
* Author: Charles Karney <karney@princeton.edu>
* Date: 1999-08-05 10:44:51 -0400
*
function random(p,r)
implicit none
DOUBLE PRECISION ulp2
parameter(ulp2=2.0d0**(-47-1))
REAL ulps
DOUBLE PRECISION mult
parameter(ulps=2.0**(-23),mult=2.0d0**23)
DOUBLE PRECISION random
REAL srandom
integer p
DOUBLE PRECISION r(0:100-1)
external rand_batch
if(p.GE.100)then
call rand_batch(p,r(0))
end if
random=r(p)+ulp2
p=p+1
return
entry srandom(p,r)
if(p.GE.100)then
call rand_batch(p,r(0))
end if
srandom=(int(mult*r(p))+0.5)*ulps
p=p+1
return
end
*
subroutine random_array(y,n,p,r)
implicit none
DOUBLE PRECISION ulp2
parameter(ulp2=2.0d0**(-47-1))
REAL ulps
DOUBLE PRECISION mult
parameter(ulps=2.0**(-23),mult=2.0d0**23)
integer n
DOUBLE PRECISION y(0:n-1)
REAL ys(0:n-1)
integer p
DOUBLE PRECISION r(0:100-1)
integer i,k,j
external rand_batch
if(n.LE.0)return
k=min(n,100-p)
do i=0,k-1
y(i)=r(i+p)+ulp2
end do
p=p+(k)
do j=k,n-1,100
call rand_batch(p,r(0))
do i=j,min(j+100,n)-1
y(i)=r(i-j+p)+ulp2
end do
p=p+(min(100,n-j))
end do
return
entry srandom_array(ys,n,p,r)
if(n.LE.0)return
k=min(n,100-p)
do i=0,k-1
ys(i)=(int(mult*r(i+p))+0.5)*ulps
end do
p=p+(k)
do j=k,n-1,100
call rand_batch(p,r(0))
do i=j,min(j+100,n)-1
ys(i)=(int(mult*r(i-j+p))+0.5)*ulps
end do
p=p+(min(100,n-j))
end do
return
end
*
subroutine rand_batch(p,r)
implicit none
integer p
DOUBLE PRECISION r(0:100-1)
integer i
DOUBLE PRECISION w(0:1009-100-1)
DOUBLE PRECISION tmp
do i=0,63-1
tmp=r(i)+r(i+100-63)
w(i)=tmp-int(tmp)
end do
do i=63,100-1
tmp=r(i)+w(i-63)
w(i)=tmp-int(tmp)
end do
do i=100,1009-100-1
tmp=w(i-100)+w(i-63)
w(i)=tmp-int(tmp)
end do
do i=1009-100,1009-100+63-1
tmp=w(i-100)+w(i-63)
r(i-1009+100)=tmp-int(tmp)
end do
do i=1009-100+63,1009-1
tmp=w(i-100)+r(i-1009+100-63)
r(i-1009+100)=tmp-int(tmp)
end do
p=0
return
end
*
subroutine random_init(seed,p,r)
implicit none
integer b
DOUBLE PRECISION del,ulp
parameter(del=2.0d0**(-14),ulp=2.0d0**(-47),b=2**14)
integer seed(0:8-1)
integer p
DOUBLE PRECISION r(0:100-1)
integer i,j,t,s(0:8-1)
logical even
integer z(0:8-1)
integer a0,a1,a2,a3,a4,a5,a6,c0
data a0,a1,a2,a3,a4,a5,a6,c0/15661,678,724,5245,13656,11852,29,1/
do i=0,8-1
s(i)=seed(i)
end do
even=.TRUE.
even=even.AND.(mod(s(7),2).EQ.0)
r(0)=(((s(7)*del+s(6))*del+s(5))*del+int(s(4)/512))*512*del
do j=1,100-1
t=c0+a0*s(0)
z(0)=mod(t,b)
t=int(t/b)+a0*s(1)+a1*s(0)
z(1)=mod(t,b)
t=int(t/b)+a0*s(2)+a1*s(1)+a2*s(0)
z(2)=mod(t,b)
t=int(t/b)+a0*s(3)+a1*s(2)+a2*s(1)+a3*s(0)
z(3)=mod(t,b)
t=int(t/b)+a0*s(4)+a1*s(3)+a2*s(2)+a3*s(1)+a4*s(0)
z(4)=mod(t,b)
t=int(t/b)+a0*s(5)+a1*s(4)+a2*s(3)+a3*s(2)+a4*s(1)+a5*s(0)
z(5)=mod(t,b)
t=int(t/b)+a0*s(6)+a1*s(5)+a2*s(4)+a3*s(3)+a4*s(2)+a5*s(1)+a6*s(0)
z(6)=mod(t,b)
t=int(t/b)+a0*s(7)+a1*s(6)+a2*s(5)+a3*s(4)+a4*s(3)+a5*s(2)+a6*s(1)
z(7)=mod(t,b)
do i=0,8-1
s(i)=z(i)
end do
even=even.AND.(mod(s(7),2).EQ.0)
r(j)=(((s(7)*del+s(6))*del+s(5))*del+int(s(4)/512))*512*del
end do
if(even)then
t=c0+a0*s(0)
z(0)=mod(t,b)
t=int(t/b)+a0*s(1)+a1*s(0)
z(1)=mod(t,b)
t=int(t/b)+a0*s(2)+a1*s(1)+a2*s(0)
z(2)=mod(t,b)
t=int(t/b)+a0*s(3)+a1*s(2)+a2*s(1)+a3*s(0)
z(3)=mod(t,b)
t=int(t/b)+a0*s(4)+a1*s(3)+a2*s(2)+a3*s(1)+a4*s(0)
z(4)=mod(t,b)
t=int(t/b)+a0*s(5)+a1*s(4)+a2*s(3)+a3*s(2)+a4*s(1)+a5*s(0)
z(5)=mod(t,b)
t=int(t/b)+a0*s(6)+a1*s(5)+a2*s(4)+a3*s(3)+a4*s(2)+a5*s(1)+a6*s(0)
z(6)=mod(t,b)
t=int(t/b)+a0*s(7)+a1*s(6)+a2*s(5)+a3*s(4)+a4*s(3)+a5*s(2)+a6*s(1)
z(7)=mod(t,b)
do i=0,8-1
s(i)=z(i)
end do
j=int((s(8-1)*100)/b)
r(j)=r(j)+(ulp)
end if
p=100
return
end
*
subroutine decimal_to_seed(decimal,seed)
implicit none
integer b
parameter(b=2**14)
character*(*)decimal
integer seed(0:8-1)
integer i,ten(0:8-1),c(0:8-1),ch
data ten/10,7*0/
external rand_axc
do i=0,8-1
seed(i)=0
c(i)=0
end do
do i=1,len(decimal)
ch=ichar(decimal(i:i))
if(ch.GE.ichar('0').AND.ch.LE.ichar('9'))then
c(0)=ch-ichar('0')
call rand_axc(ten,seed,c)
end if
end do
return
end
*
subroutine string_to_seed(string,seed)
implicit none
integer b
parameter(b=2**14)
character*(*)string
integer seed(0:8-1)
integer t,i,k,unity(0:8-1),c(0:8-1),ch
data unity/1,7*0/
external rand_axc
do i=0,8-1
seed(i)=0
c(i)=0
end do
do i=1,len(string)
ch=ichar(string(i:i))
if(ch.GT.ichar(' ').AND.ch.LT.127)then
t=mod(seed(0),2)*(b/2)
do k=0,8-1
seed(k)=int(seed(k)/2)
if(k.LT.8-1)then
seed(k)=seed(k)+(mod(seed(k+1),2)*(b/2))
else
seed(k)=seed(k)+(t)
end if
end do
c(0)=ch
call rand_axc(unity,seed,c)
end if
end do
return
end
*
subroutine set_random_seed(time,seed)
implicit none
integer time(8)
integer seed(0:8-1)
character*21 c
external decimal_to_seed
write(c(1:8),'(i4.4,2i2.2)')time(1),time(2),time(3)
write(c(9:12),'(i1.1,i3.3)') (1-sign(1,time(4)))/2,abs(time(4))
write(c(13:21),'(3i2.2,i3.3)')time(5),time(6),time(7),time(8)
call decimal_to_seed(c,seed)
return
end
*
subroutine seed_to_decimal(seed,decimal)
implicit none
integer pow,decbase,b
parameter(pow=4,decbase=10**pow,b=2**14)
character*(*)decimal
integer seed(0:8-1)
integer z(0:8-1),i,t,j,k
character*36 str
k=-1
do i=0,8-1
z(i)=seed(i)
if(z(i).GT.0)k=i
end do
str=' '
i=9
90000 continue
i=i-1
t=0
do j=k,0,-1
z(j)=z(j)+t*b
t=mod(z(j),decbase)
z(j)=int(z(j)/decbase)
end do
if(z(max(0,k)).EQ.0)k=k-1
if(k.GE.0)then
write(str(pow*i+1:pow*(i+1)),'(i4.4)')t
else
write(str(pow*i+1:pow*(i+1)),'(i4)')t
end if
if(k.GE.0)goto 90000
if(len(decimal).GT.len(str))then
decimal(:len(decimal)-len(str))=' '
decimal(len(decimal)-len(str)+1:)=str
else
decimal=str(len(str)-len(decimal)+1:)
end if
return
end
*
subroutine rand_next_seed(n,ax,cx,y)
implicit none
integer n,ax(0:8-1),cx(0:8-1)
integer y(0:8-1)
integer a(0:8-1),c(0:8-1),z(0:8-1),t(0:8-1),m,i
data z/8*0/
external rand_axc
if(n.EQ.0)return
m=n
do i=0,8-1
a(i)=ax(i)
c(i)=cx(i)
end do
90000 continue
if(mod(m,2).GT.0)then
call rand_axc(a,y,c)
end if
m=int(m/2)
if(m.EQ.0)return
do i=0,8-1
t(i)=c(i)
end do
call rand_axc(a,c,t)
do i=0,8-1
t(i)=a(i)
end do
call rand_axc(t,a,z)
goto 90000
end
*
subroutine next_seed3(n0,n1,n2,seed)
implicit none
integer n0,n1,n2
integer seed(0:8-1)
integer af0(0:8-1),cf0(0:8-1)
integer ab0(0:8-1),cb0(0:8-1)
integer af1(0:8-1),cf1(0:8-1)
integer ab1(0:8-1),cb1(0:8-1)
integer af2(0:8-1),cf2(0:8-1)
integer ab2(0:8-1),cb2(0:8-1)
data af0/15741,8689,9280,4732,12011,7130,6824,12302/
data cf0/16317,10266,1198,331,10769,8310,2779,13880/
data ab0/9173,9894,15203,15379,7981,2280,8071,429/
data cb0/8383,3616,597,12724,15663,9639,187,4866/
data af1/8405,4808,3603,6718,13766,9243,10375,12108/
data cf1/13951,7170,9039,11206,8706,14101,1864,15191/
data ab1/6269,3240,9759,7130,15320,14399,3675,1380/
data cb1/15357,5843,6205,16275,8838,12132,2198,10330/
data af2/445,10754,1869,6593,385,12498,14501,7383/
data cf2/2285,8057,3864,10235,1805,10614,9615,15522/
data ab2/405,4903,2746,1477,3263,13564,8139,2362/
data cb2/8463,575,5876,2220,4924,1701,9060,5639/
external rand_next_seed
if(n2.GT.0)then
call rand_next_seed(n2,af2,cf2,seed)
else if(n2.LT.0)then
call rand_next_seed(-n2,ab2,cb2,seed)
end if
if(n1.GT.0)then
call rand_next_seed(n1,af1,cf1,seed)
else if(n1.LT.0)then
call rand_next_seed(-n1,ab1,cb1,seed)
end if
entry next_seed(n0,seed)
if(n0.GT.0)then
call rand_next_seed(n0,af0,cf0,seed)
else if(n0.LT.0)then
call rand_next_seed(-n0,ab0,cb0,seed)
end if
return
end
*
subroutine rand_axc(a,x,c)
implicit none
integer b
parameter(b=2**14)
integer a(0:8-1),c(0:8-1)
integer x(0:8-1)
integer z(0:8-1),i,j,t
t=0
do i=0,8-1
t=c(i)+int(t/b)
do j=0,i
t=t+(a(j)*x(i-j))
end do
z(i)=mod(t,b)
end do
do i=0,8-1
x(i)=z(i)
end do
return
end
*
c4science · Help