PROGRAM gauss_rand
  IMPLICIT NONE

  INTEGER :: i,seed,index,min,max
  DOUBLE PRECISION :: gaussianrand,area
  INTEGER :: stat(-10000:10000)

  seed=11467342

  DO i=-10000,10000
     stat(i)=0
  ENDDO

  DO i=1,1000000
     index=FLOOR(gaussianrand(seed)*100+0.5)
     IF (index < min) min=index
     IF (index > max) max=index
     if (index < -10000) index=-10000
     if (index >  10000) index=10000
     stat(index)=stat(index)+1
  ENDDO

  DO i=min-1,max+1
     PRINT *,i/100.0,stat(i)/10000.0
  ENDDO

END PROGRAM gauss_rand

! ---------------------------------------------------------------

DOUBLE PRECISION FUNCTION uniformrand(seed)
    IMPLICIT NONE

    ! -----------------------------------------------------------
    ! Park-Miller "minimal" Random number generator
    ! uniform distribution ]0,1[ . See Numerical Recipes ch. 7.0
    ! -----------------------------------------------------------

    INTEGER :: seed,ia,im,iq,ir,mask
    DOUBLE PRECISION :: am
    INTEGER :: k

   PARAMETER (ia=16807,im=2147483647,am=1.0/im,iq=127773,ir=2836,mask=123459876)

    seed=IEOR(seed,mask)
    k=seed/iq
    seed=IA*(seed-k*iq)-ir*k
    IF (seed < 0) seed=seed+im
    uniformrand=am*seed
    seed=IEOR(seed,mask)

    RETURN

END FUNCTION uniformrand

! ---------------------------------------------------------------

DOUBLE PRECISION FUNCTION gaussianrand(seed)
  IMPLICIT NONE

  ! ---------------------------------------------------------
  ! Random numbers with normal (Gaussian) distribution.
  ! Mean is 0 and standard deviation is 1
  ! See W.H.Press et al., Numerical Recipes 1st ed., page 203
  ! ---------------------------------------------------------

  INTEGER :: seed

  DOUBLE PRECISION :: fac,v1,v2,r
  DOUBLE PRECISION, SAVE :: gset
  INTEGER, SAVE :: iset=0
  DOUBLE PRECISION :: uniformrand



  IF (iset.EQ.0) THEN
     DO 
	v1 = 2.*uniformrand(seed)-1.
	v2 = 2.*uniformrand(seed)-1.
	r = v1*v1+v2*v2
	IF (r.LE.1.) EXIT
     END DO
     
     if (r /= 0.0d0) then
	fac = SQRT(-2.0*LOG(r)/r)
     else
	fac=0.0d0
     endif
     gset = v1*fac
     gaussianrand = v2*fac
     iset = 1
  ELSE
     gaussianrand = gset
     iset = 0
  ENDIF

  RETURN
END FUNCTION gaussianrand