BZOJ4833 [Lydsy1704月赛] 最小公倍佩尔数

Description

\((1+\sqrt2)^n=e_n+f_n\sqrt2\),其中 \(e_n, f_n\) 都是整数。

再令 \(g_n=\mathop{\rm lcm}_{i=1}^n f_i\)

给定 \(n, p\),求 \(\left(\sum_{i=1}^n ig_i\right)\bmod p\)

数据组数 \(T\leqslant210,n\leqslant10^6,2\leqslant p\leqslant10^9+7\)\(p\) 是质数,\(f_1\dots f_n\)\(\bmod p\) 意义下不为 \(0\)

Solution

根据容斥我们可以得到

\[\mathop{\rm lcm}_{i\in S} f_i=\prod_{\emptyset\neq T\subseteq S}\left(\gcd_{i\in T}f_i\right)^{(-1)^{|T|-1}}\]

这实际上是 min-max 容斥

\[\max(S)=\sum_{\emptyset\neq T\subseteq S}(-1)^{|T|-1}\min(T)\]

的应用。

另外,可以发现 \(\gcd(f_i,f_j)=f_{\gcd(i,j)}\) (实际上,所有 \(f_i=Af_{i-1}+Bf_{i-2},f_0=0,f_1=1\)\(\gcd(A,B)=1\) 的序列 \(f\) 都满足这个性质)。

于是

\[\begin{aligned} g_n=&\mathop{\rm lcm}_{i=1}^n f_i\\ =&\prod_{\emptyset\neq T\subseteq S}f_{\gcd(T)}^{(-1)^{|T|-1}}\\ =&\prod_{i=1}^nf_i^{\sum_{i\mid t\leq n}\mu(t/i)}\\ =&\prod_{i=1}^n\left(\prod_{d\mid i}f_d^{\mu(i/d)}\right) \end{aligned}\]

所以求出所有的 \(f\) 之后莫比乌斯反演出 \(\prod_{d\mid i}f_i^{\mu(i/d)}\),然后前缀积就可以得到 \(g\) 了。

Code

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
#include <algorithm>
#include <cstdio>
#include <cstring>

typedef long long LL;

const int N = 1000050;

int n, mod, pr[N], cnt;
bool vis[N];
LL h[N], hi[N];

void Sieve() {
for (int i = 2; i < N; ++i) {
if (!vis[i]) pr[cnt++] = i;
for (int j = 0; j < cnt && i * pr[j] < N; ++j) {
vis[i * pr[j]] = 1;
if (i % pr[j] == 0) break;
}
}
}

LL pow_mod(LL a, LL b) {
if ((b %= mod - 1) < 0) b += mod - 1;
LL ans = 1;
for (a %= mod; b; b >>= 1, a = a * a % mod)
if (b & 1) ans = ans * a % mod;
return ans;
}

int main() {
Sieve();
int T;
scanf("%d", &T);
while (T--) {
scanf("%d%d", &n, &mod);
h[1] = hi[1] = 1;
for (int i = 2; i <= n; ++i)
hi[i] = pow_mod(h[i] = (h[i - 1] * 2 + h[i - 2]) % mod, mod - 2);
for (int i = 0; i < cnt && pr[i] <= n; ++i)
for (int j = n / pr[i]; j > 0; --j) {
(h[j * pr[i]] *= hi[j]) %= mod;
(hi[j * pr[i]] *= h[j]) %= mod;
}
LL ans = 0, g = 1;
for (int i = 1; i <= n; ++i) {
g = g * h[i] % mod;
ans = (ans + i * g) % mod;
}
printf("%lld\n", (ans + mod) % mod);
}
}