【模板】【系列】 杜教筛
莫比乌斯反演的前置知识：https://www.cnblogs.com/yyc-jack-0920/p/10556351.html
杜教筛可以在$O(n^{2/3})$的时间内求出积性函数的前缀和即$\sum\limits_{i=1}^n f(i)$
实现过程需要找到另一个积性函数$g$使得$f*g$为一个方便计算的函数
$\sum f*g=\sum\limits_{i=1}^n \sum\limits_{d|i} g(d)f(i/d) $
$=\sum\limits_{d=1}^n g(d) \sum\limits_{t=1}^{n/d} f(t)$
$=g(1)\sum\limits_{i=1}^n f(i) + \sum\limits_{d=2}^n g(d) \sum\limits_{i=1}^{n/d} f(i)$
设前缀和为$S(i)$ 即$g(1)S(n)=-\sum\limits_{d=2}^n g(d) S(n/d)+\sum f*g$
这样就写成了一个数论分块的递归形式，复杂度可以证明（
同时，由于杜教筛自带数论分块形式，因此在外面再套一层数论分块并不会影响复杂度
luogu 4213 【模板】杜教筛：
求$\varphi$与$\mu$的前缀和
有结论：$\varphi*1=id$、$\mu *1 =[n=1]$ 直接套用公式即可
1 #include<bits/stdc++.h> 2 #define inf 2139062143 3 #define MAXN 5001001 4 #define MOD 1000000007 5 #define ll long long 6 #define ull unsigned long long 7 #define rep(i,s,t) for(register int i=(s),i##end=(t);i<=i##end;++i) 8 #define dwn(i,s,t) for(register int i=(s),i##end=(t);i>=i##end;--i) 9 #define ren for(register int i=fst[x];i;i=nxt[i]) 10 #define Fill(a,b) memset(a,b,sizeof(a)) 11 #define pls(a,b) ((a+b)%MOD+MOD)%MOD 12 #define mns(a,b) (((a)-(b))%MOD+MOD)%MOD 13 #define mul(a,b) (1LL*(a)*(b)%MOD) 14 #define pii pair<int,int> 15 #define mp(a,b) make_pair(a,b) 16 #define pb(i,x) vec[i].push_back(x) 17 #define fi first 18 #define se second 19 using namespace std; 20 inline ll read() 21 { 22 ll x=0;bool f=1;char ch=getchar(); 23 while(!isdigit(ch)) {if(ch=='-') f=0;ch=getchar();} 24 while(isdigit(ch)) x=x*10+(ch^48),ch=getchar(); 25 return f?x:-x; 26 } 27 int lmt,ntp[MAXN],p[MAXN],tot;ll mu[MAXN],phi[MAXN]; 28 void mem(int n) 29 { 30 mu[1]=phi[1]=1;rep(i,2,n) 31 { 32 if(!ntp[i]) p[++tot]=i,mu[i]=-1,phi[i]=i-1; 33 rep(j,1,tot) if(1LL*i*p[j]>n) break; 34 else 35 { 36 ntp[i*p[j]]=1;if(i%p[j]) mu[i*p[j]]=-mu[i],phi[i*p[j]]=phi[i]*phi[p[j]]; 37 else {phi[i*p[j]]=phi[i]*p[j];break;} 38 } 39 } 40 rep(i,1,n) mu[i]+=mu[i-1],phi[i]+=phi[i-1]; 41 } 42 unordered_map<int,ll> ansm,ansp; 43 ll smu(ll n,ll res=0) 44 { 45 if(n<=lmt) return mu[n];if(ansm[n]) return ansm[n]; 46 for(ll l=2,r;l<=n;l=r+1) r=n/(n/l),res+=smu(n/l)*(r-l+1); 47 return ansm[n]=1LL-res; 48 } 49 ll sphi(ll n,ll res=0) 50 { 51 if(n<=lmt) return phi[n];if(ansp[n]) return ansp[n]; 52 for(ll l=2,r;l<=n;l=r+1) r=n/(n/l),res+=sphi(n/l)*(r-l+1); 53 return ansp[n]=1LL*n*(n+1)/2-res; 54 } 55 int main() 56 { 57 int T=read(),n;mem(lmt=5000000); 58 while(T--) {n=read();printf("%lld %lld\n",sphi(n),smu(n));} 59 }
51nod 1237 最大公约数之和V3
题目大意：
求$\sum\limits_{i=1}^n \sum\limits_{j=1}^{n} gcd(i,j)$
思路：
原式=$\sum\limits_{d=1}^{n} d \sum\limits_{i=1}^{n/d} \sum\limits_{j=1}^{n/d} [gcd(i,j)==1]$
$\sum\limits_{d=1}^n d \sum\limits_{i=1}^{n/d} \mu(i) \lfloor\frac{n}{id}\rfloor \lfloor\frac{n}{id}\rfloor$
$\sum\limits_{g=1}^n \lfloor\frac{n}{g}\rfloor \lfloor\frac{n}{g}\rfloor \sum\limits_{d|g} \mu(d)*\frac{g}{d}$
$\sum\limits_{g=1}^n \lfloor\frac{n}{g}\rfloor \lfloor\frac{n}{g}\rfloor \cdot \varphi(g)$
数论分块+杜教筛求欧拉函数前缀和即可
1 #include<bits/stdc++.h> 2 #define inf 2139062143 3 #define MAXN 5001001 4 #define MOD 1000000007 5 #define ll long long 6 #define ull unsigned long long 7 #define rep(i,s,t) for(register int i=(s),i##end=(t);i<=i##end;++i) 8 #define dwn(i,s,t) for(register int i=(s),i##end=(t);i>=i##end;--i) 9 #define ren for(register int i=fst[x];i;i=nxt[i]) 10 #define Fill(a,b) memset(a,b,sizeof(a)) 11 #define pls(a,b) (a+b)%MOD 12 #define mns(a,b) (a-(b)+MOD)%MOD 13 #define mul(a,b) (1LL*(a)*(b)%MOD) 14 #define pii pair<int,int> 15 #define mp(a,b) make_pair(a,b) 16 #define pb(i,x) vec[i].push_back(x) 17 #define fi first 18 #define se second 19 using namespace std; 20 inline ll read() 21 { 22 ll x=0;bool f=1;char ch=getchar(); 23 while(!isdigit(ch)) {if(ch=='-') f=0;ch=getchar();} 24 while(isdigit(ch)) x=x*10+(ch^48),ch=getchar(); 25 return f?x:-x; 26 } 27 int ntp[MAXN],p[MAXN],tot;ll phi[MAXN],lmt,n; 28 const int inv2=500000004; 29 void mem(int n) 30 { 31 phi[1]=1;rep(i,2,n) 32 { 33 if(!ntp[i]) p[++tot]=i,phi[i]=i-1; 34 rep(j,1,tot) if(1LL*i*p[j]>n) break; 35 else 36 { 37 ntp[i*p[j]]=1;if(i%p[j]) phi[i*p[j]]=mul(phi[i],phi[p[j]]); 38 else {phi[i*p[j]]=mul(phi[i],p[j]);break;} 39 } 40 } 41 rep(i,1,n) phi[i]=pls(phi[i-1],phi[i]); 42 } 43 map<ll,int> ansp; 44 ll sum(ll n) 45 { 46 n%=MOD;return mul(mul(n,n+1),inv2); 47 } 48 ll sphi(ll n,ll res=0) 49 { 50 if(n<=lmt) return phi[n];if(ansp[n]) return ansp[n]; 51 for(ll l=2,r;l<=n;l=r+1) r=n/(n/l),res=pls(mul(sphi(n/l),(r-l+1)%MOD),res); 52 return ansp[n]=mns(sum(n),res); 53 } 54 ll solve(ll n,ll res=0) 55 { 56 for(ll l=1,r;l<=n;l=r+1) 57 r=n/(n/l),res=pls(res,mul(mul((n/l)%MOD,(n/l)%MOD),mns(sphi(r),sphi(l-1)))); 58 return res; 59 } 60 int main() 61 { 62 n=read();mem(lmt=5000000);printf("%lld\n",solve(n)); 63 }
51nod 1238 最大公倍数之和V3
题目大意：
求$\sum\limits_{i=1}^n \sum\limits_{j=1}^{n} lcm(i,j)$
思路：
原式=$\sum\limits_ {i=1}^n \sum\limits_{j=1}^{n} \frac{ij}{gcd(i,j)}$
$\sum\limits_{d=1}^n \ \frac{1}{d}\ \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} i \cdot j [gcd(i,j)==d]$
$\sum\limits_{d=1}^n \ d\ \sum\limits_{i=1}^{n/d} \sum\limits_{j=1}^{n/d} i \cdot j [gcd(i,j)==1]$
$\sum\limits_{d=1}^n \ d \ \sum\limits_{g=1}^{n/d} \mu(g) \sum\limits_{i=1}^{n/d} i\cdot[g|i] \sum\limits_{j=1}^{n/d} j \cdot [g|j]$
$\sum\limits_{d=1}^n \ d\ \sum\limits_{g=1}^{n/d} \mu(g) g^2 \sum\limits_{i=1}^{n/gd} i \sum\limits_{j=1}^{n/gd} j$
$\sum\limits_{g=1}^n \lgroup \frac{1}{2} \lfloor \frac{n}{g} \rfloor \cdot \lgroup \lfloor \frac{n}{g} \rfloor +1\rgroup\rgroup^2 \sum\limits_{d|g} \mu(g) g^2 \frac{d}{g}$
记：$f(i)=\sum\limits_{d|i} \mu(d) d^2 \frac{i}{d}$即$f(i)=(\mu \cdot id^2) * id$ 、记前缀和$S(i)=\frac{i(i+1)}{2}$
原式=$\frac{1}{2} \sum\limits_{g=1}^{n} S(\lfloor \frac{n}{g}\rfloor)^2 f(g) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \pmb{\lbrack1\rbrack}$
在这种形式下，只需考虑求$f(i)$前缀和，不妨令$g(i)=id^2$，则有：
$f*g=((\mu \cdot id^2)*id)*id^2$
由于$dirichlet$卷积满足交换律：
$f*g=((\mu\cdot id^2)*id^2)*id$
$(f*g)(n)=\sum\limits_{g|n} (\sum\limits_{d|g} \mu(d)\cdot d^2\cdot \frac{g^2}{d^2})\cdot \frac{n}{g}$
$=\sum\limits_{g|n} g^2\cdot [g=1] \cdot \frac{n}{g}$
$=n \sum\limits_{g|n} g\cdot [g=1]$
$=n$
记$s(n)=\sum\limits_{i=1}^n f(i)$
代入杜教筛公式即有$s(n)=-\sum\limits_{i=2}^n g(i) s(n/i) +\sum\limits_{i=1}^n i$
利用平方和公式求$g(i)$的区间和即可用杜教筛求出$f(i)$的前缀和，套用式$\ \pmb{\lbrack 1 \rbrack}\ $求和即可
另：
在线性筛预处理$f(i)$时，化为$f(i)=i\sum\limits_{d|i}\mu(d)d$，由于$\mu$的良好性质，最小因子出现多次的情况可以被忽略，容易在线筛过程中推出。
1 #include<bits/stdc++.h> 2 #define ll long long 3 #define MAXN 5001001 4 #define inf 2139062143 5 #define MOD 1000000007 6 #define rep(i,s,t) for(int i=(s),i##end=(t);i<=i##end;++i) 7 #define dwn(i,s,t) for(int i=(s),i##end=(t);i>=i##end;--i) 8 #define ren for(int i=fst[x];i;i=nxt[i]) 9 #define pls(a,b) (a+b)%MOD 10 #define mns(a,b) (a-(b)+MOD)%MOD 11 #define mul(a,b) (1LL*(a)*(b))%MOD 12 using namespace std; 13 int isp[MAXN],p[MAXN],tot;ll f[MAXN],s[MAXN],lmt,n; 14 const int inv2=500000004,inv6=166666668; 15 void mem(int n) 16 { 17 f[1]=1;rep(i,2,n) 18 { 19 if(!isp[i]) p[++tot]=i,f[i]=1-i; 20 rep(j,1,tot) if(1LL*i*p[j]>n) break; 21 else {isp[i*p[j]]=1;if(i%p[j]) f[i*p[j]]=f[i]*f[p[j]];else {f[i*p[j]]=f[i];break;}} 22 } 23 rep(i,1,n) f[i]=pls(f[i-1],mul(f[i],i)); 24 } 25 map<ll,int> ansf; 26 inline ll sum(ll n) 27 { 28 n%=MOD;return mul(mul(n,n+1),mul(2*n+1,inv6)); 29 } 30 inline ll gets(ll n) 31 { 32 n%=MOD;return mul(mul(n,n+1),inv2); 33 } 34 ll sumf(ll n,ll res=0) 35 { 36 if(n<=lmt) return f[n];if(ansf[n]) return ansf[n]; 37 for(ll l=2,r;l<=n;l=r+1) r=n/(n/l),res=pls(res,mul(mns(sum(r),sum(l-1)),sumf(n/l))); 38 return ansf[n]=mns(gets(n),res); 39 } 40 ll solve(ll n,ll res=0) 41 { 42 for(ll l=1,r;l<=n;l=r+1) 43 r=n/(n/l),res=pls(res,mul(mul(gets(n/l),gets(n/l)),mns(sumf(r),sumf(l-1)))); 44 return (res+MOD)%MOD; 45 } 46 int main() 47 { 48 scanf("%lld",&n);mem(lmt=5000000);printf("%lld\n",solve(n)); 49 }
51nod 1227 平均最小公倍数
题目大意：
令$A(n)=\frac{1}{n}\sum\limits_{i=1}^n lcm(i,n)$，求$\sum_{i=a}^b A(i)$
思路：
先化简$A(n)$：
$A(n)=\frac{1}{n}\sum\limits_ {i=1}^n\frac{ni}{gcd(n,i)}$
由于$gcd$为$n$的因数，因此枚举时满足$d|n$
$A(n)=\sum\limits_{d|n} \ \frac{1}{d}\ \sum\limits_{i=1}^{n} i [gcd(i,n)==d]$
$\sum\limits_{d|n}\sum\limits_{i=1}^{n/d} i\ [gcd(i,\frac{n}{d})==1]$
$\sum\limits_{d|n}\ \sum\limits_{i=1}^{n/d}i\sum\limits_{g|gcd(i,n/d)} \mu(g)$
$\sum\limits_{d|n}\sum\limits_{g|\frac{n}{d}} \mu(g)\ \sum\limits_{i=1}^{n/g}i\cdot \lbrack g|i \rbrack $
$\sum\limits_{d|n}\sum\limits_{g|\frac{n}{d}} \mu(g)g \sum\limits_{i=1}^{n/gd}i$
$\sum\limits_{g|n} \frac{\frac{n}{g} \cdot \lgroup \frac{n}{g}+1\rgroup}{2}\sum\limits_{d|g} \mu(d)d$$
令$h(n)=\sum\limits_{d|n} \mu(d)d$即$h(n)=(\mu \cdot id)*1$
即$A(n)=\frac{1}{2} (\sum\limits_{g|n} \frac{n}{g} \cdot h(g)+\sum\limits_{g|n} (\frac{n}{g})^2 \cdot h(g))$
要求$A(n)$的前缀和，将这两项分开考虑：
前半部分$=id*((\mu\cdot id)*1)=(\mu \cdot id )*id*1=[n=1]*1=1$，计算非常简单
后半部分$id^2 *(\mu \cdot id)*1$，设为$f(i)$，考虑取$g(i)=id$，则有：
$f*g=(\mu\cdot id)*id*1*id^2=1* id^2$
$\sum\limits_{i=1}^n (f*g)(n)=\sum\limits_{i=1}^n\sum\limits_{d|i}d^2$
$=\sum\limits_{d=1}^n d^2\sum\limits_{i=1}^n [d|i]$
$=\sum\limits_{d=1}^n d^2 \cdot \lfloor \frac{n}{d}\rfloor$
则该卷积结果的前缀和可以通过数论分块求出，在杜教筛求后半部分的同时顺便求出即可
将两部分的和相加后乘$2$的逆元即可
另：
在线性筛处理$f(i)$时，将$f(i)$化为$(\varphi\cdot id)*1$，在过程中处理出最小因子最高次幂的贡献，可以顺利处理出每个质因子次数$\ge2$的情况。
1 #include<bits/stdc++.h> 2 #define inf 2139062143 3 #define MAXN 5001001 4 #define MOD 1000000007 5 #define ll long long 6 #define ull unsigned long long 7 #define rep(i,s,t) for(register int i=(s),i##end=(t);i<=i##end;++i) 8 #define dwn(i,s,t) for(register int i=(s),i##end=(t);i>=i##end;--i) 9 #define ren for(register int i=fst[x];i;i=nxt[i]) 10 #define Fill(a,b) memset(a,b,sizeof(a)) 11 #define pls(a,b) (a+b)%MOD 12 #define mns(a,b) (a-(b)+MOD)%MOD 13 #define mul(a,b) (1LL*(a)*(b)%MOD) 14 #define pii pair<int,int> 15 #define mp(a,b) make_pair(a,b) 16 #define pb(i,x) vec[i].push_back(x) 17 #define fi first 18 #define se second 19 using namespace std; 20 inline ll read() 21 { 22 ll x=0;bool f=1;char ch=getchar(); 23 while(!isdigit(ch)) {if(ch=='-') f=0;ch=getchar();} 24 while(isdigit(ch)) x=x*10+(ch^48),ch=getchar(); 25 return f?x:-x; 26 } 27 int ntp[MAXN],p[MAXN],tot;ll f[MAXN],lmt,n,m,g[MAXN]; 28 const int inv2=500000004,inv6=166666668; 29 int gcd(int a,int b){return !b?a:gcd(b,a%b);} 30 void mem(int n) 31 { 32 f[1]=1;rep(i,2,n) 33 { 34 if(!ntp[i]) p[++tot]=i,f[i]=mns(mul(i,i),i-1),g[i]=mns(f[i],1); 35 rep(j,1,tot) if(1LL*i*p[j]>n) break; 36 else 37 { 38 ntp[i*p[j]]=1;if(i%p[j]) f[i*p[j]]=mul(f[i],f[p[j]]),g[i*p[j]]=mns(f[i*p[j]],f[i]); 39 else {g[i*p[j]]=mul(g[i],mul(p[j],p[j])),f[i*p[j]]=pls(g[i*p[j]],f[i]);break;} 40 } 41 } 42 rep(i,1,n) f[i]=pls(f[i-1],f[i]); 43 } 44 map<ll,int> ansf; 45 ll gets(ll n,ll res=0) 46 { 47 n%=MOD;return mul(mul(n,n+1),mul(2*n+1,inv6)); 48 } 49 ll solve(ll n,ll sum=0,ll res=0) 50 { 51 if(n<=lmt) return f[n];if(ansf[n]) return ansf[n]; 52 sum=n%MOD;for(ll l=2,r;l<=n;l=r+1) 53 r=n/(n/l),res=pls(mul(mul(solve(n/l),inv2),mul(r-l+1,(r+l)%MOD)),res), 54 sum=pls(sum,mul((n/l)%MOD,mns(gets(r),gets(l-1)))); 55 return ansf[n]=mns(sum,res); 56 } 57 int main() 58 { 59 m=read(),n=read();mem(lmt=5000000); 60 printf("%lld\n",mul(pls(mns(solve(n),solve(m-1)),mns(n+1,m)),inv2)); 61 }
51nod 1222 最小公倍数计数
题目大意：
定义$F(i)$表示最小公倍数为$n$的二元组$(x,y)$其中$x\le y$的数量，求$\sum\limits_{i=a}^b F(i)$
思路：
题目所求显然可以通过计算不考虑$x$与$y$的大小关系的情况下的答案快速得出，即加上$x=y$的情况再除以二即可。
只需计算$F(i)$的前缀和，即：
$\sum\limits_{i=1}^n F(i)=\sum\limits_{i=1}^n \sum\limits_{j=1}^n [lcm(i,j)\le n]$
$=\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^n [\frac{ij}{d}\le n][gcd(i,j)=d]$
$=\sum\limits_{d=1}^n \sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{n/d} [ijd\le n][gcd(i,j)=1]$
$=\sum\limits_{d=1}^n \sum\limits_{k=1}^{n/d}\mu(k)\sum\limits_{i=1}^{n/dk}\sum\limits_{j=1}^{n/dk} [ijd\le \lfloor \frac{n}{k^2}\rfloor]$
考虑先枚举$k$，由于后续需要计算$ijd\le \lfloor \frac{n}{k^2}\rfloor $，故$k$只需枚举到$\sqrt n$即可，后续枚举时上界也同样可以缩小，即：
$\sum\limits_{i=1}^nF(i)=\sum\limits_{k=1}^{\sqrt n} \mu(k) \sum\limits_{d=1}^{\lfloor \frac{n}{k^2}\rfloor }\sum\limits_{i=1}^{\lfloor\frac{n}{dk^2}\rfloor} \sum\limits_{j=1}^{\lfloor \frac{n}{idk^2}\rfloor} [ijd\le \lfloor \frac{n}{k^2}\rfloor]$
令$g(n)=\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor} \sum\limits_{j=1}^{\lfloor \frac{n}{id}\rfloor} [ijd\le n]$，在枚举时不妨令$d\le i\le j$，则在枚举$d\le \sqrt[3]n,i\le \sqrt{\frac{n}{i} }$之后可以直接计算出符合条件的$j$的数量，再特判相等时候的贡献即可。
1 #include<bits/stdc++.h> 2 #define inf 2139062143 3 #define MAXN 500100 4 #define MOD 1000000007 5 #define ll long long 6 #define ull unsigned long long 7 #define rep(i,s,t) for(register int i=(s),i##end=(t);i<=i##end;++i) 8 #define dwn(i,s,t) for(register int i=(s),i##end=(t);i>=i##end;--i) 9 #define ren for(register int i=fst[x];i;i=nxt[i]) 10 #define Fill(a,b) memset(a,b,sizeof(a)) 11 #define pls(a,b) (a+b)%MOD 12 #define mns(a,b) (a-(b)+MOD)%MOD 13 #define mul(a,b) (1LL*(a)*(b)%MOD) 14 #define pii pair<int,int> 15 #define mp(a,b) make_pair(a,b) 16 #define pb(i,x) vec[i].push_back(x) 17 #define fi first 18 #define se second 19 using namespace std; 20 inline ll read() 21 { 22 ll x=0;bool f=1;char ch=getchar(); 23 while(!isdigit(ch)) {if(ch=='-') f=0;ch=getchar();} 24 while(isdigit(ch)) x=x*10+(ch^48),ch=getchar(); 25 return f?x:-x; 26 } 27 ll a,b; 28 int mu[MAXN],p[MAXN],tot,isp[MAXN]; 29 void mem(int n) 30 { 31 mu[1]=1;rep(i,2,n) 32 { 33 if(!isp[i]) p[++tot]=i,mu[i]=-1; 34 rep(j,1,tot) if(i*p[j]>n) break; 35 else {isp[i*p[j]]=1;if(i%p[j]) mu[i*p[j]]=-mu[i];else break;} 36 } 37 } 38 ll calc(ll n,ll res=0) 39 { 40 ll t;for(ll i=1;i*i*i<=n;++i) 41 for(ll j=i;i*j*j<=n;++j) 42 { 43 t=n/i/j;if(t<j) break;t-=j-1; 44 if(i==j) res+=t*3-2;else res+=t*6-3; 45 } 46 return res; 47 } 48 ll solve(ll n,ll res=0) 49 { 50 rep(i,1,sqrt(n)) res+=mu[i]*calc(n/i/i);return (res+n)>>1; 51 } 52 int main() 53 { 54 a=read(),b=read();mem(500000); 55 printf("%lld\n",solve(b)-solve(a-1)); 56 }
bzoj 4176 Lucas的数论
题目大意：
求$\sum\limits_{i=1}^n \sum\limits_{j=1}^n d(ij)$，其中$d(n)$表示$n$的约数个数
思路：
首先有结论$d(nm)=\sum\limits_{i|n}\sum\limits_{j|m} [gcd(i,j)=1]$，简略证明：
设$n=\prod\limits_{i=1} p_i^{\alpha_i}$，$m=\prod\limits_{j=1}p_j^{\beta_j}$
考虑对每个质因数$p$，符合条件的$i$,$j$一定只能有一个含有因子$p$
两种情况分别有$\beta$与$\alpha$ ，一共$\alpha + \beta$恰好表示$p$对于$nm$的贡献
代入得：
原式 $=\sum\limits_{i=1}^n \sum\limits_{j=1}^n \sum\limits_{p|i}\sum\limits_{q|j} [gcd(p,q)==1]$
$=\sum\limits_{i=1}^n \sum\limits_{j=1}^n \sum\limits_{p|i}\sum\limits_{q|j}\sum\limits_{d=1}^n \mu(d) [d|p][d|q]$
$=\sum\limits_{d=1}^n \mu(d)\sum\limits_{p=1}^{n/d} \sum\limits_{q=1}^{n/d} \sum\limits_{i=1}^{n/pd}\sum\limits_{j=1}^{n/qd} 1$
$=\sum\limits_{d=1}^n \mu(d) \lgroup \sum\limits_{p=1}^{n/d} \lfloor \frac{n}{pd}\rfloor \rgroup \lgroup \sum\limits_{q=1}^{n/d} \lfloor \frac{n}{qd}\rfloor \rgroup $
$=\sum\limits_{d=1}^n \mu(d) \lgroup \sum\limits_{i=1}^{n/d} \lfloor \frac{n}{id}\rfloor \rgroup^2 $
令$f(n)=\sum\limits_{i=1}^n \lfloor \frac{n}{i}\rfloor $，即原式$=\sum\limits_{d=1}^n \mu(d) f(\lfloor \frac{n}{d}\rfloor)^2 $
$f(i)$可以通过分块快速计算，只需预处理$\mu$的前缀和即可
1 #include<bits/stdc++.h> 2 #define inf 2139062143 3 #define MAXN 5001001 4 #define MOD 1000000007 5 #define ll long long 6 #define ull unsigned long long 7 #define rep(i,s,t) for(register int i=(s),i##end=(t);i<=i##end;++i) 8 #define dwn(i,s,t) for(register int i=(s),i##end=(t);i>=i##end;--i) 9 #define ren for(register int i=fst[x];i;i=nxt[i]) 10 #define Fill(a,b) memset(a,b,sizeof(a)) 11 #define pls(a,b) (a+b)%MOD 12 #define mns(a,b) (a-(b))%MOD 13 #define mul(a,b) (1LL*(a)*(b)%MOD) 14 #define pii pair<int,int> 15 #define mp(a,b) make_pair(a,b) 16 #define pb(i,x) vec[i].push_back(x) 17 #define fi first 18 #define se second 19 using namespace std; 20 inline int read() 21 { 22 int x=0;bool f=1;char ch=getchar(); 23 while(!isdigit(ch)) {if(ch=='-') f=0;ch=getchar();} 24 while(isdigit(ch)) x=x*10+(ch^48),ch=getchar(); 25 return f?x:-x; 26 } 27 int n,lmt,mu[MAXN],p[MAXN],tot,isp[MAXN]; 28 void mem(int n) 29 { 30 mu[1]=1;rep(i,2,n) 31 { 32 if(!isp[i]) p[++tot]=i,mu[i]=-1; 33 rep(j,1,tot) if(i*p[j]>=n) break; 34 else {isp[i*p[j]]=1;if(i%p[j]) mu[i*p[j]]=-mu[i];else break;} 35 } 36 rep(i,1,n) mu[i]=pls(mu[i-1],mu[i]); 37 } 38 map<int,int> ansm; 39 int smu(int n,int res=0) 40 { 41 if(n<=lmt) return mu[n];if(ansm[n]) return ansm[n]; 42 for(int l=2,r;l<=n;l=r+1) r=n/(n/l),res=pls(mul(r-l+1,smu(n/l)),res); 43 return ansm[n]=mns(1,res); 44 } 45 int calc(int n,int res=0) 46 { 47 for(int l=1,r;l<=n;l=r+1) r=n/(n/l),res=pls(res,mul(r-l+1,(n/l)%MOD)); 48 return mul(res,res); 49 } 50 int solve(int n,int res=0) 51 { 52 for(int l=1,r;l<=n;l=r+1) 53 r=n/(n/l),res=pls(mul(mns(smu(r),smu(l-1)),calc(n/l)),res); 54 return (res+MOD)%MOD; 55 } 56 int main() 57 { 58 n=read();mem(lmt=5000000);printf("%d\n",solve(n)); 59 }
51nod 1220 约数之和
题目大意：
求$\sum\limits_{i=1}^n \sum\limits_{j=1}^n \sigma(ij)$，其中$\sigma(n)$表示$n$的约数和
思路：
首先有结论$\sigma(nm)=\sum\limits_{i|n}\sum\limits_{j|m} \frac{im}{j}\cdot [gcd(i,j)=1]$，证明与上题类似：
设$n=\prod\limits_{i=1} p_i^{\alpha_i}$，$m=\prod\limits_{j=1}p_j^{\beta_j}$
考虑对每个质因数$p$，符合条件的$i$,$j$一定只能有一个含有因子$p$
两种情况分别对应$1\ \ p^1\cdots p^{\beta}$与$p^{\beta}\cdots p^{\alpha+\beta}$ ，一共恰好表示$p$对于$nm$的贡献
将这个式子代入：
原式 $=\sum\limits_{i=1}^n \sum\limits_{j=1}^n \sum\limits_{p|i}\sum\limits_{q|j} \frac{pj}{q} \cdot [gcd(p,q)==1]$
$=\sum\limits_{i=1}^n \sum\limits_{j=1}^n \sum\limits_{p|i}\sum\limits_{q|j}\sum\limits_{d=1}^n \mu(d) \frac{pj}{q}[d|p][d|q]$
$=\sum\limits_{d=1}^n \mu(d)\sum\limits_{p=1}^{n/d} \sum\limits_{q=1}^{n/d} \sum\limits_{i=1}^{n/pd}\sum\limits_{j=1}^{n/qd} pd \cdot j$
$=\sum\limits_{d=1}^n \mu(d)d \lgroup \sum\limits_{p=1}^{n/d} p \lfloor \frac{n}{pd}\rfloor \rgroup \lgroup \sum\limits_{q=1}^{n/d} \sum\limits_{j=1}^{n/qd} j\rgroup$
设中间两部分分别为$f(\lfloor\frac{n}{d}\rfloor)$与$g(\lfloor\frac{n}{d}\rfloor)$，其中$f(n)=\sum\limits_{i=1}^n i\lfloor\frac{n}{i}\rfloor$，$g(n)=\sum\limits_{i=1}^n\sum\limits_{j=1}^{\lfloor\frac{n}{i}\rfloor}j$
考虑$f(n)$的实际意义，即每个$i$乘以它在$\le n$范围内倍数个数，故$f(n)=\sum\limits_{i=1}^n \sigma_1(i)$
又$g(n)=\sum\limits_{j=1}^{n}j\sum\limits_{i=1}^{\lfloor\frac{n}{j}\rfloor}1=\sum\limits_{j=1}^n j \lfloor\frac{n}{j}\rfloor=f(n)$
代入原式$=\sum\limits_{d=1}^n \mu(d)d\ \lgroup \sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sigma_1(i)\rgroup^2$
对于$\mu \cdot id$的前缀和，与之前题目做法一样，卷上$id$即可解决，
对$f(n/d)$的计算，较小的情况可以用线性筛预处理，记录一下最小质因子的贡献即可；较大的情况直接采用$f(n)$的式子计算即可，复杂度分析与杜教筛类似。
1 #include<bits/stdc++.h> 2 #define inf 2139062143 3 #define MAXN 5001001 4 #define MOD 1000000007 5 #define ll long long 6 #define ull unsigned long long 7 #define rep(i,s,t) for(register int i=(s),i##end=(t);i<=i##end;++i) 8 #define dwn(i,s,t) for(register int i=(s),i##end=(t);i>=i##end;--i) 9 #define ren for(register int i=fst[x];i;i=nxt[i]) 10 #define Fill(a,b) memset(a,b,sizeof(a)) 11 #define pls(a,b) (a+b)%MOD 12 #define mns(a,b) (a-(b))%MOD 13 #define mul(a,b) (1LL*(a)*(b)%MOD) 14 #define pii pair<int,int> 15 #define mp(a,b) make_pair(a,b) 16 #define pb(i,x) vec[i].push_back(x) 17 #define fi first 18 #define se second 19 using namespace std; 20 inline int read() 21 { 22 int x=0;bool f=1;char ch=getchar(); 23 while(!isdigit(ch)) {if(ch=='-') f=0;ch=getchar();} 24 while(isdigit(ch)) x=x*10+(ch^48),ch=getchar(); 25 return f?x:-x; 26 } 27 int n,lmt; 28 int mu[MAXN],f[MAXN],g[MAXN],p[MAXN],tot,isp[MAXN]; 29 void mem(int n) 30 { 31 mu[1]=f[1]=g[1]=1;rep(i,2,n) 32 { 33 if(!isp[i]) p[++tot]=i,mu[i]=-1,f[i]=1,g[i]=1+i; 34 rep(j,1,tot) if(i*p[j]>n) break; 35 else 36 { 37 isp[i*p[j]]=1; 38 if(i%p[j]) mu[i*p[j]]=-mu[i],f[i*p[j]]=mul(f[i],g[i]),g[i*p[j]]=p[j]+1; 39 else {f[i*p[j]]=f[i],g[i*p[j]]=pls(mul(g[i],p[j]),1);break;} 40 } 41 } 42 rep(i,1,n) mu[i]=pls(mu[i-1],mul(mu[i],i)); 43 rep(i,1,n) f[i]=pls(f[i-1],mul(f[i],g[i])); 44 } 45 map<int,int> ansm; 46 int sum(int l,int r) {return (1LL*(r-l+1)*(r+l)>>1LL)%MOD;} 47 int smu(int n,int res=0) 48 { 49 if(n<=lmt) return mu[n];if(ansm[n]) return ansm[n]; 50 for(int l=2,r;l<=n;l=r+1) r=n/(n/l),res=pls(mul(sum(l,r),smu(n/l)),res); 51 return ansm[n]=mns(1,res); 52 } 53 int calc(int n,int res=0) 54 { 55 if(n<=lmt) return mul(f[n],f[n]); 56 for(int l=1,r;l<=n;l=r+1) r=n/(n/l),res=pls(res,mul(sum(l,r),n/l)); 57 return mul(res,res); 58 } 59 int solve(int n,int res=0) 60 { 61 for(int l=1,r;l<=n;l=r+1) 62 r=n/(n/l),res=pls(mul(mns(smu(r),smu(l-1)),calc(n/l)),res); 63 return (res+MOD)%MOD; 64 } 65 int main() 66 { 67 n=read();mem(lmt=5000000);printf("%d\n",solve(n)); 68 }
51nod 1584 加权约数和
题目大意：
$\sum\limits_{i=1}^n \sum\limits_{j=1}^n max(i,j)\cdot \sigma(ij)$
思路：
先将$max(i,j)$转化为另一个形式$\sum\limits_{k=1}^n [k\le i\ or\ k\ge j]$，即$n-\sum\limits_{k=1}^n[k>i][k>j]$
将式子代入：
原式$=\sum\limits_{i=1}^n\sum\limits_{j=1}^n (n-\sum\limits_{k=1}^n [k>i][k>j])\cdot \sigma(ij)$
$=n\cdot \sum\limits_{i=1}^n\sum\limits_{j=1}^n\sigma(ij)-\sum\limits_{k=1}^n\sum\limits_{i=1}^{k-1}\sum\limits_{j=1}^{k-1}\sigma(ij)$
令$f(n)=\sum\limits_{i=1}^n\sum\limits_{j=1}^n\sigma(ij)$，则原式$=n\cdot f(n)-\sum\limits_{k=1}^{n-1} f(k)$
发现$f(n)$就是上一题的式子
$f(n)=\sum\limits_{d=1}^n \mu(d)d\ \lgroup \sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sigma(i)\rgroup^2$
由于本题需要计算$f(n)$的前缀和，使用杜教筛复杂度$O(n\sqrt n)$过高
因为$n\le 10^6$且有多组数据，考虑使用$nlogn$预处理$f(n)$，以便快速查询
在线性筛出$\mu$与$\sum\limits_{i=1}^n\sigma(i)$之后，枚举$d$与$n/d$之后对一段连续的$n$的贡献相同，因此可以快速得到$f(n)$的差分数组
然后求两次前缀和即可得到$f(n)$的前缀和数组
1 #include<bits/stdc++.h> 2 #define inf 2139062143 3 #define MAXN 1001001 4 #define MOD 1000000007 5 #define ll long long 6 #define ull unsigned long long 7 #define rep(i,s,t) for(register int i=(s),i##end=(t);i<=i##end;++i) 8 #define dwn(i,s,t) for(register int i=(s),i##end=(t);i>=i##end;--i) 9 #define ren for(register int i=fst[x];i;i=nxt[i]) 10 #define Fill(a,b) memset(a,b,sizeof(a)) 11 #define pls(a,b) (a+b)%MOD 12 #define mns(a,b) (a-(b))%MOD 13 #define mul(a,b) (1LL*(a)*(b)%MOD) 14 #define pii pair<int,int> 15 #define mp(a,b) make_pair(a,b) 16 #define pb(i,x) vec[i].push_back(x) 17 #define fi first 18 #define se second 19 using namespace std; 20 inline int read() 21 { 22 int x=0;bool f=1;char ch=getchar(); 23 while(!isdigit(ch)) {if(ch=='-') f=0;ch=getchar();} 24 while(isdigit(ch)) x=x*10+(ch^48),ch=getchar(); 25 return f?x:-x; 26 } 27 int n,lmt; 28 int mu[MAXN],f[MAXN],g[MAXN],p[MAXN],tot,isp[MAXN],ans[MAXN<<1]; 29 void mem(int n) 30 { 31 mu[1]=f[1]=g[1]=1;rep(i,2,n) 32 { 33 if(!isp[i]) p[++tot]=i,mu[i]=-1,f[i]=1,g[i]=1+i; 34 rep(j,1,tot) if(i*p[j]>n) break; 35 else 36 { 37 isp[i*p[j]]=1; 38 if(i%p[j]) mu[i*p[j]]=-mu[i],f[i*p[j]]=mul(f[i],g[i]),g[i*p[j]]=p[j]+1; 39 else {f[i*p[j]]=f[i],g[i*p[j]]=pls(mul(g[i],p[j]),1);break;} 40 } 41 } 42 rep(i,1,n) mu[i]=mul(mu[i],i),f[i]=pls(mul(f[i],g[i]),f[i-1]); 43 rep(i,1,n) f[i]=mul(f[i],f[i]); 44 } 45 inline int calc(int n){return mns(mul(n,ans[n]),mul(n+1,ans[n-1]));} 46 int main() 47 { 48 mem(lmt=1000000); 49 rep(i,1,lmt) rep(j,1,lmt/i) 50 ans[i*j]=pls(ans[i*j],mul(mu[i],f[j])),ans[i*j+i]=mns(ans[i*j+i],mul(mu[i],f[j])); 51 rep(i,1,lmt) ans[i]=pls(ans[i-1],ans[i]); 52 rep(i,1,lmt) ans[i]=pls(ans[i-1],ans[i]); 53 for(int T=read(),t=1;t<=T;++t) 54 n=read(),printf("Case #%d: %d\n",t,(calc(n)+MOD)%MOD); 55 }
