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)$，通过打表观察数学分析可以得到：
扩域模拟计算，可得：
\[\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 代码
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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
//#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;
}
