BZOJ3667,4802:Rabin-Miller算法

传送门

题意：
大素数判断并分解。

题解：
Robin-Miller随机算法。
http://blog.csdn.net/thy_asdf/article/details/51347390

—————————————————————————————————————————
发现以前写得有点水啊，现在来重写一发。
Rabin-Miller检测素数：
首先，检测前9个质数$pr_i^{n-1}=1\pmod{n}$
其次，发现除了2，$n-1$为偶数（否则n不为质数）。那么可以顺便做二次检测：
$a^2=1\pmod{n}$，当n为质数时当且仅当$a=1$或$a=n-1$。

Pollard-Rho 大数分解：
设$n=p·q$，$x,y$均为随机数。
根据生日悖论，随机选取$\sqrt{n}$个$x,y$,存在$|x_i-y_j|$为$p,q$的概率为$50 \%$。
需要两两枚举$n$次。
考虑$gcd(|x_i-y_i|,n)$为$p,q$的概率，这样的差$p+q-2 = \sqrt{n}$种，那么只需要选取$n^{\frac{1}{4}}$个数，两两枚举$\sqrt{n}$次。

神奇的地方来了：
我们使用一个函数来生成伪随机。
换句话说，我们不断地 使用函数 使用函数 $f$ 来生成（看上去或者说感觉上像的）随机数。
并不是 所有的函都能够这样做，但一个神奇的函数可以。

$$f(x)=(x^2+a)\% n$$
( 我们可以自己指定$a$，也可以用 $rand()$随机函数来生成，这不是重点 )
我们从$x_1=2$ 或者其他数开始，让$x_2=f(x_1),x_3=f(x_2)...$，生成规则为$x_{n+1}=f(x_n)$

我们来看下面的伪代码:

x:2
while (.. exit criterion .. )
$\,\,\,\,$y := f(x);
$\,\,\,\,$p := GCD( | y - x | , N);
$\,\,\,\,$if ( p > 1)
$\,\,\,\,$$\,\,\,\,$ return “Found factor: p”; x := y;
end
return “Failed. :-(”

假设 $N$ = 55 , $f(x)=(x^2+2)\%55$

$x_n$	$x_{n+1}$	$gcd(x_n,x_{n+1},N)$
2	6	1
6	38	1
38	16	1
16	36	5

你可以发现对于大部分的数据这个算法能够正常运行，但是对于某些数据，它将会进入无限循环。为什么呢？这是因为存在$f$环的原因。
当它发生的时候，我们会在一个有限数集中进行无限循环
例如，我们可以构造一个伪随机函数并生成如下伪随机数：

$$2, 10, 16, 23, 29, 13, 16, 23, 29, 13 … …$$
在这个例子中，我们最终将会在$16, 23, 29, 13$这个圈中无限循环，永远找不到因子 。
那么，如何探测环的出现呢？
一种方法是记录当前产生过的所有的数$x_1, x_2, … … x_n$,并检测是否存在$x_l = x_n(l< n)$。
在实际过程中，当 $n$增长到一定大小时，可能会造成的内存不够用的情况。

另一种方法是由Floyd发明的一个算法，我们举例来说明这个聪明而又有趣的算法。
假设我们在一个很长很长的圆形轨道上行走，我们如何知道我们已经走完了一圈呢？
当然，我们可以像第一种方法那样做，但是更聪明的方法是让A 和B 按照B 的速度是
A 的速度的两倍从同一起点开始往前走，当B 第一次敢上A 时(也就是我们常说的套圈)， 我们便知道，B 已经走了至少一圈了。

a := 2;
b := 2;
while ( b != a )
$\,\,\,\,$a = f(a); // a runs once
$\,\,\,\,$b = f(f(b)); // b runs twice as fast. p = GCD( | b - a | , N);
$\,\,\,\,$p = GCD( | b - a | , N);
$\,\,\,\,$if ( p > 1)
$\,\,\,\,$$\,\,\,\,$return “Found factor: p”;
end
return “Failed. :-(“

发现以前学长写的一篇博客挺好的:关于Pollard_rho的期望复杂度证明

$n=p·q$,那么$f_1(x)=(x^2+c) \pmod{n}$的rho环比$f_2(x)=(x^1+c)\pmod{p}$的环要大，在$p$环上相遇的步数是期望$n^{\frac{1}{4}}$的，且进入$p$环时有很大可能没有进入$n$环，那么此时$x \not = y$,设$x=k_1*p+b_1$,$y=k_2*p+b_2$,$b_1=b_2$，那么此时可以求出 p <script type="math/tex" id="MathJax-Element-1005">p</script>这一因数。

还是贴4802的板吧，发现不倍长的话会被4卡T。。。

#include<bits/stdc++.h>
#include<tr1/unordered_map>
typedef long long ll;
typedef unsigned int uint;
using namespace std;
tr1::unordered_map<ll,int> S;
int a[10]={2,3,5,7,11,13,17,19,23,29};
vector<ll>fac;
inline uint unit(){
    static uint state0 = 19491001;
    state0^=(state0>>13);
    state0^=(state0<<17);
    state0^=(state0>>5);
    return state0;
}
inline ll ksc(ll a,ll b,ll mod){
    return (a*b-(ll)((long double)a/mod*b)*mod+mod)%mod;
}
inline ll power(ll a,ll b,ll mod){
    ll rs=1; a=a%mod;
    for(;b;b>>=1,a=ksc(a,a,mod)) if(b&1) rs=ksc(rs,a,mod);
    return rs;
}
inline bool check(ll x,ll n,ll t1,ll t2){
    x=power(x,t1,n); ll p=x;
    while(t2--){
        x=ksc(p,p,n); 
        if(x==1 && (p!=1 &&p!=n-1))return 1;
        p=x;
    }
    if(p!=1) return 1;
    return 0;
}
inline bool MR(ll n){
    if(n==2)return fac.push_back(n),0;
    if(n%2==0)return 1;
    if(n==3) return fac.push_back(n),0;
    if(n%6!=1&&n%6!=5)return 1;
    ll t1=n-1,t2=0;
    while(!(t1&1))t1>>=1,++t2;
    for(int i=0;i<10;i++){
        if(a[i]==n)return fac.push_back(n),0;
        if(!(n%a[i]))return 1;
        if(check(a[i],n,t1,t2))
            return 1;
    }
    return fac.push_back(n),0;
}
inline ll gcd(ll a,ll b){return b?(gcd(b,a%b)):a;}
inline ll rho(ll n,ll c){
    ll x=unit()%n+1,y=x,p=1;
    for(int k=2,i=1;p==1;i++){
        x= (ksc(x,x,n)+c)%n;
        p=(y>x)? (y-x): (x-y);
        p=gcd(p,n);
        (i==k)? (y=x,k+=k): 0;
    }
    return p;
}
inline void solve(ll n){
    if(n==1)return; 
    if(!MR(n))return;
    else{
        ll d=rho(n,unit()%n);
        while(d==n) d=rho(n,unit()%n);
        solve(n/d); solve(d);
    }
}
ll n,c;
int main(){
    scanf("%lld",&n); solve(n);
    for(int i=0;i<fac.size();i++){
        ll q=fac[i]; if(S[q])continue;
        S[q]=1; n=n/q*(q-1);
    }
    printf("%lld\n",n);
}