【问题标题】:Miller Rabin Primality test accuracy米勒拉宾素性测试精度
【发布时间】:2014-07-28 13:55:01
【问题描述】:

我知道Miller–Rabin primality test 是概率性的。但是我想将它用于programming task,这样就不会出错。

如果输入数字是 64 位整数(即 C 中的 long long),我们能否假设它以非常高的概率正确?

【问题讨论】:

  • 在我看来,long long 被定义为最小为 64 位。

标签: c algorithm cryptography primes primality-test


【解决方案1】:

对于 64 位值的 MR 测试有有效的deterministic variants - 依赖于 GRH - 已经通过利用 GPU 和其他已知结果进行了详尽的测试。

我列出了我编写的用于测试任何 64 位值的素数的 C 程序的相关部分:(n > 1),使用 Jaeschke 和 Sinclair 的确定性 MR 变体的基础。它利用 gcc 和 clang 的 __int128 扩展类型进行求幂。如果不可用,则需要显式例程。也许其他人会觉得这很有用...

#include <inttypes.h>

/******************************************************************************/

static int sprp (uint64_t n, uint64_t a)
{
    uint64_t m = n - 1, r, y;
    unsigned int s = 1, j;

    /* assert(n > 2 && (n & 0x1) != 0); */

    while ((m & (UINT64_C(1) << s)) == 0) s++;
    r = m >> s; /* r, s s.t. 2^s * r = n - 1, r in odd. */

    if ((a %= n) == 0) /* else (0 < a < n) */
        return (1);

    {
        unsigned __int128 u = 1, w = a;

        while (r != 0)
        {
            if ((r & 0x1) != 0)
                u = (u * w) % n; /* (mul-rdx) */

            if ((r >>= 1) != 0)
                w = (w * w) % n; /* (sqr-rdx) */
        }

        if ((y = (uint64_t) u) == 1)
            return (1);
    }

    for (j = 1; j < s && y != m; j++)
    {
        unsigned __int128 u = y;
        u = (u * u) % n; /* (sqr-rdx) */

        if ((y = (uint64_t) u) <= 1) /* (n) is composite: */
            return (0);
    }

    return (y == m);
}

/******************************************************************************/

static int is_prime (uint64_t n)
{
    const uint32_t sprp32_base[] = /* (Jaeschke) */ {
        2, 7, 61, 0};

    const uint32_t sprp64_base[] = /* (Sinclair) */ {
        2, 325, 9375, 28178, 450775, 9780504, 1795265022, 0};

    const uint32_t *sprp_base;

    /* assert(n > 1); */

    if ((n & 0x1) == 0) /* even: */
        return (n == 2);

    sprp_base = (n <= UINT32_MAX) ? sprp32_base : sprp64_base;

    for (; *sprp_base != 0; sprp_base++)
        if (!sprp(n, *sprp_base)) return (0);

    return (1); /* prime. */
}

/******************************************************************************/

请注意,MR (sprp) 测试稍作修改以通过值在基数是候选者的倍数的迭代中,如网站的“备注”部分所述


更新:虽然这比Niklas' 答案的基础测试更少,但重要的是要注意基础:{3, 5, 7, 11, 13, 17, 19, 23, 29} 提供了一种廉价的测试,使我们能够消除超过:29 * 29 = 841 的候选人- 只需使用 GCD。

对于(n &gt; 29 * 29),我们可以清楚地消除任何偶数作为素数。小素数的乘积:(3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 29} = 3234846615,非常适合 32 位无符号值。 gcd(n, 3234846615) 比 MR 测试便宜很多!如果结果不是(1),那么(n) &gt; 841的系数很小。

Merten (?) 定理表明,这个简单的gcd(u64, u64) 测试消除了约 68% 的奇数候选(作为复合材料)。如果您使用 M-R 搜索素数(随机或增量),而不仅仅是“一次性”测试,那么这当然值得!

【讨论】:

  • 四分之二的答案不是“最”;其他两个答案在几年前提到了“确定性 MR”。这意味着您的开场白对于仅浏览此主题的专家来说是完全误导的。但是,您的答案确实提供了一些新的东西,因为您提供了出色的 C 版本,而其他人只提供链接或 C++,因此 +1(主要是因为我喜欢干净、高效的专业代码)。
  • @DarthGizka - 我们必须同意不同意。使用像 BPSW 或 Lucas 测试这样的大锤的想法具有误导性 - 并且超出了 OP 请求的范围。 Jacobi 求和实现的隐含复杂性 - 本身并不是一项简单的练习 - 以及 SPRP-2 在任何情况下都消除了绝大多数复合材料的事实;更不用说事实错误,比如需要 9 个碱基,或者它对黎曼假设的依赖,这仅与 广义 界有关,而不是辛克莱等人的详尽检验。
  • 你是对的 - 在仔细查看 Niklas B 的帖子后,我不得不承认大多数人(Niklas 的想法是正确的,但与 user448810 不同,他并没有完全做到这一点)。再次感谢您发布您的代码 - 与大量散文相比,这样的简短而紧凑的文章可以提供更多有用的信息和见解。
  • @DarthGizka - 值得一提的是,我的开场白有点强,我的辩护可能有点迂腐。写作时总是很容易误解语气。我应该更外交一点:)
【解决方案2】:

你的电脑并不完美;它有一定的失败概率,导致计算不正确。如果 M-R 测试给出错误结果的概率大大低于其他计算机故障的概率,那么你就可以了。没有理由将 M-R 测试运行少于 64 次迭代(2 ^ 128 中的 1 次错误)。大多数示例在最初的几次迭代中都会失败,因此只有实际的素数会被彻底测试。使用 128 次迭代,错误概率为 1 in 2^256。

