RANDOM_INIT revisted.

RANDOM_INIT revisted.

[Steve Kargl]

Dear all, 

J3 defined yet another useless intrinsic subprogram with RANDOM_INIT().

Damian recently brought my attention back to this routine, and I looked
at how I implemented.  Consider the following program

   program foo
      implicit none

      integer i
      real r(5)

      call random_init(repeatable=.true., image_distinct=.false.)
      call random_number(r)
      print '(5F8.5)', r

      ! Burn some random numbers.
      do i = 1, 10
        call random_number(r)
      end do

      ! This reseeds with original seeds.
      call random_init(repeatable=.true., image_distinct=.false.)
      call random_number(r)
      print '(5F8.5)', r
   end program foo

Currently, gfortran interprets 'repeatable' to apply to the
current execution of the current image.  That is, the second
call to RANDOM_INIT() above re-initializes the PRNG with the
same set of initial seed values.

% gfortran9 -o z a.f90 
% ./z
 0.40558 0.89909 0.76782 0.10472 0.43036
 0.40558 0.89909 0.76782 0.10472 0.43036

When the executable is run again, a new set of seeds are used.

% ./z
 0.10866 0.36828 0.64442 0.00358 0.10796
 0.10866 0.36828 0.64442 0.00358 0.10796

The problem is that the language in the Standard suggests that
the same set of seeds need to be used whenever 'z' is executed
(ie., today, tomorrow, ten years from now).  Now, note that 
'repeatable' is not a compile constants, so the executable must be
prepared to have it set to .true. at run-time.  This means that
an initial set of seeds needs to be compiled into 'z'.

There are two options.  The first option would have gfortran produce
a unique processor-dependent hardcoded set of seeds everytime the
program a.f90 is compiled.  The problem with this approach is that
one could no longer do reproducible builds.  If the unique set of 
seeds changes, then of course, the resulting executable will be
different from previous builds.  The second option is to use a 
hardcoded set of seeds that do not change during compilation.  The
downside to this approach is that every executable built by gfortran
will use the exact same set of processor-dependent seeds.  With the
following patch, gfortran now generates an executable such that

% gfcx -o z a.f90
% ./z
 0.82526 0.19133 0.15550 0.98538 0.07473
 0.82526 0.19133 0.15550 0.98538 0.07473
% ./z
 0.82526 0.19133 0.15550 0.98538 0.07473
 0.82526 0.19133 0.15550 0.98538 0.07473

I can leave this executable on my system and come back next year
and it will give me the same sequence of PRN.

The patch uses an embedded simply linear congruential generator (LCG)
with a single default seed to populate the seeds for the PRNG.  If
multiple images are instantiate and image_distinct=.true., then unique
sets of seeds are drawn from the LCG to seed the PRNG.


Index: libgfortran/intrinsics/random_init.f90
===================================================================
--- libgfortran/intrinsics/random_init.f90	(revision 280157)
+++ libgfortran/intrinsics/random_init.f90	(working copy)
@@ -70,19 +70,30 @@ impure subroutine _gfortran_random_init(repeatable, im
    integer, value, intent(in) :: hidden
 
    logical, save :: once = .true.
-   integer :: nseed
+   integer :: nseed, lcg_seed
    integer, save, allocatable :: seed(:)
 
    if (once) then
       once = .false.
       call random_seed(size=nseed)
       allocate(seed(nseed))
-      call random_seed(get=seed)
       !
-      ! To guarantee that seed is distinct on multiple images, add the hidden
-      ! argument (which is the image index).
+      ! To guarantee that seed is distinct on multiple images, use the
+      ! hidden argument (which is the image index) to change the lcg_seed.
+      ! For non-coarray programs and programs compiled with -fcoarray=single,
+      ! hidden is zero, so the do-loop would not execute. Add 1.
       !
-      if (image_distinct) seed = seed + hidden
+      block
+         integer i
+         lcg_seed = 57911963
+         if (image_distinct ) then
+            do i = 1, hidden * 10 + 1
+               call _gfortran_lcg(seed)
+            end do
+         else
+            call _gfortran_lcg(seed)   
+         end if
+      end block 
    end if
 
    if (repeatable) then
@@ -90,5 +101,24 @@ impure subroutine _gfortran_random_init(repeatable, im
    else
       call random_seed();
    end if
+   contains
+      !
+      ! SK Park and KW Miller, ``Random number generators: good ones are hard
+      ! to find,'' Comm. ACM, 31(10), 1192--1201, (1988).
+      !
+      ! Implementation of a prime modulus multiplicative linear congruential
+      ! generator, which avoids overflow and provides the full period.
+      !
+      impure elemental subroutine _gfortran_lcg(i)
+         implicit none
+         integer, intent(out) :: i
+         integer, parameter :: a  = 16807    ! Multiplier
+         integer, parameter :: m = huge(a)   ! Modulus
+         integer, parameter :: q = 127773    ! Quotient to avoid overflow
+         integer, parameter :: r = 2836      ! Remainder to avoid overflow
+         lcg_seed = a * mod(lcg_seed, q) - r * (lcg_seed / q)
+         if (lcg_seed <= 0) lcg_seed = lcg_seed + m
+         i = lcg_seed
+      end subroutine _gfortran_lcg
 
 end subroutine _gfortran_random_init


-- 
Steve