各种友(e)善(xin)数论总集，从入门到绝望2---快速判断素数

因为原来的那篇已经很多了，所以在此写上第二篇。

这一章可以说是紧紧围绕的素数的主旨展开的。

前置芝士

快速乘

博主博主，平常$O(1)$都已经如此之快，难道可以$O(0)$？

不不不，都一样，只不过算的是$x\ast y\%z$，因为有时候$x\ast y$溢出了long long，但是结果并没有，所以发明了快速乘。

O(1)

$O(1)$版的非常简单。$x\ast y\%z=x\ast y-(x\ast y/z)\ast z$（在C++环境下），但是这个有什么特殊的吗？

就是用溢出对待溢出，我们先用long double(16位)得出(x*y/z)（你在比赛的时候，可以先用$sizeof(long$ $double)$得出你电脑的long double是几位的，如果不是$16$位的话，那评测机应该也是，为了稳妥还是用$log$的吧）。

然后我们乘一下，两边可能会溢出，但是我们还是能减出正确的结果。

另外还有一坨精度问题，模数大的话还是少用吧。

inline  LL  ksc(LL  x,LL  y,LL  z)
{
	LL  c=(LD)x*y/z+0.5;
	LL  ans=x*y-c*z;
	return  ans<0?ans+z:ans;
}

O(log)

有没有什么稳得一批又好用的快速乘？

当然后，假设又有个$x,y,z$。

我们把$x$拆成二进制：$c_n\ast 2^n+c_{n-1}\ast 2^{n-1}+....$，而$c$的取值只能为$0,1$
然后$x\ast y$在乘法分配率一下：$y\ast c_0\ast 2^0+y\ast c_1\ast 2^1+...$，那么我们只要边乘边模不就好起来了吗。

inline ll ksc(ll x,ll y,ll p){//计算x乘y的积
    ll res=0;//加法初始化
    while(y){
        if(y&1)res=(res+x)%p;//模仿二进制
        x=(x<<1)%p; y>>=1;//将x不断乘2达到二进制
    }return res;
}
// ll 表示 long long
//这里的代码用的是https://www.cnblogs.com/812-xiao-wen/p/10543023.html

floyd提出的判环法

判环

你以为是floyd，不是，是这样的，假设一个链表有环，怎么判环，我们这样想：$y$以$x$两倍的速度奔跑，那么当$x,y$相遇时，$y$刚好跑完几圈了，就退出。

找环

其实还有个扩展，如何找到环的起始位置。

我们设链表头走到环开始的地方步数为$m$，从链表头走到相遇地点的步数是$m+k$，然后环的长度为$n$，$x,y$分别走了$X,Y$圈。

那么$S_x=m+k+Xn,S_y=m+k+Yn$，然后又因为$S_y=2S_x$，所以$S_x=(Y-X)n$。
也就是说两人走的距离肯定是$n$的倍数。

然后我们再把$x$提到了链表开始的地方，让两个人继续开始走，当走了$m$步（$x$在环开始的地方），那么因为说了$S_y=2S_x=2(Y-X)n$，也就是说$y$走的步数应该是$n$的倍数，也就是说当第一次相遇的时候，他应该离走完这个环到环开始的地方为$m$，所以$y$也到了环开始的地方，所以$x,y$将会在环开始的地方相遇。

可惜这里不用。

生日悖论

首先我们来看看，现在我们来选数字，在1-100之间，如果我们选一个数字，那么是1的概率则是$\frac{1}{100}$，但是如果我们选两个数字，然后取差的绝对值，会怎样？我们选一个数字，然后选到他周围的两个数字的概率就变成了$\frac{1}{50}$了！（忽略第一个数字$1$和$100$的情况，那还是$\frac{1}{100}$），难道多元能增加概率！

没错，这就是生日悖论的内容。

