      program bboxm
c
c Box-Mueller algorithm for generating gaussian random numbers
c
      parameter(msamp=10000,nnum=50000)
c
      ibm=1955
      rinv=0.5/(2.**31-1.)
      do i=1,1000
      ibm=ibm*16807
      enddo
c
      pi2=8*atan(1.)
c
      avemax=0.
c
      do isamp=1,msamp
c
      xmax=0.
c
      do i=1,nnum/2
      ibm=ibm*16807
      y1=rinv*float(ibm)+0.5
      ibm=ibm*16807
      y2=rinv*float(ibm)+0.5
      x1=sqrt(-2.*log(y1))*cos(pi2*y2)
      x2=sqrt(-2.*log(y1))*sin(pi2*y2)
      xmax=max(xmax,x1)
      xmax=max(xmax,x2)
      enddo
c
      avemax=avemax+xmax
c
      enddo
c
      avemax=avemax/msamp
      print *,nnum,avemax
      end

