判断素数（Miller-Rabin、筛选）与 求素数因子（Pollard rho、试除）

Miller-Rabin算法：

用来确定一个指定的数N是否是素数。

---

费马小定理：

如果p是素数，那么p满足：（对任意a>=1）

如果我们指定一个a,如果n不满足费马定理，那么n不是素数（即合数）。（逆否）
但，在基数a下，如果n满足费马定理，n有可能是素数。
用一个样本a判断，过程会出错，但机会非常少：（不证明）
在基数2下，10000内的n值，有22个值判断错误。
对n随机选取的位数越多，出错的概率趋近于0。

---

一个推论

如果存在数x满足 x^2=1(mod n) 且 x ≠1或-1 则n是合数

证明：如果有n满足上式，不可能是奇素数
又n=2时，x=1(mod 2),n不可能是素数。

---

构造Miller-Rabin算法：

1、bool withness(a,n)
判断n是否为合数。
为求得 的结果
我们把 a^p-1拆分进行计算：
设u为 n-1的二进制数从左到右最后一个1的十进制数
设t为n-1的二进制尾部0的个数
可得 n-1 = 2^t * u

a^(n-1)=(a^u)^(2^t) (mod n)
通过先计算 a^u mod n 然后对结果进行连续平方t次来计算。（模乘法定理）

long long multi(long long x ,long long y,long long n){
    long long ans=0;
    x%=n;
    while(y){
        if(y&1){
            ans=(ans+x)%n;
        }
        x=(x<<1)%n;
        y>>=1;
    }
    return ans;
}

long long cal(long long a,long long u,long long n){
    long long ans=1;
    long long temp=a%n;
    while(u){
        if(u&1){// 末尾为1
            ans=multi(ans,temp,n);
        }
        u>>=1;
        temp=multi(temp,temp,n);
    }
    return ans;
}

bool Withness(long long a,long long n){
    long long t=0;
    long long u=(n-1);
    while(!(u&1)){
        t++;
        u>>=1;
    }
    long long x=cal(a,u,n);
    //连续平方t次
    for(int i=0;i<t;i++){
        long long pre=x;
        x=multi(x,x,n);
        if(x==1&&pre!=-1&&pre!=n-1){
            return true;
        }
    }
    if(x!=1){
        return true;
    }
     return false;
}

bool MillerRabin(long long n,long long  s){
    //检测s次
    while(s--){
        long long a=1+rand()%(n-1);//随机数a的范围为【1~n-1】
        if(Withness(a,n)){
            return false;
        }
    }
    return true;
}

令人难过的是，学习完Miller-Rabin算法做题超时orz，还是靠普通筛选法过的（苦笑）。

---

筛选法

最普通的筛选法：把自身的2~n倍标记为合数。
模板：

vector <bool> composite;
//vector <bool> prime;
void Create(int n){
    composite.resize(n+1);
    for(int i=2;i<=sqrt(n);i++){
        if(!composite[i]){//该数i没有被判断过
            for(int j=i+i;j<n;j+=i){//把该数的所有倍数判为合数
                composite[j]=true;
            }
        }
    }
}

进行普通筛选时，我们可以发现，对于有些数，被筛选了两次：比如6，被2筛选过一次，被3筛选过一次。
所以我们对此进行优化，使一个数只被筛选过一次。
线性筛选法：对于一个数由它的最小素数因子筛选掉
实现：

• 遍历n，乘以每个已经找到的素数i，把n*i 踢出。
• 如果n能够被当前素数i整除，即代表n=i*m，此时踢出的 n*i (i*i*m)不是由它的最小素因子踢出，故不判断。

for(long long i=2;i<=(maxn);i++){
        if(!composite[i]){//该数i没有被判断过
            prime[k++]=i;
        }
        for(long long j=0;j<k&&i*prime[j]<maxn;j++){//把所有素数*i判为合数
                                                                composite[i*prime[j]]=true;
                if(!(i%prime[j]))break;//是否由最小素因子踢出
        }
    }

---

求素数因子

试除法：循环1~sqrt(n),一个个试除。

