莫比乌斯反演

# 莫比乌斯反演以及数论筛法

#### cdqz xby

### 1. 数论函数

定义： 定义域为正整数集的函数

#### 积性函数

对于一个数论函数 $f$， 如果满足 $\forall x\bot y(x\bot y\rightarrow \gcd(x,y)=1),f(xy)=f(x)f(y)$， 那么 $f$ 被称为**积性函数**。

如果对于任意 $x,y$ 满足 $f(xy)=f(x)f(y)$ 那么 $f$ 被称为**完全积性函数**。

#### 常见积性函数

单位函数： $\epsilon (n)=[n=1]$

常数函数： $\operatorname{I}(n)=1$

恒等函数：$\operatorname{Id}(n)=n$

幂函数： $\operatorname{Id}_k(n)=n^k$ 

莫比乌斯函数： $\mu(n)=\begin{cases}1& n=1\\(-1)^k&n=\prod_{i=1}^kp_i\wedge \forall i,j\in[1,k],p_i\neq p_j\\ 0&\exists d>1,d^2|n  \end{cases}$

欧拉函数： $\varphi(n)=\sum_{i=1}^{n-1}[\gcd(i,n)=1]$

除数函数： $\sigma_k(n)=\sum_{d|n}d^k$

其中， 前四个函数为完全积性函数。

**大部分积性函数可以通过线性筛在 $O(n)$ 的时间内筛出。**

#### Dirichlet 卷积

定义： 如果 $h=f*g$， $h(n)=\sum_{d|n}f(n)g(\dfrac{n}{d})$。

性质：两个积性函数的狄利克雷卷积仍然是积性函数。

证明：
$$
\begin{aligned}
h(x)&=\sum_{d|x}f(d)g(\frac{n}{d}),\ h(y)=\sum_{k|y}f(k)g(\frac{n}{k})\\
h(x)h(y)&=\left(\sum_{d|x}f(d)g(\frac{n}{d}) \right)\left(\sum_{k|y}f(k)g(\frac{n}{k}) \right)\\&
=\sum_{d|x,k|y}f(d)f(k)g(\frac{x}{d})g(\frac{y}{k})\\
&=\sum_{d|x,k|y}f(dk)g(\frac{xy}{dk})=\sum_{d|xy}f(d)g(\frac{xy}{d})=h(xy)
\end{aligned}
$$

性质： 狄利克雷卷积满足结合律， 分配律和交换律。

证明从略。

##### 常用性质 & 等式

$\forall f,f*\epsilon=f$

证明显然。

$\forall f,\exists g, f*g=\epsilon$

这里的 $g$ 称为 $f$ 的**逆**， 且若 $f$ 为积性函数， 则 $g$ 也为积性函数。

证明（存在逆）：

构造 $g(n)=\dfrac{1}{f(1)}\left([n=1]-\sum_{d|n\wedge d>1}f(d)g(\dfrac{n}{d}) \right)$。

可以证明 
$$
\begin{aligned}
h(n) & =\sum_{d|n}f(d)g(\dfrac{n}{d})\\
&=f(1)g(n)+\sum_{d|n,d>1}f(d)g(\frac{n}{d})\\
&=[n=1]-\sum_{d|n,d>1}f(d)g(\frac{n}{d})+\sum_{d|n,d>1}f(d)g(\frac{n}{d})
\\
&=[n=1]=\epsilon
\end{aligned}
$$

$\mu *\operatorname{I}=\epsilon$

证明：

即证明： $\epsilon (n)=\sum_{d|n}\mu (d)$。

当 $n=1$ 时， 等式显然成立。

当 $n>1$ 时， 设 $n=\prod_i^mp_i^{k_i}$。 有 $\sum_{d|n}\mu(d)=\sum_{i=0}^m (-1)^i\tbinom{m}{i}=0$。

$\varphi*\operatorname{I} =\operatorname{Id}$

证明：

即证明 $n=\operatorname{Id}(n)=\sum_{d|n}\varphi(d)$。 将 $\sum_{d|n}\varphi(d)$ 命名为 $(1)$ 式。

考虑 $n$ 个分数： $\dfrac{1}{n},\dfrac{2}{n}\dots \dfrac{n}{n}$。

考虑第 $i$ 个分数 $\dfrac{i}{n}$， 设其化到最简分数为 $\dfrac{a}{b}$。

那么一定有： $\gcd(a,b)=1,b|n$。 所以每一个分数都会对 $(1)$ 式有 $1$ 的贡献。

对于任意 $a|n,b|n\wedge a\neq b$， 不会存在 $c,d$ 使得 $c<a,d<b,\gcd(a,c)=1,\gcd(b,d)=1$ 且 $\dfrac{c}{a}=\dfrac{d}{b}$。 证明显然。

