Dyd's Blog

He who has a strong enough why can bear almost any how.

莫比乌斯反演

Posted on Edited on In
Symbols count in article: 16k Reading time ≈ 15 mins.

数学的恶魔

莫比乌斯反演

前置知识：数论分块

基本内容

数论分块可以快速计算一些含有除法向下取整的和式（即形如 $\sum_{i=1}^{n} f(i)g( \lfloor \frac{n}{i} \rfloor)$ 的式子），只要可以在 $O(1)$ 内计算 $f(r)-f(l)$ 或已经预处理出 $f$ 的前缀和时，可以用数论分块在 $O(\sqrt{n})$ 的时间内计算答案。

它主要利用富比尼定理（Fubini’s theorem）（听起来很厉害但其实没啥）：我们注意到，有很多的 $i$ 它们的 $\lfloor \frac{n}{i} \rfloor$ 是相等的，这启发我们分段打包计算。

举个例子，如果我们要求 $\sum_{i=1}^{n} \lfloor \frac{n}{i} \rfloor$ ，我们发现很多 $i$ 都对应同一个 $\lfloor \frac{n}{i} \rfloor$ ，它们呈块状分布，在打个表仔细观察（证明），发现对于每一个块，它的最后一个数是 $n/(n/i)$ （ $c++$ 自带下取整，一般我们称其这个式子的值为 $g(i)$ ），如对于 $\sum_{i=1}^{10} \lfloor \frac{10}{i} \rfloor$ ，分成： $(1)(2)(3)(4,5)(6,7,8,9,10)$ ，其中 $1=10/(10/1),2=10/(10/2),3=10/(10/3),5=10/(10/5),10=10/(10/10)$ ，于是，求解上述问题的代码如下（常用）：

UVA11526 H(n)

1
2
3
4
for(int l=1,r;l<=n;l=r+1){
    r=n/(n/l);
    ans+=(r-l+1)*(n/l);
}

类似的，我们可以用 $O(\sqrt{n})$ 计算 $g(\lfloor \frac{n}{i} \rfloor)$ 那么只要 $f$ 好求我们就能算了！

一般而言，我们可以统计一个前缀和（积），因为每当我们使用分块跳过一个区间的时候，其所对应的函数值也跳过了一个区间。所以此时，就需要乘上那一个区间的函数值（可以自己举个例子理解一下）。

但要是不一般（比如 $O(n)$ 的前缀和都会TLE），那……去死吧！（dalao可以用杜教筛试一试）

引理一

$$
\forall a,b,c \in \mathbb{Z},\lfloor \frac{a}{bc} \rfloor=\lfloor \frac{\lfloor \frac{a}{b} \rfloor}{c} \rfloor
$$

证明：
$$
\frac{a}{b}=\lfloor\frac{a}{b}\rfloor+r(0\le r<1)
\Rightarrow\lfloor\frac{a}{bc}\rfloor=\lfloor\frac{1}{c}(\lfloor\frac{a}{b}\rfloor+r)\rfloor=\lfloor\frac{\lfloor\frac{a}{b}\rfloor}{c}+\frac{r}{c}\rfloor=\lfloor \frac{\lfloor \frac{a}{b} \rfloor}{c} \rfloor
$$

引理二

$$
\forall n\in\mathbb{N_+},{\LARGE|}{\lfloor \frac{n}{d}\rfloor|d\in\mathbb{N_+}\wedge d\le n}{\LARGE|}\le\lfloor2\sqrt{n}\rfloor
$$

其中 $|V|$ 表示集合 $V$ 的元素个数

该引理是时间复杂度的保证

证明：

对于 $d\le\lfloor\sqrt{n}\rfloor$ ，由于 $d$ 最多有 $\lfloor\sqrt{n}\rfloor$ 种不同取值， $\lfloor\frac{n}{d}\rfloor$ 也只有 $\lfloor\sqrt{n}\rfloor$ 种取值

对于 $d>\lfloor\sqrt{n}\rfloor$ ，则 $\lfloor\frac{n}{d}\rfloor\le\lfloor\sqrt{n}\rfloor$ ，也只有 $\lfloor\sqrt{n}\rfloor$ 种取值

前置知识：积性函数

