fortranでDLAを書きました。
implicit none
integer,parameter::r=1000
integer,parameter::n=100000
real*8::c
character*4::num
integer::seedsize,k,x,y,a
integer,allocatable::seed(:)
integer::p(r,r),mark(r,r),t,i,j,ii,jj
integer::kmax,check,count,anum
call random_seed(size=seedsize)
allocate(seed(seedsize))
call random_seed(get=seed)
p=0
p(r/2,r/2)=2
t=0
timeloop:do
t=t+1
write(0,*)t
do x=1,r
do y=1,r
if(x==1 .or. x==r .or. y==1 .or. y==r)then
if(p(x,y)==2) then
write(0,*)"ok"
exit timeloop
endif
endif
enddo
enddo
mark=0
check=0
do x=1,r
do y=1,r
if(p(x,y)==1)check=check+1
enddo
enddo
kmax=n-check
if(kmax==0)kmax=n+1
check=0
aloop:do k=1,kmax
if(kmax==n+1) exit aloop
do
do
call random_number(c)
if(c==0d0)cycle
a=nint(c*real(r))
if(a==0) cycle
exit
enddo
call random_number(c)
if(0d0<c .and. c<=0.25d0) then
x=1
y=a
endif
if(0.25d0<c .and. c<=0.50d0) then
x=r
y=a
endif
if(0.50d0<c .and. c<=0.75d0) then
x=a
y=1
endif
if(0.75d0<c .and. c<=1d0)then
x=a
y=r
endif
if(p(x,y)/=0)then
if(check==r*r) exit aloop
check=check+1
cycle
endif
p(x,y)=1
mark(x,y)=1
loops1:do ii=-1,1
do jj=-1,1
if(x+ii==0 .or. x+ii==r+1 .or. y+jj==0 .or. y+jj==r+1)cycle
if(p(x+ii,y+jj)==2) p(x,y)=2
exit loops1
enddo
enddo loops1
exit
enddo
enddo aloop
check=0
do y=1,r
do x=1,r
if(p(x,y)==1 .and. mark(x,y)==0)check=check+1
enddo
enddo
anum=check
count=0
do
do
call random_number(c)
x=nint(c*real(r))
if(x==0) cycle
exit
enddo
do
call random_number(c)
y=nint(c*real(r))
if(y==0) cycle
exit
enddo
if(p(x,y)==1 .and. mark(x,y)==0) then
count=count+1
check=0
do ii=-1,1
loopj:do jj=-1,1
if(x+ii==0 .or. x+ii==r+1 .or. y+jj==0 .or. y+jj==r+1)then
check=check+1
cycle
endif
if(p(x+ii,y+jj)/=0)check=check+1
enddo loopj
enddo
if(check==9)then
! write(0,*)"full"
i=0
j=0
else
do
call random_number(c)
if(c==0d0)cycle
if(0d0<c .and. c<=0.125d0 )then
i=-1
j=-1
endif
if(0.125d0<c .and. c<=0.250d0)then
i=0
j=-1
endif
if(0.250d0<c .and. c<=0.375d0)then
i=1
j=-1
endif
if(0.375d0<c .and. c<=0.500d0)then
i=-1
j=0
endif
if(0.500d0<c .and. c<=0.625d0)then
i=1
j=0
endif
if(0.625d0<c .and. c<=0.750d0)then
i=-1
j=1
endif
if(0.750d0<c .and. c<=0.875d0)then
i=0
j=1
endif
if(0.875d0<c .and. c<=1d0 )then
i=1
j=1
endif
if(x+i==0 .or. x+i==r+1 .or. y+j==0 .or. y+j==r+1)cycle
if(p(x+i,y+j)/=0) cycle
mark(x+i,y+j)=1
p(x,y)=0
p(x+i,y+j)=1
exit
enddo
endif
loopm:do ii=-1,1
do jj=-1,1
if(x+i+ii==0 .or. x+i+ii==r+1 .or. y+j+jj==0 .or. y+j+jj==r+1)cycle
if(p(x+i+ii,y+j+jj)==2) then
p(x+i,y+j)=2
exit loopm
endif
enddo
enddo loopm
end if
if(count==anum)exit
enddo
if(t/100==real(t)/real(100))then
! write(num,'(i4.4)') t
open(11,file="hoge.dat",status='replace')
do x=1,r
do y=1,r
write(11,*)x,y,p(x,y)
enddo
enddo
close(11)
endif
end do timeloop
do x=1,r
do y=1,r
print*,x,y,p(x,y)
enddo
enddo
end program main