在C和Fortran之间产生相同的随机数序列(gcc 10.3.0)

Question

我想在用C和Fortran编写的程序中产生相同的随机数序列。我使用的是运行Ubuntu20.04LTS的Windows Subsystem for Linux (Win10)中的gcc 10.3.0版。我已经尝试在C中使用相同的RNG实现，并使用相同的种子，但到目前为止，我似乎无法获得相同的输出。

根据RANDOM_NUMBER()的文档，运行时实现xoshiro256**

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

我从这里抓取了xoshiro256**：

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

这是我的C测试程序(使用高53位来计算64位结果)：

#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);
        }
    }
}

使用以下代码行编译：

gcc-10 -o test_rng_c test_rng.c

下面是C语言的输出：

0.58728832572904277053993382651242427527904510498046875000
0.99628180739434757384742624708451330661773681640625000000
0.91328887972667560646300444204825907945632934570312500000
0.47383268683575230362237107328837737441062927246093750000
0.28868422781068314719732370576821267604827880859375000000
0.90562880102745912935802152787800878286361694335937500000
0.16771219329385134155785408438532613217830657958984375000
0.27161266560374741629857453517615795135498046875000000000
0.78859890296564516543043055207817815244197845458984375000
0.95109314522241128475599225566838867962360382080078125000

这是我的Fortran测试程序：

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

使用以下代码行编译：

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

下面是Fortran的输出：

0.33898027596283408779953560951980762183666229248046875000
0.30153436367708152943123423028737306594848632812500000000
0.85599856867939150273372206356725655496120452880859375000
0.39833679489551077068654194590635597705841064453125000000
0.06316340952695409516337576860678382217884063720703125000
0.66356016219458069382852727358113043010234832763671875000
0.68412746418987169239045442736824043095111846923828125000
0.98553591281548125202505161723820492625236511230468750000
0.27351003115908101293030085798818618059158325195312500000
0.09340321315109112454422302107559517025947570800781250000

gfortran-10 -v的输出。

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)

我要么没有在C代码中正确设置种子，要么在一些关键方面稍微延迟了实现。任何洞察力都将不胜感激。

Answer 1

为了防止将种子设置为全零的常见错误，GFortran实现将提供的种子与“秘密”密钥进行异或。" secret“在引号中，因为它在任何密码意义上都不是真正的秘密，你可以在实现的源代码中找到密钥。

此外，如果您使用线程，GFortran会做一些技巧来跳过PRNG序列，以便每个线程都有自己的PRNG流。

因此，对于您的C实现，您需要复制这些特性。

可能最简单和最健壮的解决方案是为C实现创建一个ISO_C_BINDING接口，并从Fortran中调用该接口。

Answer 2

多亏了janneb发布的洞察力，我决定将其带到源代码中，并在这里找到了答案：gcc/libgfortran/intrinsics/random.c

xoshiro256**实现与问题中发布的源文件相同，并且秘密异或密钥位于名为xor_keys的常量中，并使用名为scramble_seed的方法应用。此外，当从int32用户数组读取种子数组时，种子数组被颠倒。

对xor打乱的解释直接来自源代码：

  /* 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);

下面是更新后的C示例，它产生的结果与上面的Fortran代码完全相同。从uint64也从源代码转换为浮点型：

#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);
    }
}