# Normal random number generator

How to create a subroutine based on random_number that fills an array with random numbers following a normal distribution with mean mu and standard deviation sigma?

Maybe, Box-Muller transformation

``````Zc = sqrt(-2.0 * log(U1)) * cos(2.0*pi*U2))
Zs = sqrt(-2.0 * log(U1))* sin(2.0*pi*U2))
``````

x = mu + sigma * Zc,s

example:

``````x = mu + sigma*sqrt(-2.0*log(U1))*cos(2*pi*U2)
``````

U1,2 ~Uniform[0,1)
x ~ N[mu, sigma^2]

Reference here

Here is a code adapted from Alan Millerâ€™s rnorm.f90 to use double precision:

``````module rnorm_mod
implicit none
private
public :: dp, rnorm, rnorm_vec
integer, parameter :: dp = kind(1.0d0)
contains
function rnorm_vec(n,mu,sigma) result(variates)
integer      , intent(in)           :: n
real(kind=dp), intent(in), optional :: mu, sigma
real(kind=dp)                       :: variates(n)
integer                             :: i
do i=1,n
variates(i) = rnorm()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_vec
!
FUNCTION rnorm() RESULT(fn_val)

!   Generate a random normal deviate using the polar method.
!   Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
!              normal variables', Siam Rev., vol.6, 260-264, 1964.

IMPLICIT NONE
REAL(kind=dp)  :: fn_val

! Local variables

REAL(kind=dp)   :: u, sum
REAL(kind=dp), SAVE      :: v, sln
LOGICAL, SAVE   :: second = .FALSE.
REAL(kind=dp), PARAMETER :: one = 1.0_dp, vsmall = TINY( one )

IF (second) THEN
! If second, use the second random number generated on last call

second = .false.
fn_val = v*sln

ELSE
! First call; generate a pair of random normals

second = .true.
DO
CALL RANDOM_NUMBER( u )
CALL RANDOM_NUMBER( v )
u = SCALE( u, 1 ) - one
v = SCALE( v, 1 ) - one
sum = u*u + v*v + vsmall         ! vsmall added to prevent LOG(zero) / zero
IF(sum < one) EXIT
END DO
sln = SQRT(-SCALE(LOG(sum),1)/sum)
fn_val = u*sln
END IF
END FUNCTION rnorm
end module rnorm_mod
!
program main
use rnorm_mod, only: dp, rnorm_vec
implicit none
real(kind=dp), allocatable :: x(:)
integer, parameter :: n = 10**7
integer :: ipow, iter
real(kind=dp), parameter :: mu = 2.0_dp, sigma = 3.0_dp
call random_seed()
allocate (x(n))
print "(*(a8))",   "n", "mu", "sigma"
print "(i8,2f8.4)", n , mu  ,  sigma
print "(/,'central moments',/,*(i10))",(ipow,ipow=1,4)
do iter=1,5
x = rnorm_vec(n,mu,sigma)
x = x - sum(x)/n
print "(*(f10.4))",(sum(x**ipow)/n,ipow=1,4)
end do
print*,"theoretical"
print "(*(f10.4))",0.0_dp,sigma**2,0.0_dp,3*sigma**4
end program main
``````

sample gfortran output:

``````       n      mu   sigma
10000000  2.0000  3.0000

central moments
1         2         3         4
-0.0000    8.9980    0.0211  242.9815
0.0000    9.0006   -0.0265  243.0072
0.0000    8.9938   -0.0025  242.6678
0.0000    8.9995   -0.0236  243.0038
-0.0000    8.9914    0.0200  242.3521
theoretical
0.0000    9.0000    0.0000  243.0000
``````

What if U1 = 0?

U ~ Uniform(0,1)
Ann. Math. Statist. 29(2): 610-611 (June, 1958). DOI: 10.1214/aoms/1177706645

1 Like

Still working on updating the NIST DATAPAC library, but it might be of interest

in particular

https://urbanjost.github.io/M_datapac/norran.3m_datapac.html

it still lacks a unit test, so caveat emptor. There are several dozen routines for various random distributions in the package.

1 Like

I updated the previous code to add the Box-Muller method and the ratio of uniforms method in addition to the polar method and to add timings. For Box-Muller I guard against zero uniform variates, which does not reduce speed significantly. For the code

