狄利克雷卷积&莫比乌斯反演学习笔记

数学 学习笔记

---

前置知识

莫比乌斯函数

莫比乌斯函数，即$\mu(x)=\cases{1,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x=1\\(-1)^s,\ \ \ \ \ \ \ x=p_1\cdot p_2\cdots p_s(\forall_{1\leqslant i\leqslant x,1\leqslant j\leqslant x,i\ne j}p_i\ne p_j)\\0,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ otherwise}$

狄利克雷卷积

狄利克雷卷积，可以表示为

$$设r=f\times g，则r(n)=\sum_{d\mid n}f(d)\times g(\frac{n}{d})$$

狄利克雷卷积满足交换律和结合律，且存在单位元$\epsilon(n)=[n==1]$

简证$\epsilon$为单位元：

$$有一数论函数f，则(f\times\epsilon)(n)=\sum_{d\mid n} f(d)\times\epsilon(\frac{n}{d})=f(n)\times\epsilon(1)+\sum_{d\mid n,d\ne 1}f(d)\times\epsilon(\frac{n}{d})\\=f(n)\times 1+\sum_{d\mid n,d\ne 1}f(d)\times0=f(n)+0=f(n)$$

例子：

$$\epsilon=\mu\times I，即\epsilon(n)=\sum_{d\mid n}\mu(d)\times I(\frac{n}{d})=\sum_{d\mid n}\mu(d)=\cases{1,\ &n=1\\0,\ &n>1}=[n==1]\\d=I\times I，即d(n)=\sum_{d\mid n}I(d)\times I(\frac{n}{d})=\sum_{d\mid n}1=n的因子个数$$

重要的数论函数中互相转化：

$$\mu\cdot I=\epsilon\\\varphi\cdot I=id\\id\cdot I=\sigma\\I\cdot I=d\\\mu\cdot id=\varphi$$
可以列一个表：

$\ $	$\epsilon$	$I$	$d$	$\mu$	$\varphi$	$\sigma$	$id$
$\epsilon$	$\epsilon$	$I$	$d$	$\mu$	$\varphi$	$\sigma$	$id$
$I$	$I$	$d$	$-$	$\epsilon$	$id$	$d\cdot id$	$\sigma$
$d$	$d$	$-$	$-$	$I$	$\sigma$	$-$	$I\cdot\sigma$
$\mu$	$\mu$	$\epsilon$	$I$	$-$	$-$	$id$	$\varphi$
$\varphi$	$\varphi$	$id$	$\sigma$	$-$	$-$	$id\cdot id$	$-$
$\sigma$	$\sigma$	$d\cdot id$	$-$	$id$	$id\cdot id$	$d\cdot id\cdot id$	$-$
$id$	$id$	$\sigma$	$I\cdot\sigma$	$\varphi$	$-$	$-$	$\varphi\cdot\sigma$

公式推导

1. $$对于n\in\mathbb{Z^+}，\sum_{d\mid n}\mu(d)=\cases{1,\ &n=1\\0,\ &n>1}\\证明如下：\\\sum_{d\mid n}\mu(d)=\sum_{d\mid n,p^2\mid d(p\in\mathbb{P})}\mu(d)+\sum_{d\mid n,p^2\not\mid d(p\in\mathbb{P})}\mu(d)=0+\sum_{d\mid n,p^2\not\mid d(p\in\mathbb{P})}\mu(d)\\设n=\prod_{i=1}^s p_i^{k_i}，m=\prod_{i=1}^s p_i，则\sum_{d\mid n}\mu(d)=\sum_{d\mid n,p^2\not\mid d(p\in\mathbb{P})}\mu(d)=\sum_{d\mid m}\mu(d)\\设w\mid m，则有：\\当w有0个质因子，即w=1时，有1=C^0_s种选择，贡献为\mu(w)=1=(-1)^0\\当w有1个质因子时，有C^1_s种选择，贡献为\mu(w)=(-1)^1=-1\\当w有2个质因子时，有C^2_s种选择，贡献为\mu(w)=(-1)^2=1\\\vdots\\当w有i个质因子时，有C^i_s种选择，贡献为\mu(w)=(-1)^i\\\vdots\\当w有s个质因子时，有C^s_s种选择，贡献为\mu(w)=(-1)^n\\故\sum_{d\mid n}\mu(d)=\sum_{d\mid m}\mu(d)=\sum_{i=0}^sC^i_s\cdot (-1)^i=\sum_{i=0}^s C^i_s\cdot 1^{s-i}\cdot (-1)^i\\由二项式定理得\sum_{d\mid n}\mu(d)=\sum_{i=0}^s C^i_s\cdot 1^{s-i}\cdot (-1)^i=(1+(-1))^s=0^s=0$$
2. $$设I(x)=1,\epsilon(x)=[n==1]（即狄利克雷卷积的单位元），则有\mu为I的逆元，即\mu\times I=\epsilon，证明如下：\\由狄利克雷卷积得\mu\times I=\sum_{d\mid n}\mu(d)\times I(\frac{n}{d})=\sum_{d\mid n}\mu(d)\times 1=\sum_{d\mid n}\mu(d)\\由性质1得\mu\times I=\sum_{d\mid n}\mu(d)=\cases{1,\ &n=1\\0,\ &n>1}=[n==1]=\epsilon\\即\mu\times I=\epsilon，\mu为I的逆元$$
3. 莫比乌斯反演公式 $$设I(x)=1,\epsilon(x)=[n==1]，则对于数论函数f和F=f\times I，f=F\times\mu，证明如下：\\\because F=f\times I\\\therefore F\times\epsilon=f\times I\\由莫比乌斯函数性质2得\mu\times I=\epsilon，即F\times\mu\times I=f\times I\\\Rightarrow f=F\times\mu\\可以写成：当F(x)=\sum_{d\mid n}f(d)，则f(n)=\sum_{d\mid n}\mu(d)\times F(\frac{n}{d})$$
4. 莫比乌斯函数的积性证明 $$当\gcd(n,m)=1时，\mu(n)\times\mu(m)=\mu(n\times m)，证明如下：\\当n=1时，\mu(n)=1，即n\times m=m，\mu(n)\times\mu(m)=\mu(m)\\即\mu(n)\times\mu(m)=\mu(m)=\mu(1\times m)=\mu(n\times m)（m同理）\\当存在质数p使得p^2\mid n时，有\mu(n)=0，即\mu(n)\times\mu(m)=0\\此时肯定有p^2\mid n\times m，即\mu(n\times m)=0\\整理得\mu(n)\times\mu(m)=0=\mu(n\times m)（m同理）\\否则肯定有n=\prod_{i=1}^s p_i，且\forall_{i\ne j}p_i\ne p_j，m=\prod_{i=1}^t q_i，且\forall_{i\ne j}q_i\ne q_j\\\because\gcd(n,m)=1\\\therefore \forall_{i,j} p_i\ne q_j\\\Rightarrow\mu(n\times m)=(-1)^{s+t}\\\because\mu(n)=(-1)^s，\mu(m)=(-1)^t\\\therefore\mu(n)\times \mu(m)=(-1)^{s+t}=\mu(n\times m)$$
5. $$\sum_{d\mid\gcd(i,j)}\mu(d)=[\gcd(i,j)==1]\\由莫比乌斯函数性质1得\sum_{d\mid n}\mu(d)=[n==1]，将\gcd(i,j)带入就可以得到上式$$

