!------------------!
 program randomwalk
!------------------!
 implicit none

 real(8), allocatable :: p(:)
 ! p will contain the exact probability distribution
 real(8), allocatable :: h(:,:)
 ! h will contain the generated distributions for each bit

 integer :: i,j,b,b1,b2,nwalk,nstep,x(0:63)
 real(8) :: d,pdev(0:63)
 character(8) :: fname

 call initran(1)
 
 open(10,file='read.in',status='old')
 read(10,*)nwalk,nstep
 close(10)

 allocate(p(-nstep:nstep)); p=0.d0
 allocate(h(-nstep:nstep,0:63)); h=0.d0

 call exactp(p,nstep)

 ! do nwalk random walks and collect in h the end points for each bit
 do i=1,nwalk
    call dowalk(x,nstep)
    do b=0,63
       h(x(b),b)=h(x(b),b)+1.d0
    enddo
 enddo
 h=h/dble(nwalk)

! compute the deviations from the exact random walk
 pdev=0.d0    
 do b=0,63
    d=0.d0
    do i=-nstep,nstep
       d=d+(p(i)-h(i,b))**2
    enddo
    d=sqrt(d)
    pdev(b)=pdev(b)+d
 enddo

 ! write data to files
 ! d.dat will contain the measure of deviations for each bit
 ! pxy.dat will contain the distribution for bit xy
 open(10,file='d.dat')
 fname='p00.dat'
 do b=0,63
    write(10,'(i4,2f20.8)')b,pdev(b)*sqrt(dble(nwalk))
    b1=mod(b,10)
    b2=b/10
    fname(2:2)=achar(48+b2)
    fname(3:3)=achar(48+b1)
    open(20,file=fname)
    do i=-nstep,nstep
       if (h(i,b)>1.d-12) write(20,'(i5,2f16.10)')i,h(i,b),p(i)
    enddo
    close(20)
 enddo
 close(10)

 deallocate(p)
 deallocate(h)

 end program randomwalk
!----------------------!

!-----------------------!
 subroutine dowalk(x,nt)
!-------------------------------------------------------------!
! does 64 different random walks based on the bits of integers
! generated by the random-number generator iran()
!-------------------------------------------------------------!

 implicit none

 integer :: i,b,nt,x(0:63)
 integer(8) :: r,iran

 external :: iran

 x(:)=0
 do i=1,nt
    r=iran()
    do b=0,63
       if (btest(r,b)) then
          x(b)=x(b)+1
       else
          x(b)=x(b)-1
       endif
    enddo
 enddo

 end subroutine dowalk
!---------------------!

!----------------------!
 subroutine exactp(p,n)
!------------------------------------------------!
! Calculates the exact probability distributions.
!------------------------------------------------!
 implicit none

 integer :: i,j,d,n
 real(8) :: p(-n:n)

 p=0.d0
 do j=0,n
    d=2*j-n
    do i=j+1,n
       p(d)=p(d)+log(dble(i))
    enddo
    do i=1,n-j
       p(d)=p(d)-log(dble(i))
    enddo
    p(d)=p(d)-dble(n)*log(2.d0)
    p(d)=exp(p(d))
 enddo

 end subroutine exactp
!---------------------!
 
!--------------------------!
 integer(8) function iran()
!---------------------------------------------------!
! 64-bit congruental generator (generating integers)!                !
! iran64=oran64*2862933555777941757+1013904243      !
!---------------------------------------------------!
 implicit none

 real(8)    :: dmu64
 integer(8) :: ran64,mul64,add64
 common/bran64/dmu64,ran64,mul64,add64

 ran64=ran64*mul64+add64
 iran=ran64

 end function iran
!-----------------!

!---------------------!
 subroutine initran(w)
!------------------------------------------------!
! Initializes the random-number generator        !
! - writes a new seed-integer to seed.in if w/=0 !
!------------------------------------------------!
 implicit none

 integer(8) :: irmax
 integer(4) :: w,nb,b

 real(8)    :: dmu64
 integer(8) :: ran64,mul64,add64
 common/bran64/dmu64,ran64,mul64,add64
      
 irmax=2_8**31
 irmax=2*(irmax**2-1)+1
 mul64=2862933555777941757_8
 add64=1013904243
 dmu64=0.5d0/dble(irmax)

 open(10,file='seed.in',status='old')
 read(10,*)ran64
 close(10)
 if (w.ne.0) then
    open(10,file='seed.in',status='unknown')
    write(10,*)(ran64*mul64)/5+5265361
    close(10)
 endif

 end subroutine initran
!----------------------!