``````module rnorm_mod
implicit none
private
public :: dp, rnorm_polar, rnorm_ratio_uni, rnorm_vec, rnorm_polar_vec, rnorm_ratio_uni_vec, &
rnorm_box_muller, rnorm_box_muller_vec
integer, parameter :: dp = kind(1.0d0)
real(kind=dp), parameter :: half = 0.5_dp, two_pi = 6.28318530718_dp
contains
function rnorm_vec(n,method,mu,sigma) result(variates)
integer      , intent(in)           :: n     ! # of normal variates
character (len=*), intent(in)       :: method
real(kind=dp), intent(in), optional :: mu    ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp)                       :: variates(n) ! normal variates
select case (method)
case ("polar")     ; variates = rnorm_polar_vec(n,mu,sigma)
case ("ratio_uni") ; variates = rnorm_ratio_uni_vec(n,mu,sigma)
case ("box_muller"); variates = rnorm_box_muller_vec(n,mu,sigma)
case default       ; error stop trim(method) // " is invalid"
end select
end function rnorm_vec
!
function rnorm_polar_vec(n,mu,sigma) result(variates)
integer      , intent(in)           :: n     ! # of normal variates
real(kind=dp), intent(in), optional :: mu    ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp)                       :: variates(n) ! normal variates
integer                             :: i
do i=1,n
variates(i) = rnorm_polar()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_polar_vec
!
function rnorm_ratio_uni_vec(n,mu,sigma) result(variates)
integer      , intent(in)           :: n     ! # of normal variates
real(kind=dp), intent(in), optional :: mu    ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp)                       :: variates(n) ! normal variates
integer                             :: i
do i=1,n
variates(i) = rnorm_ratio_uni()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_ratio_uni_vec
!
function rnorm_box_muller_vec(n,mu,sigma) result(variates)
integer      , intent(in)           :: n     ! # of normal variates
real(kind=dp), intent(in), optional :: mu    ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp)                       :: variates(n) ! normal variates
integer                             :: i,j
logical                             :: n_odd
n_odd = mod(n,2) /= 0
do i=1,n/2
j = 2*i - 1
variates(j:j+1) = rnorm_box_muller()
end do
if (n_odd) variates(n) = rnorm_box_muller_single_variate()
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_box_muller_vec
!
FUNCTION rnorm_ratio_uni() RESULT(fn_val) ! adapted from https://jblevins.org/mirror/amiller/ran_norm.f90
! 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.

real(kind=dp) :: fn_val
real(kind=dp), parameter :: s = 0.449871_dp, t = -0.386595_dp, a = 0.19600_dp, b = 0.25472_dp,    &
r1 = 0.27597_dp, r2 = 0.27846_dp
! local variables
real(kind=dp) :: u, v, x, y, q
!     Generate P = (u,v) uniform in rectangle enclosing acceptance region
DO
CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)
v = 1.7156_dp * (v - half)
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_dp*LOG(u)*u**2) EXIT
END DO
!     Return ratio of P's coordinates as the normal deviate
fn_val = v/u
END FUNCTION rnorm_ratio_uni
!
function rnorm_box_muller() result(variates) ! coded formulas from https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
! return two uncorrelated standard normal variates
real(kind=dp) :: variates(2)
real(kind=dp) :: u(2), factor, arg
do
call random_number(u)
if (u(1) > 0.0_dp) exit
end do
factor = sqrt(-2 * log(u(1)))
arg = two_pi*u(2)
variates = factor * [cos(arg),sin(arg)]
end function rnorm_box_muller
!
function rnorm_box_muller_single_variate() result(variate)
! return a standard normal variate
real(kind=dp) :: variate
real(kind=dp) :: u(2), factor, arg
call random_number(u)
factor = sqrt(-2 * log(u(1)))
arg = two_pi*u(2)
variate = factor * cos(arg)
end function rnorm_box_muller_single_variate
!
FUNCTION rnorm_polar() RESULT(fn_val) ! adapted from https://jblevins.org/mirror/amiller/rnorm.f90

!   Generate a random normal deviate using the polar method.
!   Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
!              normal variables', Siam Rev., vol.6, 260-264, 1964.

IMPLICIT NONE
REAL(kind=dp)  :: fn_val

! Local variables