若函数 $f(x)$ 为积性函数，则有 $\forall gcd(a,b)=1,f(ab)=f(a)f(b)$

可以证明以下函数全是积性函数：

  • $\varphi(n)$ －欧拉函数
  • $\mu(n)$ －莫比乌斯函数
  • $gcd(n,k)$ －最大公因子，当k固定的情况
  • $d(n)$ －正因子数目
  • $\sigma(n)$ －所有正因子之和
  • $\sigma k(n)$ －因子函数， $n$ 的所有正因子的 $k$ 次幂之和，当中 $k$ 可为任何复数
  • $1(n)$ －不变的函数，定义为 $1(n)=1$ （完全积性）
  • $Id(n)$ －单位函数，定义为 $Id(n)=n$ （完全积性）
  • $Idk(n)$ －幂函数，对于任何复数、实数 $k$ ，定义为 $Idk(n)=n^k$ （完全积性）
  • $\varepsilon(n)$ －定义为：若 $n=1$ ， $\varepsilon(n)=1$ ；若 $n>1$ ， $\varepsilon(n)=0$ 。别称为“对于狄利克雷卷积的乘法单位”（完全积性）
  • $\lambda(n)$ －刘维尔函数，关于能整除 $n$ 的质因子的数目
  • $\gamma(n)$ ，定义为 $\gamma(n)=(-1)^{\omega(n)}$，在此加性函数 $\omega(n)$ 是不同能整除n的质数的数目
  • 另外，所有狄利克雷特征均是完全积性的

积性函数一般可以用线性筛筛得

前置知识：狄利克雷卷积

定义

狄利克雷（ $Dirichlet$ ）卷积定义为：

对于两个数论函数 $f,g$ ，有：
$$
(f*g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})
$$

性质

狄利克雷卷积满足以下运算律：

  • 交换律： $(f* g)=(g* f)$
  • 结合律： $((f* g)* h)=(f* (g* h))$
  • 分配律： $(f* (g+h))=((f* g)+(f* h))$
  • $(f*\varepsilon)=f$ ，其中 $\varepsilon$ 为狄利克雷卷积的单位元

for example：

  • $\varepsilon=(\mu*1)\Leftrightarrow\varepsilon(n)=\sum_{d|n}\mu(d)1(\frac{n}{d})=\sum_{d|n}\mu(d)$
  • $d=(1*1)\Leftrightarrow d(n)=\sum_{d|n}1(d)1(\frac{n}{d})=\sum_{d|n}1$
  • $Id=(\varphi*1)\Leftrightarrow Id(n)=n=\sum_{d|n}\varphi(n)$
  • $\varphi=(\mu*Id)\Leftrightarrow\varphi(n)=\sum_{d|n}d\mu(\frac{n}{d})$

求法

对于卷积的第 $n$ 项可以直接枚举约数，在 $O(\sqrt{n})$ 的时间内计算

但是对于狄利克雷卷积的前 $n$ 项
$$
h[i]=\sum_{d|i}f[d]g[i/d]
$$
如果一项一项算的话就需要 $O(n\sqrt{n})$ ，但实际上可以优化

设 $x=d,y=\frac{i}{d}$ 分别枚举 $x,y$ 对于 $h[xy]+=f[x]g[y]$ 即可。
时间复杂度 $O(n\log{n})$

1
2
3
for(int i=1;i<=n;++i) h[i]=0;
for(int i=1;i<=n;++i)
    for(int j=1;j*i<=n;++j) h[i*j]=(h[i*j]+f[i]*g[j]%P)%P;

对于 $k$ 次卷积，还可以快速幂优化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
struct HanShu{
    int f[N];
    void inint(){
    	for(int i=1;i<=n;++i) f[i]=0;
    }
    HanShu operator*(const HanShu _t){
	    HanShu h;
		h.inint();
	    for(int i=1;i<=n;++i)
	        for(int j=1;j*i<=n;++j) h.f[i*j]=(h.f[i*j]+f[i]*_t.f[j]%P)%P;
	    return h;
	}
	void operator=(const HanShu _t){
	    for(int i=1;i<=n;++i) f[i]=_t.f[i];
	}
};
HanShu H_qpow(HanShu x,int y){
    HanShu res;
    int _t=0;
    while(y){
    	if(y&1){
    		++_t;
    		if(_t==1) res=x;
    		else res=res*x;
    	}
    	x=x*x;
    	y>>=1;
    }
    return res;
}

