random_normal Function

public function random_normal()

Arguments

None

Return Value real(kind=r8)

Note

Adapted from the following fortran 77 code algorithm 712, collected algorithms from acm. this work published in transactions on mathematical software, vol. 18, no. 4, december, 1992, pp. 434-435.

  • The function random_normal() returns a normally distributed pseudo-random number with zero mean and unit variance.

  • The algorithm uses the ratio of uniforms method of a.j. kinderman and j.f. monahan augmented with quadratic bounding curves.

    Author:

    • Alan Miller
    • csiro division of mathematical & information sciences
    • private bag 10, clayton south mdc
    • clayton 3169, victoria, australia

Source Code

   function random_normal()
   implicit none
   !================================================================================================
   !< @note
   !<
   !< Adapted from the following fortran 77 code
   !<      algorithm 712, collected algorithms from acm.
   !<      this work published in transactions on mathematical software,
   !<      vol. 18, no. 4, december, 1992, pp. 434-435.
   !<
   !<  + The function random_normal() returns a normally distributed pseudo-random
   !<     number with zero mean and unit variance.
   !<  + The algorithm uses the ratio of uniforms method of a.j. kinderman
   !<    and j.f. monahan augmented with quadratic bounding curves.
   !<
   !<     Author:
   !<
   !<     + Alan Miller
   !<     + csiro division of mathematical & information sciences
   !<     + private bag 10, clayton south mdc
   !<     + clayton 3169, victoria, australia
   !<
   !< @endnote
   !------------------------------------------------------------------------------------------------
      real(kind=r8) :: random_normal
      real(kind=r8) :: s, t, a, b, r1, r2, u, v, x, y, q

      s  =  0.449871_r8
      t  = -0.386595_r8
      a  =  0.19600_r8
      b  =  0.25472_r8
      r1 =  0.27597_r8
      r2 =  0.27846_r8

      !     generate p = (u,v) uniform in rectangle enclosing acceptance region

      do

        call random_number(u)
        call random_number(v)

        v = 1.7156 * (v - 0.5_r8)

      !     evaluate the quadratic form
        x = u - s
        y = abs(v) - t
        q = x**2 + y*(a*y - b*x)

      !     accept p if inside inner ellipse
        if (q < r1) exit
      !     reject p if outside outer ellipse
        if (q > r2) cycle
      !     reject p if outside acceptance region
        if (v**2 < -4.0*log(u)*u**2) exit

      enddo

      !     return ratio of p's coordinates as the normal deviate
      random_normal = v/u

   return
   endfunction random_normal