「BJOI 2019」堪破神机【数论】

题目大意

设用 1×21 \times 2 填满 2×n2 \times n 的方案数为 fnf_n,填满 3×n3 \times n 的方案数为 gng_n

给定 k,mk, m (m{2,3})(m \in \{ 2, 3 \}),求出 ansm\text{ans}_m998 244 353998\ 244\ 353 的值,其中:

ans2=1rl+1n=lr(fnk)ans3=1rl+1n=lr(gnk)\begin{aligned} \text{ans}_2 &= \frac{1}{r - l + 1} \sum_{n = l}^r \binom{f_n}{k} \\ \text{ans}_3 &= \frac{1}{r - l + 1} \sum_{n = l}^r \binom{g_n}{k} \end{aligned}

保证 1lr1018,1 \leqslant l \leqslant r \leqslant 10^{18}, 1k5011 \leqslant k \leqslant 501

思路分析

初步转化

先用第一类斯特林数,将下降幂转为 kk 次幂的形式:

n=lr(fnk)=1k!n=lrfnk=1k!n=lri=0k(1)ki[ki]fni=1k!i=0k(1)ki[ki]n=lrfni\begin{aligned} \sum_{n = l}^r \binom{f_n}{k} &= \frac{1}{k!} \sum_{n = l}^r f_n^{\underline k} \\ &= \frac{1}{k!} \sum_{n = l}^r \sum_{i = 0}^k (-1)^{k - i} \begin{bmatrix} k \\ i \end{bmatrix} f_n^i \\ &= \frac{1}{k!} \sum_{i = 0}^k (-1)^{k - i} \begin{bmatrix} k \\ i \end{bmatrix} \sum_{n = l}^r f_n^i \end{aligned}

通项公式

m=2m = 2 为例。考虑 ff 的生成函数 F(z)=z1zz2F(z) = \small \dfrac{z}{1 - z - z^2},分解为:

F(z)=15(111+52z11152z)F(z) = \frac{1}{\sqrt 5} \Bigg( \frac{1}{1 - \frac{1 + \sqrt 5}{2} z} - \frac{1}{1 - \frac{1 - \sqrt 5}{2} z} \Bigg)

二项式展开,比较 ziz^i 系数得:

fi=15((1+52)i(152)i)f_i = \frac{1}{\sqrt 5} \Bigg( \bigg( \frac{1 + \sqrt 5}{2} \bigg)^i - \bigg( \frac{1 - \sqrt 5}{2} \bigg)^i \Bigg)

展开答案

A=15,A = \small \dfrac{1}{\sqrt 5}, B=15,B = - \small \dfrac{1}{\sqrt 5}, x=1+52,x = \small \dfrac{1 + \sqrt 5}{2}, y=152y = \small \dfrac{1 - \sqrt 5}{2}

则有 fn=Axn+Bynf_n = A x^n + B y^n,代回原式得:

ans2=1k!i=0k(1)ki[ki]n=lr(Axn+Byn)i=1k!i=0k(1)ki[ki]n=lrj=0i(ij)AjBij(xjyij)n=1k!i=0k(1)ki[ki]j=0i(ij)AjBijn=lr(xjyij)n\begin{aligned} \text{ans}_2 &= \frac{1}{k!} \sum_{i = 0}^k (-1)^{k - i} \begin{bmatrix} k \\ i \end{bmatrix} \sum_{n = l}^r \big( A x^n + B y^n \big)^i \\ &= \frac{1}{k!} \sum_{i = 0}^k (-1)^{k - i} \begin{bmatrix} k \\ i \end{bmatrix} \sum_{n = l}^r \sum_{j = 0}^i \binom{i}{j} A^j B^{i - j} \big( x^j y^{i - j} \big)^n \\ &= \frac{1}{k!} \sum_{i = 0}^k (-1)^{k - i} \begin{bmatrix} k \\ i \end{bmatrix} \sum_{j = 0}^i \binom{i}{j} A^j B^{i - j} \sum_{n = l}^r \big( x^j y^{i - j} \big)^n \end{aligned}