REAL(kind=dp)   :: u, sum
REAL(kind=dp), SAVE      :: v, sln
LOGICAL, SAVE   :: second = .FALSE.
REAL(kind=dp), PARAMETER :: one = 1.0_dp, vsmall = TINY( one )

IF (second) THEN
! If second, use the second random number generated on last call

second = .false.
fn_val = v*sln

ELSE
! First call; generate a pair of random normals

second = .true.
DO
CALL RANDOM_NUMBER( u )
CALL RANDOM_NUMBER( v )
u = SCALE( u, 1 ) - one
v = SCALE( v, 1 ) - one
sum = u*u + v*v + vsmall         ! vsmall added to prevent LOG(zero) / zero
IF(sum < one) EXIT
END DO
sln = SQRT(-SCALE(LOG(sum),1)/sum)
fn_val = u*sln
END IF
END FUNCTION rnorm_polar
end module rnorm_mod
!
program main
use rnorm_mod, only: dp, rnorm_vec, rnorm_box_muller_vec
implicit none
integer          , parameter :: n = 10**8, niter = 5, nmethods = 3
real(kind=dp)    , parameter :: mu = 2.0_dp, sigma = 3.0_dp
character (len=*), parameter :: methods(nmethods) = ["polar     ", &
"ratio_uni ","box_muller"]
integer                      :: ipow, iter, imethod
real(kind=dp)                :: xmean, old_time, new_time, elapsed_time(nmethods)
real(kind=dp), allocatable   :: x(:)
call random_seed()
elapsed_time = 0.0_dp
allocate (x(n))
print "(*(a10))",    "n", "niter", "mu", "sigma"
print "(2i10,2f10.4)", n ,  niter,   mu  , sigma
do imethod=1,nmethods
print "(/,'method = ',a)",trim(methods(imethod))
print "('central moments',/,a10,*(i10))","mean",(ipow,ipow=1,4)
do iter=1,niter
call cpu_time(old_time)
x = rnorm_vec(n,methods(imethod),mu,sigma)
call cpu_time(new_time)
elapsed_time(imethod) = elapsed_time(imethod) + new_time - old_time
xmean = sum(x)/n
x = x - xmean
print "(*(f10.4))",xmean,(sum(x**ipow)/n,ipow=1,4)
end do
print*,"theoretical"
print "(*(f10.4))",mu,0.0_dp,sigma**2,0.0_dp,3*sigma**4
end do
print "(/,*(a15))", "method",(trim(methods(imethod)),imethod=1,nmethods)
print "(a15,*(f15.4))", "cpu_time",elapsed_time
print "(/,'check Box-Muller for n odd or even')"
print "(*(f10.4))",rnorm_box_muller_vec(3),rnorm_box_muller_vec(4)
end program main
``````

on WSL2 with `gfortran -O3` I get timings

``````     method          polar      ratio_uni     box_muller
cpu_time        13.0271        15.3004        12.3667
``````

and with `ifort -O3`

``````     method          polar      ratio_uni     box_muller
cpu_time        16.8491        20.3027         8.6018
``````
Full output
``````gfortran -O3 rnorm_dp.f90
(base) /mnt/c/fortran/test\$ time ./a.out
n     niter        mu     sigma
100000000         5    2.0000    3.0000

method = polar
central moments
mean         1         2         3         4
2.0002    0.0000    9.0006   -0.0003  243.0160
2.0001   -0.0000    8.9974    0.0124  242.8701
2.0000   -0.0000    9.0003   -0.0028  242.9531
2.0001    0.0000    8.9992   -0.0094  242.9653
1.9999    0.0000    9.0008   -0.0128  243.0604
theoretical
2.0000    0.0000    9.0000    0.0000  243.0000

method = ratio_uni
central moments
mean         1         2         3         4
1.9999   -0.0000    9.0011   -0.0086  243.0712
2.0003    0.0000    8.9992   -0.0098  242.9281
1.9999   -0.0000    9.0001    0.0071  243.0048
1.9998   -0.0000    8.9997   -0.0047  242.9843
2.0002   -0.0000    9.0011   -0.0098  243.0690
theoretical
2.0000    0.0000    9.0000    0.0000  243.0000

method = box_muller
central moments
mean         1         2         3         4
2.0000   -0.0000    9.0007    0.0096  243.1083
1.9992   -0.0000    9.0012    0.0070  243.0418
2.0000    0.0000    9.0015    0.0044  243.0983
2.0000   -0.0000    8.9992    0.0017  242.9084
2.0005    0.0000    9.0005   -0.0068  243.0650
theoretical
2.0000    0.0000    9.0000    0.0000  243.0000

