MillerRabin质数判定+PollardRho质因子分解

参考 matrix67，kuangbin blog

测试 POJ1811

#include<stdio.h>
#include<string.h>
#include<iostream>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#include<assert.h>
using namespace std;


typedef long long LL;
#define maxn 10000

/***********************************************************
MillerRabin素数测试
理论基础：fermat定理 和 二次探测定理
如果你每次都用前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。
*/

typedef LL LINT;
const int CHECK_CNT = 9;
int check_cnt = 9;
LINT checkNumber[9] = {2, 3, 5, 7, 11, 13, 17, 61, 24251}; //选定一些判定性强的素数

/**返回(x*y) mod modv*/
LINT modMul(LINT x, LINT y, LINT modv)
{
    x %= modv, y %= modv;
    LINT rv = 0;
    for(; y > 0; y>>= 1)
    {
        if(y & 1)
        {
            rv += x;
            if(rv >= modv) rv -= modv;
        }
        x += x;
        if(x >= modv) x -= modv;
    }
    return rv;
}

/**返回(x^y) mod modv*/
LINT modPow(LINT x, LINT y, LINT modv)
{
    x %= modv;
    LINT rv = 1;
    for(; y > 0; y>>= 1)
    {
        if(y & 1) rv = modMul(rv, x, modv);
        x = modMul(x, x, modv);
    }
    return rv;
}

/**cnt表示第几次探测*/
LINT getCheckNum(int cnt, LINT maxv)
{
    if(cnt < check_cnt && checkNumber[cnt] < maxv)
    {
        return (checkNumber[cnt]);
    }
    else
    {
        return ( rand() % (maxv - 2) + 2 );
    }
}

/**以cur为基，x-1=u*2^pt，检验x是不是合数*/
bool isOKBySecondaryProbe(LINT cur, LINT x, LINT u, LINT pt)
{
    LINT lastv, val = modPow(cur, u, x);

    //   printf("ci = %d cur = %I64d val = %I64d\n", ci, cur, val);

    for(int pi = 0; pi < pt; pi++)
    {
        lastv = val;//记录方程 x^2 % p = 1 的解x
        val = modMul(val, val, x);

        //    printf("-->pi = %d cur = %I64d val = %I64d lastv = %I64d\n", pi, cur, val, lastv);
        if(1 == val && 1 != lastv && x - 1 != lastv)
        {
            return false;
        }
    }
    return bool(1 == val);
}

bool isPrimeByMillerRabin(LINT x)
{
    assert(x > 0);
    if(2 == x) return true;
    if(x < 2 || 0 == (x & 1)) return false;

    srand((LL)time(0));
    int pt = 0;
    LINT u = x - 1;
    while(0 == (u & 1)) u >>= 1, pt++;

    //  printf("x = %I64d u = %I64d pt = %d\n", x, u, pt);

    for(int ci = 0; ci < CHECK_CNT; ci++)
    {
        LINT cur = getCheckNum(ci, x);
        if(!isOKBySecondaryProbe(cur, x, u, pt)) return false;
    }
    return true;
}


/**************************************************************/
//pollard_rho 算法进行质因数分解
int tot;
LL factor[maxn];//结果无序

LL gcd(LL a,LL b){
    if (a==0) return 1;
    if (a<0) return gcd(-a,b);
    while (b){
        LL t=a%b; a=b; b=t;
    }
    return a;
}

LL PollardRho(LL x,LL c){
    LL i=1,x0=rand()%x,y=x0,k=2;
    while (true){
        i++;
        x0=(modMul(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 FindPrimeFactor(LL n){           //递归进行质因数分解N
    if (isPrimeByMillerRabin(n)) {factor[tot++] = n; return;}

    LL p = n;
    while(p >= n) p = PollardRho(p, rand() % (n-1) +1);
    FindPrimeFactor(p), FindPrimeFactor(n/p);
}

int main()
{
   // srand(time(NULL));//POJ上G++要去掉这句话
    int T;
    scanf("%d",&T);
    long long n;
    while(T--)
    {
        scanf("%I64d",&n);
        if (isPrimeByMillerRabin(n)) {printf("Prime\n"); continue; }
        tot = 0;
        FindPrimeFactor(n);
        long long ans=factor[0];
        for(int i=1;i<tot;i++)
          if(factor[i]<ans)ans=factor[i];
        printf("%I64d\n",ans);
    }
    return 0;
}