问题 Xorshift1024 *跳转不可交换?


我一直在移动塞巴斯蒂亚诺维尼亚 xorshift1024 * PRNG与标准的C ++ 11统一随机数生成器契约兼容,并注意到了一些奇怪的行为 jump() 他提供的功能。

根据维尼亚的说法 jump() 应该相当于2 ^ 512次调用 next()。因此一系列来电 jump() 和 next() 应该是可交换的。例如,假设发电机以某种已知状态启动,

jump();
next();

应该让发电机处于相同的状态

next();
jump();

因为两者都应该相当于

for (bigint i = 0; i < (bigint(1) << 512) + 1; ++i)
    next();

假设 bigint 是一些具有极大最大值的整数类型(假设你是一个非常非常非常有耐心的人)。

不幸的是,这不适用于Vigna提供的参考实现(我将在最后包括后代;如果上面链接的实现发生变化或将来被删除)。使用以下测试代码测试前两个选项时:

memset(s, 0xFF, sizeof(s));
p = 0;

// jump() and/or next() calls...

std::cout << p << ';';
for (int i = 0; i < 16; ++i)
    std::cout << ' ' << s[i];

调用 jump() 之前 next() 输出:

1; 9726214034378009495 13187905351877324975 10033047168458208082 990371716258730972 965585206446988056 74622805968655940 11468976784638207029 3005795712504439672 6792676950637600526 9275830639065898170 6762742930827334073 16862800599087838815 13481924545051381634 16436948992084179560 6906520316916502096 12790717607058950780

一边打电话 next() 第一个结果:

1; 13187905351877324975 10033047168458208082 990371716258730972 965585206446988056 74622805968655940 11468976784638207029 3005795712504439672 6792676950637600526 9275830639065898170 6762742930827334073 16862800599087838815 13481924545051381634 16436948992084179560 6906520316916502096 12790717607058950780 9726214034378009495

显然要么我理解什么 jump() 做错了,或者有一个错误 jump() 函数,或跳转多项式数据是错误的。 Vigna声称这个跳跃函数可以计算出该期间的任何步幅,但没有详细说明如何计算它(包括在他的 关于xorshift *发电机的论文)。如何计算正确的跳转数据以验证其中某处没有拼写错误?


Xorshift1024 *参考实现; http://xorshift.di.unimi.it/xorshift1024star.c

/*  Written in 2014-2015 by Sebastiano Vigna (vigna@acm.org)

To the extent possible under law, the author has dedicated all copyright
and related and neighboring rights to this software to the public domain
worldwide. This software is distributed without any warranty.

See <http://creativecommons.org/publicdomain/zero/1.0/>. */

#include <stdint.h>
#include <string.h>

/* This is a fast, top-quality generator. If 1024 bits of state are too
   much, try a xorshift128+ generator.

   The state must be seeded so that it is not everywhere zero. If you have
   a 64-bit seed, we suggest to seed a splitmix64 generator and use its
   output to fill s. */

uint64_t s[16]; 
int p;

uint64_t next(void) {
    const uint64_t s0 = s[p];
    uint64_t s1 = s[p = (p + 1) & 15];
    s1 ^= s1 << 31; // a
    s[p] = s1 ^ s0 ^ (s1 >> 11) ^ (s0 >> 30); // b,c
    return s[p] * UINT64_C(1181783497276652981);
}


/* This is the jump function for the generator. It is equivalent
   to 2^512 calls to next(); it can be used to generate 2^512
   non-overlapping subsequences for parallel computations. */

void jump() {
    static const uint64_t JUMP[] = { 0x84242f96eca9c41dULL,
        0xa3c65b8776f96855ULL, 0x5b34a39f070b5837ULL, 0x4489affce4f31a1eULL,
        0x2ffeeb0a48316f40ULL, 0xdc2d9891fe68c022ULL, 0x3659132bb12fea70ULL,
        0xaac17d8efa43cab8ULL, 0xc4cb815590989b13ULL, 0x5ee975283d71c93bULL,
        0x691548c86c1bd540ULL, 0x7910c41d10a1e6a5ULL, 0x0b5fc64563b3e2a8ULL,
        0x047f7684e9fc949dULL, 0xb99181f2d8f685caULL, 0x284600e3f30e38c3ULL
    };

    uint64_t t[16] = { 0 };
    for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
        for(int b = 0; b < 64; b++) {
            if (JUMP[i] & 1ULL << b)
                for(int j = 0; j < 16; j++)
                    t[j] ^= s[(j + p) & 15];
            next();
        }

    memcpy(s, t, sizeof t);
}

9691
2018-01-03 08:32


起源

什么是 call()? - deviantfan
@deviantfan一个错字:)编辑。 - bcrist


