数学(2)——莫比乌斯反演做题记录


Luogu P1447 [NOI2010] 能量采集

\[\sum_{i=1}^n\sum_{j=1}^m(2\times\gcd{(i,j)-1}) \]

可以转换为

\[-(n\times m)+2\times{\color{blue}{\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)}} \]

考虑对蓝色部分莫比乌斯反演

\[\begin{aligned} f(\mathrm d)&=\sum_{i=1}^n\sum_{j=1}^m[\mathrm d=\gcd(i,j)]\\ g(\mathrm d)&=\sum_{i=1}^n\sum_{j=1}^m[\mathrm d|\gcd(i,j)]\\ g(\mathrm d)&=\lfloor\frac{n}{\mathrm d}\rfloor\lfloor\frac{m}{\mathrm d}\rfloor \end{aligned} \]

这个 \(g(\mathrm d)\) 可以帮我们做数论分块。

首先可以看出前面的部分

\[\texttt{Part}=\sum_{\mathrm d=1}^nf(\mathrm d)\times \mathrm d \]

我们可以直接枚举 \(i\) 而不是 \(i|d\)

所以有

\[f(\mathrm d)=\sum_{i=1}^{\lfloor\frac{n}{\mathrm d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{\mathrm d}\rfloor}\sum_{k|\gcd(i,j)}\mu(k)\\ f(\mathrm d)=\sum_{k=1}^{\lfloor\frac{n}{\mathrm d}\rfloor}\mu(k)\times\lfloor\frac{n}{k\times \mathrm d}\rfloor\lfloor\frac{m}{k\times \mathrm d}\rfloor \]

答案就是

\[\texttt{ans}=[2\sum_{\mathrm d=1}^n\mathrm d\times\sum_{k=1}^{\lfloor\frac{n}{\mathrm d}\rfloor}(\mu(k)\lfloor\frac{n}{k\times \mathrm d}\rfloor\lfloor\frac{m}{k\times \mathrm d}\rfloor)]-n\times m \]

代码通过数论分块和线性筛实现可以就通过这题。并且首先保证 \(\color{red}n\leq m\)

#include
#include
#include
#include
#include

#define ll long long

using namespace std;

inline int read()
{
	int x=0,w=1;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	if(ch=='-')ch=getchar(),w=-1;
	while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();
	return x*w;
}

int n,m,tot;
const int N=100005;
ll phi[N],pri[N],sum[N];
bool vis[N];

void Phi()
{
	phi[1]=1;
	for(int i=2;i<=100000;i++)
	{
		if(!vis[i])pri[++tot]=i,phi[i]=i-1;
		for(int j=1;j<=tot&&i*pri[j]<=100000;j++)
		{
			if(i*pri[j]>100000)continue;
			vis[i*pri[j]]=1;
			if(i%pri[j])phi[i*pri[j]]=phi[i]*phi[pri[j]];
			else {phi[i*pri[j]]=phi[i]*pri[j];continue;}
		}
	}
}

int main()
{
	n=read();m=read();
	Phi();
	for(int i=1;i<=100000;i++)sum[i]=sum[i-1]+phi[i];
	ll ans=0;
	ll l=1,r=0;
	while(l<=min(n,m))
	{
		r=min(n/(n/l),m/(m/l));
		ans+=(ll)(sum[r]-sum[l-1])*(n/l)*(m/l);
		l=r+1; 
	}
	printf("%lld\n",2ll*ans-(ll)n*m);
	return 0;
}

Luogu P2257 YY的GCD

这个题有个神必的巧妙做法,参考这篇题解。

这个题是

\[\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)\in \texttt{Prime}] \]

(下文质数集合用 \(\mathbb P\) 表示且首先同样保证 \(\color{red}n\leq m\))

我们考虑直接枚举质数:

\[\sum_{p\in\mathbb P}\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=p] \]

我们考虑莫比乌斯反演技巧,设:

\[f(\mathrm d)=\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=\mathrm d]\\ g(x)=\sum_{x|\mathrm d}f(\mathrm d)=\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{x}\rfloor\\ f(x)=\sum_{x|\mathrm d}\mu(\lfloor\frac{\mathrm d}{x}\rfloor)g(\mathrm d) \]

