莫比乌斯反演

以下内容中，$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;
}

这个卷积满足：

1. 交换律$f*g=g*f$
2. 结合律$f*g*h=f*(g*h)$
3. 分配率$f*(g+h)=f*g+f*h$
4. 当$f(x),g(x)$都为积性函数时，$(f*g)(x)$为积性函数
5. 当$(f*g)(x),g(x)$都为积性函数时，$f(x)$为积性函数

常见的狄利克雷卷积有：

1. $\varphi*I=id$
2. $\mu*I=\epsilon$
3. $\epsilon*f=f$
4. $id*I=\sigma$
5. $I*I=d$
6. $(\lambda*I)(n)=[\sqrt{n}\in Z]$

简单证明：（3,4,5见相关函数的定义：积性函数相关）

1. $\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$）

2. $\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$

6. $(\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})$可以通过。

再来一题

P3455 [POI2007]ZAP-Queries

先把答案用式子表示一下：

$\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;
}

P2257 YY的GCD

这题是要求：$\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;
}