整除分块

在莫比乌斯反演的应用题中，经常会有$\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor$的形式，但是$O(n)$的复杂度有时候仍然满足不了需求

为了解决这一问题，我们就要用到整除分块

很容易发现很多$\lfloor\frac{n}{i}\rfloor$的值是一样的，且所有的$\lfloor\frac{n}{i}\rfloor$为一个不下降子序列，呈块状分布

通过简单的计算可以得到，对于每个起点为$i$的块，它的值为$\lfloor\frac{n}{i}\rfloor$，终点为$\lfloor\frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor$，然后我们就可以用$O(\sqrt{n})$的算法计算$\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor$的值了

代码：

int l=1,r,ans=0;
if(n>m)
    swp(n,m);
while(l<=n){
    r=min(n/(n/l),mm/(m/l));
    ans+=(r-l+1)*(n/l);
    l=r+1;
}

莫比乌斯反演应用

1.求值：

$$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==1]$$ 证明如下：
由性质5得 $$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==1]=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid \gcd(i,j)}\mu(d)=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid i}\sum_{d\mid j}\mu(d)$$
不妨设$n\leqslant m$，则由$d\mid\gcd(i,j),i\leqslant n,j\leqslant m$显然可得$d\in(1,n)$
则有 $$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==1]=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid i}\sum_{d\mid j}\mu(d)=\sum_{i=1}^n\sum_{d\mid i}\sum_{j=1}^m\sum_{d\mid j}\mu(d)\\=\sum_{d=1}^n(\mu(d)\times\sum_{i=1}^n\sum_{j=1}^m[d\mid i,d\mid j])=\sum_{d=1}^n(\mu(d)\times\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}1)\\=\sum_{d=1}^n(\mu(d)\times\lfloor\frac{n}{d}\rfloor\times\lfloor\frac{m}{d}\rfloor)$$
然后后面的式子用整除分块可以优化到$O(\sqrt{n})$的时间复杂度
2.例题1：[POI2007]ZAP-Queries
题意：共$T$组数据，求值$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==k]$
数据范围：$T\leqslant 5\cdot 10^4，n,m\leqslant 5\cdot 10^4$
简证：原式很显然等于 $$\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i,j)==1]$$
然后由应用1得 $$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==k]=\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}(\mu(d)\cdot\lfloor\frac{\lfloor\frac{n}{k}\rfloor}{d}\rfloor\cdot\lfloor\frac{\lfloor\frac{m}{k}\rfloor}{d}\rfloor)\\=\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}(\mu(d)\cdot\lfloor\frac{n}{k\cdot d}\rfloor\cdot \lfloor\frac{m}{k\cdot d}\rfloor)$$
然后直接整除分块，时间复杂度：$O(n+T\sqrt{n})$

#include<stdio.h>
const int maxn=100005;
int i,j,k,m,n,T,cnt,l,r,ans;
int a[maxn],p[maxn],mu[maxn],sum[maxn];
inline void swp(int &a,int &b){
    a+=b,b=a-b,a-=b;
}
inline int min(int a,int b){
    return a<b? a:b;
}
int main(){
    p[1]=mu[1]=1;
    for(i=2;i<maxn;i++){
        if(p[i]==0)
            a[++cnt]=i,mu[i]=-1;
        for(j=1;j<=cnt;j++){
            if(i*a[j]>=maxn)
                break;
            p[i*a[j]]=1;
            if(i%a[j]==0)
                mu[i*a[j]]=0;
            else mu[i*a[j]]=-mu[i];
        }
    }
    for(i=1;i<maxn;i++)
        sum[i]=sum[i-1]+mu[i];
    scanf("%d",&T);
    while(T--){
        scanf("%d%d%d",&n,&m,&k);
        n/=k,m/=k;
        if(n>m)
            swp(n,m);
        l=1,ans=0;
        while(l<=n){
            r=min(n/(n/l),m/(m/l));
            ans+=(sum[r]-sum[l-1])*(n/l)*(m/l);
            l=r+1;
        }
        printf("%d\n",ans);
    }
    return 0;
}

3.例题2：[HAOI2011]Problem b
题意：共$n$组数据，求值$\sum_{i=a}^b\sum_{j=c}^d[\gcd(i,j)==k]$
数据范围：$n\leqslant 5\cdot 10^4,a\leqslant b,c\leqslant d,a,b,c,d\leqslant 5\cdot 10^4$
这道题与例题1很像，但是$i$和$j$不是从$1$开始的，因此这有点像在一个二维平面内取一个矩形，你想到了什么？容斥定理！
我们只需要写一个跟例题1一样的函数$calc(n,m,k)$求值
二维容斥：

$$ans=calc(b,d,k)-calc(b,c-1,k)-calc(a-1,d,k)+calc(a-1,c-1,k)$$

#include<stdio.h>
const int maxn=100005;
int i,j,k,m,n,a,b,c,d,T,cnt;
int z[maxn],p[maxn],mu[maxn],sum[maxn];
inline void swp(int &a,int &b){
    a+=b,b=a-b,a-=b;
}
inline int min(int a,int b){
    return a<b? a:b;
}
int calc(int n,int m,int k){
    int l=1,r,res=0;
    n/=k,m/=k;
    if(n>m)
        swp(n,m);
    while(l<=n){
        r=min(n/(n/l),m/(m/l));
        res+=(sum[r]-sum[l-1])*(n/l)*(m/l);
        l=r+1;
    }
    return res;
}
int main(){
    p[1]=mu[1]=1;
    for(i=2;i<maxn;i++){
        if(p[i]==0)
            z[++cnt]=i,mu[i]=-1;
        for(j=1;j<=cnt;j++){
            if(i*z[j]>=maxn)
                break;
            p[i*z[j]]=1;
            if(i%z[j]==0)
                mu[i*z[j]]=0;
            else mu[i*z[j]]=-mu[i];
        }
    }
    for(i=1;i<maxn;i++)
        sum[i]=sum[i-1]+mu[i];
    scanf("%d",&n);
    for(i=1;i<=n;i++){
        scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
        printf("%d\n",calc(b,d,k)-calc(b,c-1,k)-calc(a-1,d,k)+calc(a-1,c-1,k));
    }
    return 0;
}

