c
      PROGRAM ran
      IMPLICIT REAL*8(A-H,O-Z)
      dimension num(32)

      write(*,*)'Enter Initial idum and N'
      read(*,*)idum,N

      sum=0.d0
      do i=1,N
        x=gasdev(idum)
c        write(*,*)i,idum,x
        write(1,*)i,idum,x
        x1=-4.d0
        do j=1,32
	  x2=x1+0.25d0
          if(x.gt.x1.and.x.le.x2)num(j)=num(j)+1
          x1=x2
        end do
        sum=sum+x
      end do
      write(*,*)'average=',sum/N

      x1=-4.d0
      do j=1,32
        x2=x1+0.25d0
        write(*,1000)x1,x2,(x1+x2)/2.d0,num(j)
        write(2,*)x1,x2,(x1+x2)/2.d0,num(j)
        x1=x2
      end do
 1000 format('x1,x2,x12,num=',3f8.3,i5)

      stop
      end
c
c *** Randam Number with a normally distributed deviate with zero mean
c     and unit variance, using ran0 as the source of uniform deviates
c
      function gasdev(idum)
      IMPLICIT REAL*8(A-H,O-Z)
      save iset,gset
      data iset/0/

      if(iset.eq.0)then
   10   v1=2.d0*ran0(idum)-1.d0
        v2=2.d0*ran0(idum)-1.d0
        rsq=v1*v1+v2*v2
        if(rsq.ge.1.d0.or.rsq.eq.0.d0)goto 10
        fac=dsqrt(-2.d0*dlog(rsq)/rsq)
        gset=v1*fac
        gasdev=v2*fac
        iset=1
      else
        gasdev=gset
        iset=0
      end if

      return
      end

c *** Randam Function **
c
      function ran0(idum)
      IMPLICIT REAL*8(A-H,O-Z)
      parameter(IA=16807,IM=2147483647,AM=1.d0/IM,IQ=127773
     *         ,IR=2836,MASK=123459876)
c
      idum=ieor(idum,MASK)
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if(idum.lt.0)idum=idum+IM
      ran0=AM*idum
      idum=ieor(idum,MASK)
c
      return
      end