method          polar      ratio_uni     box_muller
cpu_time        13.0271        15.3004        12.3667

check Box-Muller for n odd or even
-0.4509    0.6263    0.1658   -0.5268    0.2805    1.0316   -0.5751

real	0m58.591s
user	0m58.355s
sys	0m0.230s

ulimit -s unlimited
(base) /mnt/c/fortran/test\$ ifort -O3 rnorm_dp.f90
(base) /mnt/c/fortran/test\$ time ./a.out
n     niter        mu     sigma
100000000         5    2.0000    3.0000

method = polar
central moments
mean         1         2         3         4
2.0002   -0.0000    8.9993    0.0057  242.9807
2.0002    0.0000    8.9988   -0.0029  242.8929
2.0006    0.0000    8.9994    0.0061  243.0622
2.0007   -0.0000    9.0016   -0.0039  243.1514
2.0002    0.0000    8.9988   -0.0034  242.8767
theoretical
2.0000    0.0000    9.0000    0.0000  243.0000

method = ratio_uni
central moments
mean         1         2         3         4
2.0001   -0.0000    9.0002    0.0071  242.9891
2.0000    0.0000    9.0012    0.0030  243.0431
2.0003   -0.0000    8.9997    0.0047  242.9036
1.9995    0.0000    9.0012    0.0052  243.0759
1.9997   -0.0000    9.0008    0.0046  242.9873
theoretical
2.0000    0.0000    9.0000    0.0000  243.0000

method = box_muller
central moments
mean         1         2         3         4
2.0002   -0.0000    8.9980   -0.0033  242.9163
2.0003    0.0000    8.9983   -0.0006  242.8713
1.9994    0.0000    8.9998   -0.0070  243.0213
1.9998    0.0000    9.0010    0.0029  243.1056
2.0002    0.0000    8.9998   -0.0025  243.0385
theoretical
2.0000    0.0000    9.0000    0.0000  243.0000

method          polar      ratio_uni     box_muller
cpu_time        16.8491        20.3027         8.6018

check Box-Muller for n odd or even
-0.4093   -0.5503   -0.5804    0.6342    0.7505   -1.8514   -0.9049

real	1m30.322s
user	1m30.162s
sys	0m0.141s
``````
2 Likes

One way to get normal deviates is to use the inverse error function, see here: Inverse error function - MATLAB erfinv.

This method is also discussed in the following paper:

I didnâ€™t know that Fortran has an erfinv function. Can you be more specific? I only know how to use the erfc function to get normcdf.

The suggested function is not a standard intrinsic function, but a user-provided function. One of the references provides approximations for erfinv() that can be coded fairly easily.

1 Like

The method from Giles I linked above has been designed to run fast on NVIDIA GPUâ€™s. If Iâ€™m not mistaken, the author, Mike Giles, is somehow associated with NAG. The paper also states:

These approximations were originally developed as part of a commercial maths library providing Numerical Routines for GPUs.

There exist also older approximations, that perform well on CPUâ€™s, such as

Blair, J. M., Edwards, C. A., & Johnson, J. H. (1976). Rational Chebyshev approximations for the inverse of the error function. Mathematics of computation, 30(136), 827-830.

A copy of that paper can be found here, but the tables of polynomial coefficients are missing. You can however orient yourself by the Julia implementantion: SpecialFunctions.jl/erf.jl at 49be7e6add22057776b5ed8ae32e0dae84500203 Â· JuliaMath/SpecialFunctions.jl Â· GitHub

A draft Fortran attempt can be found here: List of special mathematical functions to include Â· Issue #179 Â· fortran-lang/stdlib Â· GitHub. Unfortunately, I havenâ€™t done any proper testing of the accuracy or speed. The Julia special functions library provides both single and double precision versions.

Iâ€™ve located the routines from Giles in a third-party GitHub repository: lwafomps/erfinv.cpp at 1ff3fbbdb23d5772c1c45d201ba050e4819a23e8 Â· MersenneTwister-Lab/lwafomps Â· GitHub. A few (related) C codes and discussion can also be found in this thread: Inverse Error Function in C - Stack Overflow