前置知识：莫比乌斯函数

定义

设正整数 $N$ 分解质因数为 $N=p_1^{c_1}p_2^{c_2}…p_m^{c_m}$ ，定义函数
$$
\mu(N)=
\begin{cases}
0&\exists i\in[1,m],c_i>1\\
1&m\equiv0(mod\ 2),\forall i\in[1,m],c_i=1\\
-1&m\equiv1(mod\ 2),\forall i\in[1,m],c_i=1
\end{cases}
$$
我们称 $\mu$ 为莫比乌斯（ $M\ddot{o}bius$ ）函数。

性质和结论

  • 莫比乌斯函数是积性函数

  • $\sum_{d|n}\mu(d)=[n==1]$ ，其中 $[p]$ 表示 $p$ 为真时取 $1$ 其他时候取 $0$ ，该式也可表示为 $\sum_{d|n}\mu(d)=\varepsilon(n)$ ，也就是 $\mu*1=\varepsilon$

    证明：

    设 $n=\prod_{i=1}^{k}p_i^{c_i},n’=\prod_{i=1}^{k}p_i$

    那么 $\sum_{d|n}\mu(d)=\sum_{d|n’}\mu(d)=\sum_{i=0}^{k}C_k^i(-1)^i$

    又因为二项式定理： $(a+b)^k=\sum_{i=0}^{k}C_k^i a^i b^{k-i}$

    有： $\sum_{i=0}^{k}C_k^i(-1)^i=\sum_{i=0}^{k}C_k^i(-1)^i 1^{k-i}=((-1)+1)^k$

    当且仅当 $k=0$ ，即 $n=1$ 时原式为 $1$ ，其余时候都为 $0$

  • $[gcd(i,j)==1]=\sum_{d|gcd(i,j)}\mu(d)$ 该结论在反演中十分常用，正确性显然

筛法

若只求一项莫比乌斯函数，分解质因数即可；

若要求多项可以在线性筛素数时一并求出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int primes[N],mu[N],cnt=0;
bool v[N];
mu[1]=1; //1特判
for(int i=2;i<=n;++i){
    if(!v[i]) primes[++cnt]=i,mu[i]=-1; //质数
    for(int j=1;j<=cnt&&primes[j]*i<=n;++j){
		v[primes[j]*i]=true;
        if(i%primes[j]==0){ //若有一个质因数指数大于1
            mu[i*primes[j]]=0;
            break;
        }
        mu[i*primes[j]]=-mu[i];
    }
}

正题：莫比乌斯反演

定理

若 $F(n)$ 和 $f(n)$ 是定义在 $\mathbb{N}$ 上的函数，且有
$$
F(n)=\sum_{d|n}f(d)
$$
则：
$$
f(n)=\sum_{d|n}\mu(d)F(\frac{n}{d})
$$
我们称该定理为莫比乌斯反演定理（原来是用 $f$ 表示 $F$ ，现在反过来了）。

证明如下：
$$
\begin{align}
\sum_{d|n}\mu(d)F(\frac{n}{d})
&=\sum_{d|n}\mu(d)\sum_{i|\frac{n}{d}}f(i)\\
&=\sum_{i|n}f(i)\sum_{d|\frac{n}{i}}\mu(d)\\
&=\sum_{i|n}f(i)[\frac{n}{i}==1]\\
&=f(n)
\end{align}
$$
当然，不难发现，狄利克雷卷积也可以证明：

$$
\begin{align}
&F(n)=\sum_{d|n}f(d)\\
&\Rightarrow F=f* 1\\
&\Rightarrow F* \mu=f* 1* \mu\\
\text{又}\because&\mu* 1=\varepsilon\\
&\Rightarrow F* \mu=f* (1* \mu)=f\\
&\Rightarrow f=\mu* F
\end{align}
$$

然而，在莫比乌斯反演中更常用的是这个式子：

