BZOJ 4802: 欧拉函数（大数因数分解算法 Pollard_rho 和素数测试算法 Miller_Rabin）

Description

已知N，求phi(N)

Input

正整数N。N<=10^18

Output

输出phi(N)

Sample Input

8

Sample Output

4

Solution

• 我们知道 $$\phi (n)=n(1-\frac{1}{p_1})(1-\frac{1}{p_2})\cdot \cdot \cdot (1-\frac{1}{p_k})$$

• 其中 $n=\prod p_i^{k_i}$ ，$p_i$ 为互不相同的质数。

• 很自然，我们需要将 $n$ 分解质因数，但是 $n\le 10^{18}$ ，很大。

• 于是我们可以使用 Pollard_rho 算法。

• 其主要思想是找到两个数 $p,q$ ，使得 $n=p\ast q$ ，在不断继续分解。

• 设递归过程 $proc(x)$ ，若 $x$ 是质数（用 Miller_Rabin 算法判断），则 $x$ 是 $n$ 的一个质因数，并退出过程。

• 否则找到一个数 $y$ ，使 $gcd(x,y)>1$ ，继续递归 $proc(y)$ 、$proc(\frac{x}{y})$ 。

• 那么怎样找到这个数呢？

• 我们用到了随机化的思想和生日悖论，即找到两个数 $x1,x2$ 使得 $gcd(|x1-x2|,x)>1$ ，这样概率就会大得多。

• 设变换函数 $f(x)=(x^2+c)\%m$ ，其中 $c$ 是一个随机变量，那么令 $x1,x2$ 沿着该函数一直走。

• 可这样走是会成环走回原点的，如下图：

• 它的形状形似希腊字母 $\rho$ ，该算法名即得于此。

• 判环的话可以用 Floyd 判圈算法：$$x1=f(x1)$$$$x2=f(f(x2))$$

• 即 $x2$ 以两倍速度跑，当 $x1=x2$ 时重构即可。

• 可证单 Pollard_rho 算法的期望时间复杂度为 $O(n^{\frac{1}{4}})$ 。

• 然而其中还有一个算法 Miller_Rabin （素数测试算法）还没介绍。

• 我们当然可以 $O(n)$ 线筛出素数，但是在这种时间空间都不允许的情况下，如何快速判定单个素数呢？

• 根据费马小定理，当 $n$ 为素数时，有： $$a^{n-1}\equiv 1(modn)$$

• 但是当正整数 $n$ 满足 $a^{n-1}\equiv 1(modn)$ 时不一定 $n$ 就是素数（即逆定理不成立）。

• 但是我们想如果对于很多个 $a$ 都满足上式，是否就能说明 $n$ 就是素数了呢？

• 我们发现这样错误率依然很高，这时二次探测定理就派上用场了。

• 若一个数 $x$ 满足： $x^2\equiv 1(modp)$ ，则合法的解就只有 $x=1$ 或 $x=p-1$ 。

• 证明很简单，就是因为 $x^2-1=(x+1)(x-1)$ ，合法解显然。

• 那么 $a^{n-1}$ 就可以拆成 $a^{d\ast 2^t}$ ，其中 $d$ 是奇数。

• 那么我们先用 $a^d$ 检测，之后每次平方，若 ${(a^x)}^2\equiv 1$ 则 $a^x\equiv 1或p-1(modp)$ 。

• 如此一来，检测的成功率就大了不少，更容易检测出是否为质数了。

• 我们只需随机几个 $a$ （或拿几个质数）来判断即可，复杂度约为 $lo{g}_2^3(n)$。

• 至此，我们就成功地实现了算法 Miller_Rabin 和 Pollard_rho 。

Code

#include<cstdio>
#include<algorithm>
#include<iostream>
#include<ctime>
#include<cstdlib>
using namespace std;
typedef long long LL;
typedef long double LD;
const int N=65,P=5,prime[P]={2,3,7,61,24251};
LL n,pn;
LL f[N];
inline LL mul(LL a,LL b,LL p)
{
	a%=p,b%=p;
	LL c=(LD)a*b/p;
	c=a*b-c*p;
	if(c<0) c+=p; else
		if(c>p) c-=p;
	return c;
}
inline LL ksm(LL x,LL y,LL p)
{
	LL s=1;
	while(y)
	{
		if(y&1) s=mul(s,x,p);
		x=mul(x,x,p);
		y>>=1;
	}
	return s;
}
inline LL twice(LL a,LL p)
{
	LL d=p-1;
	int t=0;
	while(!(d&1)) d>>=1,t++;//p-1=d*2^t
	LL x,y;
	x=y=ksm(a,d,p);
	while(t--)
	{
		y=mul(x,x,p);
		if(y==1 && x^1 && x^p-1) return 0;
		x=y;
	}
	return y;
}
inline LL random(LL up)
{
	return (LL)rand()*rand()%up;
}
LL gcd(LL x,LL y)
{
	return !y?x:gcd(y,x%y);
}
inline LL abs(LL x,LL y)
{
	return x<y?y-x:x-y;
}
inline bool Miller_Rabin(LL x)
{
	for(int i=0;i<P;i++)
	{
		if(x==prime[i]) return true;
		if(twice(prime[i],x)^1) return false;
	}
	return true;
}
inline LL trans(LL x,LL y,LL z)
{
	return (mul(x,x,z)+y)%z;
}
void Pollard_rho(LL m)
{
	if(Miller_Rabin(m))
	{
		f[++f[0]]=m;
		return;
	}
	LL x1=0,x2=0,c=0,p=1;
	while(p==1 || p==m)
	{
		x1=trans(x1,c,m);
		x2=trans(trans(x2,c,m),c,m);
		while(x1==x2)
		{
			c=random(m);
			x1=x2=random(m);
			x2=trans(x2,c,m);
		}
		p=gcd(abs(x1-x2),m);
	}
	Pollard_rho(p);
	Pollard_rho(m/p);
}
inline LL getphi(LL m)
{
	sort(f+1,f+1+f[0]);
	f[0]=unique(f+1,f+1+f[0])-(f+1);
	LL sum=m;
	for(int i=1;i<=f[0];i++)
		sum=sum/f[i]*(f[i]-1);
	return sum;
}
int main()
{
	scanf("%lld",&n);
	if(n==1) return 0&puts("1");
	srand(time(NULL));
	Pollard_rho(n);
	pn=getphi(n);
	printf("%lld",pn);
	return 0;
}