I would also note `v?ErfInv` in the Intel oneAPI MKL Vector Mathematical Functions library.

Speaking of Intel MKL, the Statistical Functions include distribution generators including `v?RngGaussian`. Methods based on the Box-Muller transform or the inverse cumulative distribution function are both available.

If you donâ€™t mind linking with the C++ standard library, there you have: normal_distribution - C++ Reference

Finally, if youâ€™re using R, you could also use Râ€™s internal pseudo random number generator as described in Writing R Extensions. The interfaces needed to call these from Fortran are very simple:

``````module r_random

implicit none

interface
subroutine r_getrngstate() bind(c,name="GetRNGstate")
end subroutine
subroutine r_putrngstate() bind(c,name="PutRNGstate")
end subroutine
real(c_double) function r_unif_rand() bind(c,name="unif_rand")
use, intrinsic :: iso_c_binding, only: c_double
end function
real(c_double) function r_norm_rand() bind(c,name="norm_rand")
use, intrinsic :: iso_c_binding, only: c_double
end function
real(c_double) function r_exp_rand() bind(c,name="exp_rand")
use, intrinsic :: iso_c_binding, only: c_double
end function
end interface

end module
``````

You still have to figure out how to link your Fortran code with Râ€™s standalone math library.

1 Like

I also use the Box-Muller method in my codes, but wondering if there are better (= more efficient) methods. In addition to the methods mentioned above, is the â€śzigguratâ€ť method another option? (Iâ€™ve never used it myself, so not sure whether it is good or notâ€¦) GSL also seems to provide both routines:

https://www.gnu.org/software/gsl/doc/html/randist.html#the-gaussian-distribution

``````double gsl_ran_gaussian_ziggurat(const gsl_rng *r, double sigma)

This function computes a Gaussian random variate using
the alternative Marsaglia-Tsang ziggurat and Kinderman-Monahan-
Leva ratio methods. The Ziggurat algorithm is the fastest
available algorithm in most cases.``````

My adaptation of Alan Millerâ€™s ziggurat Fortran code is discussed at Pure random number generators and is at GitHub - Beliavsky/Ziggurat: generate uniform, normal, and exponential random deviates using algorithms from Marsaglia.

1 Like

Iâ€™ve found that both Box-Muller and Ziggurat methods are discussed in the book Random Number Generators â€“ Principles and Practices (2018) by David Johnston (Principal Engineer, Intel Corporation).

The NAG Library routine `g05skf` (`dist_normal`) uses the inverse cumulative distribution function approach with an approximation by Wichura (1988). A PDF is available here (511 kb).

Searching through ACM publications shows lots of results:

Some comments on Ziggurat can also be found in the stdlib thread: Proposal for statistical distributions Â· Issue #234 Â· fortran-lang/stdlib Â· GitHub. Another discussion is in the NumPy pull request: Patch with Ziggurat method for Normal distribution (Trac #1450) Â· Issue #2047 Â· numpy/numpy Â· GitHub

Iâ€™m surprised that weâ€™ve all forgotten Fortran stdlib already has a function for normal variates: stats_distribution_normal â€“ Fortran-lang/stdlib

``````program demo_normal_rvs
use stdlib_random, only: random_seed
use stdlib_stats_distribution_normal, only: norm => rvs_normal

implicit none
integer :: seed_put, seed_get

seed_put = 1234567
call random_seed(seed_put, seed_get)

print *, norm( )            ! single standard normal random variate
print *, norm(1.0, 2.0)     ! normal random variate mu=1.0, sigma=2.0
print *, norm(0.0, 1.0, 10) ! an array of 10 standard norml random variates

end program``````
1 Like

Most algorithms for generating normal variates rely on uniform variates. It is possible to write the module with procedures for normal variates so that toggling a single line of code determines whether the intrinsic or user-defined `random_number` is used, as shown below. A benefit of the latter is that the normal variates will be the same across compilers.

