Produce Identical Random Number Sequence between C and Fortran (gcc 10.3.0)

Question

I would like to produce the same random number sequence in both a program written in C and Fortran. I am using gcc version 10.3.0 within Windows Subsystem for Linux (Win10) running Ubuntu 20.04 LTS. I've tried using the same RNG implementation in C and using the same seeds, but thus far I can't seem to get the same output.

According to the documentation for RANDOM_NUMBER(), the runtime implements xoshiro256**

https://gcc.gnu.org/onlinedocs/gcc-10.3.0/gfortran/RANDOM_005fNUMBER.html

I grabbed xoshiro256** from here:

https://prng.di.unimi.it/xoshiro256starstar.c

This is my C test program (uses the upper 53 bits to calculate a 64-bit result):

#include <stdio.h>
#include "xoshiro256starstar.c"

int main(int argc, char** argv)
{
    size_t trials = 10u;
    int* p32state = (int*)s;

    p32state[0] = -1468754105;
    p32state[1] = -743753204;
    p32state[2] =  2022458965;
    p32state[3] = -443226243;
    p32state[4] = -23942267;
    p32state[5] = -1265286489;
    p32state[6] = -1934963269;
    p32state[7] =  1953963768;

    for( size_t i = 0u; i < trials; ++i )
    {
        uint64_t ret = next();

        if( i < 10u )
        {
            double d = ((uint64_t)(ret >> 11ull)) / (double)(1ull << 53ull);
            printf("%1.56f\n", d);
        }
    }
}

Compiled using this line:

gcc-10 -o test_rng_c test_rng.c

And here's the C output:

0.58728832572904277053993382651242427527904510498046875000
0.99628180739434757384742624708451330661773681640625000000
0.91328887972667560646300444204825907945632934570312500000
0.47383268683575230362237107328837737441062927246093750000
0.28868422781068314719732370576821267604827880859375000000
0.90562880102745912935802152787800878286361694335937500000
0.16771219329385134155785408438532613217830657958984375000
0.27161266560374741629857453517615795135498046875000000000
0.78859890296564516543043055207817815244197845458984375000
0.95109314522241128475599225566838867962360382080078125000

This is my Fortran test program:

program test_rand   

   implicit none

   real*8::v
   integer::trials,i
   integer::n
   integer, allocatable :: seed(:)
   trials = 10

   call random_seed(size=n)
   allocate(seed(n))
   if (n .eq. 8) then
     seed(1) = -1468754105
     seed(2) = -743753204
     seed(3) =  2022458965
     seed(4) = -443226243
     seed(5) = -23942267
     seed(6) = -1265286489
     seed(7) = -1934963269
     seed(8) =  1953963768
     call random_seed(put=seed)
   endif

   call random_seed(get=seed)
   do i=0,trials-1
      call random_number(v)
      Print "(f58.56)", v
   enddo

stop
end

Compiled using this line:

gfortran-10 -ffree-form -o test_rng_f test_rand.f

And here's the Fortran output:

0.33898027596283408779953560951980762183666229248046875000
0.30153436367708152943123423028737306594848632812500000000
0.85599856867939150273372206356725655496120452880859375000
0.39833679489551077068654194590635597705841064453125000000
0.06316340952695409516337576860678382217884063720703125000
0.66356016219458069382852727358113043010234832763671875000
0.68412746418987169239045442736824043095111846923828125000
0.98553591281548125202505161723820492625236511230468750000
0.27351003115908101293030085798818618059158325195312500000
0.09340321315109112454422302107559517025947570800781250000

Output of gfortran-10 -v for good measure.