当 $F(n)$ 和 $f(n)$ 满足：
$$
F(n)=\sum_{n|d}f(d)
$$
则有：
$$
f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)
$$
其实就是枚举约数变成了枚举倍数，若用卷积证明，过程都一样，下面用定义证明：
$$
\begin{align}
\sum_{n|d}\mu(\frac{d}{n})F(d)
&=\sum_{n|d}\mu(\frac{d}{n})\sum_{d|i}f(i)\\
\text{不妨设}d’=\frac{d}{n}\\
\text{原式}
&=\sum_{n|i}f(i)\sum_{d’|\frac{i}{n}}\mu(d’)\\
&=\sum_{n|i}f(i)[\frac{i}{n}==1]\\
&=f(n)
\end{align}
$$
以上就是莫比乌斯反演的常用形式（要记住）

例题一

Problem b

分析：

若我们将数对 $(x,y)$ 对应到坐标系内，其实是要求在左下角为 $(a,c)$ 右上角为 $(b,d)$ 的矩形中有多少个点 $(x,y)$ 满足 $gcd(x,y)==k$

由矩形考虑二维前缀和，设 $S(x,y)$ 表示以 $(x,y)$ 为右上角， $(0,0)$ 为左下角的矩形中满足要求的点的个数，则 $Ans=S(b,d)-S(a-1,d)-S(b,c-1)+S(a-1,c-1)$

考虑用莫比乌斯反演，我们需要定义 $F,f$ 且让 $F(n)=\sum_{n|d}f(d)$ （即满足反演条件），要注意的是， $F$ 应该是比较好计算的函数（至少要比 $f$ 好计算，不然你反演用 $F$ 表示 $f$ 干嘛），还有就是 $f$ 要与答案有关系（不然求出 $f$ 干嘛）

那么，不妨令 $f(n)=\sum_{x=1}^{N}\sum_{y=1}^{M}[gcd(x,y)==n]$ 则我们要求的 $S(N,M)$ 就是 $f(k)$

又因为 $F(n)=\sum_{n|d}f(d)$ ，故定义 $F(n)=\sum_{x=1}^{N}\sum_{y=1}^{M}[n|gcd(x,y)]$

反演得： $f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)$

又因为由 $F$ 的定义， $n|gcd(x,y)\Leftrightarrow n|x\wedge n|y$ ，故 $F(d)=\lfloor\frac{N}{d}\rfloor\lfloor\frac{M}{d}\rfloor$

所以 $f(n)=\sum_{n|d}\mu(\frac{d}{n})\lfloor\frac{N}{d}\rfloor\lfloor\frac{M}{d}\rfloor$

枚举倍数比枚举因数要简单，所以我们令 $d’=\frac{d}{n},N’=\frac{N}{n},M’=\frac{M}{n}$ ，原式化为 $f(n)=\sum_{d’}\mu(d’)\lfloor\frac{N’}{d}\rfloor\lfloor\frac{M’}{d}\rfloor$ 预处理 $\mu$ 的前缀和后用数论分块可以求

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=50000+5;
int pri[N],cnt,mu[N];
bool v[N];
LL sum[N];
void pred(){
	mu[1]=1;
	for(int i=2;i<N;++i){
		if(!v[i]) pri[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&pri[j]*i<N;++j){
			v[pri[j]*i]=true;
			if(i%pri[j]==0){
				mu[i*pri[j]]=0;
				break;
			}
			mu[i*pri[j]]=-mu[i];
		}
	}
	sum[0]=0;
	for(int i=1;i<N;++i) sum[i]=sum[i-1]+mu[i];
}
int g(int k,int x){ //得到(k/x)中x所在块的右端点
	return k/(k/x);
}
LL S(int x,int y,int k){
	LL res=0;
	x=x/k,y=y/k;
	int _t=min(x,y);
	for(int l=1,r;l<=_t;l=r+1){
		r=min(_t,min(g(x,l),g(y,l))); //每次在两个取整中找小的,注意不能超过_t
		res+=(sum[r]-sum[l-1])*(x/l)*(y/l);
	}
	return res;
}
int main(){
	int n;
	int a,b,c,d,k;
	pred();
	scanf("%d",&n);
	while(n--){
		scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
		printf("%lld\n",S(b,d,k)-S(a-1,d,k)-S(b,c-1,k)+S(a-1,c-1,k));
	}
	return 0;
}

例题二

约数个数和

分析：

