费马小定理及其应用

费马小定理：若一个数p是素数，a为整数，且gcd(a,p) == 1，则 

简单点说就是：假设a为整数，p为素数，且gcd(a,p)==1，那么a的(p-1)次方除以p的余数一定是 1 

注意：费马小定理只是素数判定的一个必要条件，素数一定满足费马小定理，满足费马小定理的数，却不一定是素数，例如Carmichael数（Carmichael数都是符合费马小定理的，但是他们都是合数）。

Euler定理：

　　我们记φ(m)为小于等于m的整数中与m互素的数的个数(这是欧拉函数的定义)。就有：

 　　

　　证明也很简单：这实际上就是费马小定理的一个推广，我们知道当m为素数时，φ(m) = m-1，这是欧拉函数的结论。将φ(m)带入就可以得到费马小定理的结论。

　　关于欧拉函数可以看这篇blog：https://www.cnblogs.com/DynastySun/p/9364673.html

逆元：

　　这个定理十分有用，我们可以证明一下费马小定理求逆元的做法，我们知道若一个数 M 是素数，则任意整数 A（要求gcd(A, M) == 1） 关于 M 的逆元就是 pow(A, M-2)，我们可以证明一下：

　　求A关于M的逆元，就需要找出满足 A * X ≡ 1 (mod M) 这个式子的 X ，因为 M 为素数，所以 A^(M-1) ≡ 1 (mod M)，我们一对比就可找到 X 的一个解 A^(M-2)。

　　  关于逆元可以看这篇blog：https://www.cnblogs.com/DynastySun/p/9362505.html

Miller_rabin算法：

　　我们可以用费马小定理来判断一个很大的数是否是素数，我们知道若一个数p是素数，则必定有 那么我们就可随机枚举出大概10个左右的数，判断这个数的p-1次方对p取模是否为1，但是这样不严谨，因为Carmichael数的存在，虽然Carmichael数的数量很少，但是这也会对算法造成影响。所以在介绍这个算法之前我们再来了解一个定理：二次探测定理，这个定理可以对算法进行改进以免将Carmichael数当做素数。

　　二次探测定理：如果p是素数，x是小于p的正整数，且, 则 x = 1 或 x = p-1。

　　有了二次探测定理我们就可以在使用费马小定理计算a^(p-1) ≡ 1(mod p)时，增加对p的二次探测，一旦发现违背二次探测条件，就可以得出p不是素数的结论。我们通过一个例子来利用这个结论，341是一个合数(341=11*31)，但是341可以通过a为2时的费马小定理2^340 ≡ 1(mod 341)，我们现在可以假设341为素数：因为(2^170)^2 ≡ 1(mod 341)，所以2^170 mod 341只能为1或者是340，我们发现确实有2^170 ≡ 1(mod 341)，我们继续往下算2^85 ≡ 32(mod 341)，这样一来就可以看出341不是素数了。

　　这就是Miller_rabin素性测试方法，不断地提取p-1中的因子2，这样我们为了方便可以把p-1表示为d*2^r（d为奇数），这样一来我们就需要不断计算a^(d*2^r)对p取余的结果，这样我们可以强化费马小定理为下面的形式：

　　尽可能的提取p-1中的因子2，把p-1表示为d*2^r，若p为素数，则a^d mod p = 1，或者存在一个i使得a^(d*2^i) % mod p = p-1 (0 <= i < r)。

　　需要说明的是，Miller_rabin素性测试也是不确定算法。我们把可以通过以a为底的Miller-Rabin测试的合数称作以a为底的强伪素数。

　　对于大数的素性判断，目前Miller-Rabin算法应用最广泛。一般底数仍然是随机选取，但当待测数不太大时，选择测试底数就有一些技巧了。比如，如果被测数小于4 759 123 141，那么只需要测试三个底数2, 7和61就足够了。当然，你测试的越多，正确的范围肯定也越大。如果你每次都用前7个素数(2, 3, 5, 7, 11, 13和17)进行测试，所有不超过341 550 071 728 320的数都是正确的。如果选用2, 3, 7, 61和24251作为底数，那么10^16内唯一的强伪素数为46 856 248 255 981。这样的一些结论使得Miller-Rabin算法在比赛中非常实用。通常认为，Miller-Rabin素性测试的正确率可以令人接受，随机选取k个底数进行测试算法的失误率大概为4^(-k)。

　　以上部分引用自：Matrix67的blog： http://www.matrix67.com/blog/archives/234

　　现在我们可以给出板子：

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
using namespace std;

typedef long long LL;
typedef unsigned long long ULL;

inline LL mut_mod(LL a, LL b, LL m){
    if(b < 0) {
        a = -a;
        b = -b;
    }

    LL ans = 0;
    a %= m;
    while(b){
        if(b&1) ans = (ans+a)%m;
        b >>= 1;
        a = (a+a)%m;
    }
    return ans;
}

inline LL pow_mod(LL a, LL b, LL m){
    LL ans = 1;
    a %= m;
    while(b){
        if(b&1) ans = mut_mod(ans, a, m);
        b >>= 1;
        a = mut_mod(a, a, m);
    }
    return ans;
}

LL primes[] = { 2, 3, 7, 61, 24251};

bool isprimes(LL n){
    int r = 0, t;
    LL x, y, tmp, d = n-1;
    while((d&1) == 0) d>>=1, r++;

    for(int i = 0; i < 5; i++){
        tmp = primes[i];
        if(tmp == n) return true;
        y = x = pow_mod(tmp, d, n);
        t = r;
        while(t--){
            y = mut_mod(x, x, n);
            if(y == 1 && x != 1 && x != n-1)
                return false;
            x = y;
        }
        if(y != 1) return false;
    }
    return true;
}

int main(){
    LL n;
    while(scanf("%lld", &n) != EOF){
        if(n > 1 && isprimes(n)) printf("YES
");
        else printf("NO
");
    }
    return 0;
}