素数与筛法相关

一.素数定义与判定.

素数：又称质数，一个大于$1$的正整数$n$定义为素数，当且仅当不存在一个数字$i\in [2,n)$满足$i|n$.

之后我们将素数集记为$P$.

合数：一个大于$1$的正整数定义为合数，当且仅当它不是个素数.

根据素数的定义，我们可以很容易的写出一个判定素数的代码：

bool Check_prime(int n){
  if (n<=1) return 0;
  for (int i=2;i<n;++i)
    if (n%i==0) return 0;
  return 1;
}

时间复杂度$O(n)$.

显然每次都枚举$n$个数来判定有些浪费，所以考虑优化这个算法.

容易注意到一个性质，若两个小于$n$的数$x,y$相乘要等于$n$，则$x,y$中必然有一个数不超过$\sqrt{n}$，所以我们只需要枚举到$\sqrt{n}$就可以了.

时间复杂度$O(\sqrt{n})$.

代码如下：

bool Check_prime(int n){
  if (n<=1) return 0;
  for (int i=2;i*i<=n;++i)
    if (n%i==0) return 0;
  return 1;
}

二.素数个数相关.

素数个数函数：定义素数个数$\pi (n)$为：
$$\pi (n)=\sum _{i=1}^n[i\in P]$$

素数定理：$\pi (n)\approx \frac{n}{\ln n}$.

素数定理的推论1：第$n$个素数的大小约为$n\ln n$.

素数定理的推论2：设第$i$个素数为$p_i$，有：
$$O\left(\sum _{p_i\le n}\frac{1}{p_i}\right)=O(\log \log n).$$

证明：
$$\begin{gathered} O\left(\sum _{p_i\le n}\frac{1}{p_i}\right)=O\left(\sum _{i\le \pi (n)}\frac{1}{i\log i}\right) \\ =O\left({\int }_1^{\pi (n)}\frac{1}{x\log x}\mathrm{d}x\right) \\ =O(\log \log \pi (n)) \\ =O(\log \log n) \end{gathered}$$

证毕.


三.素数筛法的引入与埃式筛法.

考虑如何求出$[1,n]$中的所有素数.

根据上面的判定，我们可以很容易写出一个$O(n\sqrt{n})$的算法：

bool Check_prime(int n){
  if (n<=1) return 0;
  for (int i=2;i*i<=n;++i)
    if (n%i==0) return 0;
  return 1;
}

int n,pr[N+9],cp; 

void Sieve(){
  for (int i=2;i<=n;++i)
    if (Check_prime(i)) pr[++cp]=i;
}

但是这个复杂度好像并不能满足我们的需求，有没有更快的算法呢？

考虑不要一个一个枚举后判定，而是开一个标记数组$b$.对于每个数$i$，我们把$2i,3i,4i,\cdots$都标记为$0$，最后把所有标记为$1$的数取出就是素数.

容易发现，所有非素数$i$的倍数肯定会被$i$的某个因数筛掉，所以我们就不用枚举合数的倍数了，直接枚举所有素数的倍数筛就好了.

时间复杂度$O(n\log \log n)$.

代码如下：

int n,b[N+9],pr[N+9],cp;

void Sieve(){
  for (int i=2;i<=n;++i){
  	if (!b[i]) pr[++cp]=i;
  	for (int j=1;j<=cp&&i*pr[j]<=n;++j) b[i*pr[j]]=1;
  }
}

以上筛法被称为埃式筛法.

埃式筛法还可以在$O((\sqrt{r}+r-l)\log \log r)$的复杂度内筛出区间$[l,r]$内的素数，具体参见POJ2689 Prime Distance.


四.线性筛法.

事实上，还有一种比埃式筛更快的线性筛法，学名欧拉筛，它的时间复杂度是稳定$O(n)$的.

它的具体思路是，在埃式筛法的基础上，让每一个数都只被它的最小素因子筛掉.