答案:


好的,我很抱歉,但有时会发生这种情况(我是作者)。

最初该函数有两个memcpy()。然后我意识到需要一份圆形副本。但我只更换了 第一 的memcpy()。愚蠢,愚蠢,愚蠢。网站上的所有文件都已修复。 arXiv副本正在进行更新。看到 http://xorshift.di.unimi.it/xorshift1024star.c 

顺便说一下:我没有“发布”科学意义上的任何错误,因为jump()函数不是ACM Trans的一部分。数学。柔软的。几周前,它已经在网站和arXiv / WWW版本上添加了。网络和arXiv的快速发布路径意味着,有时,一个人分发未经抛光的纸张。我只能感谢记者报告这个错误(好的,技术上StackOverflow不报告错误,但我也收到了一封电子邮件)。

不幸的是,单元测试我没有考虑情况p≠0。我主要担心的是计算多项式的正确性。如上所述,当p = 0时,该函数是正确的。

至于计算:每个发生器对应一个特征多项式P(x)。 k的跳跃多项式只是x ^ k mod P(x)。我使用fermat来计算这些功能,然后我有一些脚本生成C代码。

当然我不能测试2 ^ 512,但由于我的代码代码完美地工作在2到2 ^ 30(你可以轻松测试的范围),我相信它也可以在2 ^ 512工作。这只是fermat计算x ^ {2 ^ 512}而不是x ^ {2 ^ 30}。但是,独立验证非常受欢迎。

我的代码只适用于x ^ {2 ^ t}形式的权限。这是我需要计算有用的跳转函数。计算多项式模P(x)并不困难,因此可以想象任何值都可以具有完全通用的跳转函数,但坦率地说,我发现这完全是过度的。

如果有人对获取其他跳转多项式感兴趣,我可以提供脚本。它们将成为下一个xorshift发行版的所有其他代码的一部分,但我需要在发布之前完成文档。

