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)
I think that you are missing an operator in IX5=MOD(IX5IX3,IX1) and IX6=MOD(IX6IX4,IX2). The first argument in the MOD function likely needs a multiplication, e.g., IX6=MOD(IX6*IX4,IX2).
For a random number generator, the important statement here is the data statement. The first time the routine is called, IX1 is set to 1500419, IX2 is set to 1400159, …, IX6 is set to 3. The data statement implies the SAVE attribute, so when you call the function again instead of re-initializing the variables with the values from the data statement, the function uses whatever the values were when the returned last time.
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.
DOUBLE PRECISION FUNCTION RF()
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
end module rf_mod
use rf_mod, only: rf
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
x(i) = rf()
end program main
theory 0.50000000000000000 8.3333333333333329E-002
actual 0.50019079377094988 8.3373676712904451E-002
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.
Yes, random_number should be used in production code. Or, do you also re-implement all other standard intrinsic subprograms. After all, the implementation details of sin, cos, exp, log, etc can vary between compilers and math libraries.
I realize that it’s a foreign concept to many (here?). Reading the documentation that comes with a compiler and testing the intrinsic subprogram is always a good thing.
Would you simply use your custom PRNG without testing as it must be better that random_number? If no, and you run a few stringent tests of your custom PRNG, why not run the same tests against random_number? When you update your OS or compiler or use a different compiler, do you re-run your stringent tests?
It is in the best interest of compiler vendors to offer a robust random_number, which is not too hard to do as there are several available.
As to the quality of seeds, the entropy issue applies equally to vendor-supplied random_seed and random_init as to your custom PRNG.
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.
#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;
ix5 = (ix5*ix3) % IX1;
ix6 = (ix6*ix4) % IX2;
x = rr1*ix5+rr2*ix6;
return x < 1 ? x : x-1.0;
int main (void)
gen = unif01_CreateExternGen01 ("FLD12", rf);
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.