对于任意 $a|n$， $\forall b$ 满足 $b<a\wedge \gcd(a,b)=1$， 总存在一个数 $c$ 使得 $c<n\wedge\dfrac{b}{a}=\dfrac{c}{n}$。

综上可得 $\sum_{d|n}\varphi(d)=n=\operatorname{Id}(n)$。

$\mu*\operatorname{Id} =\varphi$

证明1：

$\mu*\operatorname{Id}=\mu*\operatorname{I}*\varphi=\epsilon *\varphi=\varphi$。

证明2：

考虑容斥原理， $\mu$ 即为容斥系数， 结论显然。

$\sigma_0(ij)=\sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]$

证明：

对于任意质数 $p$， $\sum_{x|p^a}\sum_{y|p^b}[\gcd(x,y)=1]=a+b+1=\sigma_0(p^{a+b})$。

现在任取两个正整数 $n=\prod_i p_i^{k_i},m=\prod_i p_i^{t_i}$，

则：

$$
\begin{aligned}
\sum_{x|n}\sum_{y|m}[\gcd(x,y)=1]
&=\prod_i\sigma_0(p_i^{k_i+t_i})=\sigma_0(nm)
\end{aligned}
$$

$\sigma_1(ij)=\sum_{x|i}\sum_{y|j}[\gcd(i,j)=1]\dfrac{xj}{y}$

证明同上。

$\mu(ij)=\mu(i)\mu(j)[\gcd(i,j)=1]$

证明显然。

$\varphi(ij)=\dfrac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi(\gcd(i,j))}$

证明：

将 $\varphi$ 的计算式带入计算即可。

### 2.莫比乌斯反演

如果有 $f(n)=\sum_{d|n}g(d)$， 则有 $g(n)=\sum_{d|n}\mu(d)f(\dfrac{n}{d})$。

证明：

已知 $f=g*\operatorname{I}$， 那么 $f*\mu=g*\operatorname{I}*\mu=g$， 得证。

同样的， 如果有 $f(n)=\sum_{n|d}g(d)$， 则有 $g(n)=\sum_{n|d}\mu(\dfrac{n}{d})f(d)$， 证明同上。

#### 例题

[简单题](https://www.luogu.com.cn/problem/P6156)

> 求 $\sum_{i=1}^n\sum_{j=1}^n(i+j)^k\gcd(i,j)\mu^2(\gcd(i,j))$。

推式子：
$$
\begin{aligned}
&\ \ \ \ \ \sum_{i=1}^n\sum_{j=1}^n(i+j)^k\gcd(i,j)\mu^2(\gcd(i,j))
\\
&=\sum_{i=1}^n\sum_{j=1}^n(i+j)^k\sum_{d}[\gcd(i,j)=d]d\mu^2(d)
\\
&=\sum_{d=1}^nd^{k+1}\mu^2(d)\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor\frac{n}{d} \rfloor}(i+j)^k[\gcd(i,j)=1]\\
&=\sum_{d=1}^nd^{k+1}\mu^2(d)\sum_{i=1}^{\lfloor\frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor\frac{n}{d} \rfloor}(i+j)^k\sum_{p|\gcd(i,j)}\mu(p)
\\
&=\sum_{d=1}^nd^{k+1}\mu^2(d)\sum_{p=1}^{\lfloor\frac{n}{d} \rfloor}\mu(p)p^k S(\lfloor\frac{n}{pd} \rfloor)
\\
&=\sum_{T=1}^nT^kS(\lfloor\frac{n}{T} \rfloor)\sum_{d|T}d\mu^2(d)\mu(\frac{T}{d})
\end{aligned}
$$
其中 $T=pd,S(n)=\sum_{i=1}^n\sum_{j=1}^n(i+j)^k$。

现在需要快速处理 $S(n)$ 和 $\sum_{d|T}d\mu^2(d)\mu(\frac{T}{d})$。

首先是 $S(n)$：

令 $F(n)=\sum_{i=1}^ni^k$， 则 $S(n)=\sum_{i=1}^n\sum_{j=1}^n(i+j)^k=\sum_{i=n+1}^{2n}(F(i)-F(i-n))=\sum_{i=n+1}^{2n}F(n)-\sum_{i=1}^nF(n)$。

令 $G(n)=\sum_{i=1}^nF(n)$， 那么 $S(n)=G(2n)-2G(n)$。 线性筛出自然数的 $k$ 次幂做两次前缀和即可。

接着是 $f(n)=\sum_{d|n}d\mu^2(d)\mu(\frac{n}{d})$：

