Hey, i just started learning fortran and my advisor sent me this code, it’s a pseudo random number generator but i don´t really know how to implement this or how it works, can someone explain the basics pls?
DOUBLE PRECISION FUNCTION RF()
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
data IX1,IX2,IX3,IX4,IX5,IX6/1500419,1400159,1364,1528,1,3/
RR1=1.0/FLOAT(IX1)
RR2=1.0/FLOAT(IX2)
IX5=MOD(IX5*IX3,IX1)
IX6=MOD(IX6*IX4,IX2)
RF=RR1*IX5+RR2*IX6
IF(RF.GE.1.0)RF=RF-1.0
END
Welcome to the forum. Do you really need to understand the workings of the random number generator, or do you just need to be able to use it? If it’s the latter, here is a program illustrating its use. The rf() function generates random uniform variates between 0 and 1.
module rf_mod
contains
DOUBLE PRECISION FUNCTION RF()
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
data IX1,IX2,IX3,IX4,IX5,IX6/1500419,1400159,1364,1528,1,3/
RR1=1.0/FLOAT(IX1)
RR2=1.0/FLOAT(IX2)
IX5=MOD(IX5*IX3,IX1)
IX6=MOD(IX6*IX4,IX2)
RF=RR1*IX5+RR2*IX6
IF(RF.GE.1.0)RF=RF-1.0
END FUNCTION
end module rf_mod
program main
use rf_mod, only: rf
implicit none
integer :: i
integer, parameter :: n = 10**7, dp = kind(1.0d0)
real(kind=dp) :: x(n)
! real(kind=dp) is another way of saying double precision
do i=1,n
x(i) = rf()
end do
print*,"theory",0.5_dp,1/12.0_dp
print*,"actual",sum(x)/n,sum((x-0.5)**2)/n
end program main
sample output:
theory 0.50000000000000000 8.3333333333333329E-002
actual 0.50019079377094988 8.3373676712904451E-002
just need to understand how to use it, i’m not used to this language yet so i’m still learning on how to call a function and so on. But thank you so much!
You did not say which programming language you normally use, and it is not clear what you need to “implement”. The RNG is already implemented in the Fortran 77 code that you posted.
What is your purpose/intent? To understand the code, to investigate the nature of the RNG, or simply obtain a random sequence, in which case you can use the intrinsic subroutine RANDOM_NUMBER()?
What your instructor gave you is code written in FORTRAN 77, which was great for its era but it is/should be used for backwards compatibility today. The language provides a built-in random number generator for decades now, and it has been improved over the years as well. Even if the code is meant to be for studying the general idea behind random generators, it should be written in modern Fortran instead, for better readability.
Amazing how instructors still use 77 syntax nowadays - and then we wonder why many students think Fortran is obsolete - but because they are taught the language as it was 50+ years ago, of course…
To play devil’s advocate, should random_number be used in production codes? Its implementation varies across compilers and may change from version of a compiler to the next. This limits reproducibility.
Some implementations of random_number give different sequences in repeated runs. The compiler’s runtime library may read the clock, etc., to generate and apply a different seed at the beginning of each run. As @Beliavsky noted, when debugging a program using more than one compiler, it may become necessary to use a personal RNG in the earlier stages of development.
Correct. If I were to implement a serious Monte Carlo application, I would definitely use a custom RNG. But I think the built-in one is “good enough” when you just need some random numbers and you don’t care too much about the period of the RNG used or how uniform the distribution is. At any rate I doubt the RNG of the code above is better, not to mention Fortan makes it easy to randomize the seed as well (of course, one could wonder how random the seed is that way - oh well, there is no perfect RNG).
Ha! Hence my reference to George Masaglia. His DIEHARD code is still a standard used to assess how random implementations really are. I have a partly converted (F90) version somewhere, if it would help.
i am more familiarized with python, now i am at a statistical thermodinamics research group and we work with monte carlo simulations so that’s why my advisor prefers to use his own generator instead of fortran’s. About how to implement i mean on how to call a function exatctly in fortran, what are the parameters etc, but now i think i got it
The gfortran compiler now has a very good PRNG. In fact, it is good enough, with several appealing features, that it might be time to consider standardizing that algorithm in fortran. For those who want to use other algorithms, they could continue to roll their own, but at least the fortran algorithm and its characteristics would be well defined. However, I doubt that many people need anything better than the gfortran algorithm. Also, the fortran interface to the intrinsic routine is very flexible, allowing for reproducible debugging, non reproducible production runs, and also all within sequential and parallel environments.
After failures in my attempts to build the Diehard and Dieharder RNG tests on Windows+Cygwin, I turned to the newer L’Ecuyer - Simard TestU01 suite, which I was able to build and run without problems on the RNG posted by the OP. One can read this review paper for a perspective on testing RNGS.
Here is the test program, in C, which is to be linked with the TestU01 library.
#include "unif01.h"
#include "bbattery.h"
double rf(void){
#define IX1 1500419
#define IX2 1400159
const static int ix3 = 1364, ix4 = 1528;
static int ix5 = 1, ix6 = 3;
const static double rr1 = 1.0/IX1, rr2 = 1.0/IX2;
double x;
ix5 = (ix5*ix3) % IX1;
ix6 = (ix6*ix4) % IX2;
x = rr1*ix5+rr2*ix6;
return x < 1 ? x : x-1.0;
}
int main (void)
{
unif01_Gen *gen;
gen = unif01_CreateExternGen01 ("FLD12", rf);
bbattery_SmallCrush (gen);
unif01_DeleteExternGen01 (gen);
return 0;
}
The result of the test run was that the OP’s advisor’s RNG did quite well, passing 14 out of 15 tests in the SmallCrush battery (the test that it failed was the BirthdaySpacings test). This RNG also passed all the tests in the PseudoDiehard battery.
I applied the same test batteries to the infamous RANDU RNG. RANDU failed all 15 tests in the SmallCrush battery and 8 of the 10 tests in the PseudoDiehard battery.
I am not qualified to comment about these tests except as a consumer, but I find solace in this quote of L’Ecuyer:
Bad RNGs are those that fail simple tests, whereas good RNGs fail only complicated tests that are hard to find and run.
@mecej4 Thank you for the link to the newer randomness tests. The other thing to sometimes remember about RNG is the ability to produce the same set of random numbers portably. One example might be where the objective function of an optimisation is based on a Monte Carlo simulation and we wish (have to) deploy to a cluster, GPU etc.
These two lines in the original post serve a special purpose. The first of the two maps the pair of random integers (IX5, IX6) to the real interval (0,2) . The second line takes generated values in (1,2) and maps them back to (0,1).
Naively, I tried two variants of this: (B) – Instead of subtracting 1.0, multiply only those values that are > 1 by 0.5, and (C) – multiply by 0.5, regardless of the RF value produced by the first line. Variants B and C failed 10 of the 15 TestU01 SmallCrush tests, whereas the original version failed just 1 of the 15 tests.
Agreed, but I am a bit wary of using MODULO with real type arguments. In this case, however, the second argument has an exact internal representation, so your suggestion should be fine.
OP’s post contained F77 code, so the MODULO function was not available, but MOD() was.
