笔记：杜教筛
杜教筛
逃不过ww。被模拟赛人均爆切的莫比乌斯反演T1的1e8数据范围逼着爬来学杜教筛了。
杜教筛是用来筛积性函数前缀和的筛法。复杂度 \(O(n^{\frac 2 3})\)
前置芝士
@莫比乌斯反演 笨蛋莫反还没搞清楚就被逼迫来学。
\(f(i)\) 是积性函数，计算 \(\sum_{i=1}^n f(i)\) 。
套路式：
构造两个积性函数 \(h,g\) 使得 \(h=f*g\)
令 \(s(n)=\sum_{i=1}^n f(i)\)
\[\sum_{i=1}^n h(i)\\
=\sum_{i=1}^n \sum_{d|i} g(d)*f(\frac i d)\\
=\sum_{xd=1}^n \sum_{d|i} g(d)*f(x)\\
=\sum_{x=1}^{\frac n d} \sum_{d|i} g(d)*f(x)（顺序不太对劲）\\
= \sum_{d|i} g(d) \sum_{x=1}^{\frac n d}f(x)\\
=\sum_{d=1}^n g(d) \sum_{i=1}^{\frac n d} f(i)\\
提取 \sum_{i=1}^{\frac n d} f(i) ，\\
=\sum_{d=1}^n g(d) s(\frac n d)\\
提出左边的第一项。\\
=g(1)·s(n)+\sum_{d=2}^n g(d) s(\frac n d)\\
∴ g(1)·s(n)=\sum_{i=1}^n h(i)-\sum_{d=2}^n g(d) s(\lfloor \frac n d \rfloor)
\]
只要 \(h(i)\) 的前缀和好求，那么对之后的柿子整除分块，求s(n)的时间复杂度是\(O(n^{\frac 2 3})\)
构造 \(g,h\) 只能靠经验/kk 太草了
举例
-
求 \(s(n)=\sum_{i=1}^n \mu(i)\)
\[g(1)·s(n)=\sum_{i=1}^n h(i)-\sum_{d=2}^n g(d)s(\lfloor \frac n d \rfloor) \]寻找一个积性函数 \(g\) 使得 $f·\mu =h $ 且 \(h\) 前缀和非常好求。
\(\mu*I=\epsilon\) 显然非常好求。
\[s(n)=1-\sum_{d=2}^n s(\lfloor \frac n d \rfloor)\\ \]就这样……
-
求 \(s(n)=\sum_{i=1}^n \phi(i)\)
\[g(1)·s(n)=\sum_{i=1}^n h(i) -\sum_{d=2}^n g(d)s(\lfloor \frac n d \rfloor) \]我们有：
\[\varphi * I =id \]\(g:I,h:id\)
\[s(n)=\frac {n*(n+1)} 2-\sum_{d=2}^n s(\lfloor \frac n d\rfloor) \] -
求 \(S(n)=\sum_{i=1}^{n}i\cdot \varphi(i)\)
\[g(1)·s(n)=\sum_{i=1}^n h(i) -\sum_{d=2}^n g(d)s(\lfloor \frac n d \rfloor) \]\[ f*g=h\\ h(x)=\sum_{d|x} f(x)g(\frac x d)=\sum_{d|x}( d·\phi(d))·g(\frac x d)\\ 为了把d 搞掉，尝试让g为id\\ =\sum_{d|x} \phi(d)·x\\ =x\sum_{d|x} \phi(d)\\ =x^2 \]\[g(1)·s(n)=\sum_{i=1}^n h(i) -\sum_{d=2}^n g(d)s(\lfloor \frac n d \rfloor)\\ s(n)=\frac {n·(n+1)·(2n+1)} 6-\sum_{d=2}^n d·s(\lfloor \frac n d \rfloor) \]
三个题最后都是整除分块一下求就行了。
代码实现
map存一下杜教筛的前缀和，然后递归实现的形式很有意义（？），就这样。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N=5e6+10;
bool vis[N];
int cnt,mu[N],phi[N],p[N];
ll summu[N],sumphi[N];
map<ll,ll>mpmu;
map<ll,ull>mpphi;
int bmin(ll x,ll y){ return (x>y)?y:x; }
void init(ll n){
mu[1]=1; vis[1]=1; phi[1]=1;
for(ll i=2;i<=n;i++){
if(!vis[i]) p[++cnt]=i,mu[i]=-1,phi[i]=i-1;
for(ll j=1;j<=cnt&&p[j]*i<=n;j++){
vis[p[j]*i]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
mu[i*p[j]]=-mu[i];
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
for(ll i=1;i<=n;i++){
summu[i]=summu[i-1]+mu[i];
sumphi[i]=sumphi[i-1]+phi[i];
}
return;
}
ll getsummu(ll n){
if(n<=N-10) return summu[n];
else if(mpmu.count(n)) return mpmu[n];
//杜教筛
ll ret=1;
for(ll l=2,r;l<=n;l=r+1){
r=bmin(n,n/(n/l));
ret-=getsummu(n/l)*(r-l+1);
}
return mpmu[n]=ret;
}
ull getsumphi(ll n){
if(n<=N-10) return sumphi[n];
else if(mpphi.count(n)) return mpphi[n];
//杜教筛
ull ret=1ll*n*(n+1)/2;
for(ll l=2,r;l<=n;l=r+1){
r=bmin(n,n/(n/l));
ret-=(ull)getsumphi(n/l)*(r-l+1);
}
return mpphi[n]=ret;
}
int main(){
int t; scanf("%d",&t);
while(t--){
ll n; scanf("%lld",&n);
init(N-10);
printf("%llu %lld\n",getsumphi(n),getsummu(n));
}
return 0;
}
/*
\sum_{d=1}^{n} \mu(d) \frac {(t-1)·(t-2)} 4
mu: s(n)=1-\sum_{d=2}^n s(\lfloor \frac n d \rfloor)
*/
例题
LGP3768 简单的数学题
给定 \(n\) ，求
\[({\sum_{i=1}^n \sum_{j=1}^n {i·j·gcd(i,j)}}) \bmod p
\]
\(1\le n\le 10^{10}\)
\[{\sum_{i=1}^n \sum_{j=1}^n {i·j·gcd(i,j)}}\\
\sum_{d=1}^n \sum_{d|i} \sum_{d|j} [gcd(i,j)=d] i·j·d\\
\sum_{d=1}^n \sum_{i=1}^{\frac n d} \sum_{j=1}^{\frac n d} [gcd(i,j)=1] i·j·d^3\\
\sum_{d=1}^n \sum_{i=1}^{\frac n d} \sum_{j=1}^{\frac n d} ijd^3\sum_{t|gcd(i,j)} \mu(t) \\
\sum_{d=1}^nd^3 \sum_{t=1}^{\frac n d} \mu(t) \sum_{i=1}^{\frac n {dt}} \sum_{j=1}^{\frac n {dt}} ijt^2 \\
\sum_{d=1}^nd^3 \sum_{t=1}^{\frac n d} \mu(t) t^2calc(\frac {\frac n d} {t})^2\\
先整除分块\frac n d,然后整除分块t，\mu 那里用一下杜教筛。\\
不对，你的复杂度爆炸了。
\]
\[\sum_{d=1}^nd^3 \sum_{t=1}^{\frac n d} \mu(t) t^2calc(\frac {n} {td})^2\\
令T=td,\\
\sum_{d=1}^nd^3 \sum_{\frac T d=1}^{\frac n d} \mu(\frac T d) (\frac T d)^2calc(\frac {n} {T})^2\\
\sum_{d=1}^nd^3 \sum_{T=d}^{n} \mu(\frac T d) (\frac T d)^2calc(\frac {n} {T})^2\\
把T丢出去，\\
\sum_{T=1}^n \sum_{d|T} d^3 \mu(\frac T d) (\frac T d)^2calc(\frac {n} {T})^2\\
好了，该莫比乌斯反演了。\\
f(T)= \sum_{d|T} d^3 \mu(\frac T d) (\frac T d)^2calc(\frac {n} {T})^2\\
令 F(n)=\sum_{n|T} f(T)\\
我谔谔，所以F并不好求啊搞他干什么。\sum_{d=1}^nd^3 \sum_{t=1}^{\frac n d} \mu(t) t^2calc(\frac {n} {td})^2\\
令T=td,\\
\sum_{d=1}^nd^3 \sum_{\frac T d=1}^{\frac n d} \mu(\frac T d) (\frac T d)^2calc(\frac {n} {T})^2\\
\sum_{d=1}^nd^3 \sum_{T=d}^{n} \mu(\frac T d) (\frac T d)^2calc(\frac {n} {T})^2\\
把T丢出去，\\
\sum_{T=1}^n T^2calc(\frac {n} {T})^2\sum_{d|T} \mu(\frac T d) d\\
令f(x)=\mu(x),g(x)=x\\
\sum_{d|T} \mu(\frac T d) d=f*g=\mu*id=\varphi\\
\\
所以\\
\sum_{T=1}^ncalc(\frac {n} {T})^2 T^2 \varphi(T)\\
\sum_{T=1}^n \frac {n^2*(\frac n T +1)^2} 4 \varphi(T)\\
然后，整除分块你的 \frac n T 就行了。\phi 用杜教筛。
\]
不知道为甚么代码过不去样例，气炸了。
为了抄题解代码爬去推式子。
\[\sum_{T=1}^ncalc(\frac {n} {T})^2 T^2 \varphi(T)\\
f(T)=T^2\varphi (T)\\
f*g=h\\
g(1)S(n)=\sum_{i=1}^n h(i)-\sum_{d=2}^n g(d)s(\lfloor \frac n d \rfloor)\\
h(x)=\sum_{d|x} g(d) f(\frac x d)\\
h(x)=\sum_{d|x} g(d) \phi(\frac x d) \frac {x^2} {d^2}\\
设g(d)=d^2\\
h(x)=\sum_{d|x} \phi(\frac x d) {x^2}\\
h(x)=x^2\sum_{d|x} \phi(\frac x d)\\
\phi*I=id \ 所以h(x)=x^3？\\
g(1)S(n)=\sum_{i=1}^n h(i)-\sum_{d=2}^n g(d)s(\lfloor \frac n d \rfloor)\\
S(n)=\sum_{i=1}^n i^3 -\sum_{d=2}^n d^2 s(\frac n d)\\
整除分块\frac n d\\
\sum_{i=1}^n i^3=\frac {n(n+1)^2} 4
\]
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll N=5e6+10;
int p,cnt,sum[N],phi[N],prim[N],inv2,inv6;
bool vis[N];
map<ll,int>mp;
int power(int a,int b){
int ret=1;
while(b){
if(b&1) ret=1ll*a*ret%p;
a=1ll*a*a%p; b>>=1;
}
return ret;
}
void init(ll n){
phi[1]=1; vis[1]=1;
for(ll i=2;i<=n;i++){
if(!vis[i]) vis[i]=1,phi[i]=i-1,prim[++cnt]=i;
for(int j=1;j<=cnt&&prim[j]*i<=n;j++){
vis[prim[j]*i]=1;
if(i%prim[j]==0){
phi[i*prim[j]]=1ll*phi[i]*prim[j]%p;
break;
}
phi[i*prim[j]]=1ll*phi[i]*(prim[j]-1)%p;
}
}
for(ll i=1;i<=n;i++) sum[i]=(sum[i-1]+1ll*i*i%p*phi[i]%p)%p;
return;
}
int calc(ll x){ return x%p*(x%p+1)%p*inv2%p; }
int pfh(ll x){ return x%p*(x%p+1)%p*(2*x%p+1)%p*inv6%p; }
int getphi(ll n){
if(n<=N-10) return sum[n];
if(mp.count(n)) return mp[n];
int ret=1ll*calc(n)*calc(n)%p;
for(ll l=2,r;l<=n;l=r+1){
r=min(n,n/(n/l));
ret=(ret-1ll*(pfh(r)-pfh(l-1))*getphi(n/l)%p)%p;
}
return mp[n]=ret;
}
int main(){
ll n;
scanf("%d%lld",&p,&n);
inv6=power(6,p-2);inv2=power(2,p-2);
init(N-10);
int ans=0;
for(ll l=1,r;l<=n;l=r+1){
r=min(n,n/(n/l));
int qwq=calc(n/l);
ans=(ans+1ll*qwq*qwq%p*(getphi(r)-getphi(l-1))%p)%p;
}
printf("%d\n",(ans+p)%p);
return 0;
}
/*
998244353 2000
334905957
883968974
*/
LGP1587 [NOI2016]循环之美
给定n,m,k.
\[\sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=1\&\&\frac i j 在k进制下 可以化为纯循环小数]
\]
首先考虑在什么情况下可以化为纯循环小数。
当y是999...99的因数的时候
假设有 \(m\) 个 (k-1) 时，满足要求。\(y| (k^m-1)\)
那么计算到多少个 \(k\) 的时候是可以判断不可能的呢。。。
我们观察十进制的，一个不合法的会变成\(999...9000...0\) 显然如果 \(2|y\) 或者 \(5|y\) \(y\) 就要变成这个了。
思考一下，只要 \(y\) 与 \(k\) 互质即可。
柿子变为
\[\sum_{i=1}^n \sum_{j=1}^m [gcd(i,j)=1\&\&gcd(j,k)=1]\\
人傻了，换一下枚举顺序。\\
\sum_{j=1}^m [gcd(j,k)=1]\sum_{i=1}^n [gcd(i,j)=1]\\
先考虑后面那坨？\\
\sum_{j=1}^m [gcd(j,k)=1]\sum_{i=1}^n [gcd(i,j)=1]\\
\sum_{j=1}^m [gcd(j,k)=1]\sum_{i=1}^n \sum_{d|gcd(j,i)} \mu(d)\\
\sum_{d=1}^{min(n,m)} \mu(d) \sum_{d|j}^m[gcd(j,k)=1] \sum_{d|i}^n 1\\
\sum_{d=1}^{min(n,m)} \mu(d) \sum_{d|j}^m [gcd(j,k)=1] \frac n d\\
\sum_{d=1}^{min(n,m)} \frac n d\mu(d) \sum_{d|j}^m [gcd(j,k)=1] \\
关于k：他是一个可爱的2e3，不用担心。\\
所以你刚才tm以为n,m,k同阶拆出一个不能莫反的玩意是zz行为。\\
\sum_{d=1}^{min(n,m)} \frac n d\mu(d) \sum_{j=1}^{\frac m d} [gcd(jd,k)=1] \\
为了让它们不相关，\\
\sum_{d=1}^{min(n,m)} \frac n d\mu(d) [gcd(d,k)=1]\sum_{j=1}^{\frac m d} [gcd(j,k)=1] \\
令 F(n)=\sum_{d=1}^{min(n,m)}[gcd(d,k)=1]\frac n d\mu(d) ,\\
令 G(\frac m d)=\sum_{j=1}^{\frac m d} [gcd(j,k)=1]\\
\]
\(G:\)
\[G(m)=\sum_{j=1}^{m} [gcd(j,k)=1] \\
=\lfloor \frac m k \rfloor \varphi(k)+pre[m\%k]\\
pre[m\%k]=\sum_{i=1}^{m\%k} [gcd(i,k)=1]\\
进行O(k\log k) 预处理即可 O(1) 求出答案。
\]
\(F：\)
……原来自己推出来了一个，然后以为搞不出来就删掉了。然后其实是对的。我是zz。但是我懒得再搞。
\(\frac n d\) 去掉先。
\[F(n)=\sum_{d=1}^{min(n,m)}\mu(d) [gcd(d,k)=1] \\
=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{t|gcd(d,k)}\mu(t) \\
=\sum_{d=1}^{min(n,m)}\mu(d)\sum_{t|k}\mu(t) [t|d]\\
=\sum_{t|k}\mu(t)\sum_{d=1}^{\frac {min(n,m)} t}\mu(td)\\
=\sum_{t|k}\mu(t)\sum_{d=1}^{\frac {min(n,m)} t}\mu(t)\mu(d)[gcd(t,d)=1]\\
=\sum_{t|k}\mu(t)^2\sum_{d=1}^{\frac {min(n,m)} t}\mu(d)[gcd(t,d)=1]\\
\]
将原来的F与最后的F放在一起。
\[F(n,k)=\sum_{d=1}^{min(n,m)}\mu(d) [gcd(d,k)=1] \\
F(n,k)=\sum_{t|k}\mu(t)^2\sum_{d=1}^{\frac {min(n,m)} t}\mu(d)[gcd(d,t)=1]\\
\]
可以发现，后半部分是 \(F(\frac n t,t)\) 这是可以递归的，杜教筛F。边界是 \(F(1)\)
时间复杂度
\(G: O(1)\)，\(F: O(n^{\frac 2 3})\)
预处理 \(O(k\log k)\)
\[ans=\sum_{d=1}^n \mu(d)[gcd(d,k)=1] G(\frac m d) \lfloor \frac n d \rfloor
\]
整除分块+杜教筛mu！
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e6+10,K=2e3+10;
bool vis[N];
int cnt,p[N],mu[N],sum[N],pre[K];
unordered_map<int,int>mp,f[K];
int gcd(int a,int b){ return (b==0)?a:gcd(b,a%b); }
void init(int n,int k){
for(int i=1;i<=k;i++) pre[i]=pre[i-1]+(gcd(i,k)==1);
mu[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i]) p[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&p[j]*i<=n;j++){
vis[p[j]*i]=1;
if(i%p[j]==0) break;
mu[p[j]*i]=-mu[i];
}
}
for(int i=1;i<=n;i++) sum[i]=sum[i-1]+mu[i];
return;
}
int getmu(int n){
if(n<=N-10) return sum[n];
if(mp.count(n)) return mp[n];
ll ret=1;
for(int l=2,r;l<=n;l=r+1){
r=n/(n/l);
ret-=(r-l+1)*getmu(n/l);
}
return mp[n]=ret;
}
int getf(int n,int k){
if(f[k].count(n)) return f[k][n];
if(n==0) return 0;
if(k==1) return f[k][n]=getmu(n);
ll ret=0;
for(int i=1;i*i<=k;i++){
if(k%i!=0) continue;
if(mu[i]) ret+=mu[i]*mu[i]*getf(n/i,i);
if(mu[k/i]&&i*i!=k) ret+=mu[k/i]*mu[k/i]*getf(n/(k/i),k/i);
}
return f[k][n]=ret;
}
int g(int n,int k){
return 1ll*pre[k]*(n/k)+pre[n%k];
}
int main(){
int n,m,k;
scanf("%d%d%d",&n,&m,&k);
init(N-10,k);
ll ans=0; int tmp=min(n,m);
for(int l=1,r;l<=tmp;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=1ll*(getf(r,k)-getf(l-1,k))*g(m/l,k)*(n/l);
}
printf("%lld\n",ans);
return 0;
}
QwQwQ