首先， 这显然是一个积性函数。 考虑线性筛的过程。 假设现在枚举到的 $n$ 的最小质因子为 $p$， 且 $p$ 的次数为 $k$。

则 $f(n)=f(p^k)\times f(\dfrac{n}{p^k})$。 现在需要求出 $f(p^k)$。

现在对 $k$ 的值进行分类讨论。

若 $k=0$， 显然 $f(1)=1$。

若 $k=1$， $f(p)=1+p$。

若 $k=2$， $f(p^2)=-p$。

若 $k\geqslant 3$， 则 $\mu(d)$ 与 $\mu(\dfrac{p^k}{d})$ 中至少有一个有 $p^2$  这个因子， $f(p^k)=0$。

那么现在可以对 $S$ 进行整除分块解决问题了。

``` cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int N=1e7+5,mod=998244353;
int p[700005],cnt; bool flag[N];
int F[N],f[N];
int qpow(int a,int b){
    int ans=1;
    while(b){
        if(b&1) ans=1ll*ans*a%mod;
        a=1ll*a*a%mod; b>>=1;
    }
    return ans;
}
int n; ll k;
void sieve(int n){
    F[1]=f[1]=1;
    for(int i=2; i<=n; i++){
        if(!flag[i]) p[++cnt]=i,F[i]=qpow(i,k),f[i]=i-1;
        for(int j=1,t; j<=cnt&&i*p[j]<=n; j++){
            t=i*p[j];
            flag[t]=1; F[t]=1ll*F[i]*F[p[j]]%mod;
            if(i%p[j]==0){
                int q=i/p[j];
                if(q%p[j]!=0) f[t]=1ll*(mod-p[j])*f[q]%mod;
                else f[t]=0; 
                break;
            }
            f[t]=1ll*f[i]*f[p[j]]%mod;
        }
    }
    for(int i=1; i<=n; i++) f[i]=(f[i-1]+1ll*f[i]*F[i]%mod)%mod,F[i]=(F[i-1]+F[i])%mod;
    for(int i=1; i<=n; i++) F[i]=(F[i]+F[i-1])%mod;
}
int ans;
int main(){
    scanf("%d%lld",&n,&k); k%=(mod-1);
    sieve(n*2);
    for(int l=1,r; l<=n; l=r+1){
        r=n/(n/l);
        ans=(ans+1ll*(f[r]-f[l-1]+mod)%mod*((F[n/l*2]-F[n/l]*2%mod+mod)%mod))%mod;
    }
    printf("%d\n",ans);
    return 0;
}
```

[数字表格 ](https://www.luogu.com.cn/problem/P3704)

>计算：
>
>$\prod_{i=1}^n\prod_{j=1}^mf(\gcd(i,j))\bmod(10^9+7)$
>
>多测， $T\leqslant 1000,1\leqslant n,m\leqslant 10^6$。
>
>$f(i)$ 表示斐波那契数列的第 $i$ 项， 即 $f(0)=0,f(1)=1,f(i)=f(i-1)+f(i-2)$。

连乘与连加大抵相似， 还是推式子：
$$
\begin{aligned}
\large
\prod_{i=1}^n\prod_{j=1}^mf(\gcd(i,j))&=\prod_{i=1}^n\prod_{j=1}^m\prod_{d=1}^{\min(n,m)}[\gcd(i,j)=d]f(d)\\
&=\prod_{d}^{\min(n,m)}f(d)^{\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]}\\
&=\prod_{d}^{\min(n,m)}f(d)^{\sum_{p=1}^{\lfloor\frac{\min(n,m)}{d}\rfloor}\mu(p)\lfloor\frac{n}{pd} \rfloor\lfloor\frac{m}{pd} \rfloor} \\
\end{aligned}
$$
令 $H=dp$。
$$
\begin{aligned}
\large 
\prod_{i=1}^n\prod_{j=1}^mf(\gcd(i,j))&=\prod_{H=1}^{\min(n,m)}\left(\prod_{d|H}f(d)^{\mu(\frac{H}{d})} \right)^{\lfloor\frac{n}{H} \rfloor\lfloor\frac{m}{H} \rfloor}
\end{aligned}
$$
令 $g(H)=\prod_{d|H} f(d)^{\mu(\frac{H}{d})}$， 这个可以通过枚举倍数在 $O(n\ln n)$ 的复杂度下预处理出来。

剩下的就是一个整除分块了， 总时间复杂度 $O(Tn\log mod)$。