4.例题3：YY的GCD
题意：共$T$组数据，求$\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)\in\mathbb{P}]$
数据范围：$T\leqslant 10^4，n,m\leqslant 10^7$
不妨设$n\leqslant m$，很自然可以想到枚举质数的套路：

$$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)\in\mathbb{P}]=\sum_{p=1,p\in\mathbb{P}}^n\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==p]$$
由应用2得原式等于 $$\sum_{p=1,p\in\mathbb{P}}^n\sum_{d=1}^{\lfloor\frac{n}{p}\rfloor}(\mu(d)\times\lfloor\frac{n}{p\times d}\rfloor\times\lfloor\frac{m}{p\times d}\rfloor)$$
但是复杂度还是太高了，我们都知道 $$\pi(x)=\sum_{i=1}^{x-1}[i\in\mathbb{P}]\approx\frac{x}{In(x)}$$
也就是说我们每次求解的复杂度都为$O(\frac{n}{In(n)}\sqrt{n})$，当$n=10^7$时，$\frac{n}{In(n)}\sqrt{n}=\frac{10^7}{In(10^7)}\sqrt{10^7}\approx\frac{10000000}{16}\cdot 3162\approx 3\cdot 10^9$，再加上多组数据，是肯定无法通过的
于是我们要考虑优化，令$k=p\times d$，则有： $$\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)\in\mathbb{P}]=\sum_{p=1,p\in\mathbb{P}}^n\sum_{d=1}^{\lfloor\frac{n}{p}\rfloor}(\mu(\frac{k}{p})\times\lfloor\frac{n}{k}\rfloor\times\lfloor\frac{m}{k}\rfloor)$$
然后可以枚举$k$： $$\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)\in\mathbb{P}]=\sum_{k=1}^n(\sum_{p\mid k,p\in\mathbb{P}}\mu(\frac{k}{p})\times\lfloor\frac{n}{k}\rfloor\times\lfloor\frac{m}{k}\rfloor)$$
设$f(x)=\sum_{p\mid x,p\in\mathbb{P}}\mu(\frac{x}{p})$，则原式等于 $$\sum_{k=1}^n(f(k)\times\lfloor\frac{n}{k}\rfloor\times\lfloor\frac{m}{k}\rfloor)$$
此时$\lfloor\frac{n}{k}\rfloor$和$\lfloor\frac{m}{k}\rfloor$可以用整除分块做到$O(\sqrt{n})$的复杂度
那么我们的目标就是将$f(x)=\sum_{p\mid x,p\in\mathbb{P}}\mu(\frac{x}{p})$用线性筛处理，并求出前缀和，我们可以分类讨论一下：
当$x\in\mathbb{P}$，则$f(x)=\mu(1)=1$
否则设$x=y\times z$，且$y$为$x$的最小质因子，此时有：
当$y\mid z$，即$x$的质因子$y$出现了多次时
- 当不存在质数$w$使$w^2\mid z$时，只有$\mu(\frac{x}{y})\ne 0$，其余均存在$y^2\mid \frac{x}{p}$，即$\mu(\frac{x}{p})=0$$f(x)=\sum_{p\mid x,p\ne y,p\in\mathbb{P}}\mu(\frac{x}{p})+\mu(\frac{x}{y})=\sum_{p\mid x,p\ne y,p\in\mathbb{P}}0+\mu(z)=\mu(z)$
- 当存在质数$w$使$w^2\mid z$时，$w^2\mid\frac{x}{y}=z$且其余的$p$满足$y^2\mid\frac{x}{p}$，故所有的$\mu(\frac{x}{p})=0,$，$f(x)=\sum_{p\mid x,p\in\mathbb{P}}\mu(\frac{x}{p})=\sum_{p\mid x,p\in\mathbb{P}}0=0=\mu(z)$
当$y\not\mid z$，即$\gcd(y,z)=1$，$x$的最小质因子只出现一次时
$$f(z)=\sum_{p\mid z,p\in\mathbb{P}}\mu(\frac{z}{p}),f(x)=\sum_{p\mid x,p\in\mathbb{P}}\mu(\frac{x}{p})=\sum_{p\mid x,p\in\mathbb{P}}\mu(\frac{y\times z}{p})$$
由$\gcd(y,z)=1$得对于所有的质数$p$，都有$\mu(\frac{x}{p})=\mu(\frac{y\times z}{p})=-\mu(\frac{z}{p})$
即对于满足$p\ne y$的质数$p$，$\mu(\frac{x}{p})=-\mu(\frac{z}{p})$
$$\Rightarrow f(x)=\sum_{p\mid x,p\in\mathbb{P}}\mu(\frac{x}{p})=\sum_{p\mid x,p\ne y,p\in\mathbb{P}}\mu(\frac{x}{p})+\mu(\frac{x}{y})\\=\sum_{p\mid x,p\ne y,p\in\mathbb{P}}(-\mu(\frac{z}{p}))+\mu(z)=-\sum_{p\mid x,p\ne y,p\in\mathbb{P}}\mu(\frac{z}{p})+\mu(z)=\mu(z)-f(z)$$
综上所述得 $$f(x)=\cases{1,\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x\in\mathbb{P}\\\mu(z)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x=y\times z,y\mid z\\\mu(z)-f(z)\ \ \ \ \ \ \ \ \ x=y\times z,y\not\mid z}$$
可以用线性筛在用$O(n)$的复杂度求出$\mu(x)$和$f(x)$
求出$f(x)$之后别忘了前缀和
最终复杂度：$O(n+T\sqrt{n})$
代码：

#include<stdio.h>
const int maxn=10000005;
int i,j,k,m,n,cnt,T,l,r;
long long ans;
int a[maxn],p[maxn],mu[maxn],f[maxn],sum[maxn];
inline void swp(int &a,int &b){
    a+=b,b=a-b,a-=b;
}
inline int min(int a,int b){
    return a<b? a:b; 
}
int main(){
    p[1]=mu[1]=1;
    for(i=2;i<maxn;i++){
        if(p[i]==0)
            a[++cnt]=i,mu[i]=-1,f[i]=1;
        for(j=1;j<=cnt;j++){
            if(i*a[j]>maxn)
                break;
            p[i*a[j]]=1;
            if(i%a[j]==0){
                mu[i*a[j]]=0;
                f[i*a[j]]=mu[i];
                break;
            }
            mu[i*a[j]]=-mu[i];
            f[i*a[j]]=mu[i]-f[i];
        }
        sum[i]=sum[i-1]+f[i];
    }
    scanf("%d",&T); 
    while(T--){
        scanf("%d%d",&n,&m);
        if(n>m)
            swp(n,m);
        l=1,ans=0;
        while(l<=n){
            r=min(n/(n/l),m/(m/l));
            ans+=(long long)(n/l)*(long long)(m/l)*(long long)(sum[r]-sum[l-1]);
            l=r+1;
        }
        printf("%lld\n",ans);
    }
    return 0;
}