那我们发现答案就是

\[\texttt{ans}=\sum_{p\in\mathbb P}f(p) \]

再对答案进行莫比乌斯反演

\[\begin{aligned} \texttt{ans}&=\sum_{p\in\mathbb P}\sum_{p|\mathrm d}\mu(\lfloor\frac{\mathrm d}{p}\rfloor)g(\mathrm d)\\ &=\sum_{p\in\mathbb P}\sum_{\mathrm d=1}^{\lfloor\frac{n}{p}\rfloor}\mu(\mathrm d)g(p\times \mathrm d)\\ &=\sum_{p\in\mathbb P}\sum_{\mathrm d=1}^{\lfloor\frac{n}{p}\rfloor}\mu(\mathrm d)\lfloor\frac{n}{p\times \mathrm d}\rfloor\lfloor\frac{m}{p\times \mathrm d}\rfloor\\ \end{aligned} \]

消去 \(p\times \mathrm d\) 可以考虑换元为 \(k\)

\[\begin{aligned} \texttt{ans}&=\sum_{k=1}^n\sum_{p|k,p\in\mathbb P}\mu(\lfloor\frac{k}{p}\rfloor)\lfloor\frac{n}{k}\rfloor\lfloor\frac{m}{k}\rfloor\\ &=\sum_{k=1}^n\lfloor\frac{n}{k}\rfloor\lfloor\frac{m}{k}\rfloor(\sum_{p|k}\mu(\lfloor\frac{k}{p}\rfloor)) \end{aligned} \]

这样数论分块,筛法,前缀和就可以实现了。

#include
#include
#include
#include
#include

#define ll long long

using namespace std;

inline int read()
{
	int x=0;int w=1;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	if(ch=='-')ch=getchar(),w=-1;
	while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();
	return x*w;
}

int n,m;
const int N=10000000;
const int M=10000005;

ll pri[M/10],phi[M],s[M],g[M];
bool vis[M];
ll mn;
ll ans,cnt;

void Phi()
{
	phi[1]=1;
	for(int i=2;i<=mn;i++)
	{
		if(!vis[i])phi[i]=-1,pri[++cnt]=i;
		for(int j=1;j<=cnt&&i*pri[j]<=mn;j++)
		{
			vis[i*pri[j]]=1;
			if(i%pri[j]==0)break;
			else phi[i*pri[j]]=-phi[i];
		}
	}
	for(int i=1;i<=cnt;i++)
	{
		for(int j=1;j*pri[i]<=mn;j++)
		g[j*pri[i]]+=phi[j];
	}
	for(int i=1;i<=mn;i++)
	{
		s[i]=s[i-1]+g[i];
	}
}

int main()
{
	int T;
	T=read();
	mn=10000000;
	Phi();ll ans=0;ll r;
	while(T--)
	{
		n=read();m=read();
		mn=(ll)min(n,m);ans=0;
		for(int i=1;i<=mn;i=r+1)
		{
			r=min(n/(n/i),m/(m/i));
			ans+=(ll)(s[r]-s[i-1])*(ll)(n/i)*(ll)(m/i);
		}
		printf("%lld\n",ans);
	}
	return 0;
}

Luogu P2261 [CQOI2007]余数求和

这个题好像很水,因为好像只是一个入门的除法分块练习题。

\[\because k \bmod i=k-i\times\lfloor\frac{k}{i}\rfloor\\ \therefore \sum_{i=1}^nk\bmod i=\sum_{i=1}^n(k-i\times\lfloor\frac{k}{i}\rfloor)=n\times k-\sum_{i=1}^n(i\times\lfloor\frac{k}{i}\rfloor) \]

这样就可以做除法分块

就是对 \(\displaystyle {\color{blue}{\sum_{i=1}^n(i\times\lfloor\frac{k}{i}\rfloor)}}\) 做除法分块。

代码如下

