Miller-rabin素性测试与pollard-rho快速质因数分解

~~~~

Miller-Rabin

可以在O(log n)的时间内检测一个数是否为质数。
理论上是一个概率算法，但在一定范围以内的数经过验证全部可以判出来。
所以OI范围内可以视作是一个确定算法。

费马小定理

若p是质数，且$a<p$，则有
$a^{p-1}=1(mod p)$

其逆定理不一定成立，但不符合逆定理的一定不是素数。
根据这个检测n的素性。

即：尝试若干个a（a与n不能有公约数（欧拉定理条件），即对于质数n而言a不能是n的倍数，对于合数而言没关系了）
检测$a^{n-1} mod n是否为1$

然而这还不够强，会存在判不出的情况
对于Carmichael数，第一个是561，任意与他互质的数判出来都是Yes.

二次探测

显然地，对于质数n与正整数a<=n
$a^2 = 1(mod n)$，则有$n|(a+1)(a-1)$
又因为n是质数，$a<=n$，所以a=1或a=n-1.

与费马小相同，若不符合此条件则n一定不是质数。

both are ok

两个加起来，设mr(n,a)表示以a为底数，判定n的返回结果。
将n-1中所有约数2提出来。设去除所有2之后的的数是d.
从d开始二次探测即可。最后考虑一下费马小。

int mr(ll x,ll a) {
    a%=x;
    if (a==0) return 1; //注意a不能是x的倍数。
    ll d=x-1;
    while (!(d&1)) d>>=1;
    ll t=ksm(a,d,x);
    for (; d<x; d<<=1) {
        if (t==1) return 1; //后面所有t都是1，没必要继续判了。
        if (t!=x-1) {
            t=t*t%x;
            if (t==1) return 0; //若后一个是1，且当前不是x-1与1,那么合数。
        } else t=t*t%x;
    }
    return t==1;  //要是撑到最后就要看费马小的结论了。
}

以下数据来自matrix67
1. 如果被测数小于4 759 123 141，那么只需要测试三个底数2, 7和61就足够了
2. 如果你每次都用前7个素数(2, 3, 5, 7, 11, 13和17)进行测试，所有不超过341 550 071 728 320的数都是正确的。
3. 如果选用2, 3, 7, 61和24251作为底数，那么10^16内唯一的强伪素数为46 856 248 255 981。
当然最好需要rand若干个数。

由上可知，检测一个数n的时间是logn的。

pollar rho

据说时间复杂度是$n^{0.25}$，反正比较玄学吧。

利用生日悖论的结论，作差撞因子会比直接枚举概率大（很多）（超级反直觉）。发现可以利用gcd求约数，用log的时间来做一个炒鸡优化。

大概做法

假设我要分解n，那么
首先每一次用miller rabin判一下n是否质数。

然后定义个随机函数f(x)=（x^2+c）%c

两个数a,b，a随机，b在a的基础上走一步（取一个f），他们的差是d，将d的绝对值与n求gcd，若不为1或n则找到了一个因子，分下去即可。

若为1或n，那么要对a，b进行更改，可以用一个随机函数f(x)=(x^2+c)%mo，c是rand()出来的一个常数。

如果能保证每一对a,b都是有效的，一直做下去很快就能找到因子。（生日悖论）
但会发现可能有循环节，这时候无论多少次都是是出不了解的。

用一下floyed的判圈思想，b每次走2步，a每次走1步。这样当a,b重合时（套1~2圈）就不需要继续判断了，重新找c（重构随机函数）。

为什么a,b初值要间隔一步： 若a,b一开始取相同值，则4是无论如何也判不出来的。 大概和他只有长度为2的圈有关吧。

bzoj4802，注意有longlong相乘。

#include <cstdio>
#include <iostream>
#include <cstdlib>
#include <ctime>
#define abs(a) ((a)>0?(a):-(a))
#include <algorithm>
using namespace std;
typedef long long ll;
typedef long double ld;
ll gcd(ll a,ll b) {return b==0?a:gcd(b,a%b);}

ll _(ll a,ll b,ll c) {
    if (a>=c)a%=c; if (b>=c)b%=c;
    ll d=(ld)a*b/c;
    ll r=a*b-d*c;
    return (r%c+c)%c;
}
ll n;
ll ksm(ll a,ll b,ll mo) {
    ll ret=1;
    for (ll i=1; i<=b; i*=2) {
        if (b&i) ret=_(ret,a,mo);
        a=_(a,a,mo);
    }
    return ret;
}
int mr(ll x,ll a) {
    a%=x;
    if (a==0) return 1;
    ll d=x-1;
    while (!(d&1)) d>>=1;
    ll t=ksm(a,d,x);
    for (; d<x; d<<=1) {
        if (t==1) return 1;
        if (t!=x-1) {
            t=_(t,t,x);
            if (t==1) return 0; 
        } else t=_(t,t,x);
    }
    return t==1;
}
int isprime(ll x) {
    return mr(x,2) && mr(x,7) && mr(x,61) && mr(x,rand());
}
ll fo[10000];
int cs=0;

void pr(ll x) {
    if (isprime(x)) {
        fo[++fo[0]]=x;
        return;
    }
    int flag=fo[0];
    while (flag==fo[0]) {
        ll c=rand()%x,a=rand()%x;
        ll b=(_(a,a,x)+c)%x;
        while (a!=b) {
            cs++;
            ll det=gcd(abs(a-b),x);
            if (det>1 && det<x) {
                pr(x/det);
                pr(det);
                return;
            }
            a=(_(a,a,x)+c)%x;
            b=(_(b,b,x)+c)%x;
            b=(_(b,b,x)+c)%x;
        }
    }
}
int main() {
    freopen("pol.in","r",stdin);
//  freopen("pol.out","w",stdout);
    srand(time(0));

    rand();
    ll nn=0;
    cin>>nn;
    for (ll n=nn; n<=nn; n++) {
        if (n==1) {
            cout<<"1"<<endl; return 0;
        }
        fo[0]=0;
        pr(n);
        sort(fo+1,fo+1+fo[0]);
        ll ans=1;
        for (int i=1,cs=0; i<=fo[0]+1; i++) {
            if ( i==fo[0]+1 || i!=1 && fo[i-1] != fo[i]) {
                ll d = ksm(fo[i-1],cs,n); if (d==0) d=n;
                d-=ksm(fo[i-1],cs-1,n);
                ans*=d; 
                cs=1;
            } else cs++;
        }
        cout<<ans<<endl;
    }
}