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.
References
-
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. http://doi.org/10.1214/aoms/1177706645
-
Leva, J. L. (1992). Algorithm 712; a normal random number generator. ACM Transactions on Mathematical Software, 18(4), 454–455. http://doi.org/10.1145/138351.138367