莫比乌斯反演
以下内容中，\(a|b\)表示\(a\)是\(b\)的因数，\(a\not|b\)表示\(a\)不是\(b\)的因数
前置知识：积性函数相关
狄利克雷卷积
定义：任意函数\(f(n),g(n)\)的狄利克雷卷积为\((f*g)(n)=\sum\limits_{d|n}f(d)g(\dfrac{n}{d})\)
其中\(\sum\limits_{d|n}\)是要枚举\(n\)的因数并累加
换一种写法就是\(\sum\limits_{xy=n}f(x)g(y)\)
还可以用代码写成：
int dirichlet(int n) {
int ans=0;
for(int i=1; i<=n; ++i) {
if(n%i==0) {
ans+=f[i]*g[n/i];//f[i]为f(i),g[n/i]为g(n/i)
}
}
return ans;
}
这个卷积满足：
- 交换律\(f*g=g*f\)
- 结合律\(f*g*h=f*(g*h)\)
- 分配率\(f*(g+h)=f*g+f*h\)
- 当\(f(x),g(x)\)都为积性函数时，\((f*g)(x)\)为积性函数
- 当\((f*g)(x),g(x)\)都为积性函数时，\(f(x)\)为积性函数
常见的狄利克雷卷积有：
- \(\varphi*I=id\)
- \(\mu*I=\epsilon\)
- \(\epsilon*f=f\)
- \(id*I=\sigma\)
- \(I*I=d\)
- \((\lambda*I)(n)=[\sqrt{n}\in Z]\)
简单证明：（3,4,5见相关函数的定义：积性函数相关）
- \(\varphi*I=id\)（下面以\(n=6\)为例）
将\([1,n]\)的整数都除以\(n\)：
\(1\div6=\dfrac{1}{6},5\div6=\dfrac{5}{6}\)，对应与\(6\)互质的数，个数为\(\varphi(6)\)
\(2\div6=\dfrac{1}{3},4\div6=\dfrac{2}{3}\)，对应与\(3\)互质的数，个数为\(\varphi(3)\)
\(3\div6=\dfrac{1}{2}\)，对应与\(2\)互质的数，个数为\(\varphi(2)\)
\(6\div6=\dfrac{1}{1}\)，对应与\(1\)互质的数，个数为\(\varphi(1)\)
即\(\varphi(1)+\varphi(2)+\varphi(3)+\varphi(6)=6\)（\(\sum\limits_{d|6}\varphi(d)=6\)）
- \(\mu*I=\epsilon\)
设\(n\)质因数分解后为\(\prod\limits_{i=1}^kp_i^{c_i}\)
当\(n>1\)时：
\(\sum\limits_{d|n}\mu(d)\\=\mu(1)+\mu(p_1)+\mu(p_1^2)+\cdots+\mu(p_i^{c_i})+\cdots+\mu(p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k})\\=\mu(1)+\mu(p_1)+\mu(p_1^2)+\cdots+\mu(p_i^{c_i})+\cdots+\mu(p_1^{c_1})\mu(p_2^{c_2})\cdots \mu(p_k^{c_k})\\=\prod\limits_{i=1}^k\sum\limits_{j=0}^{c_i}\mu(p_i^j)\\=\prod\limits_{i=1}^k\sum\limits_{j=0}^{1}\mu(p_i^j)\\=\prod\limits_{i=1}^k(1-1)\\=0\)
当\(n=1\)时，\(\sum\limits_{d|n}\mu(d)=\mu(1)=1\)
- \((\lambda*I)(n)=[\sqrt{n}\in Z]\)
同样设\(n\)质因数分解后为\(\prod\limits_{i=1}^kp_i^{c_i}\)
\(\sum\limits_{d|n}\lambda(d)\\=\lambda(1)+\lambda(p_1)+\lambda(p_1^2)+\cdots+\lambda(p_i^{c_i})+\cdots+\lambda(p_1^{c_1}p_2^{c_2}\cdots p_k^{c_k})\\=\lambda(1)+\lambda(p_1)+\lambda(p_1^2)+\cdots+\lambda(p_i^{c_i})+\cdots+\lambda(p_1^{c_1})\lambda(p_2^{c_2})\cdots \lambda(p_k^{c_k})\\=\prod\limits_{i=1}^k\sum\limits_{j=0}^{c_i}\lambda(p_i^j)\\=\prod\limits_{i=1}^k\sum\limits_{j=0}^{c_i}(-1)^j\\=\prod\limits_{i=1}^k[c_i\equiv0\pmod2]\)
最后那个式子只有在\(c_i\)都为偶数的时候才为\(1\)，所以\(n\)为完全平方数时才是\(1\)。
一道题
给出\(n(n\leq10^{12})\)，求\(\sum\limits_{k=1}^n\sum\limits_{i|k}\sum\limits_{j|i}\lambda(i)\lambda(j)\pmod{998244353}\)
我们可以先枚举\(i\)，再枚举满足\(i|k\)的\(k\)：
\(\sum\limits_{i=1}^n\sum\limits_{i|k}\sum\limits_{j|i}\lambda(i)\lambda(j)\)
可以提出公因式：
\(\sum\limits_{i=1}^n\lambda(i)(\sum\limits_{i|k}1)\sum\limits_{j|i}\lambda(j)\)
中间的部分就是求小于等于\(n\)的\(i\)的倍数有多少，也就是\(\lfloor\dfrac{n}{i}\rfloor\)
\(\sum\limits_{i=1}^n\lambda(i)\lfloor\dfrac{n}{i}\rfloor\sum\limits_{j|i}\lambda(j)\)
后面的部分是\(\lambda*I\)，上面说了，是\([\sqrt{n}\in Z]\)
\(\sum\limits_{i=1}^n\lambda(i)\lfloor\dfrac{n}{i}\rfloor[\sqrt{i}\in Z]\)
这样我们就只用考虑完全平方数了。
而且当\(i\)为完全平方数时\(\lambda(i)=1\)
\(\sum\limits_{i=1}^{\sqrt{n}}\lfloor\dfrac{n}{i^2}\rfloor\)
\(O(\sqrt{n})\)可以通过。
再来一题
先把答案用式子表示一下：
\(\sum\limits_{i=1}^a\sum\limits_{j=1}^b[\gcd(i,j)=d]\)
之后把\(i,j\)同时除以\(d\)：
\(\sum\limits_{i=1}^{\lfloor\frac{a}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{d}\rfloor}[\gcd(i,j)=1]\)
\([\gcd(i,j)=1]\)就是\(\epsilon(\gcd(i,j))\)，把\(\epsilon=\mu*1\)代入上式得：
\(\sum\limits_{i=1}^{\lfloor\frac{a}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{d}\rfloor}\sum\limits_{k|\gcd(i,j)}\mu(k)\)
把最后一个循环提前：
\(\sum\limits_{k=1}^{\lfloor\frac{a}{d}\rfloor}\mu(k)\sum\limits_{i=1}^{\lfloor\frac{a}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{d}\rfloor}[k|\gcd(i,j)]\)
\(k|\gcd(i,j)\)其实就是\(k|i\&k|j\)，代入上式
\(\sum\limits_{k=1}^{\lfloor\frac{a}{d}\rfloor}\mu(k)\sum\limits_{i=1}^{\lfloor\frac{a}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{d}\rfloor}[k|i\&k|j]\)
\(=\sum\limits_{k=1}^{\lfloor\frac{a}{d}\rfloor}\mu(k)\sum\limits_{i=1}^{\lfloor\frac{a}{dk}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{b}{dk}\rfloor}\)
\(=\sum\limits_{k=1}^{\lfloor\frac{a}{d}\rfloor}\mu(k)\lfloor\frac{a}{dk}\rfloor\lfloor\frac{b}{dk}\rfloor\)
线性筛\(\mu\)的前缀和并用整除分块就可以达到每次询问\(O(\sqrt{a})\)，总复杂度\(O(T\sqrt{a})\)
#include<bits/stdc++.h>
using namespace std;
const int MAXN=50005;
typedef long long ll;
bool isp[MAXN];
ll sum[MAXN];
int prime[MAXN],pcnt,mu[MAXN];
ll f(int n,int m) {
ll ans=0;
if(n>m)swap(n,m);
for(int l=1,r; l<=n; l=r+1) {
r=min(n/(n/l),m/(m/l));
ans+=(sum[r]-sum[l-1])*1ll*(n/l)*(m/l);
}
return ans;
}
int main() {
sum[1]=mu[1]=1;
for(int i=2; i<MAXN; ++i) {
if(!isp[i]) {
prime[++pcnt]=i;
mu[i]=-1;
}
for(int j=1,x; j<=pcnt&&(x=i*prime[j])<MAXN; ++j) {
isp[x]=true;
if(i%prime[j]) {
mu[x]=-mu[i];
} else {
mu[x]=0;
break;
}
}
sum[i]=sum[i-1]+mu[i];
}
int t,a,b,d;
scanf("%d",&t);
while(t--) {
scanf("%d%d%d",&a,&b,&d);
printf("%lld\n",f(a/d,b/d));
}
return 0;
}
莫比乌斯反演
第一种形式：对于\(f*I=g\)，有\(f=g*\mu\)，反之亦然
证明：
两边都与\(\mu\)做卷积
\(f*I*\mu=g*\mu\)
结合律：
\(f*(I*\mu)=g*\mu\)
由\(\epsilon=1*\mu\)得：
\(f*\epsilon=g*\mu\)
由\(f*\epsilon=f\)得：
\(f=g*\mu\)
第二种形式：\(\sum\limits_{n|d}f(d)=g(n)\)，则\(f(n)=\sum\limits_{n|d}g(d)\mu(\dfrac{d}{n})\)
证明：
\(\sum\limits_{n|d}g(d)\mu(\dfrac{d}{n})\)
将\(d\)除以\(n\)：
\(=\sum\limits_{d=1}^{+\infty}g(dn)\mu(d)\)
把\(g\)写成\(f\)
\(=\sum\limits_{d=1}^{+\infty}\mu(d)\sum\limits_{nd|t}f(t)\)
枚举\(t\)，\(d\)即是\(\dfrac{t}{n}\)的因数
\(=\sum\limits_{n|t}f(t)\sum\limits_{d|\frac{t}{n}}\mu(d)\)
后面的是\(\mu*I\)的形式，改写成\(\epsilon\)
\(=\sum\limits_{n|t}f(t)[\dfrac{t}{n}=1]=f(n)\)
来几道题：
三倍经验的题：
SP3871 GCDEX - GCD Extreme，UVA11426 拿行李（极限版） GCD - Extreme (II)，UVA11424 GCD - Extreme (I)
我们换个思路，先写成这个形式：
\(\sum\limits_{i=1}^n\sum\limits_{j=1}^{i-1}\gcd(i,j)\)
根据\(P2303\)的结论\(\sum\limits_{i=1}^n\gcd(i,n)=\sum\limits_{d|n}d\varphi(\dfrac{n}{d})\)（上面有）把式子写成：
\(\sum\limits_{i=1}^n((\sum\limits_{j=1}^{i}\gcd(i,j))-\gcd(i,i))=\sum\limits_{i=1}^n(\sum\limits_{d|i}d\varphi(\dfrac{i}{d}))-i=(\sum\limits_{i=1}^n\sum\limits_{d|i}d\varphi(\dfrac{i}{d}))-\dfrac{n(n+1)}{2}\)
设\(g(i)=\sum\limits_{d|i}d\varphi(\dfrac{i}{d})\)
我们可以试着\(O(n)\)线性筛出\(g(i)\)，同时算个前缀和就预处理完了。
下面写的\(i,prime_j\)为欧拉筛代码中的\(i,prime[j]\)
我们设\(i\)分解质因数后是\(\prod\limits_{j=1}^kp_j^{c_j}\)，\(low_i=p_1^{c_1}\)。其中\(p_1\)是最小的质因子。
而且\(prime_j\)为\(i\times prime_j\)得最小质因子（欧拉筛的性质）
则当\(prime_j|i\)时，\(p_1=prime_j\)
可以得到，\(g(1)=1,g(p)=2p-1\)，\(p\)是质数。
由于函数\(id(x)\)和\(\varphi(x)\)是积性函数
那么\(g(i)=id*\varphi\)也是积性函数
则当\(i\)与\(prime_j\)互质时，\(g(i\times prime_j)=g(i)\times g(prime_j)\)，且\(low_{i\times prime_j}=prime_j\)
当\(i\)与\(prime_j\)不互质即\(prime_j|i\)时，有\(p_1=prime_j\)（上面讲了），根据\(low\)的定义得
\(\gcd(\dfrac{i}{low_i},p_1\times low_i)=1,low_{i\times prime_j}=low_i\times prime_j\)
这样就可以用积性函数的定义求\(g(i\times prime_j)\)
\(g(i\times prime_j)=g(\dfrac{i}{low_i})\times g(p_1\times low_i)\)
但是当\(i=low_i\)时，这个式子就会变成\(g(i\times prime_j)=g(1)\times g(i\times prime_j)=g(i\times prime_j)\)，无法递推
特殊处理一下\(i=low_i\)的情况：
\(i=low_i\)即\(i=p_1^{c_1}\)，则\(i\times prime_j=p_1^{c_1+1}\)
\(g(p_1^{c_1+1})=\sum\limits_{d|p_1^{c_1+1}}d\varphi(\dfrac{p_1^{c_1+1}}{d})\)
把\(d=1\)的一项单独提出：
\(g(p_1^{c_1+1})=\sum\limits_{d|p_1^{c_1+1}\&d>1}d\varphi(\dfrac{p_1^{c_1+1}}{d})+\varphi(p_1^{c_1+1})\)
\(\sum\)内的\(d\)此时只能是\(p_1,p_1^2,\cdots,p_1^{c_1+1}\)，我们就可以把\(d\)除以\(p_1\)
\(g(p_1^{c_1+1})=\sum\limits_{d|p_1^{c_1}}dp_1\varphi(\dfrac{p_1^{c_1+1}}{dp_1})+\varphi(p_1^{c_1+1})\)
化简一下，并且提出\(p_1\)
\(g(p_1^{c_1+1})=p_1\sum\limits_{d|p_1^{c_1}}d\varphi(\dfrac{p_1^{c_1}}d)+\varphi(p_1^{c_1+1})\)
可以发现\(\sum\limits_{d|p_1^{c_1}}d\varphi(\dfrac{p_1^{c_1}}d)=g(p_1^{c_1})\)，代入
\(g(p_1^{c_1+1})=p_1g(p_1^{c_1})+\varphi(p_1^{c_1+1})\)
也就是\(g(i\times prime_j)=prime_j\times g(i)+\varphi(i\times prime_j)\)
总结一下\(g(x)=\begin{cases}1\quad x=1\\2x-1\quad x\in prime\\g(i)\times g(prime_j)\quad x=i\times prime_j,\gcd(prime_j,i)=1\\prime_j\times g(i)+\varphi(i\times prime_j)\quad x=i\times prime_j,prime_j|i,low_i=i\\g(i\times prime_j)=g(\dfrac{i}{low_i})\times g(p_1\times low_i)\quad x=i\times prime_j,prime_j|i,low_i\not=i\end{cases}\)
\(low_x=\begin{cases}0\quad i=1\\x\quad x\in prime\\prime_j\quad x=i\times prime_j,\gcd(i,prime_j)=1\\low_i\times prime_j\quad x=i\times prime_j,prime_j|i\end{cases}\)
代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN=4e6+100;
int prime[MAXN],pcnt,phi[MAXN+1],n;
ll g[MAXN+1],low[MAXN+1],sum[MAXN+1];
bool isp[MAXN+1];
int main() {
sum[1]=g[1]=phi[1]=1;
for(int i=2; i<=MAXN; ++i) {//筛出phi,g,low
if(!isp[i]) {
prime[++pcnt]=i;
phi[i]=i-1;
low[i]=i;
g[i]=(i<<1)-1;//素数
}
for(int j=1,x; j<=pcnt&&(x=i*prime[j])<=MAXN; ++j) {
isp[x]=true;
if(i%prime[j]) {
g[x]=g[i]*g[prime[j]];
phi[x]=phi[i]*phi[prime[j]];
low[x]=prime[j];//三个都要算上
} else {
low[x]=low[i]*prime[j];
phi[x]=phi[i]*prime[j];
g[x]=(low[i]==i)?g[i]*prime[j]+phi[x]:g[i/low[i]]*g[low[i]*prime[j]];
}//记得特判i=low i
}
sum[i]=sum[i-1]+g[i];//前缀和
}
while(~scanf("%d",&n)&&n) {
printf("%lld\n",sum[n]-((n+1ll)*n>>1));//输出
}
return 0;
}
这题是要求：\(\sum\limits_{p\in prime}\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=p]\)
把\(i,j\)除以\(p\)
\(\sum\limits_{p\in prime}\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor}[\gcd(i,j)=1]\)
可以对后面的部分进行莫比乌斯反演
设\(f(x)\)表示\(\gcd(i,j)\)为\(x\)的二元组的个数，\(g(x)\)表示\(\gcd(i,j)\)为\(x\)的倍数的二元组的个数，
把它俩用式子写就是：
\(f(x)=\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor}[\gcd(i,j)=x]\)
\(g(x)=\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor}[x|\gcd(i,j)]=\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor}\sum\limits_{x|d}[\gcd(i,j)=d]\)
后面的\(\sum\limits_{x|d}\)可以提前
\(g(x)=\sum\limits_{x|d}\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor}[\gcd(i,j)=d]=\sum\limits_{x|d}f(x)\)
根据莫比乌斯反演第二形式，可得：\(f(x)=\sum\limits_{x|d}g(d)\mu(\dfrac{d}{x})\)
\(f(1)=\sum\limits_{1|d}g(d)\mu(\dfrac{d}{1})=\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}g(i)\mu(i)\)
再看\(g(x)=\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{p}\rfloor}[x|\gcd(i,j)]\)这个式子，把\(i,j\)除以\(x\)可以得到
\(g(x)=\sum\limits_{i=1}^{\lfloor\frac{n}{px}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{px}\rfloor}[1|\gcd(i,j)]=\sum\limits_{i=1}^{\lfloor\frac{n}{px}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{px}\rfloor}1=\lfloor\frac{n}{px}\rfloor\lfloor\frac{m}{px}\rfloor\)
代回\(f(1)=\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}g(i)\mu(i)=\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\lfloor\frac{n}{px}\rfloor\lfloor\frac{m}{px}\rfloor\mu(i)\)
提出\(\lfloor\frac{n}{px}\rfloor\lfloor\frac{m}{px}\rfloor\)
\(f(1)=\lfloor\frac{n}{px}\rfloor\lfloor\frac{m}{px}\rfloor\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\mu(i)\)
再代回原来的式子：
\(\sum\limits_{p\in prime}\lfloor\frac{n}{px}\rfloor\lfloor\frac{m}{px}\rfloor\sum\limits_{i=1}^{\lfloor\frac{n}{p}\rfloor}\mu(i)\)
先用\(T\)枚举\(px\)，再枚举\(p|T\)
\(\sum\limits_{i=1}^T\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum\limits_{p|T\& p\in prime}\mu(i)\)
由于\(10000128\)以内的素数只有\(664583\)个，所以我们暴力计算后面的部分都不会超时。
前面的部分就可以用整除分块优化
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN=10000129;
int prime[MAXN>>3],pcnt;
bool vis[MAXN];
int mu[MAXN];
ll sum[MAXN],g[MAXN];
int read() {
int x=0,c=getchar();
bool f=1;
while(c<=47||c>=58) {
f=f&&(c^45);
c=getchar();
}
while(c>=48&&c<=57) {
x=(x<<3)+(x<<1)+(c&15);
c=getchar();
}
return f?x:-x;
}
int main() {
vis[0]=vis[1]=1;
mu[1]=1;
for(int i=2; i<MAXN; ++i) {//预处理mu
if(!vis[i]) {
prime[++pcnt]=i;
mu[i]=-1;
}
for(int j=1; i*1ll*prime[j]<MAXN&&j<=pcnt; ++j) {
vis[i*prime[j]]=1;
if(i%prime[j]==0)break;
else mu[i*prime[j]]=-mu[i];
}
}
for(int i=1; i<=pcnt; ++i) {
for(int j=1; j*1ll*prime[i]<MAXN; ++j) {
g[j*prime[i]]+=mu[j];
}
}
for(int i=1; i<MAXN; ++i) {//预处理后面部分的前缀和
sum[i]=sum[i-1]+g[i];
}
int t=read(),n,m,l,r;
ll ans;
while(t--) {
n=read(),m=read();
if(n>m)swap(n,m);
ans=0;
for(l=1; l<=n; l=r+1) {//整除分块
r=min(n/(n/l),m/(m/l));
ans+=(n/l)*(sum[r]-sum[l-1]+0ll)*(m/l);
}
printf("%lld\n",ans);
}
return 0;
}
P1829 [国家集训队]Crash的数字表格 / JZPTAB
这题是要求：\(\sum\limits_{i=1}^n\sum\limits_{j=1}^m\operatorname{lcm}(i,j)\)
根据\(\operatorname{lcm}\)的性质：
\(\sum\limits_{i=1}^n\sum\limits_{j=1}^m\dfrac{ij}{\gcd(i,j)}\)
加入一个\(d=\gcd(i,j)\)：
\(\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^m[\gcd(i,j)=d]\dfrac{ij}{d}\)
把\(i,j\)都除以\(d\)：
\(\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]ijd\)
提出\(d,i\)，把\(\epsilon=I*\mu\)代入：
\(\sum\limits_{d=1}^nd\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}i\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}j\sum\limits_{k|\gcd(i,j)}\mu(k)\)
\(k|\gcd(i,j)\)其实就是\(k|i\&k|j\)（手玩一下可以知道）
\(\sum\limits_{d=1}^nd\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}i\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}j\sum\limits_{k|i\&k|j}\mu(k)\)
直接把\(k\)提前：
\(\sum\limits_{d=1}^nd\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum\limits_{k|i\&i<\lfloor\frac{n}{d}\rfloor}i\sum\limits_{k|j\&j<\lfloor\frac{m}{d}\rfloor}j\)
把\(i,j\)都除以\(k\)
\(\sum\limits_{d=1}^nd\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum\limits_{i=1}^{\lfloor\frac{n}{dk}\rfloor}ik\sum\limits_{j=1}^{\lfloor\frac{m}{dk}\rfloor}jk\)
把两个\(k\)都提出来，把\(d\)乘进去，再设\(S(x)=1+2+\cdots+x=\dfrac{x(x+1)}{2}\)：
\(\sum\limits_{d=1}^n\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}dk^2\mu(k)S(\lfloor\frac{n}{dk}\rfloor)S(\lfloor\frac{m}{dk}\rfloor)\)
先枚举\(T=dk\)，再枚举\(k|T\)
\(\sum\limits_{T=1}^n\sum\limits_{k|T}Tk\mu(k)S(\lfloor\frac{n}{T}\rfloor)S(\lfloor\frac{m}{T}\rfloor)\)
再提出\(TS(\lfloor\frac{n}{T}\rfloor)S(\lfloor\frac{m}{T}\rfloor)\)
\(\sum\limits_{T=1}^nTS(\lfloor\frac{n}{T}\rfloor)S(\lfloor\frac{m}{T}\rfloor)\sum\limits_{k|T}k\mu(k)\)
我们设\(f(T)=\sum\limits_{k|T}k\mu(k)\)，看一下它的性质：
首先\(id,\mu\)都是积性函数，那么\(id\times\mu\)也是积性函数，而且\(I\)也是积性函数，那么\(f=(id\times\mu)*I\)也是积性函数。
那么欧拉筛这个\(f\)还需要看\(prime[j]|T\)的情况：
我们设\(T\)分解质因数后为\(p_1^{k_1}p_2^{k_2}\cdots p_c^{k_c}\)（不妨设其中\(p_1=prime[j]\)）则
\(f(Tp_1)-f(T)=\sum\limits_{k|Tp_1}k\mu(k)-\sum\limits_{k|T}k\mu(k)=\sum\limits_{k|Tp_1\&k\not|T}k\mu(k)\)
\(Tp_1=p_1^{k_1+1}p_2^{k_2}\cdots p_c^{k_c}\)，那么满足\(k|Tp_1\&k\not|T\)的数必须是\(p_1^{k_1+1}\)的倍数，而且\(k_1\ge1\)即\(k_1+1\ge2\)那么\(\mu(k)=0\)，那么后面那部分就是\(0\)，我们就得到：
\(f(Tp_1)=f(T)\)
最后总结一下：\(f(Tp_1)=\begin{cases}f(T)\qquad p_1|T\\f(T)f(p_1)\quad p_1\not|T\end{cases}\)
最后要求\(\sum\limits_{T=1}^nTS(\lfloor\frac{n}{T}\rfloor)S(\lfloor\frac{m}{T}\rfloor)f(T)\)，一边线性筛一边求就行了：
#include<bits/stdc++.h>
using namespace std;
const int mod=20101009;
const int MAXN=1e7+128;
typedef long long ll;
int n,m,ans;
bool isp[MAXN+1];
int prime[MAXN+1],f[MAXN+1],pcnt;
inline int S(int x) {
return (x*1ll*(x+1)/2)%mod;
}
int main() {
cin>>n>>m;
if(n>m)swap(n,m);
f[1]=1,isp[1]=true;
ans=(ans+S(n)*1ll*S(m))%mod;//先加上T=1时的值
for(int i=2; i<=n; ++i) {//线性筛f(T)
if(!isp[i]) {
prime[++pcnt]=i;
f[i]=mod-i+1;
}
for(int j=1,x; j<=pcnt&&(x=i*prime[j])<=MAXN; ++j) {
isp[x]=true;
if(i%prime[j]) {
f[x]=f[i]*1ll*f[prime[j]]%mod;
} else {
f[x]=f[i];
break;
}
}
ans=(ans+S(n/i)*1ll*S(m/i)%mod*i%mod*f[i])%mod;//同时计算答案
}
printf("%d\n",ans);
return 0;
}
