在 C 中实现 Miller-Rabin

Question

我正在尝试在 C99 中实现 Miller-Rabin 素性测试，但我在让它工作时遇到了一些问题。我制作了一个小型测试集来验证实现是否有效，以下是我检查素数的方法

int main() {
    int foo[11] = {0, 1, 2, 3, 4, 7, 28, 73, 125, 991, 1000};
    for (int i = 0; i < 11; i++) {
        printf("%s; ", isprime(foo[i], 5000) ? "Yes" : "No");
    }
    return 0;
}

根据列出的数字，预期输出为

不；不;是的;是的;不;是的;不;是的;不;是的;不；

但是，实现后，我得到的输出如下：

不；不;是的;是的;不;是的;不;不;不;不;不;

这是我编写算法的方式

int randint (int low, int up){
    return rand() % (++up - low)+low;
}

int modpow(int a, int b, int m) {
    int c = 1;
    while (b) {
        if (b & 1) {
            c *= a;
        }
        b >>= 1;
        a *= a;
    }
    return c % m;
}

bool witness(int a, int s, int d, int n) {
    int x = modpow(a,d,n);
    if(x == 1) return true;
    for(int i = 0; i< s-1; i++){
        if(x == n-1) return true;
        x = modpow(x,2,n);
    }
    return (x == n- 1);
}

bool isprime(int x, int j) {
    if (x == 2) {
        return true;
    }
    if (!(x & 1) || x <= 1) {
        return false;
    }
    int a = 0;
    int s = 0;
    int d = x - 1;

    while (!d&1){
        d >>=1;
        s+=1;
    }
    for(int i = 0; i < j; i++){
        a = randint(2, x-1);
        if(!witness(a,s,d,x)){
            return false;
        }
    }

    return true;
}

我做错了什么？为什么测试对于“大”素数失败，但对于非常小的素数却有效？我该如何解决这个问题？

Answer 1

使用 Visual Studio 2015 社区版我发现了两个问题。第一行：

while (!d&1){

需要是：

while (!(d&1)){

其次，正如评论中提到的，你的 modpow 函数溢出了。尝试：

int modpow(int a, int d, int m) {
    int c = a;
    for (int i = 1; i < d; i++)
        c = (c*a) % m;
    return c % m;
}

Answer 2

您的

modpow()

功能存在问题。 您可能想对参数和结果使用无符号类型（负数

m

意味着什么？）其次，它会溢出，因为它会在对

a^b

进行模减之前尝试计算

m

。您需要边做边减少

a

和

c

。

处理这个问题的最好方法是编写一些测试：

unsigned modpow(unsigned a, unsigned b, unsigned m) {
    unsigned c = 1;
    while (b) {
        if (b & 1) {
            c *= a;
        }
        b >>= 1;
        a *= a;
    }
    return c % m;
}

#include <stdio.h>
unsigned test(unsigned a, unsigned b, unsigned m, unsigned expected)
{
    unsigned actual = modpow(a, b, m);
    if (actual == expected)
        return 0;
    printf("modpow(%u, %u, %u) returned %u; expected %u\n",
           a, b, m, actual, expected);
    return 1;
}


int main()
{
    return test(0, 0, 2, 1)
        +  test(1005, 16, 100, 25)
        ;
}

第一个测试通过（但提出了一个问题 - 当

m < 2

时你想要什么结果？）；第二个失败了：

modpow(1005, 16, 100) 返回 49;预计25

让我们修改

modpow()

以减少每一步的结果：

unsigned modpow(unsigned a, unsigned b, unsigned m) {
    unsigned c = 1;
    while (b) {
        if (b & 1) {
            c *= a;
            c %= m;
        }
        b >>= 1;
        a *= a;
        a %= m;
    }
    return c;
}

现在已经过去了！ 我们可以进行另一个失败的测试：

int main()
{
    return test(0, 0, 2, 1)
        +  test(1005, 16, 100, 25)
        +  test(100000005, 16, 1000000000, 587890625)
        ;
}

modpow(100000005, 16, 1000000000) 返回 919214529;预计587890625

现在我们需要使用更大的类型来计算乘法：

unsigned modpow(unsigned a, unsigned b, unsigned m) {
    unsigned long long c = 1;
    while (b) {
        if (b & 1) {
            c = (unsigned long long)c * a % m;
        }
        b >>= 1;
        a = (unsigned long long)a * a % m;
    }
    return c;
}

一旦我们对

modpow()

函数有足够的信心，我们就可以调试算法的其余部分。

---

请注意，如果您的整数大小与我的不同，您将需要在测试中使用不同的值来复制结果。 我选择以

005

结尾的大数字，因为我们知道无论幂如何，最后两位数字都是不变的

25

。 您可能会发现

dc -e '???|p'

有助于生成测试用例（在其标准输入上提供三个参数，它将打印预期值）。

Answer 3

您可以使用此 C99 Miller-Rabin 实现：

typedef unsigned long long int ulong;

实用功能：

ulong mul_mod(ulong a, ulong b, const ulong mod) {
    ulong res = 0, c; // return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
    for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (c = b) >= mod - b ? c -= mod : 0, b += c);
    return res % mod;
}

ulong pow_mod(ulong n, ulong exp, const ulong mod) {
    ulong res = 1; // return (n ^ exp) % mod
    for (n %= mod; exp; exp & 1 ? res = mul_mod(res, n, mod) : 0, n = mul_mod(n, n, mod), exp >>= 1);
    return res;
}

米勒-拉宾算法：

int is_prime(ulong n) {
    // Perform a Miller-Rabin test, it should be a deterministic version.
    static const ulong primes[ ] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
    static const int n_primes = (int) sizeof(*primes);
    for (int i = 0; i < n_primes; ++i)
        if (n % primes[i] == 0) return n == primes[i];
    if (n < primes[n_primes - 1]) return 0;
    int res = 1, a = 0;
    ulong b, c;
    for (b = n - 1; ~b & 1; b >>= 1, ++a);
    for (int i = 0; i < n_primes && res; ++i)
        if (c = pow_mod(primes[i], b, n), c != 1) {
            for (int d = a; d-- && (res = c + 1 != n);)
                c = mul_mod(c, c, n);
            res = !res;
        }
    return res;
}