for(ll i=1;i<=n;i=r+1)
{
	if(k/i)
		r=min((ll)n,k/(k/i));
	else 
		r=n;
	sum+=(ll)(k/i)*(ll)(r-i+1)*(i+r)/2ll;
}

Luogu P1390 公约数的和

这个题是求

\[\sum_{i=1}^n\sum_{j=i+1}^n\gcd(i,j) \]

我们先考虑求出

\[\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j) \]

采用一些神必技巧,枚举一个数 \(k\) ,然后改成对 \(k\) 求和得到

\[\sum_{k=1}^n(k\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=k])\\\rightarrow\sum_{k=1}^n(k\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}[\gcd(i,j)=1])\\ \rightarrow\sum_{k=1}^n[k\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}\epsilon(\gcd(i,j))] \]

我们改为 \(\epsilon\) $\epsilon$ 是为了方便我们做 \(\texttt{Dirichlet}\) 卷积。

所以我们可以转化为

\[\sum_{k=1}^nk\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{\mathrm d|\gcd(i,j)}\mu(\mathrm d)\\ \rightarrow \sum_{k=1}^nk\sum_{\mathrm d=1}\mu(\mathrm d)\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}[\mathrm d|i]\sum_{j=1}^{\lfloor \frac{n}{k}\rfloor}[\mathrm d|j]=\sum_{k=1}^nk\sum_{\mathrm d=1}\mu(\mathrm d)\lfloor\frac{n}{k\times\mathrm d}\rfloor\lfloor\frac{n}{k\times \mathrm d}\rfloor \]

不过别忘了,这个式子只是 \(\displaystyle \sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)\)

我们推导得

\[\begin{aligned} \texttt{ans}&=\sum_{i=1}^n\sum_{j=i+1}^{n}\gcd(i,j)\\ &=\sum_{i=1}^n\sum_{j=1}^{i-1}\gcd(i,j)\\ &=\frac{\displaystyle\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)-\sum_{i=1}^n\gcd(i,i)}{2}\\ &=\frac{\displaystyle{ \color{blue}{\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)}}-\frac{(1+n)\times n}{2}}{2}\\ \end{aligned} \]

用之前那个式子对蓝色的 \(\displaystyle{\color{blue}{\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)}}\) 部分做数论分块求解。

代码:

#include
#include
#include
#include
#include
#include
#include
#include
//#include

#define ll long long
#define ull unsigned long long
#define INL inline
//Tosaka Rin Suki~
using namespace std;

const int N=2000005;

INL int read()
{
	 int x=0;int w=1;
	 char ch=getchar();
	 while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	 if(ch=='-')w=-1,ch=getchar();
	 while(ch>='0'&&ch<='9')
	 {x=(x<<1)+(x<<3)+ch-48,ch=getchar();}
	 return x*w;
}

INL int mx(int a,int b){return a>b?a:b;}
INL int mn(int a,int b){return a>1ll;
	ans>>=1;
	printf("%lld\n",ans);
	return 0;
}

LibreOJ #6539. 奇妙数论题

大家先回顾一下上一题,这题的部分分就是上一题。


题如其名,一道巧妙莫反题。

给定长为 \(n\)排列 \(a\)

\[\sum_{i=1}^n\sum_{j=1}^n\gcd(a_i,a_j)\times \gcd(i,j) \]

我们有部分分是 \(a_i=i\) 的,所以一直到 \(70\) 分的部分分就是求:

\[\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)^2 \]

就上一题稍微改一下就有了:

\[\sum_{k=1}^nk^2\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\lfloor\frac{n}{k\times d}\rfloor\lfloor\frac{n}{k\times d}\rfloor \]

正解的瓶颈其实就是排列变成无序后无从下手。

所以我们先将式子转换成:

\[\begin{aligned} \texttt{ans}&=\sum_{i=1}^n\sum_{j=1}^n\gcd(a_i,a_j)\times \gcd(i,j)\\ &=\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d]\times \gcd(a_i,a_j)\\ &=\sum_{d=1}^nd\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\sum_{i=1}^{\lfloor\frac{n}{d\times x}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d\times x}\rfloor}\gcd(a_{idx},a_{jdx}) \end{aligned} \]