```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define pr pair<int,int>
#define mk make_pair
#define lowbit(x) (x&(-x))
#define pb emplace_back
constexpr int N=1e6+5,mod=1e9+7,V=1e6;
int read(){
    int x=0,f=1; char c=getchar();
    while(('0'>c||c>'9')&&c!='-') c=getchar();
    if(c=='-') f=0,c=getchar();
    while('0'<=c&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return f?x:-x;
}
int qpow(int a,int b){
    int ans=1;
    while(b){
        if(b&1) ans=1ll*ans*a%mod;
        a=1ll*a*a%mod; b>>=1;
    }
    return ans;
}
int mu[N],n,m,T,f[N],g[N],p[N],cnt,qf[N],Inv[N];
bool flag[N];
void sieve(const int n=V){
    mu[1]=1;
    for(int i=2; i<=n; i++){
        if(!flag[i]) p[++cnt]=i,mu[i]=-1;
        for(int j=1; j<=cnt&&i*p[j]<=n; j++){
            flag[i*p[j]]=1;
            if(i%p[j]==0)
                break;
            mu[i*p[j]]=-mu[i];
        }
    }
}
int main(){
    sieve(); f[0]=0,f[1]=1;
    for(int i=2; i<=V; i++) f[i]=(f[i-1]+f[i-2])%mod;
    for(int i=1; i<=V; i++) qf[i]=qpow(f[i],mod-2);
    for(int i=0; i<=V; i++) g[i]=1;
    for(int i=1; i<=V; i++){
        if(mu[i]==0) continue;
        for(int j=i; j<=V; j+=i)
            if(mu[i]==1) g[j]=1ll*g[j]*f[j/i]%mod;
            else g[j]=1ll*g[j]*qf[j/i]%mod;
    } 
    int sum=1; for(int i=1; i<=V; i++) sum=1ll*sum*g[i]%mod;
    Inv[V]=qpow(sum,mod-2); for(int i=V; i; i--) Inv[i-1]=1ll*g[i]*Inv[i]%mod;
    for(int i=2; i<=V; i++) g[i]=1ll*g[i]*g[i-1]%mod;
    T=read();
    while(T--){
        n=read(),m=read(); int ans=1;
        for(int l=1,r; l<=min(n,m); l=r+1){
            r=min(n/(n/l),m/(m/l));
            ans=1ll*ans*qpow(1ll*g[r]*Inv[l-1]%mod,1ll*(n/l)*(m/l)%(mod-1))%mod;
        }
        printf("%d\n",ans);
    }
    return 0;
}

```