枚举 i,ji, j,后面的 n\sum_n 就是一个等比数列。快速幂计算即可。

扩域计算

还有一个问题,5\sqrt 5 在模 998 244 353998\ 244\ 353 意义下并不存在。

考虑在数域 Z[5]\mathbb{Z} \big[ \sqrt 5 \big] 下计算,将每个数表示为 a+b5a + b \sqrt 5 的形式。答案中 bb 必为 00

回到原题

m=3m = 3 也是类似的。首先当 nn 为奇数时,gn=0g_n = 0

hn=g2nh_n = g_{2n},可以求出递推式 hn=4hn1hn2h_n = 4 h_{n - 1} - h_{n - 2}

考虑 hh 的生成函数 H(z)=z14z+z2H(z) = \small \dfrac{z}{1 - 4 z + z^2},同理可解出通项公式:

hi=3+36(2+3)n+336(23)nh_i = \frac{3 + \sqrt 3}{6} \big( 2 + \sqrt 3 \big)^n + \frac{3 - \sqrt 3}{6} \big( 2 - \sqrt 3 \big)^n

在数域 Z[3]\mathbb{Z} \big[ \sqrt 3 \big] 下类似计算。总时间复杂度为 O(k2logr)O(k^2 \log r),常数较大。

代码实现

注意:在计算等比数列和时,需要判断公比是否为 11

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
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
const ll P = 998244353;
ll m, l, r, k, C[510][510], S[510][510];

ll qp(ll x, ll y) {
ll z = 1;
for (x %= P; y; y >>= 1, x = x * x % P) {
if (y & 1) z = z * x % P;
}
return z;
}

void init(ll K = 505) {
for (ll i = 0; i < K; i++) {
for (ll j = C[i][0] = 1; j <= i; j++) {
C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % P;
}
}
for (int i = S[0][0] = 1; i < K; i++) {
for (int j = 1; j <= i; j++) {
S[i][j] = (S[i - 1][j - 1] - (i - 1) * S[i - 1][j]) % P;
}
}
}

struct Z { // a + b sqrt(m)
ll a, b;
Z(ll x = 0, ll y = 0) : a(x % P), b(y % P) {}
Z conj() { return Z(a, -b); }
Z operator + (Z y) { return Z(a + y.a, b + y.b); }
Z operator - (Z y) { return Z(a - y.a, b - y.b); }
Z operator * (Z y) { return Z(a * y.a + m * b * y.b, a * y.b + b * y.a); }
Z operator / (Z y) { return *this * y.conj() * qp(y.a * y.a - m * y.b * y.b, P - 2); }
} A, B, x, y;
Z Zpow(Z x, ll y) {
Z z = 1;
for (; y; y >>= 1, x = x * x) {
if (y & 1) z = z * x;
}
return z;
}

Z solve() {
r++; Z ans = 0;
for (ll i = 0; i <= k; i++) {
for (ll j = 0; j <= i; j++) {
Z z = Zpow(x, j) * Zpow(y, i - j);
Z t = Zpow(A, j) * Zpow(B, i - j) * (S[k][i] * C[i][j] % P);
if (!(z - 1).a && !z.b) ans = ans + t * ((r - l) % P);
else ans = ans + t * (Zpow(z, r) - Zpow(z, l)) / (z - 1);
}
}
return ans;
}

int main() {
scanf("%*d %lld %lld %lld %lld", &m, &l, &r, &k);
init(); ll len = r - l + 1;
if (m == 2) {
l++, r++, m = 5;
A = Z(0, 1) / 5, x = Z(1, 1) / 2;
} else {
l = (l + 1) >> 1, r >>= 1;
if (l > r) printf("0\n"), exit(0);
A = Z(3, 1) / 6, x = Z(2, 1);
}
B = A.conj(), y = x.conj(); ll tmp = 1;
for (ll i = 1; i <= k; i++) tmp = tmp * i % P;
Z ans = solve() * qp(tmp, P - 2) * qp(len, P - 2);
printf("%lld\n", (ans.a % P + P) % P);
return 0;
}