文章目录
埃式筛法
for(int i=2;i<=n;i++) if(!v[i])
for(int j=i*i;j<=n;j+=i) v[j]=1;
复杂度: .
线性筛
for(int i=2;i<=n;i++) {
if(!v[i])prime[++tot]=i;
for(int j=1;i*prime[j]<=n;j++) {
v[i*prime[j]=1;
if(i%prime[j]==0) break;
}
}
复杂度
.
原理:从大到小枚举质因数
杜教筛
强烈建议学习之前做一些简单的莫比乌斯反演的题.(不用到杜教筛即可).
莫比乌斯反演题表
数论函数基础(前置知识):
- 简单数论函数:
- 积性函数: ,则有 ,则为积性函数. 以上函数皆为积性函数.
- 一些性质. , 显然也是积性函数.
- 积性函数的运算: 求逆.
- 整除分块
复杂度:
).
用途:求积性函数前缀和:
已知积性函数
,求
.
处理技巧:
设
为另一个积性函数
杜教筛的重点在于找到方便的 ,使得 都是简单的积性函数.
在1s内可以跑过 的数据.
复杂度证明:
有个性质:在求
时递归到的数比
的要多.
所以要算前缀和的总共有
个.
递归的求前缀和的
分别为(可能有相等的):
递归到
的时候,前面的前缀和一定被计算过了.
可以发现,直接用
的复杂度进行合并即可.
那么总复杂度为: .
我们发现对前面的一定范围的数,直接用线性筛可以均摊
,可以不用根号复杂度,
所以我们设对前
直接预处理.
此时的复杂度为 .
此时用一下均值不等式: .
Tips:答案记忆化才能保证复杂度.下面提供一种不用
的简单做法.
贴个代码:
now为读入的值.
N=now^(2/3).
ll S(ll n) {
if(n<N) return mu[n];//直接返回
ll x=now/n,&ans=s[x];
//x为存储的位置. now/n实际上是找到原整除分块内的右端点作为标志. 这样x就在n^(1/3)内啦.
if(vis[x]) return ans;
/*vis[x]=1; ans=1;
for(ll l=2,r;l<=n;l=r+1) {
r=n/(n/l);
(ans-=calc(l,r)*S(n/l)%mod) %= mod;
}
upd(ans);*/
return ans;
}
模板题
- .
代码:
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int M=1305,N=M*M;
ll now,s1[M],s2[M];
char v1[M],v2[M],cnt;
int prime[N],tot,mu[N];
ll phi[N];
void get(int x) {
phi[1]=mu[1]=1;
for(int i=2;i<=x;i++) {
if(!phi[i]) phi[i]=i-1,mu[i]=-1,prime[++tot]=i;
for(int j=1;i*prime[j]<=x;j++)
if(i%prime[j]) {
phi[i*prime[j]]=phi[i]*(prime[j]-1);
mu[i*prime[j]]=-mu[i];
}
else {
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
}
for(int i=2;i<=x;i++)
mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
ll S1(ll n) {
if(n<N) return phi[n];
ll x=now/n,&ans=s1[x];
if(v1[x]==cnt) return ans;
v1[x]=cnt; ans=(ll)n*(n+1)>>1;
for(ll l=2,r;l<=n;l=r+1) {
r=n/(n/l);
ans-=(r-l+1)*S1(n/l);
}
return ans;
}
ll S2(ll n) {
if(n<N) return mu[n];
ll x=now/n,&ans=s2[x];
if(v2[x]==cnt) return ans;
v2[x]=cnt; ans=S1(n);
for(ll l=2,r;l<=n;l=r+1) {
r=n/(n/l);
ans-=(l+r)*(r-l+1)/2*S2(n/l);
}
return ans;
}
int main() {
get(N-1);
int T; scanf("%d",&T);
while(T--) {
scanf("%lld",&now); ++cnt;
printf("%lld %lld\n",S1(now),S2(now));
}
return 0;
}
练习:
luoguP3768
51nod 1237
提示:
51nod 1238
整除分块然后用杜教筛求 的前缀和即可.(
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int M=2222,N=M*M,mod=1000000007;
ll now,s[M];
bool vis[M];
int prime[N],tot;
ll phi[N];
ll upd(ll &x) {x+=x>>63&mod;}
void get(int x) {
phi[1]=1;
for(int i=2;i<=x;i++) {
if(!phi[i]) phi[i]=i-1,prime[++tot]=i;
for(int j=1;prime[j]*i<=x;j++)
if(i%prime[j]) phi[i*prime[j]]=phi[i]*(prime[j]-1);
else {phi[i*prime[j]]=phi[i]*prime[j]; break;}
upd(phi[i]+=phi[i-1]-mod);
}
}
ll S(ll n) {
if(n<N) return phi[n];
ll x=now/n,&ans=s[x];
if(vis[x]) return ans;
vis[x]=1; ans=n%mod; ans=ans*(ans+1)/2%mod;
for(ll l=2,r;l<=n;l=r+1) {
r=n/(n/l);
upd(ans-=(r-l+1)*S(n/l)%mod);
}
return ans;
}
ll solve() {
ll res=0;
for(ll l=1,r,t;l<=now;l=r+1) {
r=now/(t=now/l); t%=mod;
upd(res+=(S(r)-S(l-1)+mod)*t%mod*t%mod-mod);
}
return res;
}
int main() {
scanf("%lld",&now);
get(N-1);
printf("%lld\n",solve());return 0;
}
51nod 1239
51nod 1220
约数相关的函数定义 表示约数的k次方和.特别地,k=0时表示约数个数.
一些性质:
-
.
证明:设 说明对于任意质因数 中至少有一个的指数为0,那么这样就一共有 种情况( 的指数为0时, 的指数有 种情况.反之亦然.去掉一个重复的 指数都为0的情况即可),它与 的 的合法指数数量一致.
-
通过上面的证明,我们不妨把 的指数 的非0指数 映射为 .
这样我们得到一个新的式子:
.
对每个质因子分开讨论,即可发现 实际上就是在完成上面所述的映射.( )
这题的化式子特别恶心.前方高能
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int M=1010,N=M*M,mod=1000000007,mod2=2*mod;
ll now,s[M];
bool vis[M];
int prime[N],tot,f[N];
ll h[N],mu[N];
void upd(ll &x) {x+=x>>63&mod;}
void get(int x) {
h[1]=mu[1]=1;
for(int i=2;i<=x;i++) {
if(!f[i]) f[i]=i,prime[++tot]=i,h[i]=1+i,mu[i]=-1;
for(int j=1,k;(k=i*prime[j])<=x;j++) {
if(i%prime[j]==0) {
f[k]=f[i]*prime[j];
if(f[k]==k) upd(h[k]=h[i]+k-mod);
else h[k]=h[k/f[k]]*h[f[k]]%mod;
break;
}
f[k]=prime[j];
h[k]=h[i]*(prime[j]+1)%mod;
mu[k]=-mu[i];
}
}
for(int i=2;i<=x;i++)
upd(h[i]+=h[i-1]-mod),mu[i]=(mu[i]*i+mu[i-1])%mod;
}
ll calc(ll x,ll y) {
return (x+y)%mod2*((y-x+1)%mod2)/2%mod;
}
ll S(ll n) {
if(n<N) return mu[n];
ll x=now/n,&ans=s[x];
if(vis[x]) return ans;
vis[x]=1; ans=1;
for(ll l=2,r;l<=n;l=r+1) {
r=n/(n/l);
(ans-=calc(l,r)*S(n/l)%mod) %= mod;
}
upd(ans);
return ans;
}
ll g(ll n) {
if(n<N) return h[n];
ll ans=0;
for(ll l=1,r;l<=n;l=++r) {
r=n/(n/l);
(ans+=calc(l,r)*(n/l)%mod) %= mod;
}
upd(ans);
return ans;
}
ll solve() {
ll ans=0,n=now;
for(ll l=1,r;l<=n;l=++r) {
r=(n/(n/l));
ll t=g(n/l);
t=t*t%mod;
(ans+=(S(r)-S(l-1)+mod)*t%mod) %=mod;
}
return ans;
}
int main() {
scanf("%lld",&now);
get(min(now,(ll)N-1));
printf("%lld\n",solve());
return 0;
}
BZOJ #3512. DZY Loves Math IV
题意简洁,推导毒瘤:
给定
.
与普通的反演不同的是,这题竟然是对每个 ,暴力求解.(其实也不是很暴力,就是和往常做法不同)
设
一波操作猛如虎,智商顿增250
所以分治处理,
记忆化配合即可.
边界:
,杜教筛求
前缀和即可.
复杂度不详,望读者给我点启发!!!
#include<map>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=1.1e6+10,mod=1e9+7;
int n,m,prime[N/10],tot,low[N],f[N],phi[N];
void get() {
low[1]=f[1]=phi[1]=1;
for(int i=2;i<N;i++) {
if(!low[i]) low[i]=i,f[i]=i-1,prime[++tot]=i;
for(int j=1,k;(k=i*prime[j])<N;j++)
if(i%prime[j]) {
low[k]=low[i]*prime[j];
f[k]=f[i]*(prime[j]-1);
}
else {
low[k]=low[i];
f[k]=f[i]*prime[j];
break;
}
phi[i]=f[i];
f[i]=(f[i]+f[i-1])%mod;
}
}
map<int,int>s;
ll S(ll n) {
if(n<N) return f[n];
if(s.count(n)) return s[n];
ll ans=n*(n+1)/2%mod;
for(ll l=2,r;l<=n;l=++r) {
r=n/(n/l);
ans -= (r-l+1) * S(n/l) % mod;
}
return s[n]=(ans%mod+mod)%mod;
}
map<int,int> ans[N];
ll C(int n,int m) {
if(!m) return 0;
if(ans[n].count(m)) return ans[n][m];
if(n==1) return ans[n][m]=S(m);
ll s=0;
for(int i=1;i*i<=n;i++)
if(n%i==0) {
s+=phi[n/i]*C(i,m/i)%mod;
if(i*i!=n) s+=phi[i]*C(n/i,m/(n/i))%mod;
}
return ans[n][m]=s%mod;
}
int main() {
scanf("%d%d",&n,&m); get();
ll sum=0;
for(int i=1;i<=n;i++)
sum+=i/low[i]*C(low[i],m)%mod;
printf("%lld\n",sum%mod);return 0;
}
BZOJ #3930. [CQOI2015]选数
题意:在 的方案数
方法1:直接套莫反公式,复杂度:
#include<map>
#include<set>
#include<queue>
#include<cmath>
#include<cstdio>
#include<vector>
#include<string>
#include<cstring>
#include<iostream>
#include<algorithm>
#define lc (x<<1)
#define rc (x<<1|1)
#define gc getchar()//(p1==p2&&(p2=(p1=buf)+fread(buf,1,size,stdin),p1==p2)?EOF:*p1++)
#define mk make_pair
#define pi pair<int,int>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N=2e6+10,size=1<<20,mod=1000000007;
//char buf[size],*p1=buf,*p2=buf;
template<class o> void qr(o &x) {
char c=gc; x=0; int f=1;
while(!isdigit(c)){if(c=='-')f=-1; c=gc;}
while(isdigit(c)) x=x*10+c-'0',c=gc;
x*=f;
}
template<class o> void qw(o x) {
if(x/10) qw(x/10);
putchar(x%10+'0');
}
template<class o> void pr1(o x) {
if(x<0)x=-x,putchar('-');
qw(x); putchar(' ');
}
template<class o> void pr2(o x) {
if(x<0)x=-x,putchar('-');
qw(x); puts("");
}
int n,m,l,r,prime[N/10],tot,mu[N]; bool v[N];
void upd(int &x) {x+=x>>31&mod;}
void upd(ll &x) {x+=x>>63&mod;}
void get(int x) {
mu[1]=1;
for(int i=2;i<=x;i++) {
if(!v[i]) prime[++tot]=i,mu[i]=-1;
for(int j=1,k;(k=i*prime[j])<=x;j++) {
v[k]=1;
if(i%prime[j]) mu[k]=-mu[i];
else break;
}
upd(mu[i]+=mu[i-1]);
}
}
map<int,int>s;
ll S(int n) {
if(n<N) return mu[n];
if(s.count(n)) return s[n];
ll ans=1;
for(ll l=2,r;l<=n;l=++r) {
r=n/(n/l);
ans-=(r-l+1)*S(n/l);
}
upd(ans%=mod);
return s[n]=ans;
}
ll power(ll a,ll b=n) {
ll c=1;
while(b) {
if(b&1) c=c*a%mod;
b /= 2; a=a*a%mod;
}
return c;
}
int main() {
qr(n); qr(m); qr(l); qr(r); l=(l-1)/m; r /= m;
get(min(r,N-1));ll ans=0;
for(int i=1,j;i<=r;i=++j) {
j=min(r/(r/i),l>=i?l/(l/i):(int)1e9);
ans+=(S(j)-S(i-1)+mod)*power(r/i-l/i)%mod;
}
upd(ans%=mod);
pr2(ans);
return 0;
}
方法2:
首先, , ,之后题目转化为选择 中的数使得 .
方法1没有利用到 这一重要性质,其实直接容斥更快.
假如
个数中有任意两个不同,由更相减损法得
.
我们一开始计算的时候先忽略掉这个相等的影响.
我们容易求得 的总方案,然后把等于 的减去即可.
#include<cstdio>
#define ll long long
using namespace std;
const int N=1e5+10,mod=1000000007;
int n,m,l,r,f[N],len;
ll power(ll a,ll b=n) {
ll c=1;
while(b) {
if(b&1) c=c*a%mod;
b /= 2; a=a*a%mod;
}
return c;
}
void upd(int &x) {x+=x>>31&mod; x-=(x>=mod?mod:0);}
int main() {
scanf("%d %d %d %d",&n,&m,&l,&r); l=(l-1)/m; r /= m;
for(int i=(len=r-l),x,y,last=-1,val; i;i--) {
x=r/i; y=l/i;
if(last!=x-y) last=x-y,val=power(last);//这一句的总共复杂度为sqrt(n)*log(n)
upd(f[i]=val-last);//减去全部相等的方案.
for(int j=i<<1;j<=len;j+=i) upd(f[i]-=f[j]);//O(len log(len))
}
printf("%d\n",(f[1]+(!l))%mod);return 0;//最后当然允许相等
}