莫比乌斯反演的前置知识：https://www.cnblogs.com/yyc-jack-0920/p/10556351.html
杜教筛可以在$O(n^{2/3})$的时间内求出积性函数的前缀和即$sumlimits_{i=1}^n f(i)$
实现过程需要找到另一个积性函数$g$使得$f*g$为一个方便计算的函数
$sum f*g=sumlimits_{i=1}^n sumlimits_{d|i} g(d)f(i/d) $
$=sumlimits_{d=1}^n g(d) sumlimits_{t=1}^{n/d} f(t)$
$=g(1)sumlimits_{i=1}^n f(i) + sumlimits_{d=2}^n g(d) sumlimits_{i=1}^{n/d} f(i)$
设前缀和为$S(i)$ 即$g(1)S(n)=-sumlimits_{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 ",sphi(n),smu(n));} 59 }
51nod 1237 最大公约数之和V3
题目大意：
求$sumlimits_{i=1}^n sumlimits_{j=1}^{n} gcd(i,j)$
思路：
原式=$sumlimits_{d=1}^{n} d sumlimits_{i=1}^{n/d} sumlimits_{j=1}^{n/d} [gcd(i,j)==1]$
$sumlimits_{d=1}^n d sumlimits_{i=1}^{n/d} mu(i) lfloorfrac{n}{id} floor lfloorfrac{n}{id} floor$
$sumlimits_{g=1}^n lfloorfrac{n}{g} floor lfloorfrac{n}{g} floor sumlimits_{d|g} mu(d)*frac{g}{d}$
$sumlimits_{g=1}^n lfloorfrac{n}{g} floor lfloorfrac{n}{g} floor 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 ",solve(n)); 63 }
51nod 1238 最大公倍数之和V3
题目大意：
求$sumlimits_{i=1}^n sumlimits_{j=1}^{n} lcm(i,j)$
思路：
原式=$sumlimits_ {i=1}^n sumlimits_{j=1}^{n} frac{ij}{gcd(i,j)}$
$sumlimits_{d=1}^n frac{1}{d} sumlimits_{i=1}^{n} sumlimits_{j=1}^{n} i cdot j [gcd(i,j)==d]$
$sumlimits_{d=1}^n d sumlimits_{i=1}^{n/d} sumlimits_{j=1}^{n/d} i cdot j [gcd(i,j)==1]$
$sumlimits_{d=1}^n d sumlimits_{g=1}^{n/d} mu(g) sumlimits_{i=1}^{n/d} icdot[g|i] sumlimits_{j=1}^{n/d} j cdot [g|j]$
$sumlimits_{d=1}^n d sumlimits_{g=1}^{n/d} mu(g) g^2 sumlimits_{i=1}^{n/gd} i sumlimits_{j=1}^{n/gd} j$
$sumlimits_{g=1}^n lgroup frac{1}{2} lfloor frac{n}{g} floor cdot lgroup lfloor frac{n}{g} floor +1 group group^2 sumlimits_{d|g} mu(g) g^2 frac{d}{g}$
记：$f(i)=sumlimits_{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} sumlimits_{g=1}^{n} S(lfloor frac{n}{g} floor)^2 f(g) pmb{lbrack1 brack}$
在这种形式下，只需考虑求$f(i)$前缀和，不妨令$g(i)=id^2$，则有：
$f*g=((mu cdot id^2)*id)*id^2$
由于$dirichlet$卷积满足交换律：
$f*g=((mucdot id^2)*id^2)*id$
$(f*g)(n)=sumlimits_{g|n} (sumlimits_{d|g} mu(d)cdot d^2cdot frac{g^2}{d^2})cdot frac{n}{g}$
$=sumlimits_{g|n} g^2cdot [g=1] cdot frac{n}{g}$
$=n sumlimits_{g|n} gcdot [g=1]$
$=n$
记$s(n)=sumlimits_{i=1}^n f(i)$
代入杜教筛公式即有$s(n)=-sumlimits_{i=2}^n g(i) s(n/i) +sumlimits_{i=1}^n i$
利用平方和公式求$g(i)$的区间和即可用杜教筛求出$f(i)$的前缀和，套用式$ pmb{lbrack 1 brack} $求和即可
另：
在线性筛预处理$f(i)$时，化为$f(i)=isumlimits_{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 ",solve(n)); 49 }
51nod 1227 平均最小公倍数
题目大意：
令$A(n)=frac{1}{n}sumlimits_{i=1}^n lcm(i,n)$，求$sum_{i=a}^b A(i)$
思路：
先化简$A(n)$：
$A(n)=frac{1}{n}sumlimits_ {i=1}^nfrac{ni}{gcd(n,i)}$
由于$gcd$为$n$的因数，因此枚举时满足$d|n$
$A(n)=sumlimits_{d|n} frac{1}{d} sumlimits_{i=1}^{n} i [gcd(i,n)==d]$
$sumlimits_{d|n}sumlimits_{i=1}^{n/d} i [gcd(i,frac{n}{d})==1]$
$sumlimits_{d|n} sumlimits_{i=1}^{n/d}isumlimits_{g|gcd(i,n/d)} mu(g)$
$sumlimits_{d|n}sumlimits_{g|frac{n}{d}} mu(g) sumlimits_{i=1}^{n/g}icdot lbrack g|i brack $
$sumlimits_{d|n}sumlimits_{g|frac{n}{d}} mu(g)g sumlimits_{i=1}^{n/gd}i$
$sumlimits_{g|n} frac{frac{n}{g} cdot lgroup frac{n}{g}+1 group}{2}sumlimits_{d|g} mu(d)d$$
令$h(n)=sumlimits_{d|n} mu(d)d$即$h(n)=(mu cdot id)*1$
即$A(n)=frac{1}{2} (sumlimits_{g|n} frac{n}{g} cdot h(g)+sumlimits_{g|n} (frac{n}{g})^2 cdot h(g))$
要求$A(n)$的前缀和，将这两项分开考虑：
前半部分$=id*((mucdot id)*1)=(mu cdot id )*id*1=[n=1]*1=1$，计算非常简单
后半部分$id^2 *(mu cdot id)*1$，设为$f(i)$，考虑取$g(i)=id$，则有：
$f*g=(mucdot id)*id*1*id^2=1* id^2$
$sumlimits_{i=1}^n (f*g)(n)=sumlimits_{i=1}^nsumlimits_{d|i}d^2$
$=sumlimits_{d=1}^n d^2sumlimits_{i=1}^n [d|i]$
$=sumlimits_{d=1}^n d^2 cdot lfloor frac{n}{d} floor$
则该卷积结果的前缀和可以通过数论分块求出，在杜教筛求后半部分的同时顺便求出即可
将两部分的和相加后乘$2$的逆元即可
另：
在线性筛处理$f(i)$时，将$f(i)$化为$(varphicdot 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 ",mul(pls(mns(solve(n),solve(m-1)),mns(n+1,m)),inv2)); 61 }
51nod 1222 最小公倍数计数
题目大意：
定义$F(i)$表示最小公倍数为$n$的二元组$(x,y)$其中$xle y$的数量，求$sumlimits_{i=a}^b F(i)$
思路：
题目所求显然可以通过计算不考虑$x$与$y$的大小关系的情况下的答案快速得出，即加上$x=y$的情况再除以二即可。
只需计算$F(i)$的前缀和，即：
$sumlimits_{i=1}^n F(i)=sumlimits_{i=1}^n sumlimits_{j=1}^n [lcm(i,j)le n]$
$=sumlimits_{d=1}^nsumlimits_{i=1}^nsumlimits_{j=1}^n [frac{ij}{d}le n][gcd(i,j)=d]$
$=sumlimits_{d=1}^n sumlimits_{i=1}^{n/d}sumlimits_{j=1}^{n/d} [ijdle n][gcd(i,j)=1]$
$=sumlimits_{d=1}^n sumlimits_{k=1}^{n/d}mu(k)sumlimits_{i=1}^{n/dk}sumlimits_{j=1}^{n/dk} [ijdle lfloor frac{n}{k^2} floor]$
考虑先枚举$k$，由于后续需要计算$ijdle lfloor frac{n}{k^2} floor $，故$k$只需枚举到$sqrt n$即可，后续枚举时上界也同样可以缩小，即：
$sumlimits_{i=1}^nF(i)=sumlimits_{k=1}^{sqrt n} mu(k) sumlimits_{d=1}^{lfloor frac{n}{k^2} floor }sumlimits_{i=1}^{lfloorfrac{n}{dk^2} floor} sumlimits_{j=1}^{lfloor frac{n}{idk^2} floor} [ijdle lfloor frac{n}{k^2} floor]$
令$g(n)=sumlimits_{d=1}^nsumlimits_{i=1}^{lfloorfrac{n}{d} floor} sumlimits_{j=1}^{lfloor frac{n}{id} floor} [ijdle n]$，在枚举时不妨令$dle ile j$，则在枚举$dle sqrt[3]n,ile 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 ",solve(b)-solve(a-1)); 56 }
bzoj 4176 Lucas的数论
题目大意：
求$sumlimits_{i=1}^n sumlimits_{j=1}^n d(ij)$，其中$d(n)$表示$n$的约数个数
思路：
首先有结论$d(nm)=sumlimits_{i|n}sumlimits_{j|m} [gcd(i,j)=1]$，简略证明：
设$n=prodlimits_{i=1} p_i^{alpha_i}$，$m=prodlimits_{j=1}p_j^{eta_j}$
考虑对每个质因数$p$，符合条件的$i$,$j$一定只能有一个含有因子$p$
两种情况分别有$eta$与$alpha$ ，一共$alpha + eta$恰好表示$p$对于$nm$的贡献
代入得：
原式 $=sumlimits_{i=1}^n sumlimits_{j=1}^n sumlimits_{p|i}sumlimits_{q|j} [gcd(p,q)==1]$
$=sumlimits_{i=1}^n sumlimits_{j=1}^n sumlimits_{p|i}sumlimits_{q|j}sumlimits_{d=1}^n mu(d) [d|p][d|q]$
$=sumlimits_{d=1}^n mu(d)sumlimits_{p=1}^{n/d} sumlimits_{q=1}^{n/d} sumlimits_{i=1}^{n/pd}sumlimits_{j=1}^{n/qd} 1$
$=sumlimits_{d=1}^n mu(d) lgroup sumlimits_{p=1}^{n/d} lfloor frac{n}{pd} floor group lgroup sumlimits_{q=1}^{n/d} lfloor frac{n}{qd} floor group $
$=sumlimits_{d=1}^n mu(d) lgroup sumlimits_{i=1}^{n/d} lfloor frac{n}{id} floor group^2 $
令$f(n)=sumlimits_{i=1}^n lfloor frac{n}{i} floor $，即原式$=sumlimits_{d=1}^n mu(d) f(lfloor frac{n}{d} floor)^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 ",solve(n)); 59 }
51nod 1220 约数之和
题目大意：
求$sumlimits_{i=1}^n sumlimits_{j=1}^n sigma(ij)$，其中$sigma(n)$表示$n$的约数和
思路：
首先有结论$sigma(nm)=sumlimits_{i|n}sumlimits_{j|m} frac{im}{j}cdot [gcd(i,j)=1]$，证明与上题类似：
设$n=prodlimits_{i=1} p_i^{alpha_i}$，$m=prodlimits_{j=1}p_j^{eta_j}$
考虑对每个质因数$p$，符合条件的$i$,$j$一定只能有一个含有因子$p$
两种情况分别对应$1 p^1cdots p^{eta}$与$p^{eta}cdots p^{alpha+eta}$ ，一共恰好表示$p$对于$nm$的贡献
将这个式子代入：
原式 $=sumlimits_{i=1}^n sumlimits_{j=1}^n sumlimits_{p|i}sumlimits_{q|j} frac{pj}{q} cdot [gcd(p,q)==1]$
$=sumlimits_{i=1}^n sumlimits_{j=1}^n sumlimits_{p|i}sumlimits_{q|j}sumlimits_{d=1}^n mu(d) frac{pj}{q}[d|p][d|q]$
$=sumlimits_{d=1}^n mu(d)sumlimits_{p=1}^{n/d} sumlimits_{q=1}^{n/d} sumlimits_{i=1}^{n/pd}sumlimits_{j=1}^{n/qd} pd cdot j$
$=sumlimits_{d=1}^n mu(d)d lgroup sumlimits_{p=1}^{n/d} p lfloor frac{n}{pd} floor group lgroup sumlimits_{q=1}^{n/d} sumlimits_{j=1}^{n/qd} j group$
设中间两部分分别为$f(lfloorfrac{n}{d} floor)$与$g(lfloorfrac{n}{d} floor)$，其中$f(n)=sumlimits_{i=1}^n ilfloorfrac{n}{i} floor$，$g(n)=sumlimits_{i=1}^nsumlimits_{j=1}^{lfloorfrac{n}{i} floor}j$
考虑$f(n)$的实际意义，即每个$i$乘以它在$le n$范围内倍数个数，故$f(n)=sumlimits_{i=1}^n sigma_1(i)$
又$g(n)=sumlimits_{j=1}^{n}jsumlimits_{i=1}^{lfloorfrac{n}{j} floor}1=sumlimits_{j=1}^n j lfloorfrac{n}{j} floor=f(n)$
代入原式$=sumlimits_{d=1}^n mu(d)d lgroup sumlimits_{i=1}^{lfloorfrac{n}{d} floor}sigma_1(i) group^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 ",solve(n)); 68 }
51nod 1584 加权约数和
题目大意：
$sumlimits_{i=1}^n sumlimits_{j=1}^n max(i,j)cdot sigma(ij)$
思路：
先将$max(i,j)$转化为另一个形式$sumlimits_{k=1}^n [kle i or kge j]$，即$n-sumlimits_{k=1}^n[k>i][k>j]$
将式子代入：
原式$=sumlimits_{i=1}^nsumlimits_{j=1}^n (n-sumlimits_{k=1}^n [k>i][k>j])cdot sigma(ij)$
$=ncdot sumlimits_{i=1}^nsumlimits_{j=1}^nsigma(ij)-sumlimits_{k=1}^nsumlimits_{i=1}^{k-1}sumlimits_{j=1}^{k-1}sigma(ij)$
令$f(n)=sumlimits_{i=1}^nsumlimits_{j=1}^nsigma(ij)$，则原式$=ncdot f(n)-sumlimits_{k=1}^{n-1} f(k)$
发现$f(n)$就是上一题的式子
$f(n)=sumlimits_{d=1}^n mu(d)d lgroup sumlimits_{i=1}^{lfloorfrac{n}{d} floor}sigma(i) group^2$
由于本题需要计算$f(n)$的前缀和，使用杜教筛复杂度$O(nsqrt n)$过高
因为$nle 10^6$且有多组数据，考虑使用$nlogn$预处理$f(n)$，以便快速查询
在线性筛出$mu$与$sumlimits_{i=1}^nsigma(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 ",t,(calc(n)+MOD)%MOD); 55 }