5.求值

$$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==k]\cdot i\cdot j$$
不妨设$n\leqslant m$，然后套路地除一个$k$，得： $$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==k]\cdot i\cdot j=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i,j)==1]\cdot (i\cdot k)\cdot(j\cdot k)$$
这里$i$和$j$乘$k$的原因是$i$和$j$的范围都缩小了$k$，故$i$和$j$都要乘$k$
然后将$[\gcd(i,j)==1]$用$\sum_{d\mid\gcd(i,j)}\mu(d)$替换，得：
$$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==k]\cdot i\cdot j=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d\mid\gcd(i,j)}\mu(d)\cdot i\cdot j\cdot k^2\\=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}(i\cdot j\cdot k^2\cdot \sum_{d\mid\gcd(i,j)}\mu(d))$$
然后仍然是套路地将$d$提前，并提出$d$，得：
$$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==k]\cdot i\cdot j=\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}(\mu(d)\cdot d^2\cdot\sum_{i=1}^{\lfloor\frac{n}{k\cdot d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k\cdot d}\rfloor}i\cdot j\cdot k^2)\\=k^2\cdot\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}(\mu(d)\cdot d^2\cdot\sum_{i=1}^{\lfloor\frac{n}{k\cdot d}\rfloor}i\sum_{j=1}^{\lfloor\frac{m}{k\cdot d}\rfloor}j)$$
然后就发现这是等差数列，因此有： $$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==k]\cdot i\cdot j=k^2\cdot\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}(\mu(d)\cdot d^2\cdot(\sum_{i=1}^{\lfloor\frac{n}{k\cdot d}\rfloor}i)(\sum_{j=1}^{\lfloor\frac{m}{k\cdot d}\rfloor}j))\\=k^2\cdot\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}(\mu(d)\cdot d^2\cdot(\frac{\lfloor\frac{n}{k\cdot d}\rfloor\cdot (\lfloor\frac{n}{k\cdot d}\rfloor+1)}{2})\cdot (\frac{\lfloor\frac{m}{k\cdot d}\rfloor\cdot (\lfloor\frac{m}{k\cdot d}\rfloor+1)}{2}))\\=k^2\cdot\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}(\mu(d)\cdot d^2\cdot\frac{\lfloor\frac{n}{k\cdot d}\rfloor\cdot (\lfloor\frac{n}{k\cdot d}\rfloor+1)\cdot \lfloor\frac{m}{k\cdot d}\rfloor\cdot (\lfloor\frac{m}{k\cdot d}\rfloor+1)}{4})$$
最后整除分块，时间复杂度$O(\sqrt{n})$
6.例题4：[SDOI2015]约数个数和
题意：共$T$组数据，求$\sum_{i=1}^n\sum_{j=1}^m d(i\cdot j)$（$d(x)=\sum_{d\mid x} 1$）
数据范围：$T\leqslant 5\cdot 10^4，n,m\leqslant 5\cdot 10^4$
先证明一个毒瘤式子 $$d(x\cdot y)=\sum_{i\mid x}\sum_{j\mid y}[\gcd(i,j)==1]$$
当$x=\prod_{i=1}^s p_i^{k_i}$因为$d(x)=\prod_{i=1}^s(k_i+1)$为$x$的因数个数，因此考虑将$x\cdot y$的每个因数对应一组互质的因数
对于每个质数$p\mid x,y$，且$x$中$p$的幂次为$a$，$y$中$p$的幂次为$b$
则$x\cdot y$中$p$的幂次为$a+b$，根据乘法原理，得$d(x\cdot y)=d(\frac{x\cdot y}{p^{a+b}})\cdot (a+b+1)$，即贡献为$a+b+1$
因为$i\mid x,j\mid y,\gcd(i,j)=1$，则$i$中$p$的幂次$u$不超过$a$，$j$中$p$的幂次$v$不超过$b$，且$u$和$v$中必有一个为$0$
- 当$u=0,v\ne 0$时，$v=1,2,\cdots b$，即有$b$种取值
- 当$v=0,u\ne 0$时，$u=1,2,\cdots a$，即有$a$种取值
- 当$u=0,v=0$时，有$1$种取值
根据加法原理得每个质数$p$会对右边的式子贡献$a+b+1$的取值，也可写作$s_i+1$
故每个质数作为因子对左边的式子贡献$a+b+1$，作为互质的数对右边的式子同样贡献$a+b+1$，即$x\cdot y$的每个因数对应一组互质的因数
然后原式可以化为$\sum_{i=1}^n\sum_{j=1}^m\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)==1]$
然后用莫比乌斯函数替换$\gcd$得： $$\sum_{i=1}^n\sum_{j=1}^md(i\cdot j)=\sum_{i=1}^n\sum_{j=1}^m\sum_{x\mid i}\sum_{y\mid j}\sum_{d\mid\gcd(x,y)}\mu(d)$$
将$x$和$y$提前，得： $$\sum_{i=1}^n\sum_{j=1}^md(i\cdot j)=\sum_{x=1}^n\sum_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum_{y=1}^m\sum_{j=1}^{\lfloor\frac{m}{y}\rfloor}\sum_{d\mid\gcd(x,y)}\mu(d)$$
然后枚举$d$得： $$\sum_{i=1}^n\sum_{j=1}^md(i\cdot j)=\sum_{d=1}^n\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{i=1}^{\lfloor\frac{n}{d\cdot x}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}}\sum_{j=1}^{\lfloor\frac{m}{d\cdot y}\rfloor}\mu(d)=\sum_{d=1}^n\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{d\cdot x}\rfloor\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{m}{d\cdot y}\rfloor\mu(d)\\=\sum_{d=1}^n\mu(d)\cdot (\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{d\cdot x}\rfloor)\cdot (\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{m}{d\cdot y}\rfloor)=\sum_{d=1}^n\mu(d)\cdot (\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{x}\rfloor)\cdot (\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{\lfloor\frac{m}{d}\rfloor}{y}\rfloor)$$
设$p=\lfloor\frac{n}{d}\rfloor$，$q=\lfloor\frac{m}{d}\rfloor$，则有： $$\sum_{i=1}^n\sum_{j=1}^md(i\cdot j)=\sum_{d=1}^n\mu(d)\cdot (\sum_{x=1}^p\lfloor\frac{p}{x}\rfloor)\cdot (\sum_{y=1}^q\lfloor\frac{q}{y}\rfloor)$$
于是可以设$f(x)=\sum_{i=1}^x\lfloor\frac{x}{i}\rfloor$，则原式可化为$\sum_{d=1}^n\mu(d)\cdot f(p)\cdot f(q)$
我们可以用$O(n)$处理$\mu$，然后再看$f$，发现$f$也是一个整除分块，需要$O(n\cdot \sqrt{n})$来预处理，然后$(\sum_{x=1}^p\lfloor\frac{p}{x}\rfloor)\cdot (\sum_{y=1}^q\lfloor\frac{q}{y}\rfloor)$的部分可以整除分块来处理
时间复杂度：$O(n\sqrt{n}+T\sqrt{n})$
代码：