先看原题给我们的函数，由于 $d$ 是积性函数（忘了的上去翻），所以有 $d(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]$ （呃呃呃，鬼想的到啊！）

简单证明：

设 $i=p_1^{\alpha_1}p_2^{\alpha_2}…p_k^{\alpha_k},j=p_1^{\beta_1}p_2^{\beta_2}…p_k^{\beta_k}$ ，则 $ij=i=p_1^{\alpha_1+\beta_1}p_2^{\alpha_2+\beta_2}…p_k^{\alpha_k+\beta_k}$ ，由约数个数定理 $d(ij)=(\alpha_1+\beta_1+1)(\alpha_2+\beta_2+1)…(\alpha_k+\beta_k+1)$

再看等式右边， $x,y$ 的因数也一定是在 $p_1,p_2,…,p_k$ 中，对于质因子 $p_1$ ，若 $x$ 含 $p_1^t,(t\ne 0)$ ，则有意义的 $y$ 一定不含 $p_1$ ，此时 $t$ 有 $\alpha_1$ 种取值，同理，若 $y$ 含 $p_1^t,(t\ne 0)$ 有 $\beta_1$ 种取值，加上一种 $x,y$ 都不含 $p_1$ ，共 $\alpha_1+\beta_1+1$ 种有意义的取值，类似的，讨论其它质因数，乘法原理得左式等于右式

QED

我们发现 $Ans=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]$ ，类似与上题，令： $f(n)=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{x|i}\sum_{y|j}[gcd(x,y)==n],F(n)=\sum_{i=1}^{N}\sum_{j=1}^{M}\sum_{x|i}\sum_{y|j}[n|gcd(x,y)]$ 那么 $Ans=f(1)$

则 $F(n)=\sum_{n|d}f(d)$ ，反演得 $f(n)=\sum_{n|d}\mu(\frac{d}{n})F(d)$

故 $f(1)=\sum_{d=1}^{N}\mu(d)F(d)$ ，为了求出 $f(1)$ ，我们来看 $F(n)$ ，由于 $[n|gcd(x,y)]$ 与 $i,j$ 无关，我们可以枚举 $x,y$ 计算有多少 $i,j$ 符合要求， $F(n)=\sum_{x=1}^{N}\sum_{y=1}^{M}\lfloor\frac{N}{x}\rfloor\lfloor\frac{M}{y}\rfloor[n|gcd(x,y)]$

同上题一样，可设 $x’=\frac{x}{n},y’=\frac{y}{n},N’=\frac{N}{n},M’=\frac{M}{n}$ ，得 $F(n)=\sum_{x’=1}^{N’}\sum_{y’=1}^{M’}\lfloor\frac{N’}{x’}\rfloor\lfloor\frac{M’}{y’}\rfloor$

调换循环顺序得 $F(n)=(\sum_{x’=1}^{N’}\lfloor\frac{N’}{x’}\rfloor)(\sum_{y’=1}^{M’}\lfloor\frac{M’}{y’}\rfloor)$

令 $h(k)=\sum_{i=1}^{k}\lfloor\frac{k}{i}\rfloor$ ，有 $F(n)=h(N’)h(M’)$ ，故 $f(1)=\sum_{d=1}^{N}\mu(d)h(\frac{N}{d})h(\frac{M}{d})$ ，预处理所有的 $h(x)$ ，数论分块解决

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=50000+5;
int pri[N],cnt,mu[N];
LL sum[N],h[N];
bool v[N];
int g(int k,int x){
	return k/(k/x);
}
void pred(){
	mu[1]=1;
	for(int i=2;i<N;++i){
		if(!v[i]) pri[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&pri[j]*i<N;++j){
			v[pri[j]*i]=true;
			if(i%pri[j]==0){
				mu[i*pri[j]]=0;
				break;
			}
			mu[i*pri[j]]=-mu[i];
		}
	}
	sum[0]=0;
	for(int i=1;i<N;++i) sum[i]=sum[i-1]+mu[i];
	for(int i=1;i<N;++i)
		for(int l=1,r;l<=i;l=r+1){
			r=min(i,g(i,l));
			h[i]+=(r-l+1)*(i/l);
		}
}
int main(){
	pred();
	int T;
	scanf("%d",&T);
	while(T--){
		int n,m;
		scanf("%d%d",&n,&m);
		LL ans=0;
		int _t=min(n,m);
		for(int l=1,r;l<=_t;l=r+1){
			r=min(_t,min(g(n,l),g(m,l)));
			ans+=(sum[r]-sum[l-1])*h[n/l]*h[m/l];
		}
		printf("%lld\n",ans);
	}
	return 0;
}