为了记录,xorshift1024 *的特征多项式是x ^ 1024 + x ^ 974 + x ^ 973 + x ^ 972 + x ^ 971 + x ^ 966 + x ^ 965 + x ^ 964 + x ^ 963 + x ^ 960 + x ^ 958 + x ^ 957 + x ^ 956 + x ^ 955 + x ^ 950 + x ^ 949 + x ^ 948 + x ^ 947 + x ^ 942 + x ^ 941 + x ^ 940 + x ^ 939 + x ^ 934 + x ^ 933 + x ^ 932 + x ^ 931 + x ^ 926 + x ^ 925 + x ^ 923 + x ^ 922 + x ^ 920 + x ^ 917 + x ^ 916 + x ^ 915 + x ^ 908 + x ^ 906 + x ^ 904 + x ^ 902 + x ^ 890 + x ^ 886 + x ^ 873 + x ^ 870 + x ^ 857 + x ^ 856 + x ^ 846 + x ^ 845 + x ^ 844 + x ^ 843 + x ^ 841 + x ^ 840 + x ^ 837 + x ^ 835 + x ^ 830 + x ^ 828 + x ^ 825 + x ^ 824 + x ^ 820 + x ^ 816 + x ^ 814 + x ^ 813 + x ^ 811 + x ^ 810 + x ^ 803 + x ^ 798 + x ^ 797 + x ^ 790 + x ^ 788 + x ^ 787 + x ^ 786 + x ^ 783 + x ^ 774 + x ^ 772 + x ^ 771 + x ^ 770 + x ^ 769 + x ^ 768 + x ^ 767 + x ^ 765 + x ^ 760 + x ^ 758 + x ^ 753 + x ^ 749 + x ^ 747 + x ^ 746 + x ^ 743 + x ^ 741 + x ^ 740 + x ^ 738 + x ^ 737 + x ^ 736 + x ^ 735 + x ^ 728 + x ^ 726 + x ^ 723 + x ^ 722 + x ^ 721 + x ^ 720 + x ^ 718 + x ^ 716 + x ^ 715 + x ^ 714 + x ^ 710 + x ^ 709 + x ^ 707 + x ^ 694 + x ^ 687 + x ^ 686 + x ^ 685 + x ^ 684 + x ^ 679 + x ^ 678 + x ^ 677 + x ^ 674 + x ^ 670 + x ^ 669 + x ^ 667 + x ^ 666 + x ^ 665 + x ^ 663 + x ^ 658 + x ^ 655 + x ^ 651 + x ^ 639 + x ^ 638 + x ^ 635 + x ^ 634 + x ^ 632 + x ^ 630 + x ^ 623 + x ^ 621 + x ^ 618 + x ^ 617 + x ^ 616 + x ^ 615 + x ^ 614 + x ^ 613 + x ^ 609 + x ^ 606 + x ^ 604 + x ^ 601 + x ^ 600 + x ^ 598 + x ^ 597 + x ^ 596 + x ^ 594 + x ^ 593 + x ^ 592 + x ^ 590 + x ^ 589 + x ^ 588 + x ^ 584 + x ^ 583 + x ^ 582 + x ^ 581 + x ^ 579 + x ^ 577 + x ^ 575 + x ^ 573 + x ^ 572 + x ^ 571 + x ^ 569 + x ^ 567 + x ^ 565 + x ^ 564 + x ^ 563 + x ^ 561 + x ^ 559 + x ^ 557 + x ^ 556 + x ^ 553 + x ^ 552 + x ^ 550 + x ^ 544 + x ^ 543 + x ^ 542 + x ^ 541 + x ^ 537 + x ^ 534 + x ^ 532 + x ^ 530 + x ^ 528 + x ^ 526 + x ^ 523 + x ^ 521 + x ^ 520 + x ^ 518 + x ^ 516 + x ^ 515 + x ^ 512 + x ^ 511 + x ^ 510 + x ^ 508 + x ^ 507 + x ^ 506 + x ^ 505 + x ^ 504 + x ^ 502 + x ^ 501 + x ^ 499 + x ^ 497 + x ^ 494 + x ^ 493 + x ^ 492 + x ^ 491 + x ^ 490 + x ^ 487 + x ^ 485 + x ^ 483 + x ^ 482 + x ^ 480 + x ^ 479 + x ^ 477 + x ^ 476 + x ^ 475 + x ^ 473 + x ^ 469 + x ^ 468 + x ^ 465 + x ^ 463 + x ^ 461 + x ^ 460 + x ^ 459 + x ^ 458 + x ^ 455 + x ^ 453 + x ^ 451 + x ^ 448 + x ^ 447 + x ^ 446 + x ^ 445 + x ^ 443 + x ^ 438 + x ^ 437 + x ^ 431 + x ^ 430 + x ^ 429 + x ^ 428 + x ^ 423 + x ^ 417 + x ^ 416 + x ^ 415 + x ^ 414 + x ^ 412 + x ^ 410 + x ^ 409 + x ^ 408 + x ^ 400 + x ^ 398 + x ^ 396 + x ^ 395 + x ^ 391 + x ^ 390 + x ^ 386 + x ^ 385 + x ^ 381 + x ^ 380 + x ^ 378 + x ^ 375 + x ^ 373 + x ^ 372 + x ^ 369 + x ^ 368 + x ^ 365 + x ^ 360 + x ^ 358 + x ^ 357 + x ^ 354 + x ^ 350 + x ^ 348 + x ^ 346 + x ^ 345 + x ^ 344 + x ^ 343 + x ^ 342 + x ^ 340 + x ^ 338 + x ^ 337 + x ^ 336 + x ^ 335 + x ^ 333 + x ^ 332 + x ^ 325 + x ^ 323 + x ^ 318 + x ^ 315 + x ^ 313 + x ^ 309 + x ^ 308 + x ^ 305 + x ^ 303 + x ^ 302 + x ^ 300 + x ^ 294 + x ^ 290 + x ^ 281 + x ^ 279 + x ^ 276 + x ^ 275 + x ^ 273 + x ^ 272 + x ^ 267 + x ^ 263 + x ^ 262 + x ^ 261 + x ^ 260 + x ^ 258 + x ^ 257 + x ^ 256 + x ^ 249 + x ^ 248 + x ^ 243 + x ^ 242 + x ^ 240 + x ^ 238 + x ^ 236 + x ^ 233 + x ^ 232 + x ^ 230 + x ^ 228 + x ^ 225 + x ^ 216 + x ^ 214 + x ^ 212 + x ^ 210 + x ^ 208 + x ^ 206 + x ^ 205 + x ^ 200 + x ^ 197 + x ^ 196 + x ^ 184 + x ^ 180 + x ^ 176 + x ^ 175 + x ^ 174 + x ^ 173 + x ^ 168 + x ^ 167 + x ^ 166 + x ^ 157 + x ^ 155 + x ^ 153 + x ^ 152 + x ^ 151 + x ^ 150 + x ^ 144 + x ^ 143 + x ^ 136 + x ^ 135 + x ^ 125 + x ^ 121 + x ^ 111 + x ^ 109 + x ^ 107 + x ^ 105 + x ^ 92 + x ^ 90 + x ^ 79 + x ^ 78 + x ^ 77 + x ^ 76 + x ^ 60 + 1


10
2018-01-03 23:51





tldr:我很确定原始代码中有一个错误:
memcpy 在 jump() 必须考虑 p 轮换也是。
在发表论文之前,作者没有进行适当的测试......