#include<stdio.h>
const int maxn=50005; 
int i,j,k,m,n,T,cnt,l,r;
long long ans;
int a[maxn],p[maxn],mu[maxn],sum[maxn];
long long f[maxn];
inline void swp(int &a,int &b){
    a+=b,b=a-b,a-=b;
}
inline int min(int a,int b){
    return a<b? a:b;
}
int main(){
    p[1]=mu[1]=1;
    for(i=2;i<maxn;i++){
        if(p[i]==0)
            a[++cnt]=i,mu[i]=-1;
        for(j=1;j<=cnt;j++){
            if(i*a[j]>=maxn)
                break;
            p[i*a[j]]=1;
            if(i%a[j]==0){
                mu[i*a[j]]=0;
                break;
            }
            mu[i*a[j]]=-mu[i];
        }
    }
    for(i=1;i<maxn;i++)
        sum[i]=sum[i-1]+mu[i];
    for(i=1;i<maxn;i++){
        l=1,ans=0;
        while(l<=i){
            r=i/(i/l);
            ans+=(r-l+1)*(i/l);
            l=r+1;
        }
        f[i]=ans;
    }
    scanf("%d",&T);
    while(T--){
        scanf("%d%d",&n,&m);
        if(n>m)
            swp(n,m);
        l=1,ans=0;
        while(l<=n){
            r=min(n/(n/l),m/(m/l));
            ans+=(long long)(sum[r]-sum[l-1])*(long long)f[n/l]*(long long)f[m/l];
            l=r+1;
        } 
        printf("%lld\n",ans);
    }
    return 0;
}

7.例题5：[国家集训队]Crash的数字表格 / JZPTAB
题意：求$\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)$
数据范围：$n,m\leqslant 10^7$
由

$$i\cdot j=\gcd(i,j)\cdot lcm(i,j)$$ 得 $$\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)=\sum_{i=1}^n\sum_{j=1}^m\frac{i\cdot j}{\gcd(i,j)}$$
设$\gcd(i,j)=d$，则等价于$d\mid i,d\mid j,\gcd(\frac{i}{d},\frac{j}{d})=1$，可得： $$\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)=\sum_{i=1}^n\sum_{j=1}^m\sum_{d\mid i,d\mid j,\gcd(\frac{d}{i},\frac{d}{j})=1}\frac{i\cdot j}{d}$$
枚举$d$可得： $$\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)==1]\cdot \frac{i\cdot j\cdot d^2}{d}\\=\sum_{d=1}^n(d\cdot \sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)==1]\cdot i\cdot j)$$
这里$i$和$j$乘$d^2$的原因是$i$和$j$的范围都缩小了$d$
把$d$替换为$k$得 $$\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)=\sum_{k=1}^n(k\cdot \sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i,j)==1]\cdot i\cdot j)$$
设$f(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)==1]\cdot i\cdot j)$，则原式为$\sum_{k=1}^n(k\cdot f(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor))$
用$(n=\lfloor\frac{n}{k}\rfloor,m=\lfloor\frac{m}{k}\rfloor,k=1)$带入应用5则有： $$f(n,m)=\sum_{d=1}^{\lfloor\frac{n}{1}\rfloor}(\mu(d)\cdot d^2\cdot\frac{\lfloor\frac{n}{1\cdot d}\rfloor\cdot (\lfloor\frac{n}{1\cdot d}\rfloor+1)\cdot \lfloor\frac{m}{1\cdot d}\rfloor\cdot (\lfloor\frac{m}{1\cdot d}\rfloor+1)}{4})\\=\sum_{d=1}^n(\mu(d)\cdot d^2\cdot\frac{\lfloor\frac{n}{d}\rfloor\cdot (\lfloor\frac{n}{d}\rfloor+1)\cdot \lfloor\frac{m}{d}\rfloor\cdot (\lfloor\frac{m}{d}\rfloor+1)}{4})$$
我们求$f$就只需要处理$\mu(i)\cdot i^2$的前缀和然后数论分块就可以了！
原式$\sum_{k=1}^n(k\cdot f(\lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor))$也是一次数论分块
代码：

#include<stdio.h>
const int maxn=10000005;
const long long mod=20101009;
long long i,j,k,m,n,cnt,l,r,ans;
long long a[maxn],p[maxn],mu[maxn],sum[maxn];
inline void swp(long long &a,long long &b){
    a+=b,b=a-b,a-=b;
}
inline long long min(long long a,long long b){
    return a<b? a:b;
}
inline long long calc(long long x,long long y){
    return (x*(x+1)/2%mod)*(y*(y+1)/2%mod)%mod;
}
long long f(long long n,long long m){
    if(n>m)
        swp(n,m);
    long long l=1,r,res=0;
    while(l<=n){
        r=min(n/(n/l),m/(m/l));
        res=(res+(sum[r]-sum[l-1]+mod)%mod*calc(n/l,m/l))%mod;
        l=r+1;
    }
    return res;
}
int main(){
    p[1]=mu[1]=1;
    for(i=2;i<maxn;i++){
        if(p[i]==0)
            a[++cnt]=i,mu[i]=-1;
        for(j=1;j<=cnt;j++){
            if(i*a[j]>=maxn)
                break;
            p[i*a[j]]=1;
            if(i%a[j]==0){ 
                mu[i*a[j]]=0;
                break; 
            }
            mu[i*a[j]]=-mu[i];
        }
    }
    for(i=1;i<maxn;i++)
        sum[i]=(sum[i-1]+mu[i]*i*i%mod)%mod;
    scanf("%lld%lld",&n,&m);
    if(n>m)
        swp(n,m);
    l=1;
    while(l<=n){
        r=min(n/(n/l),m/(m/l));
        ans=(ans+((l+r)*(r-l+1)/2)%mod*f(n/l,m/l)%mod)%mod;
        l=r+1;
    }
    printf("%lld\n",ans); 
    return 0;
}

8.例题6：[SDOI2017]数字表格
题意：求$\prod_{i=1}^n\prod_{j=1}^m f_{\gcd(i,j)},f_x=\cases{0,\ &n=0\\1\ &n=1\\f_{n-2}+f_{n-1}\ &n>1}$
不妨设$n\leqslant m$则：
首先套路地枚举$\gcd$：