例题三

龙哥的问题

本题比较简单，甚至不用莫反

考虑答案：
$$
\begin{align}
Ans
&=\sum_{i=1}^{n}gcd(i,n)\\
&=\sum_{d|n}(d\times\sum_{i=1}^{n}[gcd(i,n)==d])\\
&=\sum_{d|n}(d\times\sum_{i=1}^{n}[gcd(\frac{i}{d},\frac{n}{d})==1])\\
&=\sum_{d|n}d\times\varphi(\frac{n}{d})
\end{align}
$$
对于每一个 $n$ 枚举约数（ $int$ 范围内最多的一个数约数只有 $1600$ 多个） ，再 $\sqrt{n}$ 求 $\varphi$ ，本题的范围可以通过

那可否继续优化呢？令 $k=\frac{n}{d}$ ，再继续化一下：
$$
\begin{align}
Ans
&=\sum_{d|n}d\times\varphi(\frac{n}{d})\\
&=\sum_{k|n}\frac{n}{k}\times\varphi(k)\\
&=\sum_{k|n}\frac{n}{k}\times k\prod(1-\frac{1}{p_i})\\
&=n\sum_{k|n}\prod(1-\frac{1}{p_i})
\end{align}
$$
其中 $k_t=\prod p_i^{c_{i,t}}$ ，又因为 $k|n$ ，不妨设 $n=\prod p_i^{\alpha_i}$ ，下面的变化比较常用，要记住：
$$
\begin{align}
&k|d\\
\Rightarrow&\forall k_t,0\le c_{i,t}\le\alpha_i\\
\Rightarrow&\sum_{k|n}k=\sum_{t}\prod_{i} p_i^{c_{i,t}}=\prod_{t}\sum_{i} p_i^{c_{i,t}}\\
&=(p_1^0+p_1^2+…+p_1^{\alpha_1})(p_2^0+p_2^2+…+p_2^{\alpha_2})…
\end{align}
$$

如果不理解的话可以这样想：每个括号里取一项作 $p_i^{c_{i,t}}$ 乘起来刚好对应一个 $k_t$

那么考虑上式对答案的贡献，我们发现对于所有的 $k$ ，若取第 $i$ 个质因子的指数为 $0$ （即取 $p_i^0$ ）对答案贡献为 $1$ ，其他情况（即 $p_i^c,c\ne0$ ）对答案的贡献都是 $1-\frac{1}{p_i}$ ，故答案可化为：
$$
\begin{align}
Ans
&=n\prod_{i}(1+(1-\frac{1}{p_i})+(1-\frac{1}{p_i})+…+(1-\frac{1}{p_i}))\\
&=n\prod_{i}(1+a(1-\frac{1}{p_i}))
\end{align}
$$
分解 $n$ 的质因子，可以在 $O(\sqrt{n})$ 时间求解，但注意，本题对 $long\ long$ 卡的非常死，甚至在计算 $ans$ 时先乘再除都会炸 $long\ long$

代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include<bits/stdc++.h>
#define LL long long
using namespace std;

int main(){
	LL n;
	scanf("%lld",&n);
	LL ans=n; //最初的n 
	for(LL i=2;i*i<=n;++i)
		if(n%i==0){
			LL a=0;
			while(n%i==0) ++a,n/=i;
			ans/=i,ans*=i+a*i-a; //一定要先除 
		}
	if(n>1) ans/=n,ans*=n+n-1;
	printf("%lld\n",ans);
	return 0;
}

例题四

来看一道难题：

数字表格

这道题要点很杂很多，我们慢慢看：

首先明确
$$
Ans=\prod_{i=1}^{n}\prod_{j=1}^{m}fb(gcd(i,j))
$$
其中 $fb(x)$ 表示Fibonacci数列第 $x$ 项