具体来说，我们枚举$j$并筛掉$i{p}_j$的时候，若$p_j|i$，那么$p_{j+1}$必然不是$i{p}_{j+1}$的最小素因子，此时就可以直接退出了.

代码如下：

int n,b[N+9],pr[N+9],cp;

void Sieve(){
  for (int i=2;i<=n;++i){
  	if (!b[i]) pr[++cp]=i;
  	for (int j=1;j<=cp&&i*pr[j]<=n;++j){
  	  b[i*pr[j]]=1;
  	  if (i%pr[j]==0) break;
  	} 
  }
}

五.素因子与唯一分解定理.

素因子：定义一个数$a$是$n$的素因子，当且仅当$a\in P$且$a|n$.

唯一分解定理：对于任意一个正整数$n$，必然可以被唯一分解为：
$$n=\prod _{i=1}^k{p}_i^{c_i}$$

其中$p_i\in P$.

结论1：一个数$n$的素因子个数不超过$\log n$个，即：
$$\sum _{i=1}^k{c}_i\le \log n$$

结论2：一个数$n$本质不同的素因子个数不超过$\log \log n$个，即：
$$k\le \log \log n$$

现在我们又有了一个新问题，如何分解出$n$的所有$p_i$与对应的$c_i$.

首先一个显然的性质，必然只有最多一个素因子大于$\sqrt{n}$，所以我们可以只从$2$枚举到$\sqrt{n}$.

此时我们存一个暂存变量$now$初始为$n$，每当我们枚举到一个数$i|now$，就让$now$不停除$i$到$now$没有$i$这个因数，此时把$i$放入$p$中并把除的此时作为对应的$c$即可.

时间复杂度$O(\sqrt{n})$.

代码如下：

int d[2][C+9],cd;

void Get_d(int n){
  int now=n;
  for (int i=2;i*i<=n;++i)
    if (now%i==0)
      for (d[0][++cd]=i;now%i==0;now/=i) ++d[1][cd];
  if (now>1) d[0][++cd]=now,d[1][cd]=1;
}

六.求1到n的唯一分解.

现在我们有一个新的问题，如何求$[1,n]$中所有数的唯一分解.

我们当然可以用上面单个数的分解方法做到$O(n\sqrt{n})$分解，但是还可以用筛法做到$O(n\log n)$.

具体来说是这样的，由于在线性筛的时候每次都是用最小的素因子筛掉某个数，所以我们可以顺便记录每个数的最小素因子.有了每个数的最小素因子之后，我们就可以用一个递归或者循环来求解了.

时间复杂度$O(n\log n)$.

代码如下：

int n,b[N+9],pr[N+9],cp;
int ml[N+9];

void Sieve(){
  for (int i=2;i<=n;++i){
  	if (!b[i]) pr[++cp]=i,ml[i]=i;
  	for (int j=1;j<=cp&&i*pr[j]<=n;++j){
  	  b[i*pr[j]]=1;
  	  ml[i*pr[j]]=pr[j];
  	  if (i%pr[j]==0) break;
  	}
  }
}

int d[2][N+9][C+9],cd[N+9];

void Get_d(){ 
  for (int i=2;i<=n;++i){
  	int now=i;
  	for (;now>1;now/=ml[now])
  	  if (ml[now]==d[0][i][cd[i]]) ++d[1][i][cd[i]];
  	  else d[0][i][++cd[i]]=ml[now],d[1][i][cd[i]]=1;
  }
}

//主程序中
Sieve();
Get_d();

七.阶乘的素因子个数统计.

我们有方法在$O(\log n)$的时间复杂度内统计$n!$的某个素因子$p$的出现次数.

首先，我们先把$[1,n]$中有因数$p$的数的个数，显然为$\lfloor \frac{n}{p}\rfloor$.