$$\prod_{i=1}^n\prod_{j=1}^m f_{\gcd(i,j)}=\prod_{d=1}^n\prod_{i=1}^n\prod_{j=1}^m[\gcd(i,j)==d]\cdot f_d=\prod_{d=1}^n f_d^{\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)==d]}$$
然后仍然套路地用莫比乌斯函数替换$\gcd$： $$\prod_{i=1}^n\prod_{j=1}^m f_{\gcd(i,j)}=\prod_{d=1}^n f_d^{\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\\gcd(i,j)==1]}\\=\prod_{d=1}^n f_d^{\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{k\mid \gcd(i,j)}\mu(k)}\\=\prod_{d=1}^n f_d^{\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{i=1}^{\lfloor\frac{n}{k\cdot d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k\cdot d}\rfloor}\mu(k)}\\=\prod_{d=1}^n f_d^{\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\cdot\lfloor\frac{n}{k\cdot d}\rfloor\cdot\lfloor\frac{m}{k\cdot d}\rfloor}$$
然后套路地令$t=k\cdot d$，则原式为：$\prod_{t=1}^n\prod_{d\mid t} f_d^{\mu(\frac{t}{d})\cdot \lfloor\frac{n}{t}\rfloor\cdot\lfloor\frac{m}{t}\rfloor}=\prod_{t=1}^n(\prod_{d\mid t} f_d^{\mu(\frac{t}{d})})^{\lfloor\frac{n}{t}\rfloor\cdot\lfloor\frac{m}{t}\rfloor}$
设$g(x)=\prod_{d\mid x} f_d^{\mu(\frac{x}{d})}$，就可以在外面直接预处理掉$f(x)$，因为复杂度太高所以考虑每个点更新它的倍数
最后原式为$\prod_{t=1}^n g(t)^{\lfloor\frac{n}{t}\rfloor\cdot \lfloor\frac{m}{t}\rfloor}$，直接数论分块搞掉
代码：

#include<stdio.h>
const int maxn=1000005;
const long long mod=1000000007;
int i,j,k,m,n,cnt,T,l,r;
long long ans;
int a[maxn],p[maxn],mu[maxn];
long long f[maxn],g[maxn],mul[maxn];
inline void swp(int &a,int &b){
    a+=b,b=a-b,a-=b;
}
inline int min(int a,int b){
    return a<b? a:b;
} 
long long ksm(long long a,int b){
    long long res=1;
    while(b>0){
        if(b&1)
            res=(res*a)%mod;
        a=(a*a)%mod,b>>=1;
    }
    return res;
}
long long inv(long long a,long long b){
    return ksm(a,b-2);
}
int main(){
    p[1]=mu[1]=1;
    f[0]=0,f[1]=g[1]=1;
    for(i=2;i<maxn;i++){
        f[i]=(f[i-1]+f[i-2])%mod;
        g[i]=inv(f[i],mod);
        if(p[i]==0)
            a[++cnt]=i,mu[i]=-1;
        for(j=1;j<=cnt;j++){
            if(i*a[j]>=maxn)
                break;
            p[i*a[j]]=1;
            if(i%a[j]==0){
                mu[i*a[j]]=0;
                break;
            }
            mu[i*a[j]]=-mu[i];
        }
    }
    for(i=0;i<maxn;i++)
        mul[i]=1;
    for(i=1;i<maxn;i++){
        if(mu[i]==0)
            continue;
        for(j=i;j<maxn;j+=i){
            if(mu[i]==1)
                mul[j]=mul[j]*f[j/i]%mod;
            if(mu[i]==-1)
                mul[j]=mul[j]*g[j/i]%mod;
        }
    }
    for(i=2;i<maxn;i++)
        mul[i]=mul[i]*mul[i-1]%mod;
    scanf("%d",&T);
    while(T--){
        scanf("%d%d",&n,&m);
        if(n>m)
            swp(n,m);
        l=1,ans=1;
        while(l<=n){
            r=min(n/(n/l),m/(m/l));
            ans=(ans*ksm(mul[r]*inv(mul[l-1],mod)%mod,(long long)(n/l)*(long long)(m/l)%(mod-1)))%mod;
            l=r+1;
        }
        printf("%lld\n",(ans+mod)%mod);
    }
    return 0;
}

9.例题7：[MtOI2019]幽灵乐团
我有一个绝妙的想法，但是这里空白太小，写不下
题解：[MtOI2019]幽灵乐团解题报告