[【UR #5】怎样跑得更快](https://uoj.ac/problem/62)

> 解方程：
> $\sum_{j=1}^n\gcd(i,j)^c\operatorname{lcm(i,j)^d}x_j\equiv b_i\pmod {998244353}$

首先化简：

$\sum_{j=1}^n \gcd(i,j)^{c-d}i^dj^dx_j\equiv b_i$

令 $y_j=x_j\times j^d,c_i=\dfrac{b_i}{i^d},f(x)=x^{c-d}$。

则有： $\sum_{i=1}^nf(\gcd(i,j))y_j\equiv c_i$。

接下来运用到一个莫比乌斯反演题经常用到的技巧： 把函数拆开。

令 $f=\operatorname{I}*f'$， 即 $f(n)=\sum_{d|n}f'(d)$。

根据莫比乌斯反演， $f'(n)=\sum_{d|n}f(n)$。 那么计算 $f'$ 可以通过对每个数枚举倍数做到 $O(n\ln n)$ 的复杂度。

事实上， $f'(n)=f(n)-\sum_{d>1\wedge d|n}f'(d)$， 也可以在 $O(n \ln n)$ 的复杂度下计算出来。

接着推式子：
$$
\begin{aligned}
\sum_{i=1}^nf(\gcd(i,j))y_j &=\sum_{i=1}^n\sum_{d|\gcd(i,j)}f'(d) y_j\\
&=\sum_{d|i}f'(d) \sum_{j>d\wedge d|j}y_j

\end{aligned}
$$
令 $t_d=\sum_{j>d\wedge d|j}y_j$， 那么 $\sum_{d|i}f'(d)t_d\equiv c_i$。

那么可以用求 $f'$ 的方法求出 $f'(d)t_d$。 接着就可以求出 $t_d$。

通过与求 $f'$ 相近的方法可以推出求 $y$ 的方法：

$y_i=t_i-\sum_{j>i\wedge i|j}y_j$。

下标从大往小枚举， 减去倍数即可。

最后是判无解的问题。 只有一种情况出现无解： 不存在逆元。

显然 $i^d,j^d$ 一定存在逆元， 只可能 $f’$ 不存在逆元。 假如 $f'(d)$ 不存在逆元， 分两种情况：

若 $f'(d)t_d=0$ 则有无穷多组解。

若 $f'(d)t_d\neq 0$， 则无解。

总时间复杂度 $O(qn\ln n)$。

```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define pr pair<int,int>
#define mk make_pair
#define lowbit(x) (x&(-x))
#define pb emplace_back
constexpr int N=1e5+5,mod=998244353;
inline void dec(int &a,const int b){a-=b; if(a<0) a+=mod;}
int read(){
    int x=0,f=1; char c=getchar();
    while(('0'>c||c>'9')&&c!='-') c=getchar();
    if(c=='-') f=0,c=getchar();
    while('0'<=c&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return f?x:-x;
}
int qpow(int a,int b){
    int ans=1;
    while(b){
        if(b&1) ans=1ll*ans*a%mod;
        a=1ll*a*a%mod; b>>=1;
    }
    return ans;
}
int n,q,c,d,f[N],g[N],ans[N],t[N],Inv[N];
bool check(){
    for(int i=1; i<=n; i++) if(1ll*ans[i]*f[i]%mod!=t[i]) return 0;
    return 1;
}
int main(){
    n=read(),c=read(),d=read(),q=read(); c-=d; c=(c%(mod-1)+mod-1)%(mod-1);
    for(int i=1; i<=n; i++) g[i]=qpow(i,c);
    for(int i=1; i<=n; i++) for(int j=i+i; j<=n; j+=i) dec(g[j],g[i]);
    for(int i=1; i<=n; i++) Inv[i]=qpow(qpow(i,d),mod-2);
    for(int i=1; i<=n; i++) f[i]=g[i],g[i]=qpow(g[i],mod-2);
    while(q--){
        for(int i=1; i<=n; i++) t[i]=1ll*read()*Inv[i]%mod;
        for(int i=1; i<=n; i++) for(int j=i+i; j<=n; j+=i) dec(t[j],t[i]);
        for(int i=1; i<=n; i++) ans[i]=1ll*t[i]*g[i]%mod;
        if(!check()){puts("-1");continue;}
        for(int i=n; i>=1; i--) for(int j=i+i; j<=n; j+=i) dec(ans[i],ans[j]);
        for(int i=1; i<=n; i++) ans[i]=1ll*ans[i]*Inv[i]%mod;
        for(int i=1; i<=n; i++) printf("%d ",ans[i]);
        puts("");
    }
    return 0;
}

```

事实上存在 $O(qn\log\log n)$ 的解法。 注意到定义 $f'$ 等几个函数时是 Dirichlet 前缀和， 相当于一个高位前缀和。 只需要做一个高维差分即可（每个质数是一个维）。

### 3.杜教筛

杜教筛解决的是积性函数求前缀和的问题。 它可以在 $O(n^{\frac{2}{3}})$ 的时间求出部分积性函数的前缀和。

令 $S_f(n)=\sum_{i=1}^nf(i)$。

假如现在要求积性函数 $f$ 的前缀和 $S_f$， 我们构造另外两个积性函数 $g,h$， 使得 $h=f*g$。

那么：
$$
\begin{aligned}
S_h(n)&=\sum_{i=1}^nh(n)\\
&=\sum_{i=1}^n\sum_{j|i}g(j)f(\frac{i}{j})\\
&=\sum_{j=1}^n\sum_{j|i\wedge i\leqslant n}g(j)f(\frac{i}{j})\\
&=\sum_{j=1}^n\sum_{k\times j\leqslant n}g(j)f(k)\\
&=\sum_{j=1}^ng(j)S_f(\lfloor\frac{n}{j}\rfloor)
\end{aligned}
$$
就有:
$$
S_f(n)=\frac{S_h(n)-\sum_{j=2}^ng(j) S_f(\lfloor\frac{n}{j}\rfloor)}{g(1)}
$$
积性函数在 $1$ 处的值都为 $1$。 所以上式可以简化为：
$$
S_f(n)=S_h(n)-\sum_{j=2}^ng(j) S_f(\lfloor\frac{n}{j}\rfloor)
$$
可以发现， 右边那部分可以通过整除分块递归处理。 这要求 $h,g$ 两个函数可以快速求前缀和。

假设 $h,g$ 都可以 $O(1)$ 求前缀和，那么复杂度就是 $O(n^{\frac{3}{4}})$。

当然， 复杂度是可以继续优化的。 用线性筛预处理出前 $n^\frac{2}{3}$ 的答案， 那么复杂度就降到了 $O(n^\frac{2}{3})$。 具体证明需要一些微积分知识。

如果 $h,g$ 也能通过杜教筛求前缀和， 因为是并列的， 所以复杂度不变。

还有一个问题： 如果记忆化时使用 $map$， 复杂度就会多一个 $\log$。 当然大部分情况都是没什么问题的。 除了使用哈希表之外， 还有一种解决方法。

注意到： 
$$
\lfloor\frac{\lfloor\frac{x}{y} \rfloor}{z} \rfloor=\lfloor\frac{x}{yz} \rfloor
$$

所以， 每次要求的前缀和 $S_f(m)$ 中的 $m$ 一定是某个 $\lfloor\dfrac{n}{x} \rfloor$， 并且当 $x>n^\frac{1}{3}$ 时会直接查表， 所以只用开一个数组记录 $x<n^\frac{1}{3}$ 时的答案即可（一般用 $map$ 就够了）。

#### 例题

[模板题](https://www.luogu.com.cn/problem/P4213)

要求求出 $\mu,\varphi$ 的前缀和。

根据前面数论函数的知识， $\mu *\operatorname{I}=\epsilon,\varphi *\operatorname{I}=\operatorname{Id}$， 直接杜教筛就可以了。

```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define pr pair<int,int>
#define mk make_pair
#define lowbit(x) (x&(-x))
#define pb emplace_back
constexpr int M=5e6,N=1305;
int read(){
    int x=0,f=1; char c=getchar();
    while(('0'>c||c>'9')&&c!='-') c=getchar();
    if(c=='-') f=0,c=getchar();
    while('0'<=c&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return f?x:-x;
}
int n,T,mu[M+10],p[M+10],cnt,ansmu[N];
ll phi[M+10],ansphi[N];
bool vis1[N],vis2[N],flag[M+10];
void sieve(const int n=M){
    mu[1]=phi[1]=1;
    for(int i=2; i<=n; i++){
        if(!flag[i]) p[++cnt]=i,mu[i]=-1,phi[i]=i-1;
        for(int j=1; j<=cnt&&i*p[j]<=n; j++){
            flag[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=p[j]*phi[i];
                break;
            }
            mu[i*p[j]]=-mu[i]; phi[i*p[j]]=phi[i]*(p[j]-1);
        }
    }
    for(int i=2; i<=n; i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
ll solve1(int x){
    if(x<=M) return phi[x];
    const int v=n/x;
    if(vis1[v]) return ansphi[v];
    vis1[v]=1; ll &e=ansphi[v];
    e=x*(x+1ll)/2;
    for(unsigned l=2,r; l<=(unsigned)x; l=r+1){
        r=(unsigned)x/(x/l);
        e-=solve1(x/l)*(r-l+1);
    }
    return e;
}
int solve2(int x){
    if(x<=M) return mu[x];
    const int v=n/x;
    if(vis2[v]) return ansmu[v];
    vis2[v]=1; int &e=ansmu[v];
    e=1;
    for(unsigned l=2,r; l<=x; l=r+1){
        r=x/(x/l);
        e-=solve2(x/l)*(r-l+1);
    }
    return e;
}
int main(){
    T=read();sieve();
    while(T--){
        n=read();
        memset(vis1,0,sizeof vis1);
        memset(vis2,0,sizeof vis2);
        printf("%lld %d\n",solve1(n),solve2(n));
    }
    return 0;
}

```

[简单的数学题](https://www.luogu.com.cn/problem/P3768)

>  求 
> $$
> \left(\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\right)\bmod p
> $$
>  $n\leqslant 10^{10}$。

用 $\operatorname{Id}=\varphi*\operatorname{I}$ 将 $\gcd$ 拆开可以得到原式的等价形式：
$$
\sum_{d=1}^n\varphi(d)d^2S(\lfloor\frac{n}{d}  \rfloor)
$$
其中 $S(n)=\dfrac{n\times (n+1)}{2}$。

对 $S$ 进行数论分块， 现在的问题时如何求 $\varphi(i)i^2$ 的前缀和。 由于数据范围较大， 使用杜教筛解决。

现在已知 $f=\varphi\operatorname{Id}^2$， 如何确定 $g,h$ 呢？

**引理**

已知完全积性函数 $f$ 和数论函数 $g,h$， 那么 $(f\cdot g)*(f\cdot h)=f\cdot(g*h)$。

证明：
$$
\begin{aligned}
((f\cdot g)*(f\cdot h))(n)&=\sum_{d|n}f(n)g(n)f(\frac{n}{d})h(\frac{n}{d})\\
&=f(n)\sum_{d|n}g(n)h(\frac{n}{d})\\
&=f\cdot(g*h)
\end{aligned}
$$

现在已知 $\operatorname{Id}^2$ 为完全积性函数， $\varphi*\operatorname{I}=\operatorname{Id}$， 可以得知 $g=\operatorname{Id}^2\cdot\operatorname{I},h=\operatorname{Id}^2\cdot(\varphi*\operatorname{I})=\operatorname{Id}^3$。

那么问题就解决了。

``` cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define pr pair<int,int>
#define mk make_pair
#define lowbit(x) (x&(-x))
#define pb emplace_back
constexpr int N=1e7+5,V=1e7;
int read(){
    int x=0,f=1; char c=getchar();
    while(('0'>c||c>'9')&&c!='-') c=getchar();
    if(c=='-') f=0,c=getchar();
    while('0'<=c&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return f?x:-x;
}
ll n,s[N];
int p[N],cnt,mod,phi[N],inv4,inv6;
bool flag[N];
int qpow(int a,int b){
    int ans=1;
    while(b){
        if(b&1) ans=1ll*ans*a%mod;
        a=1ll*a*a%mod; b>>=1;
    }
    return ans;
}
void sieve(const int n=V){
    s[1]=1; phi[1]=1;
    for(int i=2; i<=n; i++){
        if(!flag[i]) p[++cnt]=i,phi[i]=i-1;
        for(int j=1; j<=cnt&&i*p[j]<=n; j++){
            flag[i*p[j]]=1;
            if(i%p[j]==0){
                phi[i*p[j]]=phi[i]*p[j];
                break;
            }
            phi[i*p[j]]=phi[i]*(p[j]-1);
        }
    }
    for(int i=1; i<=n; i++) s[i]=(1ll*i*i%mod*phi[i]%mod+s[i-1])%mod;
}
map <int,int> Q;
int calc3(ll n){
    n%=mod;
    return n*n%mod*(n+1)%mod*(n+1)%mod*inv4%mod;
}
int calc2(ll n){
    n%=mod;
    return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
int solve(ll n){
    if(n<=V) return s[n];
    if(Q.count(n)) return Q[n];
    int ans=calc3(n);
    for(ll l=2,r; l<=n; l=r+1){
        r=n/(n/l);
        ans-=1ll*(calc2(r)-calc2(l-1)+mod)%mod*solve(n/l)%mod; ans%=mod;
    }
    ans=(ans+mod)%mod;
    return Q[n]=ans;
}
int ans;
int main(){
    scanf("%d%lld",&mod,&n);sieve();
    inv4=qpow(4,mod-2); inv6=qpow(6,mod-2);
    for(ll l=1,r; l<=n; l=r+1){
        r=n/(n/l);
        ans+=1ll*(solve(r)-solve(l-1)+mod)%mod*calc3(n/l)%mod; ans%=mod;
    }
    printf("%d\n",ans);
    return 0;
}

```

### 4. Min_25 筛

一些约定：

- $P$ 表示全体素数集合。
- $p$ 表示某一个素数。
- $p_k$ 表示从小往大第 $k$ 个素数， 定义 $p_0=1$。

- $s_i$ 表示数 $i$ 的最小质因子。

Min_25 筛能在 $O(\dfrac{n^{\frac{3}{4}}}{\log n})$ 的时间复杂度下解决一类积性函数的前缀和问题。 有如下要求：

- $f(p)$ 能拆成多个完全积性函数的和（一般是多项式）。
- $f(p^k)$ 可以快速求出。

Min_25 筛一般会去求 $2\sim n$ 的函数值的和， 最后再加上 $1$。

##### 第一步： 求解所有质数函数值的和。

首先将 $f(p)$ 拆成若干完全积性函数。 对于每一个完全积性函数 $h(x)$：

令 $G(m,x)=\sum_{i=2}^m[i\in P\vee s_i>p_x]h(i)$， 即在 $[2,m]$ 范围内所有质数以及最小质因数大于 $p_x$ 的合数的 $h$ 的和。

对于在范围内的任意合数， 显然其最小质因子 $\leqslant \sqrt m$。

所以当 $p_x>\sqrt m$ 时， 显然有 $G(m,x)=G(m,x-1)$。 所以我们只用考虑 $x\leqslant \sqrt m$ 的情况。

考虑如何递推。 显然 $G(m,x-1)$ 比 $G(m,x)$ 多了 $p_x$ 的贡献。 考虑减去。 那么就有：
$$
G(m,x)=G(m,x-1)-h(p_x)\times\left(G(\lfloor\frac{m}{p_x}\rfloor,x-1)-\sum_{i=1}^{x-1}h(p_i) \right)
$$
显然， 所有 $G(m,|P|)$ 的和就是 $\leqslant m$ 的所有质数的函数值和。 为了方便， 记为 $G(m)$。

##### 第二步： 求出答案

令 $F(m,x)=\sum_{i=2}^m[s_i>p_x]f(i)$。 显然最后答案为 $F(n,0)+1$。

考虑如何递推
$$
F(m,x)=G(m)-\sum_{i=1}^xf(p_i)+\sum_{j>x} \sum_{k\geqslant 1\wedge p_j^k\leqslant m}f(p_j^k)\times\left(F(\lfloor\frac{m}{p_j^k}\rfloor)+[k\neq 1] \right)
$$
然后就可以递归求解了。  可以证明， 每种状态只会访问一次， 所以不用记忆化。

##### 关于实现

可以发现， $F$ 和 $G$ 的第一维其实只有 $O(\sqrt n)$ 个会用到的取值， 做一次整除分块就可以知道这些取值。

##### 复杂度分析

已知 $\leqslant n$ 的质数数量级为 $\dfrac{n}{\log n}$。

求 $G$ 的复杂度为 $\sum_{i=1}^{\sqrt {n}} O(\dfrac{\sqrt{i}}{\log {\sqrt {i}}})$。

求 $F$ 的复杂度为 $\sum_{i=1}^{\sqrt {n}} O(\dfrac{\sqrt {\frac{n}{i}}}{\log {\sqrt {\frac{n}{i}}}})$。

显然后者更大。 将分母提出来：
$$
\sum_{x=1}^{\sqrt n}O\sqrt\frac{n}{x} =\int_{1}^{\sqrt n}\sqrt \frac{n}{x}dx=O(n^\frac{3}{4})
$$
所以总复杂度为 $O(\dfrac{n^{\frac{3}{4}}}{\log n})$。

[模板题](https://www.luogu.com.cn/problem/P5325)代码：

```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
constexpr int N=1e5+5,M=1e5,mod=1e9+7;
int qpow(int a,int b){
    int ans=1;
    while(b){
        if(b&1) ans=1ll*ans*a%mod;
        a=1ll*a*a%mod; b>>=1;
    }
    return ans;
}
const int Inv6=qpow(6,mod-2);
ll n,w[N<<1];
int p[N],cnt,t[N],r[N],pret[N],prer[N];
bool flag[N];
inline int add(int x,int y){x+=y; if(x>=mod) x-=mod; return x;}
void sieve(const int n=M){
    for(int i=2; i<=n; i++){
        if(!flag[i]) p[++cnt]=i,pret[cnt]=add(pret[cnt-1],1ll*i*i%mod),prer[cnt]=add(prer[cnt-1],i);
        for(int j=1; j<=cnt&&i*p[j]<=n; j++){
            flag[i*p[j]]=1;
            if(i%p[j]==0) break;
        } 
    }
}
int G[N<<1],g1[N<<1],g2[N<<1],tot,pos1[N],pos2[N];
int F(ll x,int y){
    if(x<p[y]) return 0;
    int pos=(x<=M?pos1[x]:pos2[n/x]),ans=(G[pos]-(pret[y]-prer[y]+mod)%mod+mod)%mod;
    for(int k=y+1; k<=cnt&&1ll*p[k]*p[k]<=x; k++){
        ll sum=p[k];
        for(int j=1; sum<=x; j++,sum*=p[k])
            ans=(ans+1ll*sum%mod*((sum-1)%mod)%mod*(F(x/sum,k)+(j!=1))%mod)%mod;
    }
    return ans;
}
int main(){
    scanf("%lld",&n); sieve();
    for(ll l=1,r; l<=n; l=r+1){
        r=n/(n/l); ll v=n/l;
        if(v<=M) pos1[v]=++tot;
        else pos2[n/v]=++tot;
        w[tot]=v; v%=mod;
        g1[tot]=v*(v+1)%mod*(2*v+1)%mod*Inv6%mod;
        g2[tot]=v*(v+1)/2%mod;
        g1[tot]--,g2[tot]--;
    }
    for(int j=1; j<=cnt; j++) for(int i=1; i<=tot&&1ll*p[j]*p[j]<=w[i]; i++){
        ll nxt=w[i]/p[j];
        int pos=(nxt<=M?pos1[nxt]:pos2[n/nxt]);
        g1[i]-=1ll*p[j]*p[j]%mod*(g1[pos]-pret[j-1])%mod; g1[i]=(g1[i]%mod+mod)%mod;
        g2[i]-=1ll*p[j]*(g2[pos]-prer[j-1])%mod; g2[i]=(g2[i]%mod+mod)%mod;
    }
    for(int i=1; i<=tot; i++) G[i]=(g1[i]-g2[i]+mod)%mod;
    printf("%d\n",(F(n,0)+1)%mod);
    return 0;
}

```