3 Random number generators

3.2 A Normal Random Number Generator in Fortran

c**   Subroutine to return normally distributed pseudorandom numbers
c**   with zero mean and unit standard deviation, with the permission
c**   of the association for computing machines (acm transactions on
c**   mathematical software, vol 18, no 4, december 1992, pp. 454-455,
c**   by J. Leva.
c**   This implementation uses Tausworthe random number generator for
c**   uniformly distributed random numbers.

      function randn(i1,i2)
      implicit real (a-h,o-z)
      integer*4 i1,i2

      data s,t,a,b / 0.449871, -0.386595, 0.19600, 0.25472 /
      data r1,r2 / 0.27597, 0.27846 /

50    u=combtaus(i1,i2)
      v=combtaus(i1,i2)
      v=1.7156*(v-0.5)
      x=u-s
      y=abs(v)-t
      q=x**2 + y*(a*y - b*x)
      if (q .lt. r1) go to 100
      if (q .gt. r2) go to 50
      if (v**2 .gt. -4.0*alog(u)*u**2) go to 50
100   randn=v/u
      return

      end