由于 $\prod$ 不好处理，考虑枚举 $gcd$ ：
$$
Ans=\prod_{d=1}^{N}[fb(d)]^{f(d)}
$$
其中 $N=\min(n,m)$ ， $f(d)=\sum_{i=1}^{m}\sum_{j=1}^{n}[gcd(i,j)==d]$

发现这个 $f$ 的定义我们非常熟悉，于是一通莫反：
$$
\begin{align}
\text{令}&F(d)=\sum_{i=1}^{m}\sum_{j=1}^{n}[d|gcd(i,j)]\\
\text{则有}&F(d)=\sum_{r|d}f(r)\\
\text{反演得}&f(d)=\sum_{d|r}\mu(\frac{r}{d})F(r)\\
\text{又有}&d|gcd(i,j)\Leftrightarrow d|i,d|j\\
\text{故}&F(d)=\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor\\
\text{带入得}&f(d)=\sum_{d|r}\mu(\frac{r}{d})\lfloor\frac{n}{r}\rfloor\lfloor\frac{m}{r}\rfloor\\
\text{令}w=\frac{r}{d},\text{有}&f(d)=\sum_{w=1}^{\lfloor\frac{N}{d}\rfloor}\mu(w)\lfloor\frac{n}{wd}\rfloor\lfloor\frac{m}{wd}\rfloor
\end{align}
$$
此时我们发现在莫反后， $f(d)$ 可以在 $O(\sqrt{n})$ 内计算，但是，我们观察 $Ans$ ，对于每个 $d$ ， $f(d)$ 都必须重新计算（因为 $f$ 与 $n,m,d$ 都有关），这样的时间负责度为 $O(Tn(\sqrt{n}+\log{n}))$ 直接爆炸好吧

必须考虑优化，但 $fb$ 和 $\prod$ 都太不好搞了， $f$ 也不好在变形（我在这里卡了好久），既然两个都不好搞，我们不妨把它们一起来看看：
$$
Ans=\prod_{d=1}^{N}[fb(d)]^{\sum_{w=1}^{\lfloor\frac{N}{d}\rfloor}\mu(w)\lfloor\frac{n}{wd}\rfloor\lfloor\frac{m}{wd}\rfloor}
$$
我们发现一个关键点：每次询问只修改 $n,m$ 而 $d,w$ 对于每次询问来说是一样的，这启发我们把几个与 $n,m$ 无关的项拿出来预处理：

$$
\begin{align}
Ans=\prod_{d=1}^{N}[fb(d)]^{\sum_{w=1}^{\lfloor\frac{N}{d}\rfloor}\mu(w)\lfloor\frac{n}{wd}\rfloor\lfloor\frac{m}{wd}\rfloor}\\
=\prod_{d=1}^{N}{[fb(d)]^{\sum_{w=1}^{\lfloor\frac{N}{d}\rfloor}\mu(w)}}^{\lfloor\frac{n}{wd}\rfloor\lfloor\frac{m}{wd}\rfloor}\\
\text{令}p=wd,\text{有}Ans
=\prod_{p=1}^{N}{[\prod_{d|p}fb(d)]^{\mu(\frac{p}{d})}}^{\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor}\\
\text{令}h(p)=\prod_{d|p}fb(d)^{\mu(\frac{p}{d})},\text{有}Ans
=\prod_{p=1}^{N}h(p)^{\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor}
\end{align}
$$

发现指数位置的 $\lfloor\frac{n}{p}\rfloor\lfloor\frac{m}{p}\rfloor$ 可以数论分块，而 $h$ 可以提前预处理出前缀积，如此，回答询问的总时间变为 $O(T\sqrt{n}\log{n})$

