      program ex093
      implicit real*8(a-h,o-z)
      integer*4 iseed
      real*8 hist(-50:50)  ! hist(-50), hist(-49), ....., hist(50)
c
      iseed = 123456789
      n = 100000
      dx = 3.0d0/50
      pi = 4.0d0*atan(1.0d0)
c
      do j = -50,50
         hist(j) = 0.0d0
      end do
c
      do i = 1,n
         x1 = rnd(iseed)
         x2 = rnd(iseed)
           if(x2.gt.0) then
              y1 = cos(2*pi*x1)*sqrt(-2*log(x2))
              y2 = sin(2*pi*x1)*sqrt(-2*log(x2))
              j = nint(y1/dx)
                if(j.ge.-50.and.j.le.50) then
                   hist(j) = hist(j) + 1.0d0
                endif
           endif
      end do
c
      do j = -50,50
         x = dx*j
         write(*,100) x, hist(j)/n/dx, exp(-x**2/2)/sqrt(2*pi)
      end do
c
 100  format(3(1x,f10.4))
      end
c **********************************************************************
      function rnd(iseed)
c **********************************************************************
      implicit real*8(a-h,o-z)
c...32 bit (=4byte) machine
      integer*4 iseed,il,ic,ih
      il = 843314861
      ic = 453816693
      ih = 2**30
      xmax = 2.0d0**31
c
      iseed = iseed*il + ic
      if(iseed.lt.0) iseed = (iseed + ih) + ih
      rnd = iseed/xmax
c
      end