【讨论】:

  • 不幸的是,这个答案没有告知这个范围内确定性 M-R 测试的现实。概率逻辑大致正确,但与 2^64 范围内的确定性测试并不真正相关
【解决方案3】:

Miller–Rabin 确实是概率性的,但您可以任意用准确性换取计算时间。如果你测试的数字是素数,它总是会给出正确的答案。有问题的情况是当一个数字是复合的,但据报道是素数。我们可以使用formula on Wikipedia来限制这个错误的概率:如果你随机选择k不同的碱基并测试它们,错误概率小于4-k。所以即使使用k = 9,你也只有百万分之三的机会出错。而使用k = 40 左右,它变得非常不可能。

也就是说,有一个deterministic version of Miller–Rabin,依赖于广义黎曼假设的正确性。对于范围 u 最多264,检查a = 2, 3, 5, 7, 11, 13, 17, 19, 23就够了。 I have a C++ implementation online 在许多编程竞赛中经过现场测试。这是无符号 64 位整数模板的实例化:

bool isprime(uint64_t n) { //determines if n is a prime number
    const int pn = 9, p[] = { 2, 3, 5, 7, 11, 13, 17, 19, 23 };
    for (int i = 0; i < pn; ++i)
        if (n % p[i] == 0) return n == p[i];
    if (n < p[pn - 1]) return 0;
    uint64_t s = 0, t = n - 1;
    while (~t & 1)
        t >>= 1, ++s;
    for (int i = 0; i < pn; ++i) {
        uint64_t pt = PowerMod(p[i], t, n);
        if (pt == 1) continue;
        bool ok = 0;
        for (int j = 0; j < s && !ok; ++j) {
            if (pt == n - 1) ok = 1;
            pt = MultiplyMod(pt, pt, n);
        }
        if (!ok) return 0;
    }
    return 1;
}

PowerModMultiplyMod 只是在给定模数下使用平方和-{multiply,add} 进行乘法和求幂的基元。

【讨论】:

  • 模板有什么用?该问题被标记为 c
  • 参考/推理?
  • @CristianCiupitu 只是将其视为伪代码。 Niklas 几乎不能在这里发布一个 biginteger 库。 PowerModMultiplyMod 是什么很明显,但是实现它们很烦人,并且无助于理解算法。
  • @CristianCiupitu 这是从工作实现中复制的 C++ 代码,但您应该能够通过将每个出现的 T 替换为 uint64_tunsigned long long 来手动实例化模板。
  • MultiplyModPowerMod 的实现也可以在我链接的完整源文件中查找
【解决方案4】:

对于n n;见http://miller-rabin.appspot.com/

更快的素数检验对基数 2 执行强伪素检验,然后是卢卡斯伪素检验。它需要的时间大约是单个强伪素测试的 3 倍,因此比 7 碱基米勒拉宾测试快两倍多。代码更复杂,但并不令人生畏。

如果您有兴趣,我可以发布代码;在 cmets 中告诉我。

【讨论】:

  • 对我来说,这组基数适用于所有 n 除了,它对 n=4033 和 n=4681 失败,因为这两个数字是基数 2 和 235。据我了解,使用 Miller-Rabin,您只测试小于 n 的基数,因此在 n=4033 和 n=4681 的情况下,这意味着仅测试基数 2 和 235。怎么可能这逃脱了 Jim Sinclair 的注意(被引用为这个 7 基组的发现者)?还是我在这里遗漏了什么?
  • 是 325,而不是 235。你测试了所有七个碱基。
  • 啊,是的,235,不是 325,谢谢。上面的手指很滑。我的代码确实说 325(我使用剪切和粘贴来插入集合)。其他人向我指出,即使你应该测试所有 7 个碱基,你实际上应该忽略 n 倍数的碱基(例如,当测试 n = 5,您应该忽略基数 235、9375 和 450775,我没有这样做)。我刚刚修改了我的程序以在 a ≡ 0 (mod n) 时忽略 a,它现在似乎在做正确的事情。实际上,如果 a ≡ 0 (mod n) 而不是忽略 a,我会返回 true(可能是素数)。
  • 在任何情况下,绝大多数综合候选人都将无法通过 2-SPRP 测试。例如,所有 24 位奇数组合中只有 63 个通过 2-SPRP 测试。我认为你低估了卢卡斯伪素测试的复杂性,雅各比和。如果我们谈论的是加密整数大小,我可能会同意,但我认为这个范围内 平均而言“更快”的想法是误导性的。
  • 见我的blog
【解决方案5】:

在 Miller-Rabin 的每次迭代中,您需要选择一个随机数。如果你不走运,这个随机数不会显示某些组合。一个小例子就是2^341 mod 341 = 2,通过了测试

但测试保证它只让复合材料以

你应该看看Baillie–PSW primality test。虽然它可能有误报,但没有已知的例子,根据维基百科,verified 没有低于 2^64 的合数通过测试。所以它应该符合您的要求。

【讨论】:

  • 你能在这里引用算法吗?此外,维基百科的链接会很好。那么我会认为这是一个很好的答案。
  • @JanDvorak Markdown 不接受维基百科链接中的破折号。
  • 谢谢@CodesInChaos,至少你回答了!!
猜你喜欢
  • 2011-12-02
  • 1970-01-01
  • 1970-01-01
  • 2017-04-24
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2018-11-05
  • 2011-04-13
相关资源
最近更新 更多