gfortran how do I increment random_seed by 2^128

Question

The gfortran page on random_seed says that when using OMP threads, each thread increments its seed by 2^128. I am wondering how I increment the seed by 2^128 manually. I wrote a little test program to set the master seed at all 0, and then see what the seeds were, but I don't understand what I'm seeing. What I'd like to know is for example what I put in the subroutine increment_by_2_tothe_128

program main
  implicit none

  character(len=32)    :: arg
  integer              :: n
  integer              :: i
  integer              :: nthreads
  integer, allocatable :: seed(:, :)
  integer, allocatable :: master_seed(:)
  real,    allocatable :: rn(:)

  call get_command_argument(1, arg)
  read(arg, *) nthreads

  call random_seed(size=n)
  allocate(seed(n, nthreads))
  allocate(master_seed(n))
  allocate(rn(nthreads))
  master_seed = 0
  seed = 0
  call random_seed(put=master_seed)

  ! call increment_by_2_tothe_128(n)

  call omp_set_num_threads(nthreads)

  !$OMP PARALLEL DO
  do i=1,nthreads
     call random_number(rn(i))
     call random_seed(get=seed(:,i))
  end do

  do i=1,nthreads
     print *, i
     print *, rn(i)
     print *, seed(:,i)
  end do

end program main


subroutine increment_by_2_tothe_128(n)
  implicit none
  integer, intent(in) :: n
  integer             :: current_seed(n)
  integer             :: increment_seed(n)
  call random_seed(get=current_seed)
  ! what goes here:
  ! incrememt_seed = current_seed + 2**128
  call random_seed(put=increment_seed)
end subroutine increment_by_2_tothe_128

Answer 1

You cannot do that manually. You need the access to the random number generator to be able to do that, but the internals are not exposed to Fortran programmers. And you obviously cannot call the generator 2^128 times.

If you need to do the shift, you need to use some pseud-random number generator that does expose the internals and at the same time allows this kind of shift. That can be, for example, the xoroshiro PRNG family that is used internally by gfortran. These generators have a specialized function for this shift:

All generators, being based on linear recurrences, provide jump functions that make it possible to simulate any number of calls to the next-state function in constant time, once a suitable jump polynomial has been computed. We provide ready-made jump functions for a number of calls equal to the square root of the period, to make it easy generating non-overlapping sequences for parallel computations, and equal to the cube of the fourth root of the period, to make it possible to generate independent sequences on different parallel processors.

These generators are most often implemented in C, but Fortran implementations also exist (subroutine rng_jump is the jump function, disclaimer: the link goes to my repository, no guarantees for the quality).