What does this code mean?

Help
25

Why, would you expect a different behaviour between mod and modulo for positive numbers? I’m missing something?

I will also suggest not to use float as its result will be single precision instead of double, even though the difference will be small.
As other modernisation I wont use data statements as it is hiding which variable will be updated and which one is a parameter. I will write something like:

module rfb_mod
    implicit none
    integer, parameter :: dp = kind(1.0d0)

contains
    real(dp) function rfb()
        integer, parameter :: ix1 = 1500419
        integer, parameter :: ix2 = 1400159
        integer, parameter :: ix3 = 1364
        integer, parameter :: ix4 = 1528

        integer :: ix5 = 1
        integer :: ix6 = 3

        real(dp) :: rr1, rr2 

        rr1 = 1.0/real(ix1, dp)
        rr2 = 1.0/real(ix2, dp)

        ix5 = MODULO(ix5*ix3, ix1)
        ix6 = MODULO(ix6*ix4, ix2)
        rfb = rr1*ix5 + rr2*ix6

        rfb = modulo(rfb, 1.0_dp)
        ! or 
        ! if (rfb >= 1.0) rfb = rfb -1.0
        ! which is the fastest
    end function
end module

by the way the OP is also saying that he started learning Fortran so it may be better to show him a more modern syntax.
Of course there are other possible improvements.

26

In a simple try

if (rfb >= 1.0) rfb = rfb -1.0

is faster than

rfb = modulo(rfb, 1.0_dp)

At least with gfortran

28

off-topic but, why did you macro two and not the other two?

29

I am not familiar with the history of this random number generator, but am surprised at it’s good behaviour.
The addition of two random sequences can cause problems.
RF = RR1*IX5 + RR2*IX6
This combines two random number sequences; IX5 and IX6, so I thought the addition could be a triangular distribution, but the “MODULO” approach appears to retain the uniform distribution.
Omitting this “modulo” test does produce a triangular distribution in [0:2]

The use of small initial values for IX5 and IX6 could imply poor initial values, which is reported in the testing in my code below for small “i”. The histogram shows good uniform distribution for large “i”.
I prefer the use of DBLE instead of FLOAT, and parameters could help.

For Monte-Carlo analysis, the use of your own random number generator does have an advantage, as reproducible results can be more easily controlled.

  module rf_mod
    contains
      DOUBLE PRECISION FUNCTION RF ()
        implicit none
!        IMPLICIT DOUBLE PRECISION (A-H, O-Z)
        integer, parameter :: IX1 = 1500419, IX3 = 1364
        integer, parameter :: IX2 = 1400159, IX4 = 1528
        real*8,  parameter :: RR1 = 1.0 / dble (IX1)
        real*8,  parameter :: RR2 = 1.0 / dble (IX2)
        integer :: IX5 = 1, IX6 = 3
!
        IX5 = MOD ( IX5*IX3, IX1 )
        IX6 = MOD ( IX6*IX4, IX2 )
        RF  = RR1*IX5 + RR2*IX6
        IF (RF >= 1.0) RF = RF - 1.0
      END FUNCTION RF
  end module rf_mod

  program main
    use rf_mod, only: rf
    implicit none
    integer, parameter :: n = 10**8, dp = kind(1.0d0)
    integer :: i, k, hist(0:10), ih
    real(kind=dp) :: xi, s, ss  ! , x(n)
    ! real(kind=dp) is another way of saying double precision
      hist = 0
      k = 4
      s = 0
      ss = 0
      do i=1,n
         xi = rf()
         s  = s + xi
         ss = ss + (xi-0.5)**2
         ih = xi*10.
         hist(ih) = hist(ih)+1
         if ( i < k ) cycle
         print*,"actual",i, s/i, ss / i
         k = k*2
      end do
      print*,"actual",n, s/n, ss / n
    !  print*,"actual",n,sum(x)/n,sum((x-0.5)**2)/n
      print*,"theory",0, 0.5_dp,1/12.0_dp
      do i = 0,10
        write (*,*) i,hist(i)
      end do
  end program main
30

Another MC RNG with similar properties is the Wichmann-Hill RNG. This RNG combines three Lehman RNGs, each of which uses constants that are five or fewer decimal digits in size. The period of the combined RNG has been established to be close to 7 trillion.

Your observation regarding starting the generator with small seeds is correct, but a common remedy is to “warm up” the RNG with sufficient iterations to push the zeros out of the big end(s) before using the output of the RNG.

The Wichmann-Hill RNG also passed 14 out of the 15 tests in the SmallCrush TestU01 battery.