长版:

next() 呼叫只改变16中的一个 s 数组元素,带索引的元素 pp 从0开始,每个都增加 next() 呼叫,15之后再次变为0。我们打电话吧 s[p] “当前”数组元素。实施的另一种(较慢)可能性 next() 将是当前元素始终是第一个,没有 p,而不是递增 p 整体 s 旋转数组(即第一个元素移动到最后一个位置,前一个第二个元素成为第一个元素)。

独立于当前 p 价值,16次来电 next() 应该导致相同的 p 价值与以前一样,即。整个循环完成,当前元素与16个调用之前的位置相同。 jump() 应该做2 ^ 512 next(),2 ^ 512是16的倍数,所以有一个跳跃, p 它之前和之后的值应该是相同的。

你可能已经注意到你的两个不同的结果只旋转了一次,即。一个解决方案是“9726214034378009495 somethingelse”,一个是“somethingelse 9726214034378009495”

......因为你做到了   next() 之前/之后 jump() 和 jump()  无法处理 p 除了0

如果你用16测试它 next() (或32或0或......)之前/之后 jump() 而不是一个,两个结果是相等的。原因是,在跳跃中,而对于 s 数组当前元素/ p 按原样处理 next()t 数组在语义上旋转,以便当前元素始终是第一个(t[j] ^= s[(j + p) & 15];)。然后,在函数终止之前, memcpy(s, t, sizeof t); 从中复制新值 t 回到 s 根本不考虑轮换。只需更换 memcpy 适当的循环包括 p 抵消,那应该没问题。

(好吧,但这并不意味着 jump() 真是一样的 2 ^ 512  next()。但至少它 可以 是。)


4
2018-01-03 09:45



我决定定下来 p = 0 之后 memcpy 相反,因为正如你所提到的,它只是一个循环缓冲区。因为你回答了隐含的问题,我现在暂时不接受,但不是明确的问题。 - bcrist
@bcrist你的意思是,你如何验证是否 jump() 实际上相当于2 ^ 512 next() ? - deviantfan
理想情况下,算法生成魔术位向量 JUMP 任意的 ñ。 - bcrist
@bcrist本文第8章究竟有什么不清楚(假设前面的7章也被阅读过了)? - deviantfan


正如维尼亚自己说的那,这实际上是一个错误。

在开发Java实现时,我发现,如果没有弄错的话,对正确的实现有一点小的改进:

如果你从p到p-1循环更新t数组,那么你可以将它回忆到状态并且它将正常工作。

此外,循环更新t变得更紧,因为您不需要每次都添加p + j。例如:

    int j = p;
    do {
        t[j] ^= s[j];
        ++j;
        j &= 15;
    } while (j != p);

好吧,正如bcrist正确指出的那样,之前的代码是错误的,因为p会改变JUMP数组中的每个位。我提出的最佳选择如下:

void jump() {
    static const uint64_t JUMP[] = { 0x84242f96eca9c41dULL,
        0xa3c65b8776f96855ULL, 0x5b34a39f070b5837ULL, 0x4489affce4f31a1eULL,
        0x2ffeeb0a48316f40ULL, 0xdc2d9891fe68c022ULL, 0x3659132bb12fea70ULL,
        0xaac17d8efa43cab8ULL, 0xc4cb815590989b13ULL, 0x5ee975283d71c93bULL,
        0x691548c86c1bd540ULL, 0x7910c41d10a1e6a5ULL, 0x0b5fc64563b3e2a8ULL,
        0x047f7684e9fc949dULL, 0xb99181f2d8f685caULL, 0x284600e3f30e38c3ULL
    };

    uint64_t t[16] = { 0 };
    const int base = p;
    int j = base;
    for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
        for(int b = 0; b < 64; b++) {
            if (JUMP[i] & 1ULL << b) {
                int k = p;
                do {
                    t[j++] ^= s[k++];
                    j &= 15;
                    k &= 15;
                } while (j != base);
            }
            next();
        }

    memcpy(s, t, sizeof t);
}

由于p最终将具有其原始值,因此应该可行。

不太确定它是否实际上是性能的提高,因为我正在为增量和按位AND交换一个附加项。
我认为它不会慢,即使增量和添加一样昂贵,因为j和k更新之间缺乏数据依赖性。希望它可能稍快一点。

意见/更正非常受欢迎。


0
2017-11-25 10:52



你测试过这个吗?除非我误解了你的建议,否则我觉得它不会起作用。 next() 增量台/套 p 所以为了使它正常工作,你需要旋转元素 t 每一次也是如此。 - bcrist
你是对的,@ bcrist!对我感到羞耻...我用另一个(希望是正确的)提议来编辑我的答案,虽然这次没有明显更快...... - Juan