下面考虑如何计算 $h$ ，由于因数不好枚举，我们改为枚举倍数：
$$
\begin{align}
h(p)
&=\prod_{d|p}fb(d)^{\mu(\frac{p}{d})}\\
&=\prod_{i=1}^{N}fb(\frac{p}{i})^{\mu(i)}\\
\text{变形得}h(ij)
&=\prod_{i=1}^{N}\prod_{j=1}^{\lfloor\frac{N}{i}\rfloor}fb(j)^{\mu(i)}
\end{align}
$$
预处理时间为 $O(n\log{n})$

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=1000000+5,P=1000000000+7;
LL pri[N],cnt,mu[N],fb[N],h[N],inv_fb[N],mul[N];
bool v[N];
int g(int k,int x){
	return k/(k/x);
}
LL qpow(LL x,LL y){
	LL res=1;
	while(y){
		if(y&1) res=(res*x)%P;
		x=(x*x)%P;
		y>>=1;
	}
	return res;
}
void pred(){
	mu[1]=1;
	for(int i=2;i<N;++i){
		if(!v[i]) pri[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&pri[j]*i<N;++j){
			v[pri[j]*i]=true;
			if(i%pri[j]==0){
				mu[i*pri[j]]=0;
				break;
			}
			mu[i*pri[j]]=-mu[i];
		}
	}
	fb[0]=0;
	fb[1]=1;
	inv_fb[1]=1; //特殊处理
	for(int i=2;i<N;++i) fb[i]=(fb[i-1]+fb[i-2])%P,inv_fb[i]=qpow(fb[i],P-2); //求逆元
	for(int i=1;i<N;++i) h[i]=1;
    for(int i=1;i<N;++i){
        if(!mu[i])continue;
        for(int j=1;j*i<N;++j)
        	h[j*i]=h[j*i]*(mu[i]==1?fb[j]:inv_fb[j])%P; //注意mu可能为负一
    }
	mul[0]=1;
	for(int i=1;i<N;++i) mul[i]=(mul[i-1]*h[i])%P;
}
LL get_ans(int x,int y){
	LL res=1,as;
	int _t=min(x,y);
	for(int l=1,r;l<=_t;l=r+1){
		r=min(_t,min(g(x,l),g(y,l)));
		as=mul[r]*qpow(mul[l-1],P-2)%P; //前缀积得底数
		res=res*qpow(as,1ll*(x/l)*(y/l)%(P-1))%P;
	}
	return res;
}
int main(){
	pred();
	int T;
	scanf("%d",&T);
	int n,m;
	while(T--){
		scanf("%d%d",&n,&m);
		printf("%lld\n",get_ans(n,m));
	}
	return 0;
}

补充：变换方法

下面记录一些常见的变换方法：

  • $$
    gcd(i,j)==1\Leftrightarrow\sum_{d|i\wedge d|j}\mu(d)
    $$

  • 像例题四一样的完全展开预处理无关项

  • 像例题二一样 $d(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]$

  • 改变未知数使得两个 $\sum$ 分离，预处理其中一个（也可以看作将两数积设出，继续化简）

    如YY的GCD，莫反后得
    $$
    \sum_{k\in prime}\sum_{d=1}^{\lfloor\frac{min(N,M)}{k}\rfloor}\mu(d)\lfloor\frac{N}{kd}\rfloor\lfloor\frac{M}{kd}\rfloor
    $$
    若直接计算，对于每个询问，最坏时间复杂度为 $O(n)$ 会TLE，可以变化如下：

    $$
    \begin{align}
    \text{设}T
    &=kd\\
    \text{则原式}
    &=\sum_{k\in primes}\sum_{d=1}^{\lfloor\frac{min(N,M)}{k}\rfloor}\mu(d)\lfloor\frac{N}{T}\rfloor\lfloor\frac{M}{T}\rfloor\\
    &=\sum_{T=1}^{min(N,M)}\sum_{k|T \wedge k\in prime}\mu(\frac{T}{k})\lfloor\frac{N}{T}\rfloor\lfloor\frac{M}{T}\rfloor\\
    &=\sum_{T=1}^{min(N,M)}\lfloor\frac{N}{T}\rfloor\lfloor\frac{M}{T}\rfloor\sum_{k|T\wedge k\in prime}\mu(\frac{T}{k})\\
    \end{align}
    $$

    后面的 $\sum$ 与询问无关，可以预处理

  • 注意卷积中的一些常用等式，如 $Id*\mu=\varphi$ ：

    $$
    \begin{align}
    Ans
    &=\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)\\
    \text{莫反得}Ans
    &=\sum_{r=1}^{n}\lfloor\frac{n}{r}\rfloor\lfloor\frac{n}{r}\rfloor\sum_{d|r}d*\mu(\frac{r}{d})\\
    &=\sum_{r=1}^{n}\lfloor\frac{n}{r}\rfloor\lfloor\frac{n}{r}\rfloor\varphi(r)
    \end{align}
    $$