void Factor(int n){
    vector <int> res;
    for(int i=2;i<=sqrt(n);i++){
        while(!n%i){//找到因子
            res.push_back(i);
            n/=i;
        }
    }
}

---

Pollard rho 求大素数因子

上述的Miller-Rabin 用来判断一个大素数，Pollard rho用也会运用到：
步骤：

1. 判断数n是否为素数（Miller-Rabin）
2. 找出n的一个因子i（可以不是素因子）
3. 此时你获得了两个因子i,n/i，分别对这两个因子递归再求其因子
   pollard rho 的思想和试除的思想类似，区别在于找因子过程不同。

找因子是一个数学过程：

1. 随机选取一个数 X1

2. 在X1的情况下找到一个X2，满足：gcd(x1-x2,n)>1 ;
   (说明：X1,X2之差 与 n 存在最大公约数)

3. 找到 n 的一个因子 ，即 gcd(x1-x2,n)

如何找X2：x[i]=(x[i-1]*x[i-1]%n+c)%n

求最大公约数：

long long gcd(long long a, long long b) { //(a>b)
    return b ? gcd(b, a%b) : a;
}

long long Pollard_Rho(long long n, int c) {
    long long i = 1, k = 2, x = rand()%(n-1)+1, y = x;
    while(true) {
        i++;
        x = (multi(x, x, n) + c)%n;
        long long p = gcd((y-x+n)%n, n);
        if(p != 1 && p != n) return p;
        if(y == x) return n;
        if(i == k) {
            y = x;
            k <<= 1;
        }
    }
}

void find(long long n, int c) {
    if(n == 1) return;
    if(Miller_Rabin(n)) {
        f.push_back(n);
        return;
    }
    long long p = n, k = c;
    while(p >= n) p = Pollard_Rho(p, c--);
    find(p, k);
    find(n/p, k);
}

Miller-rabin 与 Pollard rho 总和：

#define S 20

using namespace std;
typedef unsigned long long LL;
vector <long long> f;
long long factor;

long long gcd(long long a, long long b) { //(a>b)
    return b ? gcd(b, a%b) : a;
}

//计算 (a*b)%c.   a,b都是long long的数，直接相乘可能溢出的
//  a,b,c <2^63
long long mult_mod(long long a,long long b,long long 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;
}



//计算  x^n %c
long long pow_mod(long long x,long long n,long long mod)//x^n%c
{
    if(n==1)return x%mod;
    x%=mod;
    long long tmp=x;
    long long ret=1;
    while(n)
    {
        if(n&1) ret=mult_mod(ret,tmp,mod);
        tmp=mult_mod(tmp,tmp,mod);
        n>>=1;
    }
    return ret;
}


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_mod(a,u,n);
        for(int j=0;j<t;j++)
        {
            y=mult_mod(x,x,n);
            if (y==1&&x!=1&&x!=n-1)
                return false;
            ///其中用到定理，如果对模n存在1的非平凡平方根，则n是合数。
            ///如果一个数x满足方程x^2≡1 (mod n),但x不等于对模n来说1的两个‘平凡’平方根：1或-1，则x是对模n来说1的非平凡平方根
            x=y;
        }
        if (x!=1)///根据费马小定理,若n是素数，有a^(n-1)≡1(mod n).因此n不可能是素数
            return false;
    }
    return true;
}
long long Pollard_Rho(long long n, int c) {
    long long i = 1, k = 2, x = rand()%(n-1)+1, y = x;
    while(true) {
        i++;
        x = (mult_mod(x, x, n) + c)%n;
        long long p = gcd((y-x+n)%n, n);
        if(p != 1 && p != n) return p;
        if(y == x) return n;
        if(i == k) {
            y = x;
            k <<= 1;
        }
    }
}

void Find(long long n) {
    if(n == 1) return;
    if(miller_rabbin(n)) {
        f.push_back(n);
        return;
    }
    long long p = n;
    while(p >= n) p = Pollard_Rho(p, rand()%(n-1)+1);
    Find(p);
    Find(n/p);
}