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

Question

我想在用 C 和 Fortran 编写的程序中生成相同的随机数序列。我在运行 Ubuntu 20.04 LTS 的 Windows 子系统 (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 实现 XOR 使用“秘密”密钥提供的种子。引号中的“秘密”，因为它在任何加密意义上都不是真正的秘密，您可以在实现的源代码中找到密钥。

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

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

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

Answer 2

感谢janneb 发表的见解，我决定将其带到源代码并在这里找到答案：gcc/libgfortran/intrinsics/random.c

xoshiro256** 实现与问题中发布的源文件相同，秘密 XOR 密钥位于名为 xor_keys 的 const 中，并使用名为 scramble_seed 的方法应用。此外，种子数组在从 int32 用户数组中读取时会反转。

xor scramble的解释直接出自源码：

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