Sergey Budaev

Jun 13, 2018

Gaussian random numbers in Fortran

The HEDTOOLS tools library has a module for working with random numbers BASE_RANDOM​. There is, in particular, a set of procedures for generating Gaussian random values: ​RNORM and RNORM_ARRAY. These are based on the Kinderman & Monahan, augmented with quadratic bounding curves method (Leva, 1992: algorithm 712, Trans. Math. Software, 18, 4, 434-435​).

I have made a quick comparison of the quality of the Gaussian random numbers generated by the simple Box-Muller method (Box & Muller, 1958​)

Classical (ancient) Fortran code:

normrand_number = dsqrt(-2.*dlog(drand(0)))*dcos(2.*pi*drand(0))

that has been used in TEG codes so far...

and the algorithm 712 as implemented in HEDTOOLS using this test program (see attachment).

Fortran code for the test program:

program test_bm
  use csv_io
  use base_random, rand_x => rand     ! Alias rand() as rand_x() for ifort.
  !use IFPORT, only : rand_x => rand  ! This is the Intel Fortran tweak.

  integer, parameter :: prec = 8, arrsize=100000
  character(len=255), parameter :: filename1="file_01.csv", filename2="file_02.csv"
  real(kind=prec), dimension(arrsize) :: norand1, norand2
  real :: timer_start, timer_end
  ! Generating Box-Muller random numbers
  call cpu_time(timer_start)  ! START
  do i=1, arrsize
    norand1(i) = sqrt(-2.*log(rand_x(0)))*cos(2.*pi*rand_x(0))
  end do
  call cpu_time(timer_end)    ! END
  print *, "Box-Muller took: ", timer_end - timer_start
  ! Write random normal data to CSV
  call CSV_MATRIX_WRITE(norand1, filename1)

  ! Generating based on algorithm 712
  call cpu_time(timer_start)  ! START
  call RNORM_ARRAY(norand2)
  call cpu_time(timer_end)    ! END
  print *, "Alg. 712 took: ", timer_end - timer_start
  ! Write random normal data to CSV
  call CSV_MATRIX_WRITE(norand2, filename2)

end program test_bm

Comparison of Box-Muller and A712

The alg. 712 looks slightly faster than the simple Box-Muller transform.

alg. 712 is much better, as the Box-Muller significantly deviates from the normal distribution, alg. 712 does not (using the Anderson-Darling test from the nortest R package).

  # Gaussian random numbers by Box-Muller deviate from the Normal distribution:
  > ad.test(data_bm$X1)
      Anderson-Darling normality test
  data:  data_bm$X1
  A = 581.7, p-value < 2.2e-16
  # Gaussian random numbers by Kinderman & Monahan's A712 do not deviate from the Normal distribution:
  > ad.test(data_a712$X1)
      Anderson-Darling normality test
  data:  data_a712$X1
  A = 0.46975, p-value = 0.2474

So, the alg. 712 procedure implemented in HEDTOOLS should be used instead of the Box-Muller method.


  • Box, G. E. P., & Muller, M. E. (1958). A note on the generation of random normal deviates. The Annals of Mathematical Statistics, 29(2), 610–611. ​

  • Leva, J. L. (1992). Algorithm 712; a normal random number generator. ACM Transactions on Mathematical Software, 18(4), 454–455. ​