生日悖论的重要思想是什么，$1-x$的范围，如果有$\sqrt{x}$个数字的话，重复的概率就会高达$50\%$，恐不恐怖。

至于证明，在这里贴上大佬的证明。

Miller_rabin

相信在第一章里面，你们已经学会了费马小定律了，那就不讲了QMQ。

前言

我们都知道判断一个数字是不是素数，有一种方法就是试除法，直接从$2$枚举到$\sqrt{p}$，但是有没有一种方法，能比$O(\sqrt{p})$还快有准确无误呢？

答案是并没有，但是如果你要求的是很大概率的话，打我可以告诉你的是，Miller_rabin就是这么一种算法，基本上准确无误，就连强伪素数都能跑过去，是什么呢？

二次探测

我们都知道，选取一个$p$数字，然后用费马小定理判断一下，如果不是$1$，那就不是素数，但是存在这么一种数字，能满足费马小定理但是不是素数的一类数字，我们又要怎么判断呢？

这里就要引入一个定理了，这个定理可以很大概率的判断是不是素数，加上费马小定理。

如果$p$是质数且$a^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)(a<p)$，那么$a=1,p-1$。
我们可以来证明一下：
$$a^2\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
$$a^2-1\equiv (\mathrm{m}\mathrm{o}\mathrm{d}p)$$
$$(a+1)(a-1)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}p)$$
那么，因为$a<p$且$p$是质数，所以$a=1,p-1$。

那么我们就可以把$p-1$分成$2^t\ast k$，然后随机选取一个值$x$，然后计算$x^k$，然后继续不断取平方：$x^{2^i\ast k}$，然后不断的用二次探测来检测，更重要的是我们最后还可以用用费马小定理，当然，$x$我们可以手动取几个素数来多判几次，不知道为什么，素数成功概率大一点QMQ。

很明显是$log$的。

代码

inline  LL  ksc(LL  x,LL  y,LL  z)
{
	LL  c=(LD)x*y/z+0.5;
	LL  ans=x*y-c*z;
	return  ans<0?ans+z:ans;
}
inline  LL  ksm(LL  x,LL  m,LL  mod)
{
	if(m==0)return  1%mod;
	LL  ans=1;
	while(m>1)
	{
		m&1?ans=ksc(ans,x,mod):0;
		x=ksc(x,x,mod);m>>=1;
	}
	return  ksc(ans,x,mod);
}
inline  int  log2(LL  &x)
{
	int  ans=0;
	while(x%2==0)ans++,x>>=1;
	return  ans;
}
int  su[]={2,3,5,7,11,23,29,61};
inline  bool  pd(LL  x)//判断一个素数 
{
	for(int  i=0;i<=7;i++)
	{
		if(x<=su[i])return  1;
		LL  y=x-1;int  tt=log2(y);
		y=ksm(su[i],y,x);
		while(tt--)
		{
			LL  z=ksc(y,y,x);
			if(z==1  &&  y!=1  &&  y!=x-1)return  0;
			y=z;
		}
		if(y!=1)return  0;
	}
	return  1;
}

Pollard-Rho

前言

你是否想快速的分解一个素数？

想吗？少年。

—来自SaDiao博主的一席话。

例题

都看到了，就是想咯QMQ。

例题

假设我们要分的数字是$p$

构建随机数列

我们构建一个随机函数$x_i=x_{i-1}\ast x_{i-1}+C(\mathrm{m}\mathrm{o}\mathrm{d}p)$（$C$为我们自己定的常数）,这么强？

同时$x_1=2$，我们发现这个序列每个都跟前面的数字有关，那么不就是类似链表了吗？而且因为是模了后序列，所以会有环（根据生日悖论，出现相同的数字概率是$O(\sqrt{p})$），就可以用$floyd$判环法了。