#include<stdio.h>
const int maxn=100005;
long long i,j,k,m,n,T,p,A,B,C,l,r,cnt;
long long a[maxn],c[maxn],mu[maxn],varphi[maxn],fac1[maxn],fac2[maxn],w1[maxn],f1[maxn],g1[maxn],w2[maxn],f2[maxn],g2[maxn],w3[maxn],f3[maxn],g3[maxn],sum[maxn];
//a[i]={Prime},c[i]=[i\in{Prime}?]
//mu[i]=(i==1?)1:((i=\prod_{j=1}^s p_i,(p_x!=p_y(x!=y)))? 0:(-1)^s)
//fac1[i]=i!,fac2[i]=\prod_{j=1}^i i^i
//w1[i]=\prod_{d|i}i^{\mu(d/i)},f1[i]=\prod_{j=1}^i w1[i],g1[i]=f1[i]^{-1}
//w2[i]=w1[i]*(i^2),f2[i]=\prod_{j=1}^i w2[i],g2[i]=f2[i]^{-1}
//w3[i]=i^{\varphi(i)},f3[i]=\prod_{j=1]^i w3[i],g3[i]=f3[i]^{-1}
//sum[i]=\sum_{j=1}^i varphi[i]
inline long long min(long long a,long long b){
    return a<b? a:b;
}
long long ksm(long long a,long long b,long long p){
    long long res=1;
    while(b>0){
        if(b&1)
            res=(res*a)%p;
        a=(a*a)%p,b>>=1;
    }
    return res;
}
long long inv(long long a,long long b){
    return ksm(a,b-2,b);
}
void init(){
    c[1]=mu[1]=varphi[1]=1;
    for(long long i=2;i<maxn;i++){
        if(c[i]==0)
            a[++cnt]=i,mu[i]=-1,varphi[i]=i-1;
        for(long long j=1;j<=cnt;j++){
            if(i*a[j]>=maxn)
                break;
            c[i*a[j]]=1;
            if(i%a[j]==0){
                mu[i*a[j]]=0;
                varphi[i*a[j]]=varphi[i]*a[j];
                break;
            }
            mu[i*a[j]]=-mu[i];
            varphi[i*a[j]]=varphi[i]*(a[j]-1);
        }
    }
    fac1[0]=1;
    for(long long i=1;i<maxn;i++)
        fac1[i]=fac1[i-1]*i%p; 
    fac2[0]=1;
    for(long long i=1;i<maxn;i++)
        fac2[i]=fac2[i-1]*ksm(i,i%(p-1),p)%p; 
    for(long long i=0;i<maxn;i++)
        w1[i]=w2[i]=1;
    for(long long i=1;i<maxn;i++){
        if(mu[i]==0)
            continue;
        for(long long j=i;j<maxn;j+=i){
            if(mu[i]==1)
                w1[j]=(w1[j]*(j/i))%p;
            if(mu[i]==-1)
                w1[j]=(w1[j]*inv(j/i,p))%p;
        }
    }
    f1[0]=g1[0]=1;
    for(long long i=1;i<maxn;i++)
        f1[i]=f1[i-1]*w1[i]%p,g1[i]=inv(f1[i],p);
    for(long long i=1;i<maxn;i++)
        w2[i]=ksm(w1[i],i*i%(p-1),p);
    f2[0]=g2[0]=1;
    for(long long i=1;i<maxn;i++)
        f2[i]=f2[i-1]*w2[i]%p,g2[i]=inv(f2[i],p);
    for(long long i=1;i<maxn;i++)
        w3[i]=ksm(i,varphi[i],p);
    f3[0]=g3[0]=1;
    for(long long i=1;i<maxn;i++)
        f3[i]=f3[i-1]*w3[i]%p,g3[i]=inv(f3[i],p);
    for(long long i=1;i<maxn;i++)
        sum[i]=sum[i-1]+varphi[i];
}
inline long long getsum(long long x,long long p){
    return x*(x+1)/2%p;
}
long long getnmt(long long A,long long B,long long C,long long t){
    if(t==0){
        long long res=ksm(fac1[A],B*C%(p-1),p);
        return res;
    } 
    if(t==1){
        long long res=ksm(fac2[A],getsum(B,p-1)*getsum(C,p-1)%(p-1),p);
        return res;
    }
    if(t==2){
        long long l=1,r,res=1;
        while(l<=min(min(A,B),C)){
            r=min(min(A/(A/l),B/(B/l)),C/(C/l));
            res=res*ksm(f3[r]*g3[l-1]%p,(A/l)*(B/l)%(p-1)*(C/l)%(p-1),p)%p*ksm(fac1[A/l],(sum[r]-sum[l-1]+(p-1))%(p-1)*(B/l)%(p-1)*(C/l)%(p-1),p)%p;
            l=r+1;
        }
        return res;
    }
}
long long getres(long long A,long long B){
    long long l=1,r,res=1;
    while(l<=min(A,B)){
        r=min(A/(A/l),B/(B/l));
        res=res*ksm(f1[r]*g1[l-1]%p,(A/l)*(B/l)%(p-1),p)%p;
        l=r+1;
    }
    return res;
}
long long getdmt(long long A,long long B,long long C,long long t){
    if(t==0){ 
        long long l=1,r,res=1;
        while(l<=min(A,B)){
            r=min(A/(A/l),B/(B/l));
            res=(res*ksm(f1[r]*g1[l-1]%p,(A/l)*(B/l)%(p-1),p))%p;
            l=r+1;
        }
        return ksm(res,C,p);
    }
    if(t==1){
        long long l=1,r,res=1;
        while(l<=min(A,B)){
            r=min(A/(A/l),B/(B/l));
            res=(res*ksm(f2[r]*g2[l-1]%p,getsum(A/l,p-1)*getsum(B/l,p-1)%(p-1),p))%p;
            l=r+1;
        }
        return ksm(res,getsum(C,p-1),p);
    }
    if(t==2){
        long long l=1,r,res=1;
        while(l<=min(min(A,B),C)){
            r=min(min(A/(A/l),B/(B/l)),C/(C/l));
            res=res*ksm(getres(A/l,B/l),(sum[r]-sum[l-1]+(p-1))%(p-1)*(C/l)%(p-1),p)%p*ksm(f3[r]*g3[l-1]%p,(A/l)*(B/l)%(p-1)*(C/l)%(p-1),p)%p;
            l=r+1;
        }
        return res;
    }
}
long long getans(long long A,long long B,long long C,long long t){
    long long nmt=getnmt(A,B,C,t)*getnmt(B,A,C,t)%p,dmt=getdmt(A,B,C,t)*getdmt(A,C,B,t)%p;
    return nmt*inv(dmt,p)%p;
}
int main(){
    scanf("%lld%lld",&T,&p);
    init();
    while(T--){
        scanf("%lld%lld%lld",&A,&B,&C);
        printf("%lld %lld %lld\n",getans(A,B,C,0),getans(A,B,C,1),getans(A,B,C,2)); 
    }
    return 0;
}

10.例题8：[SDOI2018]旧试题
题意：求$\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C d(i\cdot j\cdot k)\bmod (10^9+7)$其中$d(x)=\sum_{w\mid x}1$即$x$的因数个数
我有一个绝妙的想法，但是这里空白太小，写不下
题解：[SDOI2018]旧试题解题报告
代码：