这里就是根据上一题的套路,\(\gcd(i,j)\) 看成 \(\gcd(i,j)\times 1\) 于是可以在数论分块的部分直接带入对应下标 \(a_i\)\(\gcd\)

但是求 \(\gcd(a_i,a_j)\) 也就是 \(\gcd(a_{idx},a_{jdx})\) 的复杂度降不下来,考虑怎么做。

设:

\[f(x)=\sum_{i=1}^{\lfloor\frac{n}{x}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{x}\rfloor}\gcd(a_{ix},a_{jx}) \]

我们观察这个函数,其实就是要求在 \(a\) 这个排列中下标分别为 \(x\) 的倍数的元素彼此之间 \(\gcd\) 的和。

我们设 \(\mathbb S_x\)大小为 \(x\) 的倍数的下标的集合

于是有:

\[f(x)=\sum_{i\in\mathbb S_x}\sum_{j\in\mathbb S_x}\gcd(a_i,a_j) \]

我们进一步抽象,我们将 \(i\in\mathbb S_x\) 映射到每一个 \(a_i\)\(\forall i\in \mathbb S_x,\)\(a_i\) 的集合为 \(\mathbf S_x\)

现在有:

\[f(x)=\sum_{i\in\mathbf S_x}\sum_{j\in\mathbf S_x}\gcd(i,j)\Rightarrow \sum_d d\sum_{i\in\mathbf S_x}\sum_{j\in \mathbf S_x}[\gcd(i,j)=d] \]

\(i,j\) 的集合与范围发生了变化,显然我们不能再将 \([\gcd(i,j)=d]\) 化为 \([\gcd(i,j)=1]\) 去做。

我们有另一种方法求 \([x=d]\) 这一类的式子。

我们观察 \([x=d]\) 长得很像 \(\epsilon(x)=[x=1]\)

我们有 \(\mu*1=\epsilon\) 所以有 \(\displaystyle \epsilon(x)=\sum_{d|x}\mu(d)\)

那么我们的 \([x=d]\) 可以化为 \(\displaystyle[\frac{T}{d}=1\space\mathrm{and}\space T|d]\)

类似的我们展开得 \(\displaystyle [x=d]=\sum_{d|T,T|x}\mu(\frac{T}{d})\)

得到(这个地方看官方题解看了很久才懂的,感觉很 tricky):

\[\begin{aligned} f(x)&=\sum_d d\sum_{i\in\mathbf S_x}\sum_{j\in \mathbf S_x}\sum_{d|T,T|\gcd(i,j)}\mu(\frac{T}{d})\\ &=\sum_d d\sum_{d|T}\mu(\frac{T}{d})(\sum_{i\in \mathbf S_x}[T|i])^2\\ &=\sum_T({\color{blue}{\sum_{d|T}d\times\mu(\frac{T}{d})}})(\sum_{i\in \mathbf S_x}[T|i])^2 \end{aligned} \]

我们知道有 \(\mu * \mathrm{id}=\varphi\)

观察上面蓝色部分的 \(\displaystyle{\color{blue}{\sum_{d|T}d\times\mu(\frac{T}{d})}}\) 就是 \(\varphi(T)\)

最后得到:

\[f(x)=\sum_T \varphi(T)(\sum_{i\in \mathbf S_x}[T|i])^2 \]

我们回到 \(\texttt{ans}\) 的部分。

有:

\[\texttt{ans}=\sum_{d=1}^nd\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)f(dx) \]

具体实现方法:

  • 首先我们筛出得到 \(\forall x\in[1,n]\) 每个 \(\mathbf S_x\)
  • 依次统计对于每个 \(\mathbf S_x\) 每个元素的的因子的贡献。
  • 得到 \(f(x)\) 并筛出 \(\mu(x)\) 后直接求和即可。

复杂度分析的话,看到 loj 官方题解的评论区中有这样的结论:

我们设对于 \(x\)\(x\) 的因子个数为 \(d(x)\)

