素数&筛法
素数
素数与合数
•若一个大于一的正整数P，其约数只有1和P本身，称其为素数(质数)。若其有超过两个约数，则为合数。
素数无限定理
•正整数集中包含无限个素数。
•反证：假设素数有限，设其为\(p_1~p_n\)，构造\(s=1+π(pi)\)，若s是素数，矛盾；若s是合数，则\((p_1-p_n)\)都不是\(s\)的约数，与算术基本定理矛盾。
•特别的：1~n中大约有\(O(n/ln(n))\)个素数
互质
•若两个数A,B其最大公约数为1，则称A,B互质。
•此时k*A(0<=k<B)会遍历整个mod B剩余系。
•若A,B最大公约数为g，则k*A只能遍历mod B剩余系中g的倍数部分。
筛法
朴素筛——\(O(\sqrt n)\)
注意循环里千万不能有sqrt
for(int i=2;i*i<=n;i++)
if(n%i==0) {prime[++prime[0]]=i;while(n%i==0)n/=i;}
if(n>1)prime[++prime[0]]=n;
Eratosthenes 埃氏筛
这个不常用
void Era(int n) {
for(int i=0;i<=n;i++) is_prime[i]=1;
is_prime[0]=is_prime[1]=0;
for(int i=2;i<=n;i++) {
if(is_prime[i]) {
p[++p[0]]=i;
if((ll)i*i<=n) {
for(int j=i*i;j<=n;j+=i)
is_prime[j]=0;
}
}
}
}
欧拉线性筛——\(O(n)\)
•维护每个数的最小素因子mindiv，从2开始，若其mindiv未知，则为素数，将其小于本身的素数倍的mindiv设置成素数。
•每个数只会被其最小素因子筛一次。
#include <iostream>
#include <cstdio>
using namespace std;
const int N=100000001;
int n;
int mindiv[N],prime[N];
int main(){
scanf("%d",&n);
for(int i=2;i<=n;i++){
if(!mindiv[i]) mindiv[i]=prime[++prime[0]]=i;
for(int j=1;j<=prime[0]&&i*prime[j]<=n&&prime[j]<=mindiv[i];j++)
mindiv[i*prime[j]]=prime[j];
}
for(int i=1;i<=prime[0];i++)
printf("%d\n",prime[i]);
return 0;
}
线性筛还可以筛 φ ，d (约数个数), sd(约数和)
void PHI(int x){
phi[1]=1;
for(int i=2;i<=x;i++){
if(!mindiv[i]) mindiv[i]=p[++p[0]]=i,phi[i]=i-1;
for(int j=1;j<=p[0]&&i*p[j]<=x;j++){
mindiv[i*p[j]]=p[j];
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}else phi[i*p[j]]=phi[i]*phi[p[j]];
}
}
}
d和sd的证明见这篇博客
#include<cstdio>
#include<algorithm>
typedef long long ll;
inline int gi(){
int x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9')f=ch=='-'?-1:f,ch=getchar();
while(ch>='0'&&ch<='9')x=x*10+ch-'0',ch=getchar();
return x*f;
}
const int N=1000001;
int prime[N],d[N],f[N];
bool mindiv[N];
int main(){
int n=gi();
for(int i=1;i<=n;++i) d[i]=1;
for(int i=2;i<=n;++i){
if(!mindiv[i]) prime[++prime[0]]=i,d[i]=2,f[i]=1;
for(int j=1;i*prime[j]<=n&&j<=prime[0];++j){
mindiv[i*prime[j]]=1;
if(i%prime[j]==0){
f[i*prime[j]]=f[i]+1;
d[i*prime[j]]=d[i]/(f[i]+1)*(f[i*prime[j]]+1);
break;
}
f[i*prime[j]]=1;
d[i*prime[j]]=d[i]*d[prime[j]];
}
}
ll ans=0;
for(int i=1;i<=n;++i) ans+=d[i];
printf("%lld\n",ans);
return 0;
}
sd[1]=1;
for(int i=2;i<=n;++i){
if(!mindiv[i]) prime[++prime[0]]=i,sd[i]=sp[i]=i+1;
for(int j=1;i*prime[j]<=n&&j<=prime[0];++j){
mindiv[i*prime[j]]=1;
if(i%prime[j]==0){
sp[i*prime[j]]=sp[i]*prime[j]+1;
sd[i*prime[j]]=sd[i]/sp[i]*sp[i*prime[j]];
break;
}
sp[i*prime[j]]=1+prime[j];
sd[i*prime[j]]=sd[i]*sd[prime[j]];
}
}
例1：计算区间中素数的个数
(L≤R≤2147483647,R-L≤1000000)
上面复杂度最优的欧拉筛也只能40。。。还是只用了一遍加二分查找
正解：
筛出1-50000中的所有质数,并且对合数打上标记.
在L-R的范围呢用我们已求出的质数筛出其中的合数(设p为质数,则i*p一定不为质数),并对其打上标记(埃氏筛的思想)
遍历L-R,没有打标记的元素即为我们所求的素数
#include <iostream>
#include <cstdio>
using namespace std;
const int maxn=1e7+1;
#define int long long
int L,R;
int prime[maxn],cnt,ans;
int mindiv[maxn];
bool vis[maxn];
void Gprime(){
for(int i=2;i<=50000;i++){
if(!mindiv[i]) mindiv[i]=prime[++prime[0]]=i;
for(int j=1;j<=prime[0]&&i*prime[j]<=50000&&prime[j]<=mindiv[i];j++)
mindiv[i*prime[j]]=prime[j];
}
}
signed main(){
scanf("%lld%lld",&L,&R);
L=L==1?2:L;
Gprime();
for(int i=1;i<=prime[0];++i){
int p=prime[i],start=(L+p-1)/p*p > 2*p ? (L+p-1)/p*p : 2*p;//注意最小为2p
for(int j=start;j<=R;j+=p)vis[j-L+1]=1;
}
for(int i=1;i<=R-L+1;++i)if(!vis[i])ans++;
printf("%lld\n",ans);
}
例2：Prime Distance
Description
•给定L,R<=2147483647,L-R<1e6 求区间\([L,R]\)中所有相邻的素数中差最小与最大的素数对
Solution
(和上一题完全一样的思路....）要找到 \([l, r]\) 中相邻最近和最远的素数对肯定是需要找出 \([l, r]\) 内所有素数 . 但是无论是直接线性打表还是暴力都处理不了int的数据.素数 p 只对 \(>=p^2\)的数有贡献，可以先给 \(\sqrt r\) 内的素数打个表, 再用 \(\sqrt r\) 内的素数去筛选 \([l, r]\) 内的合数, （模拟埃氏筛）然后再遍历一次 [l, r], 记录答案即可；复杂度 : $\sum_{p<=\sqrt(r)} \frac{R-L}{p} $ p是素数，
例3：质数生成器
Description
要求生成质数的范围是\([m,n]\) (\(1 <= m <= n <= 10^9, n-m<=1000000\))
Solution
与 1 思路完全一样
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#define N 10000003
using namespace std;
int pd[N],prime[N],T,l,r,vis[N];//prime[0]
void init(){
pd[1]=1;
for(int i=2;i<=100000;i++){
if(!pd[i])prime[++prime[0]]=i;
for(int j=1;j<=prime[0];j++){
if(i*prime[j]>100000)break;
pd[prime[j]*i]=1;
if(i%prime[j]==0)break;
}
}
}
int main(){
init();
scanf("%d",&T);
for(int i=1;i<=T;i++){
scanf("%d%d",&l,&r);
for(int j=1;j<=prime[0];j++){
int t=(l-1)/prime[j]+1;
if(t<=1)t=2;
for(int k=t*prime[j];k<=r;k+=prime[j])
vis[k-l+1]=i;
}
for(int j=1;j<=r-l+1;j++)
if(vis[j]!=i && j+l-1!=1)printf("%d\n",j+l-1);
puts("");
}
return 0;
}
Miller-rabin素数测试
https://www.cnblogs.com/clockwhite/p/12148877.html
Pollard-Rho算法
https://www.luogu.com.cn/blog/JasonZhou200337/pollard-rho-suan-fa-jian-jie
dingding 10-1 18:30
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
#define il inline
typedef long long LL;
inline LL read() {
LL x=0;char ch=getchar();
while(!isdigit(ch))ch=getchar();
while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
return x;
}
il LL qmul(LL a,LL b,LL md) {
LL c=(long double)a*b/md+0.5;
c=a*b-c*md;
return c<0?c+md:c;
}
il LL qpow(LL a,LL b,LL md) {
LL r=1;
for(;b;b>>=1,a=qmul(a,a,md)) if(b&1) r=qmul(r,a,md);
return r;
}
il bool check(LL a,LL p,LL md) {
LL t=qpow(a,p,md);
if(t==md-1) return 1;
if(t==1) return p&1?1:check(a,p/2,md);
return 0;
}
il bool Miller_Rabin(LL n) {
if(n<=1) return 0;
if(n<=3) return 1;
if((n&1)==0||n%3==0) return 0;
static int p[]={2,3,5,7,11,13,17,19};
for(int i=0;i<8;++i)
if(n==p[i]) return 1;
else if(!check(p[i],n-1,n)) return 0;
return 1;
}
il LL gcd(LL a,LL b) {
return b?gcd(b,a%b):a;
}
LL rho(LL n) {
LL a=((LL)rand()<<30^rand())%n;
LL x=1,y=0,p,t;
for(int k=1;;y=x,k<<=1) {
p=1;
for(int i=0;i<k;++i) {
x=qmul(x,x,n)+a,x>=n?x-=n:0;
p=qmul(p,abs(x-y),n);
}
if((t=gcd(p,n))>1) return t;
}
}
LL mxdiv;
void divi(LL n) {
if(n==1) return;
if(Miller_Rabin(n)) {
mxdiv=max(mxdiv,n);
return;
}
LL t;
do t=rho(n); while(t==n);
while(n%t==0) n/=t;
divi(n),divi(t);
}
int main() {
int T; LL n; T=read();
while(T--) {
n=read();
if(Miller_Rabin(n)) {puts("Prime"); continue;}
mxdiv=0;
divi(n);
printf("%lld\n",mxdiv);
}
return 0;
}