#include<stdio.h>
#include<vector>
#include<string.h>
using namespace std;
const long long maxn=1000005,maxm=1000005,mod=1000000007;
long long i,j,k,a,b,c,e,mx,mn,T,cnt,tot,ans;
long long lst[maxn],d[maxn],p[maxn],mu[maxn],ok[maxn],ord[maxn],deg[maxn],f[maxn],from[maxm],to[maxm],lcm[maxm],mrk[maxn];
vector<long long>v[maxn],w[maxn];
inline long long min(long long a,long long b){
    return a<b? a:b;
}
inline long long max(long long a,long long b){
    return a>b? a:b;
}
inline long long gcd(long long a,long long b){
    return b==0? a:gcd(b,a%b);
}
int main(){
    p[1]=mu[1]=1;
    for(i=2;i<maxn;i++){
        if(p[i]==0)
            lst[++cnt]=i,mu[i]=-1;
        for(j=1;j<=cnt;j++){
            if(i*lst[j]>=maxn)
                break;
            p[i*lst[j]]=1;
            if(i%lst[j]==0){
                mu[i*lst[j]]=0;
                break;
            }
            mu[i*lst[j]]=-mu[i];
        }
    }
    for(i=1;i<maxn;i++)
        if(mu[i]!=0)
            ok[++tot]=i,ord[i]=tot;
    for(i=1;i<maxn;i++){
        for(j=1;i*j<maxn;j++)
            d[i*j]++;
        f[i]=(f[i-1]+d[i])%mod;
    }
    scanf("%lld",&T);
    while(T--){
        memset(deg,0,sizeof(deg));
        memset(v,0,sizeof(v));
        memset(w,0,sizeof(w));
        scanf("%lld%lld%lld",&a,&b,&c);
        mn=min(min(a,b),c),mx=max(max(a,b),c);
        e=ans=0;
        for(i=1;i<=tot;i++){
            if(ok[i]>mx)
                break;
            for(j=1;j<=tot;j++){
                if(ok[i]*ok[j]>mx)
                    break;
                if(mu[ok[i]*ok[j]]==0)
                    continue;
                for(k=j+1;k<=tot;k++){
                    if(ok[i]*ok[j]*ok[k]>mx)
                        break;
                    if(mu[ok[i]*ok[k]]==0||gcd(ok[j],ok[k])>1)
                        continue;
                    from[++e]=ord[ok[i]*ok[j]],to[e]=ord[ok[i]*ok[k]],lcm[e]=ok[i]*ok[j]*ok[k];
                    deg[from[e]]++,deg[to[e]]++;
                } 
            }
        }
        for(i=1;i<=tot;i++){
            if(ok[i]>mn)
                break;
            ans+=mu[ok[i]]*mu[ok[i]]*mu[ok[i]]*f[a/ok[i]]*f[b/ok[i]]*f[c/ok[i]];
        }
        for(i=1;i<=e;i++){
            v[from[i]].push_back(to[i]),w[from[i]].push_back(lcm[i]);
            v[to[i]].push_back(from[i]),w[to[i]].push_back(lcm[i]);
        }
        for(i=1;i<=tot;i++){
            if(ok[i]>min(a,b))
                break;
            for(j=0;j<v[i].size();j++){
                long long x=ok[i],y=ok[v[i][j]],z=w[i][j];
                ans=(ans+mu[x]*mu[y]*mu[y]*f[a/z]*f[b/z]*f[c/y]%mod+mod)%mod;
                ans=(ans+mu[x]*mu[x]*mu[y]*f[a/x]*f[b/z]*f[c/z]%mod+mod)%mod;
                ans=(ans+mu[x]*mu[x]*mu[y]*f[a/z]*f[b/x]*f[c/z]%mod+mod)%mod;
            }
        }
        memset(v,0,sizeof(v));
        memset(w,0,sizeof(w));
        for(i=1;i<=e;i++){
            if(deg[from[i]]>=deg[to[i]])
                v[from[i]].push_back(to[i]),w[from[i]].push_back(lcm[i]);
            else v[to[i]].push_back(from[i]),w[to[i]].push_back(lcm[i]);
        }
        for(i=1;i<=tot;i++){
            if(ok[i]>mx)
                break;
            for(j=0;j<v[i].size();j++)
                mrk[v[i][j]]=w[i][j];
            for(j=0;j<v[i].size();j++){
                long long x=v[i][j];
                for(k=0;k<v[x].size();k++){
                    long long y=v[x][k],p=mrk[y],q=w[i][j],r=w[x][k];
                    if(mrk[y]==0)
                        continue;
                    long long st1,st2,st3,st4,st5,st6;
                    st1=f[a/p]*f[b/q]*f[c/r]%mod;
                    st2=f[a/p]*f[b/r]*f[c/q]%mod;
                    st3=f[a/q]*f[b/p]*f[c/r]%mod;
                    st4=f[a/q]*f[b/r]*f[c/p]%mod;
                    st5=f[a/r]*f[b/p]*f[c/q]%mod;
                    st6=f[a/r]*f[b/q]*f[c/p]%mod;
                    ans=(ans+mu[ok[i]]*mu[ok[x]]*mu[ok[y]]*(st1+st2+st3+st4+st5+st6)%mod+mod)%mod;
                }
            }
            for(j=0;j<v[i].size();j++)
                mrk[v[i][j]]=0;
        }
        printf("%lld\n",(ans+mod)%mod);
    }
    return 0;
}

11.例题9：P4240 毒瘤之神的考验
题意：求$\sum_{i=1}^n\sum_{j=1}^m\varphi(i\cdot j)\bmod 998244353$
我有一个绝妙的想法，但是这里空白太小，写不下
题解：P4240 毒瘤之神的考验解题报告

#include<stdio.h>
#include<vector>
#define int long long
using namespace std;
const int maxn=100005,mod=998244353,maxt=55;
int T,n,m,cnt,t;
int a[maxn],p[maxn],miu[maxn],phi[maxn],nphi[maxn],f[maxn];
vector<int>g[maxn],S[maxt][maxt];
int ksm(int a,int b){
    int res=1;
    while(b){
        if(b&1)
            res=res*a%mod;
        a=a*a%mod,b>>=1;
    }
    return res;
}
void init(){
    t=50;
    p[1]=miu[1]=phi[1]=1;
    for(int i=2;i<maxn;i++){
        if(p[i]==0)
            a[++cnt]=i,miu[i]=-1,phi[i]=i-1;
        for(int j=1;j<=cnt;j++){
            if(i*a[j]>=maxn)
                break;
            p[i*a[j]]=1;
            if(i%a[j]==0){
                miu[i*a[j]]=0;
                phi[i*a[j]]=phi[i]*a[j];
                break;
            }
            miu[i*a[j]]=-miu[i];
            phi[i*a[j]]=phi[i]*(a[j]-1);
        }
    }
    for(int i=1;i<maxn;i++)
        nphi[i]=ksm(phi[i],mod-2);
    for(int i=1;i<maxn;i++)
        for(int j=1;i*j<maxn;j++)
            f[i*j]=(f[i*j]+i*nphi[i]%mod*miu[j])%mod;
    for(int i=1;i<maxn;i++){
        g[i].push_back(0);
        for(int j=1;i*j<maxn;j++)
            g[i].push_back((g[i][j-1]+phi[i*j])%mod);
    }
    for(int i=1;i<=t;i++)
        for(int j=1;j<=t;j++){
            S[i][j].push_back(0);
            for(int k=1;k*max(i,j)<maxn;k++)
                S[i][j].push_back((S[i][j][k-1]+f[k]*g[k][i]%mod*g[k][j]%mod)%mod);
        }
}
int solve(int n,int m){
    int res=0,l,r;
    for(int i=1;i<=max(n,m)/t;i++)
        res=(res+f[i]*g[i][n/i]%mod*g[i][m/i]%mod)%mod;
    l=max(n,m)/t+1;
    while(l<=min(n,m)){
        r=min(n/(n/l),m/(m/l));
        res=(res+(S[n/l][m/l][r]-S[n/l][m/l][l-1]+mod)%mod)%mod;
        l=r+1;
    }
    return res;
}
signed main(){
    init();
    scanf("%lld",&T);
    while(T--){
        scanf("%lld%lld",&n,&m);
        printf("%lld\n",solve(n,m));
    }
    return 0;
}

完结撒花！

续集：杜教筛学习笔记