Using built-in specs.
COLLECT_GCC=gfortran-10
COLLECT_LTO_WRAPPER=/usr/lib/gcc/x86_64-linux-gnu/10/lto-wrapper
OFFLOAD_TARGET_NAMES=nvptx-none:amdgcn-amdhsa:hsa
OFFLOAD_TARGET_DEFAULT=1
Target: x86_64-linux-gnu
Configured with: ../src/configure -v --with-pkgversion='Ubuntu 10.3.0-1ubuntu1~20.04' --with-bugurl=file:///usr/share/doc/gcc-10/README.Bugs --enable-languages=c,ada,c++,go,brig,d,fortran,objc,obj-c++,m2 --prefix=/usr --with-gcc-major-version-only --program-suffix=-10 --program-prefix=x86_64-linux-gnu- --enable-shared --enable-linker-build-id --libexecdir=/usr/lib --without-included-gettext --enable-threads=posix --libdir=/usr/lib --enable-nls --enable-bootstrap --enable-clocale=gnu --enable-libstdcxx-debug --enable-libstdcxx-time=yes --with-default-libstdcxx-abi=new --enable-gnu-unique-object --disable-vtable-verify --enable-plugin --enable-default-pie --with-system-zlib --enable-libphobos-checking=release --with-target-system-zlib=auto --enable-objc-gc=auto --enable-multiarch --disable-werror --with-arch-32=i686 --with-abi=m64 --with-multilib-list=m32,m64,mx32 --enable-multilib --with-tune=generic --enable-offload-targets=nvptx-none=/build/gcc-10-S4I5Pr/gcc-10-10.3.0/debian/tmp-nvptx/usr,amdgcn-amdhsa=/build/gcc-10-S4I5Pr/gcc-10-10.3.0/debian/tmp-gcn/usr,hsa --without-cuda-driver --enable-checking=release --build=x86_64-linux-gnu --host=x86_64-linux-gnu --target=x86_64-linux-gnu --with-build-config=bootstrap-lto-lean --enable-link-mutex
Thread model: posix
Supported LTO compression algorithms: zlib zstd
gcc version 10.3.0 (Ubuntu 10.3.0-1ubuntu1~20.04)

I'm either not setting the seeds correct in the C code, or the implementations defer slightly in some key way. Any insight would be appreciated.

Answer 1

In order to guard against the somewhat common error of setting the seed to all zeros, the GFortran implementation XOR's the provided seed with a "secret" key. "Secret" in quotes because it's not really secret in any kind of cryptographic sense, you can find the key right there in the source code of the implementation.

Also if you're using threads, the GFortran does some tricks to jump the PRNG sequence so that each thread gets its own PRNG stream.

So for your C implementation you'd need to replicate these features.

Probably the simplest and most robust solution is to create an ISO_C_BINDING interface to your C implementation and call that one from Fortran.

Answer 2

Thanks to the insight posted by janneb, I decided to take it to the source code and found the answer here: gcc/libgfortran/intrinsics/random.c

The xoshiro256** implementation is identical to the source file posted in the question, and the secret XOR keys are in const called xor_keys and applied with a method called scramble_seed. In addition, the seed array is reversed as it is read from the int32 user array.

The explanation of the xor scramble directly from the source code:

  /* We put it after scrambling the bytes, to paper around users who
 provide seeds with quality only in the lower or upper part.  */
  scramble_seed (master_state.s, seed);

Here is the updated C example that produces identical results to the above Fortran code. Conversion to float from uint64 also from source:

#include <stdio.h>
#include "xoshiro256starstar.c"

const uint64_t xor_keys[] = {
  0xbd0c5b6e50c2df49ULL, 0xd46061cd46e1df38ULL, 0xbb4f4d4ed6103544ULL, 0x114a583d0756ad39ULL
};

void scramble_seed(uint64_t *dest, const uint64_t *src) {
    for (size_t i = 0u; i < 4u; i++)
    {
        dest[i] = src[i] ^ xor_keys[i];
    }
}

void reverse_array(int *a, size_t n) {
    if( n < 2 ) return;
    for(size_t l = 0u, r = n-1u; l<r; ++l,--r)
    {
        int temp = a[l];
        a[l] = a[r];
        a[r] = temp;
    }
}

int main(int argc, char** argv) {
    size_t trials = 10u;
    int* p32state = (int*)s;

    p32state[0] = -1468754105;
    p32state[1] = -743753204;
    p32state[2] =  2022458965;
    p32state[3] = -443226243;
    p32state[4] = -23942267;
    p32state[5] = -1265286489;
    p32state[6] = -1934963269;
    p32state[7] =  1953963768;

    reverse_array(p32state,8u);
    scramble_seed(s,s);

    const uint32_t mask = ~ (uint32_t) 0u << 8u;
    for( size_t i = 0u; i < trials; ++i )
    {
        uint64_t ret = (uint32_t)(next() >> 32u) & mask;
        return (float)ret * 0x1.p-32f;
        printf("%1.56f\n", d);
    }
}