质数

质数（又称素数）：只能被1和自身整除的大于1的正整数。

质数是初等数论的重要研究对象，有关质数的常用算法有如下一些：

线性筛法：O(n)时间筛出n以内的所有质数

质数判定：判断给定数n是否是质数。大质数判定有Miller Rabin测试法

合数分解：找到给定合数n的一个非平凡因子。大数分解有Pollard Rho分解法

BTW：一些和整除、因子、质数相关的题目，先做质因数分解会后豁然开朗，int范围内的数至多有1300多个因子，先分解再枚举所有因子，甚至对所有因子双层循环的时间都是可以接受的。

下面给出一个包含线性筛、Miller Rabin测试、Pollard Rho分解、质因数分解的质数工具类

#include <iostream>
#include <cmath>
#include <cassert>
#include <map>
#include <cstring>
#include <vector>
using namespace std;

typedef long long ll;

template<typename T>
T gcd(T a, T b) { while(b) { T r = a%b; a = b; b = r;} return a;}
template<typename T>
T lcm(T a, T b) { return a/gcd(a,b)*b; }

// 打表大小，[0, UPBOUND)查表，[UPBOUND, UPBOUND*UPBOUND)根据表测试，[UPBOUND*UPBOUND, ll_max)米勒罗宾测试
const ll UPBOUND = 5e6;
class CPrime {
public:
    ll modMul(ll a,ll b,ll m) {
        ll t=0;
        a=(a%m+m)%m;
        b=(b%m+m)%m;
        while(b){
            if(b&1){
                t=(t+a)%m;
            }
            a=(a+a)%m;
            b>>=1;
        }
        return t;
    }

    ll modExp(ll a,ll b,ll m) {
        ll t=1,y=a%m;
        while(b){
            if(b&1){
                t=modMul(t,y,m);
            }
            y=modMul(y,y,m);
            b=(b>>1);
        }
        return t;
    }

    bool miller_rabin(ll n,ll b) {
        ll m=n-1;
        int j=0;
        while(!(m&1)){
            m>>=1;
            j++;
        }
        ll v=modExp(b,m,n);
        if(v==1 || v==n-1)return 1;
        for(int i=1;i<j;i++){
            v=modMul(v,v,n);
            if(v==n-1)return 1;
        }
        return 0;
    }

    bool _isprime(ll n)
    {
        const int K=10;//Miller_Rabin的偏真正确率为75%， isprime的正确率为1-(1/4)^K
        if(n<2)return 0;
        if(n==2 || n==3)return 1;
        for(int i=0;i<K;i++){
            if(!miller_rabin(n,rand()%(n-2)+2))return 0;
        }
        return 1;
    }

    //pollard_rho分解，给出N的一个非1因数，返回N时为一次没有找到,C为一个[1,N]的随机数 
    ll pollard_rho(ll C, ll N)
    {
        ll I, X, Y, K, D;
        I = 1;
        X = rand() % N;
        Y = X;
        K = 2;
        do{
            I++;
            D = gcd(N + Y - X, N);
            if (D > 1 && D < N) return D;
            if (I == K) Y = X, K *= 2;
            X = (modMul(X, X, N) + N - C) % N;
        }while (Y != X);
        return N;
    }

    //找出N的最小(质)因数
    ll _min_divisor(ll N)
    {
        if (isprime(N)) return N;
        while(1){
            ll T = pollard_rho(rand() % (N - 1) + 1, N);
            if (T < N){
                  ll A, B;
                  A = _min_divisor(T);
                  B = _min_divisor(N / T);
                  return A < B ? A : B;
            }
        }
    }

    int a[UPBOUND];
    int p[UPBOUND];
    int pn = 0;

    //号称线性的筛素数算法，实际性能确实不错 
    //p[]={2,3,5,7,...},pn为小于UPBOUND的素数个数 
    //若i是合数a[i]为i的最小因子，若i是素数a[i]=0 
    void primefilter(){
        int i, j;
        for(i = 2; i < UPBOUND; ++i){
            if(!(a[i])) p[pn++] = i;
            for(j = 0; (j<pn && i*p[j]<UPBOUND && (p[j]<=a[i]||a[i]==0)); ++j) {
                a[i*p[j]] = p[j];
            }
        }
    }
public:
    CPrime(){
        memset(a, 0, sizeof(a));
        memset(p, 0, sizeof(p));
        primefilter();
    }

    bool isprime(ll n){
        if(n<UPBOUND) return (n>0?!a[n]:0);
        int c = min(3500, pn-1); // 3500个足够支持到1e9
        if(n<p[c]*p[c]) {
            for(int i=0; i<c; i++) {
                if(n%p[i]==0) return false;
            }
            return true;
        } else {
            int c = min(1000, pn);
            for(int i=0; i<c; i++) {
                if(n%p[i]==0) return false;
            }
            return _isprime(n);
        }
    }

    // 返回质数表第n个质数，p[0]=2, ..., p[n]<UPBOUND
    int nth_prime(int n) {
        return n<pn?p[n]:-32768;
    }

    // 返回n的最小因子
    ll min_divisor(ll n) {
        if(n<UPBOUND) return a[n]?a[n]:n;
        int c = min(3500, pn-1); // 3500个足够支持到1e9
        if(n<p[c]*p[c]) {
            for(int i=0; i<c; i++) {
                if(n%p[i]==0) return p[i];
            }
            return n;
        } else {
            // miller_rabin 和 pollard_rho 需要较多的运算
            // 如果有一个小因子不妨先找出来
            int c = min(1000, pn);
            for(int i=0; i<c; i++) {
                if(n%p[i]==0) return p[i];
            }
            return _min_divisor(n);
        }
    }

    // 分解质因数
    vector<pair<ll, int>> factorize(ll n) {
        vector<pair<ll, int>> result;
        while(n!=1) {
            ll p = min_divisor(n);
            int c = 0;
            do {
                n/=p;
                ++c;
            } while(n%p==0);
            result.emplace_back(make_pair(p, c));
        }
        return result;
    }

    // 欧拉函数
    ll phi(ll n) {
        auto factors = factorize(n);
        for(auto [p, c] : factors) {
            n = n / p * (p-1);
        }
        return n;
    }
} Prime;