通过平方运算得到模块求幂的溢出可能性

Question

我期待实施fermat's little theorem进行主要测试。下面是我写的代码：

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)值作为2。

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

我测试了这个，它运行良好，直到10^9。我正在使用C++

Answer 1

如果模量接近您可以使用的最大整数类型的极限，事情会变得有些复杂。如果你不能使用实现biginteger算术的库，你可以通过分解低阶和高阶部分的因子来自动进行模乘。

如果模数m太大以至于2*(m-1)溢出，事情变得非常繁琐，但是如果2*(m-1)没有溢出，它是可以忍受的。

让我们假设你有和使用一个64位无符号整数类型。

可以通过分割因素低和高的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是相同的，对于p4有64次迭代。然后

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的并且需要利用2 模m相乘的p4当前值，

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; 

Answer 2

您可以在几个阶段执行乘法运算。例如，假设你想计算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。请注意，在这种情况下，您使用的数字是允许的最大值的一半。

如果拆分两部分不足（我怀疑是这种情况），使用相同的模式分成三部分。三分裂应该工作。

一个更简单的方法是繁殖学校的方式。它与先前的方法相对应，但是将其编号与它所具有的数字一样多地写入一个数字。

祝你好运！