我们再设一个函数：$y_i=x_{i-1}\ast x_{i-1}+C(\mathrm{m}\mathrm{o}\mathrm{d}q)$（$q$为$p$的一个质因数），$y_1=2$，那么这个序列也会出现循环节的，我们再在两个循环节上找到两个位置$i,j$，使得$i<j,l(i)=l(j)$，然后我们就会发现$|x_i-x_j|$含有$q$，也就是$gcd(|x_i-x_j|,p)\ne 1$，那么我们就可以找到一个素数了。

但是我们并不知道$q$，我们又怎么找到$i,j$呢，我们会发现其实$i,j$就是$y$数列循环节上的对应位置，而上面也只是提供了一种可行性，也就是说我们可以用$floyd$找环法来在$x$数列中找，如果$gcd$为$1$，那么继续找环，如果$gcd$为$p$（差为0就会这样），说明我们找到了模数为$q$的环，可惜也是模数为$p$的环，那么我们就退出，然后$C++$，如果两个都不是，那么我们不就找到了一个质因数了吗？

再看看概率有多大，原本我们找两个数字的差找到质因数的概率应该比较小，但是如果我们是$gcd\ne 1$的话，那概率不就大了吗？而且期望的循环节大小为$\sqrt{p}$，不就好起来了吗？

而且判断一个数字是不是素数就靠Miller_rabin了。

优(ka)化(chang)

我们要发现一个事情：$x\%z$，可以等同于$z\ast ((x/gcd(x,z))\%(z/gcd(x,z)))$。

那么这有什么用呢？我们可以把几个差乘起来，然后模一下（如果出现了质因子的话不会因为模了而消失掉，上面写了），至于乘几次，我选择的是$127$，当然，$pow(p,0.1)$也有人用，不要太大就可以了，我们后面则叫乘了$step$次。

那么我们也是乘完后GCD。

1. 如果为$1$继续。
2. 如果为$n$的话，说不定中间有质因子呢，我们也回到$step$次一起，不乘起来，一次次慢慢来。
3. 如果两个都不是，恭喜，喜提质因子一枚。

然后自我感觉良好，优化了不少的常数。

注意事项

我们会发现，打完代码有事后还是会卡住的。

为什么，因为$4$是个神奇的数字，我们可以把$C$从$1$到$4$枚举一遍，会发现差统统不会涉及到$2$，而其他$2$的次方，比如$2^i$，在$C=2^i-4$的时候，肯定跳一次就能得到结果，$x=0,y=C$，然后就可以把$2$筛出来，为什么$4$不行，因为$C=0$了，所以我们就只会一遍遍的得到$4$然后重来。

那我们特判$4$不就行了？不不不，特判要讲究艺术。

你想想，我们要是特判$\%2==0$不是一样的吗，而且还造福了其他的数字，尤其是$2$的次方，不加这个很有可能就老是到$2^i-4$才跳出来。

时间复杂度证明

一次的Pollard-Rho的复杂度是多少？（这道题目得用几次Pollard-Rho）

设$N=A\ast B(A<=B)$，那么$A<=\sqrt{(}N)$，我们是在$x$找$A$的循环节，期望复杂度为$O(\sqrt{A})$，那么不就是$O(N^{\frac{1}{4}})$吗？

但是其实不然，因为加上GCD什么乱起八糟的，正宗的应该是：
$O(\frac{N^{\frac{1}{4}}logN}{step}+logN)$，但原本就是玄学算法你加这么多干嘛，而且在long long范围内我们的$step$肯定大于$logN$，毕竟我们的$step$原本就是$127$吗。

所以我们还是写$O(N^{\frac{1}{4}})$，好看又好写。

代码

开了O2在luogu跑了600+ms，快的飞起，不开也有1.44s了，这不就快的飞起了吗。

可以试试这个跑不跑得出结果46856248255981。

强伪素数呀，都跑过去了。

