implicit none
      REAL ran1,x,d,V,N
      double precision xNs,xN2s,srxNs,srxN2s,Nx,srxNs2,sigma
      integer K,i,s
      parameter (N=100000,K=1000)
c
      do s=1,K
      x=0
      do i=1,N
      if(ran1(d)<0.5) then
      x=x+1
      else
      x=x-1
      endif
      enddo
c
      xNs=xNs+x
      xN2s=xN2s+x*x
c
      enddo
c
      srxNs= xNs/K
      srxN2s= xN2s/K
      srxNs2=srxNs**2
c
      V=srxN2s-srxNs2
      sigma=sqrt(V)
      write(11,'(10F12.2)')
     & log10(N),log10(sigma),xN2s,xNs,srxNs2,srxN2s,V
c
      end
c
c
c
c
      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      REAL ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     * NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
      idum=max(-idum,1)
      do 11 j=NTAB+8,1,-1
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      if (j.le.NTAB) iv(j)=idum
 11     continue
      iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
c
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END