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

Question

我想在用C和Fortran编写的程序中生成相同的随机数序列。我在运行Ubuntu 20.04 LTS的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.56fn", 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实现将提供的种子与"；秘密"；键"秘密"在引号中，因为它在任何加密意义上都不是真正的秘密，所以您可以在实现的源代码中找到密钥。

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

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

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

Answer 2

多亏了janneb发布的见解，我决定把它带到源代码中，并在这里找到了答案：gcc/libgfortran/intrinsics/random.c

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

直接从源代码解释异或加扰：

/* 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.56fn", d);
}
}