Magic的Miller Rabin素数测定和Pollard's Rho质因子分解法

boshi大佬帮助了litble，使愚蠢的litble（可能？）学会了这种神奇的算法。
orz boshi。

例题:poj1811

Miller Rabin

又称米勒挝饼， 一种神奇的快速判定一个数是否是素数的方法。如果该方法判定出一个数是合数，那么该数一定是合数。如果判定出是素数，那该数仍有极小的概率是合数。

费马小定理：若$p$是一个质数，则有$a^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$
二次探测原理：若$p$是一个质数，有$a^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$，那么$a$只能是$1$或者$p-1$。

$p$不是质数情况下，同时满足这两个条件的值大概有$\frac{1}{4}$，所以多取几个随机值检验，出错的概率就很低了。

算法过程：

typedef long long LL;
LL qmul(LL x,LL y,LL p) {//防止取模爆炸的快速乘
	x%=p,y%=p;LL re=0;
	for(LL i=y;i;i>>=1,x=(x+x)%p) if(i&1) re=(re+x)%p;
	return re;
}
LL qpow(LL x,LL y,LL p) {//快速幂
	x%=p;LL re=1;
	for(LL i=y;i;i>>=1,x=qmul(x,x,p)) if(i&1) re=qmul(re,x,p);
	return re;
}
int Miller_Rabin(LL n) {
	if(n==2||n==3||n==5||n==7||n==11||n==13) return 1;
	if(n==1||n%2==0||n%3==0||n%5==0||n%7==0||n%11==0||n%13==0) return 0;//先用几个小质因子测试
	LL t=n-1,k=0;
	while(!(t&1)) t>>=1,++k;//n-1=t*2^k
	for(RI kas=1;kas<=10;++kas) {
		LL x=qpow(rand()%(n-2)+2,t,n),y;
		for(LL i=1;i<=k;++i) {
			y=x,x=qmul(x,x,n);//每次平方
			if(x==1&&y!=1&&y!=n-1) return 0;//用二次探测原理检验
		}
		if(x!=1) return 0;//用费马小定理检验
	}
	return 1;
}

Pollard’s Rho

又称泼辣的肉， 一种神奇地将极大数进行质因数分解的算法。

我们知道，朴素的质因数分解是$O(\sqrt{n})$的，这样如果$n$出到$10^{18}$就令人束手无策了。

于是我们希望利用随机化，更快地寻找所有质因数。

说道随机化，我们脑海里可能会有一个最简单最暴力的思路：每次随机一个数，判断该数是不是$n$的质因子。

这样随机到的概率很小，但是根据生日悖论，$gcd(a-b,n)$是$n$的质因子的概率就大得多了。

我们利用一种比较高效的伪随机数来实现这个随机的$a$和$b$，一般这个伪随机函数为$f(x)=x^2+c$，$c$是一个根据你的兴趣决定的数，然后每次生成下一个随机数，不断计算。如果$d=gcd(a-b,n)$不等于$1$或$n$，说明找到了一个$n$的约数。如果这个约数是质因数，那么皆大欢喜，否则递归处理$\frac{n}{d}$和$d$即可。

判断是否是质因数就用米勒挝饼好了，肉夹馍？ 整个算法已经比较完整，唯一的问题是该伪随机数可能陷入死循环。伪随机数因为每次下一步都是固定的，所以构成基环内向树。

判环的思想也比较简单，设想一只兔子和一只乌龟在漫漫长路上前进，兔子如果追上了乌龟，就说明存在一个环。那么每次兔子走几步，乌龟走一步即可。

不过呢，实际检测，这个“兔子走几步”的“几”在每次移动乌龟后倍增，效率比较高。
嗯，代码如下：

int js;LL pri[100];
LL gcd(LL a,LL b) {return b?gcd(b,a%b):a;}
LL Pollard_Rho(LL n) {
	LL x=rand()%n,y=x,c=rand()%n,QvQ=2;int i=1;
	while(1) {
		++i,x=(qmul(x,x,n)+c)%n;//生成伪随机数
		if(y==x) return 1;//兔子追上了乌龟
		LL d=gcd((y-x+n)%n,n);
		if(d>1) return d;//注意，有可能在某个c意义下，d一直都是n
		if(i==QvQ) y=x,QvQ<<=1;//移动乌龟
	}
}
void getpri(LL n) {
	if(n==1) return;
	if(Miller_Rabin(n)) {pri[++js]=n;return;}//如果是质因子
	LL p=n;while(p==n) p=Pollard_Rho(n);
	getpri(n/p),getpri(p);//递归处理
}