2

我正在寻找实施费马小定理进行质数测试。这是我写的代码：

lld expo(lld n, lld p) //2^p mod n
{
    if(p==0)
        return 1;
    lld exp=expo(n,p/2);
    if(p%2==0)
        return (exp*exp)%n;
    else
        return (((exp*exp)%n)*2)%n;
}

bool ifPseudoPrime(lld n)
{
    if(expo(n,n)==2)
        return true;
    else
        return false;
}

注意：我取了a(<=n-1)as的值2。

现在，数字 n 可以大到10^18。这意味着变量exp可以达到 附近的值10^18。这进一步意味着表达式(exp*exp)可以达到10^36导致溢出的高度。我该如何避免这种情况。

我对此进行了测试，它运行良好，直到10^9. 我正在使用C++

4

2 回答 2

8

如果模数接近您可以使用的最大整数类型的限制，事情就会变得有些复杂。如果您不能使用实现大整数运算的库，您可以通过将因子拆分为低阶和高阶部分来自己滚动模乘法。

如果模数m大到2*(m-1)溢出，事情就变得很麻烦，但如果2*(m-1)不溢出，那是可以忍受的。

让我们假设您拥有并使用 64 位无符号整数类型。

您可以通过将因子拆分为低 32 位和高 32 位来计算模积，然后将积拆分为

a = a1 + (a2 << 32)    // 0 <= a1, a2 < (1 << 32)
b = b1 + (b2 << 32)    // 0 <= b1, b2 < (1 << 32)
a*b = a1*b1 + (a1*b2 << 32) + (a2*b1 << 32) + (a2*b2 << 64)

计算a*b (mod m)，m <= (1 << 63)减少四个产品中的每一个模m，

p1 = (a1*b1) % m;
p2 = (a1*b2) % m;
p3 = (a2*b1) % m;
p4 = (a2*b2) % m;

合并这些变化的最简单方法是

for(i = 0; i < 32; ++i) {
    p2 *= 2;
    if (p2 >= m) p2 -= m;
}

p3和 64 次迭代相同p4。然后

s = p1+p2;
if (s >= m) s -= m;
s += p3;
if (s >= m) s -= m;
s += p4;
if (s >= m) s -= m;
return s;

这种方式不是很快，但是对于这里需要的少数乘法，它可能已经足够快了。应通过减少班次数来获得小幅加速；首先计算(p4 << 32) % m，

for(i = 0; i < 32; ++i) {
    p4 *= 2;
    if (p4 >= m) p4 -= m;
}

那么所有的p2,p3的当前值p4需要乘以 2 32模m,

p4 += p3;
if (p4 >= m) p4 -= m;
p4 += p2;
if (p4 >= m) p4 -= m;
for(i = 0; i < 32; ++i) {
    p4 *= 2;
    if (p4 >= m) p4 -= m;
}
s = p4+p1;
if (s >= m) s -= m;
return s;
于 2013-05-07T19:49:43.497 回答
0

您可以分几个阶段执行乘法。例如，假设您要计算 X*Y mod n。取 X 和 Y 并将它们写为 X = 10^9*X_1 + X_0, Y = 10^9*Y_1 + Y_0。然后计算所有四个乘积 X_i*Y_j mod n，最后计算 X = 10^18*(X_1*Y_1 mod n) + 10^9*( X_0*Y_1 + X_1*Y_0 mod n) + X_0*Y_0。请注意，在这种情况下，您使用的数字是允许的最大值的一半。

如果分成两部分还不够（我怀疑是这种情况），请使用相同的模式分成三部分。一分为三应该可行。

一个更简单的方法是乘以学校的方式。它对应于前一种方法，但在与数字一样多的部分中写入一个数字。

祝你好运！

于 2013-05-07T19:48:09.980 回答