转载请注明
目录
1.前言
这一篇博客是以“数论函数(一)”为基础的。本文主要研究:
1.讨论求解因子问题时使用的筛法,包括了普通筛法,Eratosthenes筛法,欧拉筛法等,主要讨论的是线性筛法
2.利用杜教筛在低于线性时间复杂度快速求解积性函数前缀和一类的问题
2.筛法(以筛素数为例子)
2.1普通筛法一O()
#include<bits/stdc++.h>
const int N = 1e5+200;
using namespace std;
int prime[N];
void init()
{
int cnt=0;
for(int i=1;i<=N;++i)
{
int flag=1;
for(int j=2;j<=i;++j)
{
if(i%j==0)
{
flag=0;
break;
}
}
if(flag) prime[++cnt]=i;
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<prime[n]<<endl;
}
2.2.普通筛法而O()
#include<bits/stdc++.h>
const int N = 1e5+200;
using namespace std;
int prime[N];
void init()
{
int cnt=0;
for(int i=2;i<=N;++i)
{
int flag=1;
for(int j=2;j*j<=i;++j)
{
if(i%j==0) {flag==0;break;}
}
if(flag) prime[++cnt]=i;
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<prime[n]<<endl;
}
2.3Eratosthenes筛法 O()
#include<bits/stdc++.h>
const int N = 1e5+200;
using namespace std;
int prime[N];
bool vis[N];
void init()
{
int cnt=0;
for(int i=2;i<=N;++i)
{
if(!vis[i]) prime[++cnt]=i;
for(int j=2;j*i<=N;++j)
{
vis[j]=1;
}
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<prime[n]<<endl;
}
2.4欧拉筛(线性筛筛素数)O( )
#include<bits/stdc++.h>
const int N = 1e5+200;
using namespace std;
int prime[N];
bool vis[N];
void init()
{
int cnt=0;
for(int i=2;i<=N;++i)
{
if(!vis[i]) prime[++cnt]=i;
for(int j=1;i*prime[j]<=N&&j<=cnt;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
break;
}
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<prime[n]<<endl;
}
解释:
在欧拉筛中,我们必须要知道那个关键步骤:“if(i%prime[j]==0) break;”
在上面的表格中,我们发现,如果按照最小因子*剩余因子 = 合数的形式的话,就能保证无冗余。
为什么呢?
因为,若i%prime[j]==0,等价于i = k* prime[j],所以 i * prime[j+1] = k' * prime[j] ,就等于重复筛了。
由打表发现,每次用欧拉筛的时候,到n/2的地方,就已经全部筛完了,剩下的(n/2+1,n)之间都没有,所以欧拉筛的时间复杂度的下限是 3/2 n
说明:线性筛法更加快速,没有重复筛(没有冗余),那么就有线性筛法的时间复杂度趋近于线性时间复杂度。
3.线性筛
3.1筛素数(也叫做欧拉筛)
上面提到过,这里不再赘述
3.2筛欧拉函数
#include<bits/stdc++.h>
const int N = 2e7;
using namespace std;
int prime[N+20],phi[N+20],cnt;
bool vis[N+20];
void init()
{
phi[1]=1;
cnt=0;
memset(phi,0,sizeof(phi));
for(int i=2;i<=N;++i)
{
if(!vis[i])//if(!phi[i])
{
prime[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
else phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<phi[n]<<endl;
}
这里,使用了两个性质:(p是质数)
性质一:i % p ==0: phi(i * p) = phi(i) * p
性质二:i % p !=0 phi(i * p) = phi(i) * (p-1)
3.3筛莫比乌斯函数
#include<bits/stdc++.h>
const int N = 2e7;
using namespace std;
int prime[N+20],mu[N+20],cnt;
bool vis[N+20];
void init()
{
cnt=0;
mu[1]=1;
memset(vis,0,sizeof(vis));
for(int i=2;i<=N;++i)
{
if(!vis[i])
{
prime[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
else mu[i*prime[j]]=-mu[i];
}
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<mu[n]<<endl;
}
在这里,用到了三个东西来筛莫比乌斯函数:(p是质数)
一、平方数,比如2^2,3^2,...n^2,一定满足i%p=0,因此,mu(平方数)=0;
二、若 i%p==0,那么有 i=k*p,那么 i*p= p^2 * q,那么 i * p就有平方因子,所以,mu(i*p) = 0,
三、若i % p != 0, 我们分三种情况,第一:i是质数,又因为p是质数,所以mu(i * p) = (-1)^2 = 1 = - mu(p)
第二:i 含有平方根因子,又因为平方数已经被筛过了,也满足mu(i * p) = 0 = -mu(i),
第三,i 不含有平方根因子,那么再加上p这个因子,也满足mu(i * p) = -mu(p) = -mu(i)
因此,有mu(i * p) = -mu(i)
3.4筛最小质因子
#include<bits/stdc++.h>
const int N = 2e7;
using namespace std;
int prime[N+20],mnfact[N+20],cnt;
bool vis[N+20];
void init()
{
cnt=0;
mnfact[1]=1;
memset(vis,0,sizeof(vis));
for(int i=2;i<=N;++i)
{
mnfact[i]=min(mnfact[i],i);
if(!vis[i])
{
prime[++cnt]=i;
mnfact[i]=i;
}
for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
{
vis[i*prime[j]]=1;
mnfact[i*prime[j]]=prime[j];
if(i%prime[j]==0)
break;
}
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<mnfact[n]<<endl;
}
由于线性筛使用的本来就是使用最小因子这个概念,所以,很容易就写成了。
3.5筛最大质因子
#include<bits/stdc++.h>
const int N = 2e7;
using namespace std;
int prime[N+20],mxfact[N+20],cnt;
bool vis[N+20];
void init()
{
cnt=0;
memset(mxfact,0,sizeof(mxfact));
mxfact[1]=1;
for(int i=2;i<=N;++i)
{
if(!vis[i])
{
prime[++cnt]=i;
mxfact[i]=i;
}
for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
{
vis[i*prime[j]]=1;
mxfact[i*prime[j]]=max(mxfact[i],prime[j]);
if(i%prime[j]==0)
{
break;
}
}
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<mxfact[n]<<endl;
}
这里,分成两种情况:(p为质数)
一、i是质数,那么 mxfact[ i ] = i;
二、i不是质数,就有 i = k* p, mxfact[ i* p ] = max( p,mxfact[ i ] )
3.6筛不同种类质因子数目
#include<bits/stdc++.h>
using namespace std;
const int N= 2e7;
int prime[N+20],difffact[N+20],cnt;
bool vis[N+20];
void init()
{
cnt=0;
memset(difffact,0,sizeof(difffact));
difffact[1]=1;
for(int i=2;i<=N;++i)
{
if(!vis[i])
{
prime[++cnt]=i;
difffact[i]=2;
}
for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
difffact[i*prime[j]]=difffact[i];
break;
}
else
difffact[i*prime[j]]=difffact[i]+1;
}
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<difffact[n]<<endl;
}
这里,分成了三种情况:(p是质数)
一、i是质数,那么difffact[ i ] = 2
二、i不是质数,i % p = 0,那么就有i = k*p, difffact[i * p] = difffact[i](已经包含了p这个因子)
三、i不是质数,i % p != 0 ,那么就有 difffact[i * p] =difffact[ i ]+1( 加上了p这个因子)
3.7筛所有因子数目
#include<bits/stdc++.h>
using namespace std;
const int N= 2e7;
int prime[N+20],allfact[N+20],cnt,mnfact_num[N+20];
bool vis[N+20];
void init()
{
cnt=0;
memset(allfact,0,sizeof(allfact));
memset(mnfact_num,0,sizeof(mnfact_num));
allfact[1]=1;
mnfact_num[1]=1;
for(int i=2;i<=N;++i)
{
if(!vis[i])
{
prime[++cnt]=i;
allfact[i]=2;
mnfact_num[i]=1;
}
for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
allfact[i*prime[j]]=allfact[i]/(1+mnfact_num[i])*(2+mnfact_num[i]);
mnfact_num[i*prime[j]]=mnfact_num[i]+1;
break;
}
else
{
allfact[i*prime[j]]=allfact[i]*2;
mnfact_num[i*prime[j]]=1;
}
}
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<allfact[n]<<endl;
}
其中,allfact 【x】表示x的约数个数,mnfact_num【x】表示x的最小质因数的个数
这里,分成了三种情况:(p是质数)
一、i是质数,显然,allfact【i】=2,mnfact_num【i】=1
二、i不是质数,i % p = 0, 那么 有 i = k* p,假设 i = p * q, 那么这个时候就有 i * p =p^2 *q,根据 约数个数计算公式,就有
allfact【i * p】 = allfact【i】/ (1 + mnfact_num【i】) * (2 + mnfact_num【i】);(质因子为p的个数增加了一个)
而且,mnfact_num【i * p】 = mnfact_num【i * p】+1;(理由同上)
三、i不是质数,i % p != 0, 那么有 allfact【i * p】 = allfact【i 】* 2(因为增加了一个新的因子 p,就要乘上(1+p的个数)=》(1+1)=》* 2)
而且,mnfact_num【i * p】=1;
3.8筛所有因子和
#include<bits/stdc++.h>
#include<tr1/unordered_map>
#define register int rint
#define INF 0x3f3f3f3f
#define MOD 1000000007
#define mem(a,b) memset(a,b,sizeof(a))
#define PI 3.141592653589793
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int>PII;
const int N= 2e7;
const int lim=2e9;
const double esp=1e-6;
tr1::unordered_map<ll,ll>w;
template<typename T>inline void read(T &x)
{
x=0;
static int p;p=1;
static char c;c=getchar();
while(!isdigit(c)){if(c=='-')p=-1;c=getchar();}
while(isdigit(c)) {x=(x<<1)+(x<<3)+(c-48);c=getchar();}
x*=p;
}
int prime[N+20],factsum[N+20],nil_sum[N+20],cnt;
bool vis[N+20];
void init()
{
cnt=0;
factsum[1]=1;
memset(vis,0,sizeof(vis));
memset(nil_sum,0,sizeof(nil_sum));
for(int i=2;i<=N;++i)
{
if(!vis[i])
{
prime[++cnt]=i;
factsum[i]=1+i;
nil_sum[i]=1;
}
for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
factsum[i*prime[j]]=factsum[i]*prime[j]+nil_sum[i];
nil_sum[i*prime[j]]=nil_sum[i];
break;
}
else
{
factsum[i*prime[j]]=factsum[i]*(1+prime[j]);
nil_sum[i*prime[j]]=factsum[i];
}
}
}
}
int main()
{
init();
int n;
while(cin>>n) cout<<factsum[n]<<endl;
}
其中,factsum【x】表示x的约数和,nil_sum【x】表示x的约数中不能被x的最小素因子整除的约数和
这里,分成了三种情况:(p是质数)
一、i是质数,显然,factsum【i】=1+x,nil_sum【i * p】=1;
二、i不是质数,i % p =0, factsum【i】 可以写成 x + x*p +...+x*p^k , factsum【i * p】 可以表示成 x+ x *p+...+x*p^k+x*p^(k+1),
那么 factsum【i * p】 = factsum【i】* p + (factsum【i】中不包含p因子的那一部分的和),在这里,i的最小因子就是p,所以,那一部分就是 nil_sum【i】
所以 factsum【i * p】 = factsum【i】* p + nil_sum【i】
而且,nil_sum【i * p】 = nil_sum【i】(因为 i 中 已经包含了因子p,所以没有改变)
三、i不是质数,i%p != 0 ,因此,p与i互质,所以,相当于增加了一个多项式(1+p),
所以,factsum【i * p】= factsum【i】*(1+p)
而且,有 nil_sum【i * p】= factsum【i】(因为 p 是 i * p 的最小素因子,意味着要除去i * p 中含有p因子的那一个部分的约数和,就等价于factsum【i】)
3.9发现
由上这些例子,我们发现,只要是与因子有关的函数,都可以利用线性筛。
4.杜教筛
4.1线性筛和杜教筛
在上面,我们见识到了线性筛的魅力和威力,似乎有了一个很好的武器,可为什么又出来一个杜教筛呢?
首先,观察上面的筛法,发现n不会超过2e7(其实可以大一点点。那么考虑一种情况,当n<=1e10或者更大呢?
这个时候,线性时间复杂度也凉了,所以要使用低于线性时间复杂度的方法。
于是,杜教筛诞生了。
4.2杜教筛的主要功能
在低于线性时间复杂度下,求解一类积性函数前缀和。
时间复杂度的最好情况是 O()
其中,有关积性函数的一些基本概念已经在“数论函数(一)”中涉及到了,这里不再赘述
4.3杜教筛的基本套路
假设 是积性函数,设 , 设 (S(n)是我们要求的那一部分)
那么就有
=>
=>
=>
=> (最后就是求S(n),也就是到了求这一部分)
为了保证好求S(n):
1.保证 好求 (f和g的卷积和好求)(注:单位元e的前缀和为1)
2.分块+记忆化递归,以O()求出 S(n)
4.4两个注意点
4.4.1第一个注意点
若在使用杜教筛的过程中,等好的右边得到了 。那么,i 从1~n中, 有不超过 个值,前个值 分别一 一对应后 个值。
比如:设n=20,那么有:
前: 20(20/1), 10(20/2), 6(20/3), 5(20/4),
后 : 1(20/20) , 2(20/10), 3(20/6), 4(20/5),
而且,对于每一个固定的 ,i 的浮动范围是
4.4.2第二个注意点
, 这也是一个常用的形式
4.5杜教筛的应用
4.5.1求前n数的约数和(第一个注意点的情况)
#include<bits/stdc++.h>
#define register int rint
#define INF 0x3f3f3f3f3f
#define MOD 1000000007
#define mem(a,b) memset(a,b,sizeof(a))
#define PI 3.141592653589793
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int>PII;
const int N= 2e5+500;
const int Max=300;
const double esp=1e-6;
inline int rd() {
char c = getchar(); int x = 0, f = 1;
while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
return x * f;
}
ll n;
ll ans=0;
int main()
{
while(cin>>n)
{
ans=0;
ll i,q,a,b;
for(i=1;i*i<=n;++i)
{
ans=(ans+i*(n/i))%MOD;
a=n/(i+1)+1;
b=n/i;
ans=(ans+i*(a+b)*(b-a+1)/2)%MOD;
}
i--;
if(i==n/i)
ans-=i*(n/i);
cout<<ans<<endl;
}
}
4.5.2求前n个数的约数个数(第一个注意点的情况)
#include<bits/stdc++.h>
#define register int rint
#define INF 0x3f3f3f3f3f
#define MOD 1000000007
#define mem(a,b) memset(a,b,sizeof(a))
#define PI 3.141592653589793
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int>PII;
const int N= 2e5+500;
const int Max=300;
const double esp=1e-6;
inline int rd() {
char c = getchar(); int x = 0, f = 1;
while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
return x * f;
}
ll n;
ll ans=0;
int main()
{
while(cin>>n)
{
ans=0;
ll i,a,b;
for(i=1;i*i<=n;++i)
{
a=n/(1+i)+1;
b=n/i;
ans=(ans+(b-a+1)*i)%MOD;
a=n/(n/i+1)+1;
b=n/(n/i);
ans=(ans+(b-a+1)*(n/i))%MOD;
}
i--;
if(i==n/i)
{
a=n/(1+i)+1;
b=n/i;
ans-=(b-a+1)*i;
}
cout<<ans<<endl;
}
}
4.5.3求莫比乌斯函数的前缀和
#include<bits/stdc++.h>
#include<tr1/unordered_map>
typedef long long ll;
const int N = 2e7;
using namespace std;
bool vis[N+20];
int prime[N+20],mu[N+20],cnt;
tr1::unordered_map<int,int>w;
void init()
{
cnt=0;
mu[1]=1;
memset(vis,0,sizeof(vis));
for(int i=2;i<=N;++i)
{
if(!vis[i])
{
prime[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
else mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=N;++i) mu[i]+=mu[i-1];
}
int djs(ll n)
{
if(n<=N) return mu[n];
if(w[n]) return w[n];
int ans;
ll l,r;
for(l=2;l<=n;l=r+1)
{
r=n/(n/l);
ans-=(r-l+1)*djs(n/l);
}
w[n]=ans;
return ans;
}
int main()
{
init();
ll n;
while(cin>>n)
{
cout<<djs(n)<<endl;
}
}
在这里,我们根据公式
得到,
4.5.4求欧拉函数的前缀和
#include<bits/stdc++.h>
#include<tr1/unordered_map>
const int MOD=1e9+7;
typedef long long ll;
const int N = 2e7;
using namespace std;
bool vis[N+20];
int prime[N+20],phi[N+20],cnt;
int inv2=5e8 + 4;
tr1::unordered_map<int,int>w;
void init()
{
cnt=0;
phi[1]=1;
memset(vis,0,sizeof(vis));
for(int i=2;i<=N;++i)
{
if(!vis[i])
{
prime[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt&&i*prime[j]<=N;++j)
{
vis[i*prime[j]]=1;
if(i%prime[j]==0)
{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
else phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
for(int i=1;i<=N;++i) phi[i]=(phi[i]+phi[i-1])%MOD;
}
ll djs(ll n)
{
if(n<=N) return phi[n];
if(w[n]) return w[n];
ll ans=((n%MOD+1)*(n%MOD))%MOD*inv2%MOD;
ll l,r;
for(l=2;l<=n;l=r+1)
{
r=n/(n/l);
ans=ans-((r-l+1)%MOD)*djs(n/l)%MOD;
}
w[n]=ans;
return ans;
}
int main()
{
init();
ll n;
while(cin>>n)
{
cout<<djs(n)<<endl;
}
}
在这里,我们根据公式
4.5.5杜教筛的套路总结
1.找到一个好的积性函数,从而得到一个好的S(n)公式,以及便于求卷积前缀和(f 和 g的卷积前缀和)
2.采用分块 + unordered_map 进行记忆化递归
5.参考文献
1.https://blog.csdn.net/Ruger008/article/details/80245687
2.https://www.cnblogs.com/TheRoadToTheGold/p/8228969.html#_label1
3.https://blog.csdn.net/skywalkert/article/details/50500009
4.https://www.cnblogs.com/peng-ym/p/9446555.html
在此,像以上原创作者表示感谢