``````module random_number_mod
implicit none
private
public :: dp, set_seed, random_number
integer, parameter :: dp = kind(1.0d0)
real(kind=dp), save :: ds = 123.0_dp
contains
!
subroutine set_seed(seed)
real(kind=dp), intent(in) :: seed
ds = seed
end subroutine set_seed
!
impure elemental subroutine random_number(harvest)
! adapted from function ggl in http://www.netlib.org/random/2drwtest.f
real(kind=dp), intent(out) :: harvest
real(kind=dp), parameter :: d2 = 2147483647.d0
ds  = dmod(16807.d0*ds,d2)
!     Generate U(0,1] distributed random numbers:
harvest = ds/d2
end subroutine random_number
end module random_number_mod
!
module random_normal_mod
! comment the line below to use the intrinsic random_number
use random_number_mod, only: random_number
implicit none
private
public :: dp, rnorm_polar, rnorm_ratio_uni, rnorm_vec, rnorm_polar_vec, rnorm_ratio_uni_vec, &
rnorm_box_muller, rnorm_box_muller_vec
integer, parameter :: dp = kind(1.0d0)
real(kind=dp), parameter :: half = 0.5_dp, two_pi = 6.28318530718_dp
contains
function rnorm_vec(n,method,mu,sigma) result(variates)
integer      , intent(in)           :: n     ! # of normal variates
character (len=*), intent(in)       :: method
real(kind=dp), intent(in), optional :: mu    ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp)                       :: variates(n) ! normal variates
select case (method)
case ("polar")     ; variates = rnorm_polar_vec(n,mu,sigma)
case ("ratio_uni") ; variates = rnorm_ratio_uni_vec(n,mu,sigma)
case ("box_muller"); variates = rnorm_box_muller_vec(n,mu,sigma)
case default       ; error stop trim(method) // " is invalid"
end select
end function rnorm_vec
!
function rnorm_polar_vec(n,mu,sigma) result(variates)
integer      , intent(in)           :: n     ! # of normal variates
real(kind=dp), intent(in), optional :: mu    ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp)                       :: variates(n) ! normal variates
integer                             :: i
do i=1,n
variates(i) = rnorm_polar()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_polar_vec
!
function rnorm_ratio_uni_vec(n,mu,sigma) result(variates)
integer      , intent(in)           :: n     ! # of normal variates
real(kind=dp), intent(in), optional :: mu    ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp)                       :: variates(n) ! normal variates
integer                             :: i
do i=1,n
variates(i) = rnorm_ratio_uni()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_ratio_uni_vec
!
function rnorm_box_muller_vec(n,mu,sigma) result(variates)
integer      , intent(in)           :: n     ! # of normal variates
real(kind=dp), intent(in), optional :: mu    ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp)                       :: variates(n) ! normal variates
integer                             :: i,j
logical                             :: n_odd
n_odd = mod(n,2) /= 0
do i=1,n/2
j = 2*i - 1
variates(j:j+1) = rnorm_box_muller()
end do
if (n_odd) variates(n) = rnorm_box_muller_single_variate()
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_box_muller_vec
!
FUNCTION rnorm_ratio_uni() RESULT(fn_val) ! adapted from https://jblevins.org/mirror/amiller/ran_norm.f90
! 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.

real(kind=dp) :: fn_val
real(kind=dp), parameter :: s = 0.449871_dp, t = -0.386595_dp, a = 0.19600_dp, b = 0.25472_dp,    &
r1 = 0.27597_dp, r2 = 0.27846_dp
! local variables
real(kind=dp) :: u, v, x, y, q
!     Generate P = (u,v) uniform in rectangle enclosing acceptance region
DO
CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)
v = 1.7156_dp * (v - half)
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_dp*LOG(u)*u**2) EXIT
END DO
!     Return ratio of P's coordinates as the normal deviate
fn_val = v/u
END FUNCTION rnorm_ratio_uni
!
function rnorm_box_muller() result(variates) ! coded formulas from https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
! return two uncorrelated standard normal variates
real(kind=dp) :: variates(2)
real(kind=dp) :: u(2), factor, arg
do
call random_number(u)
if (u(1) > 0.0_dp) exit
end do
factor = sqrt(-2 * log(u(1)))
arg = two_pi*u(2)
variates = factor * [cos(arg),sin(arg)]
end function rnorm_box_muller
!
function rnorm_box_muller_single_variate() result(variate)
! return a standard normal variate
real(kind=dp) :: variate
real(kind=dp) :: u(2), factor, arg
do
call random_number(u)
if (u(1) > 0.0_dp) exit
end do
factor = sqrt(-2 * log(u(1)))
arg = two_pi*u(2)
variate = factor * cos(arg)
end function rnorm_box_muller_single_variate
!
FUNCTION rnorm_polar() RESULT(fn_val) ! adapted from https://jblevins.org/mirror/amiller/rnorm.f90

