What does this code mean?

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
1 Like

Welcome.

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.

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
2 Likes

The original message did not have code formatted as code, so * was suppressed. Now that the code is formatted it is compilable.

1 Like

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!

Tell your advisor that Fortran has had a built-in random number generator for 30 years. :slight_smile:

call random_number( r )
5 Likes

Do not tell your advisor, unless she has no connection to George Masaglia, and welcome from me too.

2 Likes

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()?

The RNG in the Fortran code combines two multiplicative congruential random number generators, presumably to increase the period of the RNG.

Here, for example, is the equivalent C code.

#include <stdio.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;
printf("%06X %06X  ",ix5,ix6);
x = rr1*ix5+rr2*ix6;
return x < 1 ? x : x-1.0;
}

main(){
int i;
for(i=0; i<20; i++)printf("%2d  %9.6f\n",i,rf());
}

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…

4 Likes

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.

2 Likes

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.

1 Like

Good point, but we have the stdlib statistical uniform distribution module and, the more basic, Pseudorandom Number Generator Module.

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

1 Like

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.

1 Like

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.

2 Likes

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

1 Like

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.

The previous internal implementation of gfortran’s random_number was also very good.

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.

[Professor Pierre L’Ecuyer is famed as an expert on RNGs.]

2 Likes