#include<cstdio>
#include<cstring>
#include<cmath>
using  namespace  std;
typedef  long  double  LD;
typedef  long  long  LL;
inline  LL  zabs(LL  x){return  x<0?-x:x;}
inline  LL  mymin(LL  x,LL  y){return  x<y?x:y;}
inline  LL  mymax(LL  x,LL  y){return  x>y?x:y;}
inline  LL  ksc(LL  x,LL  y,LL  z)
{
	LL  c=(LD)x*y/z+0.5;
	LL  ans=x*y-c*z;
	return  ans<0?ans+z:ans;
}
inline  LL  ksm(LL  x,LL  m,LL  mod)
{
	if(m==0)return  1%mod;
	LL  ans=1;
	while(m>1)
	{
		m&1?ans=ksc(ans,x,mod):0;
		x=ksc(x,x,mod);m>>=1;
	}
	return  ksc(ans,x,mod);
}
inline  int  log2(LL  &x)
{
	int  ans=0;
	while(x%2==0)ans++,x>>=1;
	return  ans;
}
int  su[]={2,3,5,7,11,23,29,61};
inline  bool  pd(LL  x)//判断一个素数 
{
	for(int  i=0;i<=7;i++)
	{
		if(x<=su[i])return  1;
		LL  y=x-1;int  tt=log2(y);
		y=ksm(su[i],y,x);
		while(tt--)
		{
			LL  z=ksc(y,y,x);
			if(z==1  &&  y!=1  &&  y!=x-1)return  0;
			y=z;
		}
		if(y!=1)return  0;
	}
	return  1;
}
inline  LL  gcd(LL  x,LL  y)//实测二进制版GCD只比原来的快了20+ms，估计是因为优化减少了GCD的调用次数，凸显不出优势。
{
	int  ans=0;
	while(x  &&  y)
	{
		if(x&1  &&  y&1)
		{
			y>x?x^=y^=x^=y:0;
			x=(x-y)>>1;
		}
		else  if(x&1)y>>=1;
		else  if(y&1)x>>=1;
		else  x>>=1,y>>=1,ans++;
	}
	return  (x+y)<<ans;
}
inline  LL  Pol(LL  now,LL  step,LL  add)
{
	if(now%2==0)return  2;//防止毒瘤的4的情况 
	LL  x=2,y=2,d=1;
	while(1)
	{
		LL  tx=x,ty=y;
		for(int  i=1;i<=step;i++)
		{
			x=ksc(x,x,now)+add;x>=now?x-=now:0;
			y=ksc(y,y,now)+add;y>=now?y-=now:0;
			y=ksc(y,y,now)+add;y>=now?y-=now:0;
			d=ksc(d,zabs(x-y),now);
		}
		d=gcd(d,now);
		if(d==1)continue;
		else  if(d!=now)return  d;
		x=tx;y=ty;
		for(int  i=1;i<=step;i++)
		{
			x=ksc(x,x,now)+add;x>=now?x-=now:0;
			y=ksc(y,y,now)+add;y>=now?y-=now:0;
			y=ksc(y,y,now)+add;y>=now?y-=now:0;
			d=gcd(zabs(x-y),now);
			if(d!=1)return  d%now;
		}
	}
}
inline  LL  work(LL  n)
{
	if(pd(n)  ||  n==1)return  n;
	LL  tmp=0,step=127/*玄学步数*/,add=1;
	while(!tmp)tmp=Pol(n,step,add++);
    //
	if(n/tmp<tmp)tmp=n/tmp;//使得n/tmp>=tmp
	LL  ans=work(n/tmp);
	if(ans>=tmp)return  ans;
	return  mymax(ans,work(tmp));
    //实现上的一个优化，优化空间小，但是能优化，而且不会耗多少空间，基本正优化
}
LL  n;
int  main()
{
	int  T;
	scanf("%d",&T);
	while(T--)
	{
		scanf("%lld",&n);
		LL  ans=work(n);
		if(ans==n)printf("Prime\n");
		else  printf("%lld\n",ans);
	} 
	return  0;
}