素数测试

作者：jostree 转载请注明出处 https://www.cnblogs.com/jostree/p/10274890.html

费马小定理

如果p是素数，a是小于p的正整数，那么a^(p-1) mod p = 1。

首先我们证明这样一个结论：如果p是一个素数的话，那么对任意一个小于p的正整数a，a, 2a, 3a, ..., (p-1)a除以p的余数正好是一个1到p-1的排列。例如，5是素数，3, 6, 9, 12除以5的余数分别为3, 1, 4, 2，正好就是1到4这四个数。

反证法，假如结论不成立的话，那么就是说有两个小于p的正整数m和n使得na和ma除以p的余数相同。不妨假设n>m，则p可以整除a(n-m)。但p是素数，那么a和n-m中至少有一个含有因子p。这显然是不可能的，因为a和n-m都比p小。

用同余式表述，我们证明了： (p-1)! ≡ a * 2a * 3a * … * (p-1)a (mod p)
也即： (p-1)! ≡ (p-1)! * a^(p-1) (mod p)
两边同时除以(p-1)!，就得到了我们的最终结论： 1 ≡ a^(p-1) (mod p)

Miller-Rabbin 素数测试

如果p是素数，x是小于p的正整数，且x^2 mod p = 1，那么要么x=1，要么x=p-1。这是显然的，因为x^2 mod p = 1相当于p能整除x^2-1，也即p能整除(x+1)(x-1)。由于p是素数，那么只可能是x-1能被p整除(此时x=1)或x+1能被p整除(此时x=p-1)。

从而我们选取不同的x，对p其进行进行素数测试，即可。代码如下：

#include <iostream>
#include <string>
#include <cstring>
#include <cstdlib>
#include  <ctime>

using namespace std;
typedef long long LL;
const int S = 100;

LL mul(LL a, LL b, LL n)
{
LL res = 0;
while( b )
{
if( b & 1 )
{
res += a;
res %= n;
}
a = (a + a) % n;
b = b >> 1;
}
return res;
}

LL pow(LL a, LL b, LL n)
{
LL res = 1;
while( b )
{
if( b & 1 )
{
res = mul(res, a, n);
}
a = mul(a, a, n);
b = b >> 1;
}
return res;
}

bool miller_rabbin(LL n)
{
if( n == 2 ) return true;
if(n < 2 || !(n & 1)) return false;
int t = 0;
LL a, x, y, u = n - 1;
while((u & 1) == 0) t++, u>>=1;
for( int i = 0 ; i < S ; i++ )
{
a = rand() % (n - 1) + 1;
x = pow(a, u, n);
for( int j = 0 ; j < t ; j++ )
{
y = mul(x, x, n);
if( y == 1 && x != 1 && x != n - 1 )
{
return false;
}
x = y;
}
if( x != 1 ) return false;
}
return true;
}

int main(int argc, char *argv[])
{
int T;
scanf("%d", &T);
LL N;
srand(time(0));
while( T-- )
{
cin>>N;
if( miller_rabbin(N) )
{
cout<<"Yes"<<endl;
}
else
{
cout<<"No"<<endl;
}
}
}