素数测试

作者：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;
        }
    }
}