然后，我们发现所有$p$的都有了$1$的贡献，但是$p^2$的倍数会有$2$的贡献，所以$p^2$的倍数的贡献应该再加上$1$，也就是$\lfloor \frac{n}{p^2}\rfloor$.

以此类推，我们可以得到答案为：
$$ans=\sum _{i=1}^{+\infty }\lfloor \frac{n}{p^i}\rfloor$$

显然计算这个式子可以做到$O(\log n)$.


八.Miller-Rabbin素性测试.

如何用快于$O(\sqrt{n})$的时间复杂度确定$n$是否是一个素数呢？

一个不幸的事实，保证正确性的做法中，最快的是$O(\sqrt{n})$的朴素和一个被称为AKS素性测试的上界至少为$O({\log }^{6}n)$的算法.

但是我们为什么一定要让这个算法一定正确？有没有什么虽然不能确保正确但正确率很高的做法呢？

考虑费马小定理的逆命题，即：
$$a^{n-1}\equiv 1(modn)\Rightarrow a\in P$$

这个命题虽然不是真命题，但它的出错概率非常小，通常随机个位数个$a$就可以保证很高的正确率.

同时有一个二次探测定理：对于一个素数$n$，若有$x$满足$x^2\equiv 1(modn)$，则必然有$x=1$或$x=n-1$.

证明：
首先$n=2$时定理显然成立.
$$x^2\equiv 1(modn)\iff n|(x^2-1)\iff n|(x+1)(x-1)\iff n|(x+1)\vee n|(x-1)$$

由于$n\in P$，所以只能取$x=1$或$x=n-1$.
证毕.

这个时候我们考虑在判定一个数$n$是否为素数的时候，先把$2$和$n$为偶数的情况判掉.之后把$n-1$中所有$2$提出来后的值记为$t$，并把$2$的个数记为$cnt$，即$n-1=2^{cnt}t$.

然后对于每一个随机值$a$，若$a^t\equiv 1(modn)$或$a^t\equiv n-1(modn)$，则判定这是个素数.

否则，我们把$cnt$个$2$乘回去，记当前值为$now$，若$no{w}^2\equiv 1(modn)$且$now\not\equiv 1(modn)$且$now\not\equiv n-1(modn)$，说明不满足二次探测定理，判定这不是个素数.否则把$now$平方继续下一次操作.

注意通常要这个算法来处理的时候，$n^2$的大小已经可以爆long long了，所以需要用到快速乘.

时间复杂度$O(k\log n)$，其中$k$表示取$a$的次数.

代码如下：

const int pr[9]={2,3,5,7,11,13,17,19,23};

LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
LL add(LL a,LL b,LL p){return a+b>=p?a+b-p:a+b;}
LL sub(LL a,LL b,LL p){return a-b<0?a-b+p:a-b;}

LL mul(LL a,LL b,LL p){
  LL t=(long double)a/p*b+1e-8,res=a*b-t*p;
  return (res<0?res+p:res)%p;
}

void sadd(LL &a,LL b,LL p){a=add(a,b,p);}
void ssub(LL &a,LL b,LL p){a=sub(a,b,p);}
void smul(LL &a,LL b,LL p){a=mul(a,b,p);}
LL Power(LL a,LL k,LL p){LL res=1;for (;k;k>>=1,smul(a,a,p)) if (k&1) smul(res,a,p);return res;}

bool Miller_rabbin(LL n){
  if (n==2) return 1;
  if (n&1^1||n<=1) return 0;
  LL t=n-1;
  int cnt=0;
  for (;t&1^1;t>>=1) ++cnt;
  for (int i=0;i<9;++i){
  	if (n==(LL)pr[i]) return 1;
  	LL a=Power(pr[i]%n,t,n);
	if (a==1||a==n-1) continue;
	for (int j=1;j<=cnt;++j){
	  LL t=mul(a,a,n);
	  if (a^1&&a^n-1&&t==1) return 0;
	  a=t;
	}
	if (a^1) return 0;
  }
  return 1;
}

