Miller_Rabin

#include <iostream>
#include <string>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#define range(i, a, b) for (int i = a; i <= b; ++i)
#define rerange(i, a, b) for (int i = a; i >= b; --i)
#define LL long long
#define fill(arr, tmp) memset(arr, tmp, sizeof(arr))
using namespace std;
const int S = 20;
LL mult_mod(LL a, LL b, LL c)
{
    a %= c;
    b %= c;
    long long ret = 0;
    while (b)
    {
        if (b & 1)
        {
            ret += a;
            ret %= c;
        }
        a <<= 1;
        if (a >= c)
            a %= c;
        b >>= 1;
    }
    return ret;
}
LL pow_mod(LL x, LL n, LL mod)
{
    if (n == 1)
        return x % mod;
    x %= mod;
    LL tmp = x;
    LL ret = 1;
    while (n)
    {
        if (n & 1)
            ret = mult_mod(ret, tmp, mod);
        tmp = mult_mod(tmp, tmp, mod);
        n >>= 1;
    }
    return ret;
}
bool check(LL a, LL n, LL x, LL t)
{
    LL ret = pow_mod(a, x, n);
    LL last = ret;
    range(i, 1, t)
    {
        ret = mult_mod(ret, ret, n);
        if (ret == 1 && last != 1 && last != n - 1)
            return true;
        last = ret;
    }
    if (ret != 1)
        return true;
    return false;
}
bool Miller_Rabin(LL n)
{
    if (n < 2)
        return false;
    if (n == 2)
        return true;
    if ((n & 1) == 0)
        return false;
    LL x = n - 1;
    LL t = 0;
    while ((x & 1) == 0)
    {
        x >>= 1;
        t++;
    }
    range(i, 0, S - 1)
    {
        LL a = rand() % (n - 1) + 1;
        if (check(a, n, x, t))
            return false;
    }
    return true;
}
LL factor[100];
int tol;
LL gcd(LL a, LL b)
{
    if (a == 0)
        return 1;
    if (a < 0)
        return gcd(-a, b);
    while (b)
    {
        long long t = a % b;
        a = b;
        b = t;
    }
    return a;
}
LL Pollard_rho(LL x, LL c)
{
    LL i = 1, k = 2;
    LL x0 = rand() % x;
    LL y = x0;
    while (1)
    {
        i++;
        x0 = (mult_mod(x0, x0, x) + c) % x;
        LL d = gcd(y - x0, x);
        if (d != 1 && d != x)
            return d;
        if (y == x0)
            return x;
        if (i == k)
        {
            y = x0;
            k += k;
        }
    }
}
void findfac(LL n)
{
    if (Miller_Rabin(n))
    {
        factor[tol++] = n;
        return;
    }
    LL p = n;
    while (p >= n)
        p = Pollard_rho(p, rand() % (n - 1) + 1);
    findfac(p);
    findfac(n / p);
}
int main()
{
    long long n;
    while (scanf("%lld", &n) != EOF)
    {
        tol = 0;
        /*
        findfac(n);
        for(int i=0;i<tol;++i)cout<<factor[i]<<" ";
        printf("
");
        */
        if (Miller_Rabin(n))
            printf("Yes
");
        else
            printf("No
");
    }
    return 0;
}
原文地址:https://www.cnblogs.com/HexQwQ/p/13438675.html