一．预备知识
积性函数的内容在另外一篇随笔中。
二．重要知识
关于莫比乌斯反演，有很多很好的博客，这里推荐几个
莫比乌斯反演入门 对于这篇博客中的内容，只需要知道 (mu) 其实是一种容斥系数即可，具体的乱七八糟的过程无需理会
莫比乌斯反演简要笔记
（一）莫比乌斯函数
１．莫比乌斯函数 (mu(n)) 的性质
- 设 (n=p_{1}^{k_{1}}*p_{2}^{k_{2}}*...*p_{m}^{k_{m}})，其中 (p) 为素数，那么
[mu(n)=egin{cases}
1 quadquadquadquad n=1 \\
(-1)^{m} quad prod_{i=1}^{m}k_{i}=1 \\
0 quadquadquadquad otherwise
end{cases}]
-
莫比乌斯函数是积性函数，即(mu(a)mu(b)=mu(a cdot b))
-
(sum_{d|n} mu(d)=[n=1]) ,这一点根据二项式定理即可证明
２．线性筛莫比乌斯函数
需要用到的性质
根据 (mu) 的性质可以得到，如果 (n=p*m) ，其中 (p) 为 (n) 的质因子，如果 (p|m) 那么显然 (mu(n)=0) ，否则 (mu(n)=-mu(m))
miu[1]=1;
For (i,2,n) {
if (!vis[i]) {
p[++N]=i;
miu[i]=-1;
}
For (j,1,N) {
if (i*p[j]>n) break;
vis[i*p[j]]=1;
if (i%p[j]==0) {
miu[i*p[j]]=0;
break;
}
else
miu[i*p[j]]=-miu[i];
}
}
（二）莫比乌斯反演
如果 (f(n)) , (g(n)) 为数论函数，并且满足$$f(n)=sum_{d|n}g(d)$$
则有莫比乌斯反演$$g(n)=sum_{d|n}mu(d)*f(frac{n}{d})$$
证明在上面的博客中写得很清楚。
但是那个变形式的证明我一直没有弄得很明白，到底如何交换变量？
１．例题一
(bzoj2190)
求
[sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j)=k]
]
[=sum_{i=1}^{leftlfloor frac{n}{k}
ight
floor}sum_{j=1}^{leftlfloor frac{m}{k}
ight
floor}[gcd(i,j)=1]
]
因为
[sum_{d|n}^{n}mu(d)=[n=1]
]
所以，可以化为
[=sum_{i=1}^{leftlfloor frac{n}{k}
ight
floor}sum_{j=1}^{leftlfloor frac{m}{k}
ight
floor}sum_{d|gcd(i,j)}mu(d)
]
改变枚举量，可以化为
[sum_{d=1}^{leftlfloorfrac{n}{k}
ight
floor} mu(d) leftlfloorfrac{n}{kd}
ight
floor leftlfloorfrac{m}{kd}
ight
floor
]
如果是暴力就是(O(n))的，可以用分块优化到(sqrt{n})
因为有很多取值是连续的，对于相等的段，我们求取(mu)的前缀和，即可批量计算这一个段的答案。
ll solve(int n,int m,int k) {
if (n>m) swap(n,m);
n/=k; m/=k;
ll ans=0;
for (int i=1,nxt=1;i<=n;i=nxt+1) {
nxt=min(n/(n/i),m/(m/i));
ans+=1ll*(sum[nxt]-sum[i-1])*(n/i)*(m/i);
}
return ans;
}
２．例题二
求
[sum_{i=1}^{n}sum_{j=1}^{m} gcd(i,j)^k
]
一个重要的方法枚举 (dx) 的乘积
[sum_{T=1}^{n} leftlfloorfrac{n}{T}
ight
floor leftlfloorfrac{m}{T}
ight
floor sum_{d|T}d^{k}mu(frac{T}{d})
]
然后设
[f(T)=sum_{d|T}d^{k}mu(frac{T}{d})
]
再根据积性函数的性质化简即可
具体的过程，上面博客里有，下面对其中的几点做一下说明
- 分块套分块的方法首先是分块(leftlfloor frac{n}{d} ight floor)，因为有很多是一样的，然后再分块 (mu) 后面的那两个
- 后面化简(f(n))的部分因为只有在 (d) 取到 (p_{i}^{x^{i}}) 或者 (p_{i}^{x^{i}-1}) 的时候 ((frac{T}{d})) 取到 (1) 或者 (p_{i}) 的时候 (mu) 值不为０。然后线筛积性函数。
注意，这道题目的解答过程对后面的题目很有帮助，一定要全部弄明白。
他推理的过程中有很多的变量都写错了，但是大体公式还是没问题的。
三．题目
（一）(bzoj2820)
求
[sum_{i=1}^{n}sum_{j=1}^{m}[gcd(i,j) isprime]
]
分析：
枚举质数
[=sum_psum_{i=1}^{leftlfloor frac{n}{p}
ight
floor}sum_{j=1}^{leftlfloor frac{m}{p}
ight
floor}[gcd(i,j)=1]
]
[=sum_{p}sum_{d=1}^{leftlfloor frac{n}{p}
ight
floor}mu(d)leftlfloor frac{n}{pd}
ight
floorleftlfloor frac{m}{pd}
ight
floor
]
然后枚举乘积(pd)
[=sum_{T=1}leftlfloor frac{n}{T}
ight
floor leftlfloor frac{n}{T}
ight
floorsum_{p|T}mu(frac{T}{p})
]
但是，如果设
[f(T)=sum_{p|T}mu(frac{T}{p})
]
(f(T))不是积性函数
所以需要用到一个新的方法，在预处理的时候，枚举质数 (p) ，枚举 (T) ，然后将 (p) 的每一个倍数 (T) 都加上 (mu(frac{T}{p})) ，然后求 (f(n)) 就可以做到 (O(1))
再加上前面的分块求 (T) ，因此总时间复杂度为 (O(Tsqrt{n}+maxn))
当然，还有一种方法，详情看here
最后再上代码
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
typedef long long ll;
typedef double dd;
#define For(i,j,k) for (int i=j;i<=k;++i)
#define Forr(i,j,k) for (int i=j;i>=k;--i)
#define Set(a,p) memset(a,p,sizeof(a))
using namespace std;
template<typename T>bool chkmax(T &a,T b) { return a<b?a=b,1:0; }
template<typename T>bool chkmin(T &a,T b) { return a>b?a=b,1:0; }
const int maxn=1e7+1e2;
const int maxx=1e7;
int n,m,N;
int p[maxn],miu[maxn],vis[maxn],sum[maxn],t[maxn];
inline void read(int &x) {
x=0;
int p=1;
char c=getchar();
while (!isdigit(c)) {if (c=='-') p=-1; c=getchar();}
while (isdigit(c)) {x=(x<<1)+(x<<3)+(c-'0'); c=getchar();}
x*=p;
}
inline void init() {
miu[1]=1;
For (i,2,maxx) {
if (!vis[i]) {
p[++N]=i; miu[i]=-1;
}
For (j,1,N) {
if (1ll*i*p[j]>maxx) break;
vis[i*p[j]]=1;
if (i%p[j]==0) {
miu[i*p[j]]=0;
break;
}
else miu[i*p[j]]=-miu[i];
}
}
For (i,1,N) {
for (int j=p[i];j<=maxx;j+=p[i]) t[j]+=miu[j/p[i]];
}
For (i,1,maxx) t[i]=t[i-1]+t[i];
}
int main() {
init();
int tt; read(tt);
while (tt--) {
read(n); read(m);
if (n>m) swap(n,m);
ll ans=0;
for (int i=1,nxt=1;i<=n;i=nxt+1) {
nxt=min(n/(n/i),m/(m/i));
ans+=1ll*(n/i)*(m/i)*(t[nxt]-t[i-1]);
}
printf("%lld
",ans);
}
return 0;
}
（二）(P3768)
求
[sum_{i=1}^{n}sum_{j=1}^{n}ijgcd(i,j)
]
[=sum_{d=1}^{n}dsum_{i=1}^{n}sum_{j=1}^{n}ij[gcd(i,j)=d]
]
[=sum_{d=1}^{n}d^3sum_{i=1}^{leftlfloor frac{n}{d}
ight
floor}sum_{j=1}^{leftlfloor frac{n}{d}
ight
floor}ij[gcd(i,j)=1]
]
[=sum_{d=1}^{n}d^3sum_{i=1}^{leftlfloor frac{n}{d}
ight
floor}sum_{j=1}^{leftlfloor frac{n}{d}
ight
floor}ijsum_{x|gcd(i,j)}mu(x)
]
[=sum_{d=1}^{n}d^3sum_{x=1}^{leftlfloor frac{n}{d}
ight
floor}mu(x)x^2sum_{i=1}^{leftlfloor frac{n}{dx}
ight
floor}sum_{j=1}^{leftlfloor frac{n}{dx}
ight
floor}ij
]
令(F(n)=sum_{i=1}^{n}sum_{j=1}^{n}ij=sum_{i=1}^{n}i^3=(sum_{i=1}^{n})^2)，所以
[=sum_{d=1}^{n}d^3sum_{x=1}^{leftlfloor frac{n}{d}
ight
floor}mu(x)x^2F(leftlfloor frac{n}{dx}
ight
floor)
]
令(T=dx)，改变枚举变量（干脆把向下取整去掉，反正程序中是这样的，不然写起来好麻烦）
[=sum_{T=1}^{n}sum_{d|T}d^3mu(frac{T}{d})(frac{T}{d})^2F(frac{n}{T})
]
[=sum_{T=1}^{n}F(frac{n}{T})T^2sum_{d|T}dmu(frac{T}{d})
]
然后，根据(mu)函数的性质，我们可以得到
[id*mu=1*varphi*mu=varepsilon*varphi=varphi
]
所以原式可化为
[=sum_{T=1}^{n}F(frac{n}{T})T^2varphi(T)
]
这看着就清爽多了！前半部分分块处理，后面杜教筛
设(f(i)=i^2*varphi(i))，(S(n)=sum_{i=1}^{n}f(i))
杜教筛公式：
[g(1)*S(n)=sum_{i=1}^{n}(f*g)(i)-sum_{i=2}^{n}g(i)*S(frac{n}{i})
]
设(g(x)=x^2)，那么
[(f*g)(i)=sum_{d|i}f(d)*g(frac{i}{d})=sum_{d|i}d^2*varphi(d)*(frac{i}{d})^2=i^2*sum_{d|i}varphi(d)=i^3
]
所以
[S(n)=sum_{i=1}^{n}i^3-sum_{i=2}^{n}i^2*S(frac{n}{i})
]
再然后，就可以开始打程序了！
当然，事实上是，用 (varphi) 更简单here
[=sum_{d=1}^{n}varphi(d)sum_{d|i}sum_{d|j}ij ][=sum_{d=1}^{n}varphi(d) d^2 sum_{i=1}^{leftlfloor frac{n}{d} ight floor}i^3 ]
我这篇代码被卡常了，对着别人的程序调了好久，但是还是很慢，不知道为什么......
在计算(sum)和(sq)是时候，(n)一进去就要取模！
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<map>
typedef double dd;
#define For(i,j,k) for (register ll i=j;i<=k;++i)
#define Forr(i,j,k) for (register ll i=j;i>=k;--i)
#define REP(i,j,k) for (register int i=j;i<=k;++i)
#define Set(a,p) memset(a,p,sizeof(a))
#define ll long long
using namespace std;
template<typename T>bool chkmax(T& a,T b) {return a<b?a=b,1:0;}
template<typename T>bool chkmin(T& a,T b) {return a>b?a=b,1:0;}
#define maxn 8000000+100
int maxx=8000000;
ll n,modd,Inv2,Inv6,ans;
int N,p[maxn];
ll phi[maxn];
bool vis[maxn];
map<ll,ll>Phi;
inline void init() {
phi[1]=1;
for (register int i=2;i<=maxx;++i) {
if (!vis[i]) {
p[++N]=i; phi[i]=i-1;
}
for (register int j=1;j<=N;++j) {
if (i*p[j]>maxx) break;
int k=i*p[j];
vis[k]=1;
if (i%p[j]==0) {
phi[k]=phi[i]*p[j]%modd;
break;
}
else phi[k]=phi[i]*(p[j]-1)%modd;
}
}
for (register int i=1;i<=maxx;++i) phi[i]=(phi[i-1]+phi[i]*i%modd*i%modd)%modd;
}
inline ll quick(ll a,ll b) {
ll s=1;
while (b) {
if (b&1) s=s*a%modd;
a=a*a%modd; b>>=1;
}
return s;
}
inline ll sum(ll n) {
n%=modd;
return n*(n+1)%modd*Inv2%modd;
}
inline ll sq(ll n) {
n%=modd;
return n*(n+1)%modd*(n+n+1)%modd*Inv6%modd;
}
ll calc(ll n) {
if (n<=N) return phi[n];
if (Phi[n]) return Phi[n];
ll s=sum(n); s=s*s%modd;
for (ll i=2,nxt;i<=n;i=nxt+1) {
nxt=n/(n/i);
s-=(sq(nxt)-sq(i-1))*calc(n/i)%modd; s%=modd;
}
Phi[n]=(s+modd)%modd;
return (s+modd)%modd;
}
int main() {
scanf("%lld%lld",&modd,&n);
Inv2=quick(2,modd-2); Inv6=quick(6,modd-2);
init();
for (ll i=1,nxt=1;i<=n;i=nxt+1) {
nxt=n/(n/i);
ll ss=sum(n/i); ss=ss*ss%modd;
ans+=ss*(calc(nxt)-calc(i-1))%modd;
ans=(ans+modd)%modd;
}
printf("%lld
",(ans+modd)%modd);
return 0;
}
（三）(bzoj2693)，(P4313)
求
[sum_{i=1}^{n}sum_{i=1}^{m}lcm(i,j)
]
ymy大佬的博客中推导过程和我一样的，所以我就懒得写了
[=sumlimits_{i=1}^nsumlimits_{j=1}^m frac{ij}{gcd(i,j)}
]
[=sumlimits_{d=1}^n frac{1}{d}sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} ij[gcd(i,j)=d]
=sumlimits_{d=1}^n dsumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor} ij[gcd(i,j)=1]]
[=sumlimits_{d=1}^n dsumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor} ijsumlimits_{x|gcd(i,j)}mu(x)
=sumlimits_{d=1}^n dsumlimits_{x=1}^{lfloorfrac{n}{d}
floor}mu(x)x^2sumlimits_{i=1}^{lfloorfrac{n}{dx}
floor}sumlimits_{j=1}^{lfloorfrac{m}{dx}
floor} ij]
设(f(T)=sumlimits_{i=1}^{lfloorfrac{n}{T} floor}sumlimits_{j=1}^{lfloorfrac{m}{T} floor} ij=sumlimits_{i=1}^{lfloorfrac{n}{T} floor}isumlimits_{j=1}^{lfloorfrac{m}{T} floor}j)，所以
[=sumlimits_{T=1}^n f(T)sumlimits_{d|T}mu(d)*d^2frac{T}{d}
=sumlimits_{T=1}^n f(T)*Tsumlimits_{d|T}mu(d)*d]
设(g(n)=sumlimits_{d|n}mu(d)*d) ,则(g) 为积性函数，可以线筛。
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
typedef long long ll;
typedef double dd;
#define For(i,j,k) for (int i=j;i<=k;++i)
#define Forr(i,j,k) for (int i=j;i>=k;--i)
#define Set(a,p) memset(a,p,sizeof(a))
using namespace std;
template<typename T>bool chkmax(T& a,T b) {return a<b?a=b,1:0;}
template<typename T>bool chkmin(T& a,T b) {return a>b?a=b,1:0;}
const int maxn=10000000+100;
const int maxx=10000000;
const ll modd=100000009;
int n,m,N;
int vis[maxn],p[maxn];
ll inv;
ll f[maxn],sum[maxn];
inline void read(int &x) {
x=0;
int p=1;
char c=getchar();
while (!isdigit(c)) {if (c=='-') p=-1; c=getchar();}
while (isdigit(c)) {x=(x<<1)+(x<<3)+(c-'0'); c=getchar();}
x*=p;
}
inline void init() {
f[1]=1;
For (i,2,maxx) {
if (!vis[i]) {
p[++N]=i; f[i]=1-i;
}
For (j,1,N) {
int k=i*p[j];
if (k>maxx) break;
vis[k]=1;
if (i%p[j]==0) {
f[k]=f[i]; break;
}
else f[k]=f[i]*(1-p[j])%modd;
}
}
For (i,1,maxx) f[i]=(f[i-1]+f[i]*i%modd+modd)%modd;
For (i,1,maxx) {
sum[i]=sum[i-1]+i;
if (sum[i]>modd) sum[i]-=modd;
}
}
int main() {
init();
int tt; read(tt);
while (tt--) {
ll ans=0;
read(n); read(m);
if (n>m) swap(n,m);
for (int i=1,nxt;i<=n;i=nxt+1) {
nxt=min(n/(n/i),m/(m/i));
ans+=sum[n/i]*sum[m/i]%modd*(f[nxt]-f[i-1])%modd;
ans%=modd;
}
printf("%lld
",(ans+modd)%modd);
}
return 0;
}
（四）(bzoj3994)
设 (d(n)) 为 (n) 的约数个数，求
[sum_{i=1}^{n}sum_{j=1}^{m}d(ij)
]
分析：
[=sum_{i=1}^{n}sum_{j=1}^{m}sum_{p|i}sum_{q|j}[gcd(p,q)=1]
]
[=sum_{p=1}^{n}sum_{q=1}^{m}leftlfloor frac{n}{p}
ight
floor leftlfloor frac{m}{q}
ight
floor[gcd(p,q)=1]
]
[=sum_{i=1}^{n}sum_{j=1}^{m}leftlfloor frac{n}{i}
ight
floor leftlfloor frac{m}{j}
ight
floor[gcd(i,j)=1]
]
[=sum_{d=1}^{n}mu(d)sum_{i=1}^{leftlfloor frac{n}{d}
ight
floor}sum_{j=1}^{leftlfloor frac{m}{d}
ight
floor}leftlfloor frac{n}{id}
ight
floorleftlfloor frac{m}{id}
ight
floor
]
然后，一个智障的变化，我一下子还没想出来
[=sum_{d=1}^{n}mu(d)sum_{i=1}^{leftlfloor frac{n}{d}
ight
floor}leftlfloor frac{n}{id}
ight
floorsum_{j=1}^{leftlfloor frac{m}{d}
ight
floor}leftlfloor frac{m}{id}
ight
floor
]
设
[f(n)=sum_{i=1}^{n}leftlfloor frac{n}{i}
ight
floor
]
这个显然可以分块，然后就可以化为
[=sum_{d=1}^{n}mu(d)f(leftlfloor frac{n}{d}
ight
floor)f(leftlfloor frac{m}{d}
ight
floor)
]
后面的在 (O(nsqrt{n})) 预处理之后可以 (O(1)) 计算，前面的分块 (O(sqrt{n}))
总时间复杂度 (O(nsqrt{n}+Tsqrt{n}))
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
typedef long long ll;
typedef double dd;
#define For(i,j,k) for (register int i=j;i<=k;++i)
#define Forr(i,j,k) for (register int i=j;i>=k;--i)
#define Set(a,p) memset(a,p,sizeof(a))
using namespace std;
template<typename T>bool chkmax(T &a,T b) { return a<b?a=b,1:0; }
template<typename T>bool chkmin(T &a,T b) { return a>b?a=b,1:0; }
const int maxn=50000+100;
const int maxx=50000;
int n,m,N;
int p[maxn],vis[maxn],mu[maxn],sum[maxn],s[maxn];
inline void read(int &x) {
x=0;
int p=1;
char c=getchar();
while (!isdigit(c)) {if (c=='-') p=-1; c=getchar();}
while (isdigit(c)) {x=(x<<1)+(x<<3)+(c-'0'); c=getchar();}
x*=p;
}
inline void init() {
mu[1]=1;
For (i,2,maxx) {
if (!vis[i]) {
p[++N]=i; mu[i]=-1;
}
For (j,1,N) {
if (i*p[j]>maxx) break;
int k=i*p[j];
vis[k]=1;
if (i%p[j]==0) {
mu[k]=0; break;
}
else mu[k]=-mu[i];
}
}
For (i,1,maxx) {
sum[i]=sum[i-1]+mu[i];
for (int p=1,nxt=1;p<=i;p=nxt+1) {
nxt=i/(i/p);
s[i]+=(nxt-p+1)*(i/p);
}
}
}
inline void solve() {
if (n>m) swap(n,m);
ll ans=0;
for (int i=1,nxt=1;i<=n;i=nxt+1) {
nxt=min(n/(n/i),m/(m/i));
ans+=1ll*s[n/i]*s[m/i]*(sum[nxt]-sum[i-1]);
}
printf("%lld
",ans);
}
int main() {
init();
int tt; read(tt);
while (tt--) {
read(n); read(m);
solve();
}
return 0;
}