题解：最小公倍佩尔数

题目传送门

gcd-lcm 容斥

所谓 gcd-lcm 容斥，即对于集合 \(T\)，有：

\[\operatorname*{lcm}_{i\in T}i=\prod_{S\subseteq T}\left(\gcd_{i\in S}i\right)^{(-1)^{\vert S\vert-1}} \]

证明

不妨考虑 $i$ 的质因子。发现对于质数 $p$，有：

$$ \begin{aligned} \operatorname{lcm}\left(p^a,p^b\right)&=p^{\max(a,b)}\\ \gcd\left(p^a,p^b\right)&=p^{\min(a,b)}\\ \end{aligned} $$

类似的套用 Min-max 容斥可以得到指数间关系，化简即可得到上式。

考虑 \(g(n)=\operatorname{lcm}(f(1),f(2),\cdots,f(n))\) 并不好处理，因为 \(\operatorname{lcm}\) 不可模。

因此考虑转化，记 \(T=\set{1,2,3,\cdots,n}\)，有：

\[\begin{aligned} g(n)&=\operatorname{lcm}(f(1),f(2),\cdots,f(n-1),f(n))\\ &=\operatorname*{lcm}_{i\in T}f(i)\\ &=\prod_{S\subseteq T}\left(\gcd_{i\in S}f(i)\right)^{(-1)^{\vert S\vert-1}}\\ \end{aligned} \]

递推数列性质

考虑 \(f(i)\)。显然，\(f(0)=0,f(1)=1\)。

因为 \(\left(1+\sqrt2\right)^n=e(n)+f(n)\)，通过打表观察数学分析可以得到：

\[f(n)=2f(n-1)+f(n-2)\\ \]

扩域模拟计算，可得：

\[\left(e(n)+f(n)\sqrt2\right)\left(1+\sqrt2\right)=(e(n)+2f(n))+(e(n)+f(n))\sqrt2 \]

进而可推出：

\[e(n)=e(n-1)+2f(n-1)\\ f(n)=e(n-1)+f(n-1) \]

递归代入，可得：

\[\begin{aligned} f(n)&=f(n-1)+2\sum_{i=1}^{n-2}f(i)\\ &=f(n-1)+2f(n-2)+\sum_{i=1}^{n-3}f(i)\\ &=2f(n-1)+(n-2) \end{aligned} \]

可以发现，这是一个类似于斐波那契数列的线性递推数列，同样具有一些广义斐波那契数列的美妙性质。

有：

\[\gcd(f(i),f(j))=f(\gcd(i,j)) \]

证明请参见我的 Blog。因为本人认为 OI 范围内写数学相关题目，常用性质证明过程并不重要。

莫比乌斯反演

根据我们得到的 \(f\) 的性质，有：

\[\begin{aligned} g(n)&=\prod_{S\subseteq T}\left(\gcd_{i\in S}f(i)\right)^{(-1)^{\vert S\vert-1}}\\ &=\prod_{S\subseteq T}\left(f\left(\gcd_{i\in S}i\right)\right)^{(-1)^{\vert S\vert-1}} \end{aligned} \]

\(\gcd\) 相关信息，考虑莫比乌斯反演。

先套路地枚举 \(\gcd\limits_{i\in S}i\) 算贡献：

\[\begin{aligned} g(n)&=\prod_{S\subseteq T}\left(f\left(\gcd_{i\in S}i\right)\right)^{(-1)^{\vert S\vert-1}}\\ &=\prod_{d=1}^nf(d)^{\sum\limits_{S\subseteq T}\left[\gcd\limits_{i\in S}i=d\right](-1)^{(\vert S\vert-1)}}\\ \end{aligned} \]

记：

\[\begin{aligned} h(d)&=\sum\limits_{S\subseteq T}\left[\gcd\limits_{i\in S}i=d\right](-1)^{\vert S\vert-1}\\ H(i)&=\sum_{i\mid d}h(d) \end{aligned} \]

则根据莫比乌斯变换，有：

\[h(i)=\sum_{i\mid d}\mu\left(\dfrac di\right)H(i) \]

同理，通过打表观察数学分析，我们可以惊喜地发现 \(H(i)=1\) 恒成立。

取集合 \(L=\left\lbrace1,2,3,\cdots,\left\lfloor\dfrac nd\right\rfloor\right\rbrace\)。

考虑到：

\[\begin{aligned} H(i)&=\sum_{i\mid d}h(d)\\ &=\sum_{A\subseteq L,A\neq\varnothing}(-1)^{\vert A\vert-1}\\ &=-\left(\sum_{A\subseteq L}(-1)^{\vert A\vert}-1\right)\\ &=-(0-1)\\ &=1 \end{aligned} \]

因此有：

\[h(i)=\sum_{i\mid d}\mu\left(\dfrac di\right) \]

最终推导

我们终于得到了 \(g(n)\)：

\[\begin{aligned} g(n)&=\prod_{d=1}^n\prod_{i\mid d}f^{\mu\left(\frac di\right)}(d)\\ &=g(n-1)\prod_{i\mid n}f^{\mu\left(\frac ni\right)}(n) \end{aligned} \]

预处理可以做到 \(\mathcal O(n\log n)\)。

AC 代码

//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
typedef __int128 lll;
typedef long long ll;
constexpr const int N=1e6;
int n,P;
int g[N+1],f[N+1],mu[N+1],s[N+1];
//a+b\sqrt2
struct special{
	lll a,b;
};
special operator *(const special &x,const special &y){
	return {(1ll*x.a*y.a+2ll*x.b*y.b)%P,(1ll*x.a*y.b+1ll*x.b*y.a)%P};
}
special& operator *=(special &x,const special &y){
	return x=x*y;
}
int qpow(int base,int n){
	int ans=1;
	while(n){
		if(n&1){
			ans=1ll*ans*base%P;
		}
		base=1ll*base*base%P;
		n>>=1;
	}
	return ans;
} 
void pre(){
	static int vis[N+1],prime[N+1],size;
	mu[1]=1;
	for(int i=2;i<=N;i++){
		if(!vis[i]){
			vis[i]=i;
			prime[++size]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=size&&i*prime[j]<=N;j++){
			vis[i*prime[j]]=prime[j];
			if(i%prime[j]==0){
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
}
int main(){
	/*freopen("test.in","r",stdin);
	freopen("test.out","w",stdout);*/
	
	ios::sync_with_stdio(false);
	cin.tie(0);cout.tie(0);
	
	pre();
	int T;
	cin>>T;
	while(T--){
		cin>>n>>P;
		special x={1,0};
		for(int i=1;i<=n;i++){
			x*={1,1};
			f[i]=x.b;
		}
		int ans=0;
		g[0]=1;
		fill(s+1,s+n+1,1);
		for(int d=1;d<=n;d++){
			int inv=qpow(f[d],P-2);
			for(int t=1;d*t<=n;t++){
				switch(mu[t]){
					case 1:
						s[d*t]=1ll*s[d*t]*f[d]%P;
						break;
					case -1:
						s[d*t]=1ll*s[d*t]*inv%P;
				}
			}
		}
		for(int i=1;i<=n;i++){
			g[i]=1ll*g[i-1]*s[i]%P;
			ans=(ans+1ll*i*g[i])%P;
		}
		cout<<ans<<'\n';
	}
	
	cout.flush();
	
	/*fclose(stdin);
	fclose(stdout);*/
	return 0;
}