As Vaclav Kotesovec said Aug 30 2018, \(\displaystyle \sum_{k\leq n}d^2(k)\)is asymptotic to \(\Theta(\frac{1}{\pi^2}n\log^3 n+n\log^2 n)\).

代码(其实就是出题人的写法,好像多数采用了埃筛):

#include

#define ll long long
#define ull unsigned long long
#define INL inline
#define Re register

//Tosaka Rin Suki~

INL int read()
{
	int x=0,w=1;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();if(ch=='-')w=-1,ch=getchar();
	while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();return x*w;
}

const int N=1e5+5,MOD=1e9+7;

std::vector  factors[N];
std::vector  pack;
int mu[N],g[N],f[N],d[N],n,a[N],ans;


int main()
{
	//freopen(".in","r",stdin);
	//freopen(".out","w",stdout);
	n=read();for(int i=1;i<=n;i++)a[i]=read();
	mu[1]=1;
	for(int i=1;i<=n;i++)
	{
		factors[i].push_back(i);
		for(int j=i*2;j<=n;j+=i)
			mu[j]-=mu[i],factors[j].push_back(i);
	}
	for(int i=1;i<=n;i++)
		for(int j=i;j<=n;j+=i)
			g[j]=(g[j]+mu[j/i]*i)%MOD;
	for(int i=1;i<=n;i++)
	{
		for(int j=i;j<=n;j+=i)
			for(int k=0;k

Luogu P3704 [SDOI2017]数字表格

在机房复盘 SDOI 2017 的模拟赛第一次场 A 黑题写题解纪念一下。

这个题要你求:

\[\prod_{i=1}^n\prod_{j=1}^m \mathrm{fib_{\gcd(i,j)}} \]

其中 \(\mathrm{fib_{i}}\) 表示斐波拉契数列第 \(i\) 项。

\(\prod\) 连乘运算做不了 \(\sum\) 的一些卷积公式(尤其是 \(\texttt{Dirichlet}\) 卷积)。我们考虑把这个式子换成可以用 \(\sum\) 推的情况。

首先一般性地,我们保证 \({\color{red}{n\leq m}}\)

我们枚举 \(\displaystyle\sum_i\sum_j f(i,j)\) 这一类式子的时候会习惯性地枚举一个 \(k\) 让其满足和式内部的条件。

这题即使是 \(\prod\) 也可以套路化地使用这个方法。

所以原式可以转换为

\[\prod_{k=1}^n\prod_{i=1}^n\prod_{j=1}^m \mathrm{fib}_k\times[\gcd(i,j)=k] \]

我们发现每个 \(\mathrm{fib}_d\) 都会与自己相乘很多次,于是我们可以将每个 \(d\) 的贡献考虑写成幂的形式。

显然每个 \(\mathrm{fib}_d\) 的幂的指数就是其等于 \(\gcd(i,j)\) 对应的数对 \((i,j)\) 的个数。

现在我们有

\[\prod_{k=1}^n\mathrm {fib}_k^{\displaystyle\bigg({\color{blue}{\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k]}}\bigg)} \]

我们把指数特地标记出来。因为我们发现他是一个老朋友。

我们在做形似 \(\displaystyle\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)\) 的一类莫反题是会将式子转换成 \(\displaystyle \sum_{k=1}^n\bigg(k{\color{blue}{\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=k]}}\bigg)\)

我们考虑怎么将这个式子转化成可做的形式(以下皆是莫反题老套路,跟不上请先去做一些基础莫反题,不是幽灵乐团,别信 CYJian)。

\[\begin{aligned} {\color{blue}{\texttt{blue part}}}&=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[\gcd(i,j)=1]\\ &=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\epsilon(\gcd(i,j))\\ &=\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d|\gcd(i,j)}\mu(d)\\ &=\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}[ d|i]\sum_{j=1}^{\lfloor \frac{m}{k}\rfloor}[d|j]\\ &=\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\lfloor\frac{n}{k\times d}\rfloor\lfloor\frac{m}{k\times d}\rfloor \end{aligned} \]

