# Normal random number generator

**URL:** <https://fortran-lang.discourse.group/t/normal-random-number-generator/3724>\
**Category:** Help\
**Created:** [June 11, 2022, 2:10am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724 "2022-06-11T02:10:20Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![fortran4r](https://avatars.discourse-cdn.com/v4/letter/f/9de0a6/32.png) [@fortran4r](https://fortran-lang.discourse.group/u/fortran4r)\
**Post date:** [June 11, 2022, 2:10am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/1 "2022-06-11T02:10:20Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![alozada](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/alozada/32/298_2.png) [@alozada](https://fortran-lang.discourse.group/u/alozada)\
**Post date:** [June 11, 2022, 2:51am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/2 "2022-06-11T02:51:44Z")

</div>

Maybe, Box-Muller transformation

```auto
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:

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

```

_U1,2 ~Uniform[0,1)_  
_x ~ N[mu, sigma^2]_

Reference [here](http://www.giuseppeforte.me/uploads/7/9/7/6/79768858/boxmul.pdf)

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [June 11, 2022, 2:53am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/3 "2022-06-11T02:53:33Z")

</div>

Here is a code adapted from Alan Miller’s [rnorm.f90](https://jblevins.org/mirror/amiller/rnorm.f90) to use double precision:

```auto
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:

```auto
       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

```

---

<div class="post-metadata">

**Author:** ![fortran4r](https://avatars.discourse-cdn.com/v4/letter/f/9de0a6/32.png) [@fortran4r](https://fortran-lang.discourse.group/u/fortran4r)\
**Post date:** [June 11, 2022, 4:08am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/4 "2022-06-11T04:08:47Z")

</div>

What if U1 = 0?

---

<div class="post-metadata">

**Author:** ![alozada](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/alozada/32/298_2.png) [@alozada](https://fortran-lang.discourse.group/u/alozada)\
**Post date:** [June 11, 2022, 4:41am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/5 "2022-06-11T04:41:23Z")

</div>

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

![box_muller](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/c/c6100a4af0cd569d2edb475f5ee5c90e309f28ea.png)

---

<div class="post-metadata">

**Author:** ![urbanjost](https://avatars.discourse-cdn.com/v4/letter/u/0ea827/32.png) [@urbanjost](https://fortran-lang.discourse.group/u/urbanjost)\
**Post date:** [June 11, 2022, 4:55am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/6 "2022-06-11T04:55:15Z")

</div>

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

> **[GitHub - urbanjost/M\_datapac: The NIST DATAPAC package modularized and made...](https://github.com/urbanjost/M_datapac)**
>
> The NIST DATAPAC package modularized and made available as an fpm(1) package - GitHub - urbanjost/M\_datapac: The NIST DATAPAC package modularized and made available as an fpm(1) package

in particular

[https://urbanjost.github.io/M\_datapac/norran.3m\_datapac.html](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.

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [June 11, 2022, 1:55pm UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/7 "2022-06-11T13:55:05Z")

</div>

I updated the previous code to add the Box-Muller method and the [ratio of uniforms method](https://jblevins.org/mirror/amiller/ran_norm.f90) 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

```auto
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)
! 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_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**
>
> ```auto
> 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
> 
> ```

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [June 11, 2022, 10:47pm UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/8 "2022-06-11T22:47:57Z")

</div>

One way to get normal deviates is to use the inverse error function, see here: [Inverse error function - MATLAB erfinv](https://mathworks.com/help/matlab/ref/erfinv.html#bup2lvj).

This method is also discussed in the following paper:

- [Approximating the `erfinv` function | Mike Giles](https://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf) (preprint)
- [Approximating the erfinv Function - ScienceDirect](https://www.sciencedirect.com/science/article/pii/B9780123859631000101)

---

<div class="post-metadata">

**Author:** ![fortran4r](https://avatars.discourse-cdn.com/v4/letter/f/9de0a6/32.png) [@fortran4r](https://fortran-lang.discourse.group/u/fortran4r)\
**Post date:** [June 12, 2022, 2:59am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/9 "2022-06-12T02:59:17Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [June 12, 2022, 3:30am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/10 "2022-06-12T03:30:33Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [June 12, 2022, 10:48am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/11 "2022-06-12T10:48:41Z")

</div>

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](https://www.nag.com/people/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](https://www.ams.org/journals/mcom/1976-30-136/S0025-5718-1976-0421040-7/S0025-5718-1976-0421040-7.pdf), but the tables of polynomial coefficients are missing. You can however orient yourself by the Julia implementantion: [https://github.com/JuliaMath/SpecialFunctions.jl/blob/49be7e6add22057776b5ed8ae32e0dae84500203/src/erf.jl#L242](https://github.com/JuliaMath/SpecialFunctions.jl/blob/49be7e6add22057776b5ed8ae32e0dae84500203/src/erf.jl#L242)

A draft Fortran attempt can be found here: [List of special mathematical functions to include · Issue #179 · fortran-lang/stdlib · GitHub](https://github.com/fortran-lang/stdlib/issues/179#issuecomment-647178168). 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.

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [June 12, 2022, 12:01pm UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/12 "2022-06-12T12:01:29Z")

</div>

I’ve located the routines from Giles in a third-party GitHub repository: [lwafomps/erfinv.cpp at 1ff3fbbdb23d5772c1c45d201ba050e4819a23e8 · MersenneTwister-Lab/lwafomps · GitHub](https://github.com/MersenneTwister-Lab/lwafomps/blob/1ff3fbbdb23d5772c1c45d201ba050e4819a23e8/src/erfinv.cpp#L136). A few (related) C codes and discussion can also be found in this thread: [Inverse Error Function in C - Stack Overflow](https://stackoverflow.com/questions/27229371/inverse-error-function-in-c)

I would also note [`v?ErfInv`](https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/vector-mathematical-functions/vm-mathematical-functions/special-functions/v-erfinv.html) in the Intel oneAPI MKL Vector Mathematical Functions library.

Speaking of Intel MKL, the Statistical Functions include distribution generators including [`v?RngGaussian`](https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/statistical-functions/random-number-generators/distribution-generators/continuous-distributions/vrnggaussian.html). Methods based on the [Box-Muller transform](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform) or the [inverse cumulative distribution function](https://en.wikipedia.org/wiki/Quantile_function) are both available.

If you don’t mind linking with the C++ standard library, there you have: [normal\_distribution - C++ Reference](https://cplusplus.com/reference/random/normal_distribution/)

Finally, if you’re using R, you could also use R’s internal pseudo random number generator as described in [Writing R Extensions](https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Random-numbers). The interfaces needed to call these from Fortran are very simple:

```fortran
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](https://colinfay.me/r-installation-administration/the-standalone-rmath-library.html).

---

<div class="post-metadata">

**Author:** ![septc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/septc/32/77_2.png) [@septc](https://fortran-lang.discourse.group/u/septc)\
**Post date:** [June 12, 2022, 2:17pm UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/13 "2022-06-12T14:17:05Z")

</div>

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](https://www.gnu.org/software/gsl/doc/html/randist.html#the-gaussian-distribution)

```c
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.
```

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [June 12, 2022, 3:18pm UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/14 "2022-06-12T15:18:42Z")

</div>

My adaptation of Alan Miller’s ziggurat Fortran code is discussed at [Pure random number generators](https://fortran-lang.discourse.group/t/pure-random-number-generators/2170) and is at [GitHub - Beliavsky/Ziggurat: generate uniform, normal, and exponential random deviates using algorithms from Marsaglia](https://github.com/Beliavsky/Ziggurat).

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [June 12, 2022, 3:42pm UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/15 "2022-06-12T15:42:52Z")

</div>

I’ve found that both Box-Muller and Ziggurat methods are discussed in the book [_Random Number Generators – Principles and Practices (2018)_](https://doi.org/10.1515/9781501506062) by David Johnston (Principal Engineer, Intel Corporation).

The NAG Library routine [`g05skf` (`dist_normal`)](https://www.nag.com/numeric/nl/nagdoc_latest/flhtml/g05/g05skf.html) uses the inverse cumulative distribution function approach with an approximation by [Wichura](https://doi.org/10.2307/2347330) (1988). A PDF is available [here (511 kb)](https://csg.sph.umich.edu/abecasis/gas_power_calculator/algorithm-as-241-the-percentage-points-of-the-normal-distribution.pdf).

Searching through ACM publications shows lots of results:

- [Algorithm 488: A Gaussian pseudo-random number generator](https://doi.org/10.1145/361604.361629) (1974)
- [Algorithm 712: A normal random number generator](https://doi.org/10.1145/138351.138367) (1992)
- [Fast pseudorandom generators for normal and exponential variates](https://dl.acm.org/doi/abs/10.1145/225545.225554) (1996)
- [An Improved Ziggurat Method to Generate Normal Random Samples](https://www.doornik.com/research/ziggurat.pdf) (2005)
- [Gaussian random number generators](https://doi.org/10.1145/1287620.1287622) (2007)

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

I’m surprised that we’ve all forgotten Fortran stdlib already has a function for normal variates: [stats\_distribution\_normal – Fortran-lang/stdlib](https://stdlib.fortran-lang.org/page/specs/stdlib_stats_distribution_normal.html)

```fortran
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
```

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [June 12, 2022, 9:49pm UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/16 "2022-06-12T21:49:14Z")

</div>

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.

```auto
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)
! 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_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

```

---

<div class="post-metadata">

**Author:** ![fortran4r](https://avatars.discourse-cdn.com/v4/letter/f/9de0a6/32.png) [@fortran4r](https://fortran-lang.discourse.group/u/fortran4r)\
**Post date:** [June 13, 2022, 12:57am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/17 "2022-06-13T00:57:00Z")

</div>

```auto
 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?

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [June 13, 2022, 1:31am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/18 "2022-06-13T01:31:40Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [June 13, 2022, 2:07am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/19 "2022-06-13T02:07:14Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![septc](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/septc/32/77_2.png) [@septc](https://fortran-lang.discourse.group/u/septc)\
**Post date:** [June 14, 2022, 1:58am UTC](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724/20 "2022-06-14T01:58:36Z")

</div>

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…

[Next page](https://fortran-lang.discourse.group/t/normal-random-number-generator/3724.md?page=2)