!   Generate a random normal deviate using the polar method.
!   Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
!              normal variables', Siam Rev., vol.6, 260-264, 1964.

IMPLICIT NONE
REAL(kind=dp)  :: fn_val

! Local variables

REAL(kind=dp)   :: u, sum
REAL(kind=dp), SAVE      :: v, sln
LOGICAL, SAVE   :: second = .FALSE.
REAL(kind=dp), PARAMETER :: one = 1.0_dp, vsmall = TINY( one )

IF (second) THEN
! If second, use the second random number generated on last call

second = .false.
fn_val = v*sln

ELSE
! First call; generate a pair of random normals

second = .true.
DO
CALL RANDOM_NUMBER( u )
CALL RANDOM_NUMBER( v )
u = SCALE( u, 1 ) - one
v = SCALE( v, 1 ) - one
sum = u*u + v*v + vsmall         ! vsmall added to prevent LOG(zero) / zero
IF(sum < one) EXIT
END DO
sln = SQRT(-SCALE(LOG(sum),1)/sum)
fn_val = u*sln
END IF
END FUNCTION rnorm_polar
end module random_normal_mod
!
program main
use random_normal_mod, only: dp, rnorm_vec, rnorm_box_muller_vec
implicit none
integer          , parameter :: n = 10**6, niter = 5, nmethods = 3
real(kind=dp)    , parameter :: mu = 2.0_dp, sigma = 3.0_dp
character (len=*), parameter :: methods(nmethods) = ["polar     ", &
"ratio_uni ","box_muller"]
integer                      :: ipow, iter, imethod
real(kind=dp)                :: xmean, old_time, new_time, elapsed_time(nmethods)
real(kind=dp), allocatable   :: x(:)
call random_seed()
elapsed_time = 0.0_dp
allocate (x(n))
print "(*(a10))",    "n", "niter", "mu", "sigma"
print "(2i10,2f10.4)", n ,  niter,   mu  , sigma
do imethod=1,nmethods
print "(/,'method = ',a)",trim(methods(imethod))
print "('central moments',/,a10,*(i10))","mean",(ipow,ipow=1,4)
do iter=1,niter
call cpu_time(old_time)
x = rnorm_vec(n,methods(imethod),mu,sigma)
call cpu_time(new_time)
elapsed_time(imethod) = elapsed_time(imethod) + new_time - old_time
xmean = sum(x)/n
x = x - xmean
print "(*(f10.4))",xmean,(sum(x**ipow)/n,ipow=1,4)
end do
print*,"theoretical"
print "(*(f10.4))",mu,0.0_dp,sigma**2,0.0_dp,3*sigma**4
end do
print "(/,*(a15))", "method",(trim(methods(imethod)),imethod=1,nmethods)
print "(a15,*(f15.4))", "cpu_time",elapsed_time
print "(/,'check Box-Muller for n odd or even')"
print "(*(f10.4))",rnorm_box_muller_vec(3),rnorm_box_muller_vec(4)
end program main
``````
`````` subroutine rnorm(x, mu, sigma)
integer(8) :: n
real(8) :: x(:), mu, sigma
real(8), allocatable :: u1(:), u2(:)
n = size(x, kind = 8)
allocate(u1(n), u2(n))
call random_number(u1)
call random_number(u2)
u1 = 1d0 - u1
u2 = 1d0 - u2
x = mu + sigma*sqrt(-2d0*log(u1))*cos(2d0*pi*u2)
deallocate(u1, u2)
end subroutine
end
``````

I followed your reference to write this subroutine, but it is slow. Any advice on how to make it faster?

My `function rnorm_box_muller` generates two normal variates rather than one without doing twice as many computations, by computing a `sin` as well as a `cos`.

Many math libraries have a subroutine that computes both sine and cosine values of the argument. In addition to the range reduction, some of the polynomial evaluation steps can be shared, making the total effort less than the sum of the effort for the two separate evaluations.

I wonder what is the typical value of `n` (the size of array `x`)? If it is very large (such that it requires 64-bit integer for the size), I guess allocating the two temporary arrays `u1` and `u2` may be very time-consumingâ€¦