这个东西已经可以 \(\Theta(\sqrt n)\) 级别的数论分块做了,但是我们要把其代入原式继续看。

现在原式变成了:

\[\prod_{k=1}^n\mathrm {fib}_k^{\bigg (\displaystyle\sum_{d=1}^{\lfloor\frac{n}{k}\rfloor}\mu(d)\lfloor\frac{n}{k\times d}\rfloor\lfloor\frac{m}{k\times d}\rfloor\bigg)} \]

我们设 \(T=k\times d\)

我们发现这样方便让我们将 \(\displaystyle\sum_{d=1}^{\lfloor \frac{n}{k}\rfloor}\)\(\mathrm{fib}_k\) 以幂的形式结合。所有变量用 \(T\) 表示出来后贡献到一个 \(\prod\) 中。

写出来就是:

\[\prod_{T=1}^n\bigg({\color{blue}{\prod_{k|T}\mathrm{fib}_k^{\mu(\frac{T}{k})}}}\bigg)^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor} \]

蓝色部分我们可以通过类似筛法预处理的方法求出:

  • 线筛 \(\mu\) 函数。
  • 对于每一个 \(k\) 枚举其倍数 \(T\) ,乘上这个 \(\mathrm{fib}_k\) 算入这个 \(T\) 的贡献。
  • 因为 \(\mu\) 只有三种取值,有影响的只有 \(\{1,-1\}\) 两种。
  • 所以问题就变成了乘上 \(\mathrm {fib}_{k}\) 或除以 \(\mathrm {fib}_k\)。这里就变成 \(\mathrm{inv fib}\)(逆元)和 \(\mathrm{fib}\) 两种情况。

外层我们做数论分块的时候,每次求内部只需要一次 \(\Theta (\log p)\)\(p\) 是模数)和 \(\Theta(\log(\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor))\)的快速幂分别求逆元和与外层指数的幂。

预处理复杂度是 \(\Theta(n\log n+n\log p)\),数论分块复杂度是 \(\Theta(t\sqrt n (\log p+\log \lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor))\) 小写 \(t\) 是数据组数。

实现看代码(考试时忘了写逆元卡了好久导致 T2 LCT 都没写完):

#include

#define ll long long
#define ull unsigned long long
#define INL inline
#define Re register

//Tosaka Rin Suki~

INL int read()
{
	int x=0,w=1;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();if(ch=='-')w=-1,ch=getchar();
	while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+ch-48,ch=getchar();return x*w;
}

const int N=1e6+5;
const ll MOD=1e9+7;

ll f[N],finv[N],prime[N],mu[N],cnt,prod[N];
bool vis[N];

INL ll fpow(ll x,ll p)
{
	ll res=1;
	while(p)
	{
		if(p&1)res=(1ll*res*x)%MOD;
		x=(1ll*x*x)%MOD;p>>=1;
	}
	return res;
}

INL void get_mu(ll n)
{
	f[0]=0,f[1]=1;
	for(int i=2;i<=n;i++)f[i]=(f[i-1]+f[i-2])%MOD;
	for(int i=0;i<=n;i++)finv[i]=fpow(f[i],MOD-2),prod[i]=1;
	vis[1]=1,vis[0]=1,mu[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!vis[i])prime[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
		{
			vis[prime[j]*i]=1;
			if(i%prime[j])mu[i*prime[j]]=-mu[i];
			else {mu[i*prime[j]]=0;break;}
		}
	}
	for(int i=1;i<=n;i++)
	{
		for(int j=i;j<=n;j+=i)
		{
			int now=j/i;
			if(mu[now]==-1)prod[j]=1ll*prod[j]*finv[i]%MOD;
			else if(mu[now]==1)prod[j]=1ll*prod[j]*f[i]%MOD;
		}
	}
	for(int i=1;i<=n;i++)prod[i]=1ll*prod[i-1]*prod[i]%MOD;
}

int main()
{
	//freopen("product.in","r",stdin);
	//freopen("product.out","w",stdout);
	int T=read();
	get_mu(1000000);
	while(T--)
	{
		ll n=read(),m=read();
		if(n