九.Pollard-rho算法.

Pollard-rho算法是一种用于寻找一个数$n$的不是$1$和$n$的因数的算法.

Pollard-rho算法的基本思想是，生成一个随机数列$f_i$，然后把随机数列中相邻两项做差并取绝对值，用这些值去和$n$取gcd，只要取出的gcd大于$1$且小于$n$就可以分解出这个因数.

有一个较为优秀的伪随机函数$f_i=(f_{i-1}^2+c)modn$，其中$c$是一个随机的常数，通过这个伪随机函数有较大概率得到一个优秀的因数.

同时，由于数列是伪随机数列，所以还是会不可避免出现循环.此时我们有一个Floyd判环法，从起点出发，记录两个值$f_i,f_{2i}$，一个用一倍速跑，一个用两倍速跑，若它们相同说明出现循环，随即跳出重新选择一个随机的$c$.

不过具体实现的时候，我们可以不用真的记录两个值跑，而是在跑了$2^i$的距离后记录一下当前值$t$，直到跑到$2^{i+1}$前都与$t$比较是否出现循环.

同时通过Pollard-rho算法配合分治得到一个数$n$的所有素因子的时间复杂度是期望$O(n^{\frac{1}{4}}\log n)$的，不过是在随机数列是真随机数列的时候.

若随机数列为$f_i=(f_{i-1}^2+c)modn$，则这样实现的Pollard-rho算法复杂度还是未知的，不过一般情况下我们仍然将其看为$O(n^{\frac{1}{4}}\log n)$.

代码如下：

const int pr[9]={2,3,5,7,11,13,17,19,23};

LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
LL add(LL a,LL b,LL p){return a+b>=p?a+b-p:a+b;}
LL sub(LL a,LL b,LL p){return a-b<0?a-b+p:a-b;}

LL mul(LL a,LL b,LL p){
  LL t=(long double)a/p*b+1e-8,res=a*b-t*p;
  return (res<0?res+p:res)%p;
}

void sadd(LL &a,LL b,LL p){a=add(a,b,p);}
void ssub(LL &a,LL b,LL p){a=sub(a,b,p);}
void smul(LL &a,LL b,LL p){a=mul(a,b,p);}
LL Power(LL a,LL k,LL p){LL res=1;for (;k;k>>=1,smul(a,a,p)) if (k&1) smul(res,a,p);return res;}

bool Miller_rabbin(LL n){
  if (n==2) return 1;
  if (n&1^1||n<=1) return 0;
  LL t=n-1;
  int cnt=0;
  for (;t&1^1;t>>=1) ++cnt;
  for (int i=0;i<9;++i){
  	if (n==(LL)pr[i]) return 1;
  	LL a=Power(pr[i]%n,t,n);
	if (a==1||a==n-1) continue;
	for (int j=1;j<=cnt;++j){
	  LL t=mul(a,a,n);
	  if (a^1&&a^n-1&&t==1) return 0;
	  a=t;
	}
	if (a^1) return 0;
  }
  return 1;
}

LL Pollard_rho(LL n){
  if (Miller_rabbin(n)) return n;
  for (;2333;){
  	LL x=rand()%(n-1)+1,y=x;
  	int c=rand()%(n-1)+1,len=2;
  	for (int i=0;2333;){
  	  x=add(mul(x,x,n),c,n);
  	  LL g=gcd(abs(x-y),n);
	  if (g>1&&g<n) return g;
	  if (x==y) break;
  	  if (++i==len) y=x,len<<=1;
  	}
  }
}

set<LL>s;

void Divide_s(LL n){
  if (n==1||Miller_rabbin(n)) {if (n^1) s.insert(n);return;}
  LL t=Pollard_rho(n);
  for (;t==n;t=Pollard_rho(n)); 
  for (;n%t==0;n/=t);
  Divide_s(t);
  Divide_s(n);
}