[Hard] P5320 [BJOI2019] 勘破神机

简要题意

给定 $m$，记 $f(n)$ 表示用 $1\times 2$ 的骨牌覆盖 $m\times n$ 的网格的方案数。给定 $l,r,k$，你需要求：

$$\frac{1}{r-l+1}\sum_{i=l}^{r}\binom{f(i)}{k}$$

答案对 $p=998,244,353$ 取模。

$1\leq l\leq r\leq 6\times 10^{18},p\nmid r-l+1,m\in\{2,3\},1\leq k\leq 501$。

思路

由于 $m\in\{2,3\}$，所以可以对于 $m$ 的不同取值分类讨论，求出此时 $f(n)$ 的表达式。

对于 $m=2$，我们发现骨牌覆盖一定由下面两种基本型组成：

所以很容易列出递推式：

$$f(i)=\begin{cases} f(i-1)+f(i-2)& i\geq 2\\ 1& i\leq 1 \end{cases}$$

那么就可以注意到：

$$f(n)=\mathrm{fib}(n+1)=\frac{\sqrt{5}}{5}\left(-\left(\frac{1-\sqrt{5}}{2}\right)^{n+1}+\left(\frac{1+\sqrt{5}}{2}\right)^{n+1}\right)$$

然后是 $m=3$，首先由三种平凡的基本型，它们都是 $3\times 2$ 的：

然后其实还有两种方法，类似这样，它们是 $3\times 2n$ 的（其中 $n\in[2,\infty)\cap\mathbb{N}$）：

即第一行 / 最后一行全是 $1\times 2$ 的骨牌覆盖，两侧是 $2\times 1$ 的骨牌覆盖，其余部分都是 $1\times 2$​ 的骨牌覆盖。

那么可以列出递推式：

$$g(n)=3g(n-2)+2\sum_{k=2}^{\lfloor\frac{n}{2}\rfloor} g(n-2k)$$

边界 $g(0)=1,g(1)=0$​。

注意到 $g(n)$ 只会被 $g(m)$（其中 $n\equiv m\pmod{2}$ 且 $m<n$）贡献，又因为 $g(1)=0$，所以我们发现 $g(n)=0$ 当且仅当 $2\nmid n$。

不妨令 $g'(n)=g(2n)$，则可以得到：

$$g'(n)=3g'(n-1)+2\sum_{k=2}^{n}g'(n-k)=3g'(n-1)+2\sum_{k=0}^{n-2}g'(k)$$

边界 $g'(0)=1$，现在需要找到它的通项公式，不妨构造它的生成函数：

$$\begin{aligned} \mathbf{G}'(z)&=3z\mathbf{G}'(z)+2\frac{z^2}{1-z}\mathbf{G}'(z)+1\\ (1-z)(1-3z)\mathbf{G}'(z)&=2z^2\mathbf{G}'(z)+1-z\\ ((1-z)(1-3z)-2z^2)\mathbf{G}'(z)&=1-z\\ \mathbf{G}'(z)&=\frac{1-z}{(1-z)(1-3z)-2z^2}=\frac{1-z}{1-4z+z^2} \end{aligned}$$

这也是一个有理式，因此它是一个（二阶）线性递推，尝试推导它的通项公式：

$$\mathbf{G}'(z)=\frac{1-z}{1-4z+z^2}=\frac{1-z}{(z-2+\sqrt{3})(z-2-\sqrt{3})}=\frac{A}{z-2+\sqrt{3}}+\frac{B}{z-2-\sqrt{3}}$$

其中：

$$\begin{aligned} &A(z-2-\sqrt{3})+B(z-2+\sqrt{3})=(A+B)z+(-2-\sqrt{3})A+(-2+\sqrt{3})=1-z\\ &\implies \begin{cases} A+B=-1\\ (-2-\sqrt{3})A+(-2+\sqrt{3})B=1 \end{cases} \implies \begin{cases} A=-\frac{3-\sqrt{3}}{6}\\ B=-\frac{3+\sqrt{3}}{6} \end{cases} \end{aligned}$$

由于：

$$\begin{aligned} &\frac{A}{z-2+\sqrt{3}}=A\sum_{i\geq 0}(-1)^iz^i(\sqrt{3}-2)^{-1-i}=\sum_{i\geq 0}A(-1)^i(-2-\sqrt{3})^{i+1}z^i=\sum_{i\geq 0} -A(2+\sqrt{3})^{i+1}z^i\\ &\frac{B}{z-2-\sqrt{3}}=B\sum_{i\geq 0}(-1)^iz^i(-\sqrt{3}-2)^{-1-i}=\sum_{i\geq 0}B(-1)^i(\sqrt{3}-2)^{i+1}z^i=\sum_{i\geq 0}-B(2-\sqrt{3})^{i+1}z^i \end{aligned}$$

故：

$$g'(n)=[z^n]\mathbf{G}'(z)=\frac{3-\sqrt{3}}{6}(2+\sqrt{3})^{n+1}+\frac{3+\sqrt{3}}{6}(2-\sqrt{3})^{n+1}$$

因此：

$$g(n)=\begin{cases} \frac{3-\sqrt{3}}{6}(2+\sqrt{3})^{\frac{n}{2}+1}+\frac{3+\sqrt{3}}{6}(2-\sqrt{3})^{\frac{n}{2}+1}& 2\mid n\\ 0& 2\nmid n \end{cases}$$

根据通项公式的形式，可以发现最终我们求解的问题可以转换为：

给定 $n,k,\alpha,\beta,x,y$，你需要求：

$$\sum_{i=0}^{n}\binom{\alpha x^i+\beta y^i}{k}$$

尝试用下降幂转普通幂公式推式子：

$$\begin{aligned} \sum_{i=0}^{n}\binom{\alpha x^i+\beta y^i}{k}&=\frac{1}{k!}\sum_{i=0}^{n}(\alpha x^i+\beta y^i)^{\underline{k}}=\frac{1}{k!}\sum_{i=0}^{n}\sum_{j\geq 0}(-1)^{k-j}{k\brack j}(ax^i+\beta y^i)^j\\ &=\frac{1}{k!}\sum_{j\geq 0}(-1)^{k-j}{k\brack j}\sum_{i=0}^{n}(\alpha x^i+\beta y^i)^j \end{aligned}$$

考察内侧的和式（这一部分的 $k$ 被用来当求和指标了，另外我们认为 $x^ky^{j-k}\neq 1$）：

$$\begin{aligned} \sum_{i=0}^{n}(\alpha x^i+\beta y^i)^j&=\sum_{i=0}^{n}\sum_{k\geq 0}\binom{j}{k}(\alpha x^i)^k(\beta y^i)^{j-k}=\sum_{k\geq 0}\binom{j}{k}\sum_{i=0}^{n}(\alpha x^i)^k (\beta y^i)^{j-k}\\ &=\sum_{k\geq 0}\binom{j}{k}\alpha^k\beta^{j-k}\sum_{i=0}^{n}x^{ik}y^{ij-ik}=\sum_{k\geq 0}\binom{j}{k}\alpha^k\beta^{j-k}\sum_{i=0}^{n}(x^ky^{j-k})^i\\ &=\sum_{k\geq 0}\binom{j}{k}\alpha^k\beta^{j-k}\cdot \frac{(x^ky^{j-k})^{n+1}-1}{x^ky^{j-k}-1} \end{aligned}$$

计算后面这个东西，时间复杂度 $O(j\log n)$，所以总时间复杂度 $O(k^2\log n)$。

你以为这就完了吗？实际上是没有的。因为 $3,5$ 在 $\mathbb{Z}/p\mathbb{Z}$ 种都是非二次剩余的！换句话说我们在 $\mathbb{Z}/p\mathbb{Z}$ 中找不到和式的 $\alpha,\beta,x,y$！

我们可以在 Mathematica 检验这一点：

{JacobiSymbol[5, 998244353], JacobiSymbol[3, 998244353]} (* 输出 {-1, -1} *)

不过这也是很容易解决的，下面以 $\sqrt{3}$ 为例。

我们定义 $\mathbf{j}=\sqrt{3}$，那么 $\mathbf{j}$ 与 $1$ 在有限域 $\mathbb{Z}/p\mathbb{Z}$ 上线性无关，所以可以构造一组次数为 $2$ 的域扩张 $K$，选取基底 $\{1,\mathbf{j}\}$（即 $\forall i\in K,\exists a,b\in\mathbb{Z}/p\mathbb{Z},i=a+b\mathbf{j}$。那么 $K$ 亦是域，因此可以进行所有原本 $\mathbb{Z}/p\mathbb{Z}$ 上的运算。这样就解决了不存在 $\sqrt{3}$ 的问题。

时间复杂度 $O(k^2\log n)$。

代码

#include <bits/stdc++.h>
using namespace std;

using i64 = long long;
const int N = 505;
constexpr int mod = 998244353, inv2 = (mod + 1) >> 1;
int Add(int x, int y){ return (x + y) >= mod ? (x + y - mod) : (x + y); }
int Sub(int x, int y){ return (x - y) < 0 ? (x - y + mod) : (x - y); }
int Mul(int x, int y){ return 1ll * x * y % mod; }

int fastpow(int x, int y){
    int ret = 1;
    for(;y;y>>=1,x=Mul(x, x)){ if(y & 1) ret = Mul(ret, x); }
    return ret;
}

template<int S>
struct Quadratic{
    int a, b;
    using T = Quadratic<S>;
    Quadratic(int a = 0, int b = 0): a(a), b(b) {}
    T operator+(const T& rhs) const { return {Add(a, rhs.a), Add(b, rhs.b)}; }
    T operator-(const T& rhs) const { return {Sub(a, rhs.a), Sub(b, rhs.b)}; }
    T operator*(const T& rhs) const {
        return {Add(Mul(a, rhs.a), Mul(Mul(b, rhs.b), S)), Add(Mul(a, rhs.b), Mul(b, rhs.a))};
    }
    T operator~() const {
        int inv = fastpow(Sub(Mul(a, a), Mul(Mul(b, b), S)), mod - 2);
        return {Mul(a, inv), Mul(mod - b, inv)};
    }
    T operator/(const T& rhs) const { return *this * ~rhs; }
    template<class U>
    T operator^(U y) const {
        T a = *this, ret = 1;
        for(;y;y>>=1,a=a*a){
            if(y & 1) ret = ret * a;
        }
        return ret;
    }
    int extract() const { return assert(!b), a; }
    T conj() const { return {a, mod - b}; }
};

template<int S>
using Q = Quadratic<S>;
using Q5 = Q<5>;
using Q3  = Q<3>;
int binom[N][N], stirling[N][N], m;

void init(int n){
    for(int i=0;i<=n;i++){
        binom[i][0] = 1;
        stirling[i][0] = i == 0;
        for(int j=1;j<=i;j++){
            binom[i][j] = Add(binom[i - 1][j], binom[i - 1][j - 1]);
            stirling[i][j] = Add(Mul(i - 1, stirling[i - 1][j]), stirling[i - 1][j - 1]);
        }
    }
}

template<int S>
Q<S> solve(i64 n, int k, Q<S> a, Q<S> b, Q<S> x, Q<S> y){
    Q<S> ans;
    for(int j=0;j<=k;j++){
        Q<S> ret = 0;
        for(int i=0;i<=j;i++){
            Q<S> t = (x ^ i) * (y ^ (j - i)), tmp;
            if(t.a == 1 && t.b == 0) tmp = (n + 1) % mod;
            else tmp = ((t ^ (n + 1)) - 1) / (t - 1);
            ret = ret + tmp * binom[j][i] * (a ^ i) * (b ^ (j - i));
        }
        if((k - j) & 1) ans = ans - ret * stirling[k][j];
        else ans = ans + ret * stirling[k][j];
    }
    for(int i=1;i<=k;i++) ans = ans * fastpow(i, mod - 2);
    return ans;
}

int solve(i64 n, int k){
    if(m == 2){
        Q5 a = Q5(0, mod - 1) / 5, b = a.conj();
        Q5 x = Q5(inv2, mod - inv2), y = x.conj();
        a = a * x, b = b * y;
        return solve(n, k, a, b, x, y).extract();
    }
    else{
        Q3 a = Q3(inv2, mod - fastpow(6, mod - 2)), b = a.conj();
        Q3 x = Q3(2, 1), y = x.conj();
        a = a * x, b = b * y;
        return solve(n, k, a, b, x, y).extract();
    }
}

void solve(){
    i64 L, R; int k;
    cin >> L >> R >> k;
    int ans = fastpow((R - L + 1) % mod, mod - 2);
    if(m == 3) L = (L + 1) >> 1, R >>= 1;
    cout << Mul(ans, Sub(solve(R, k), solve(L - 1, k))) << '\n';
}

signed main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    int t; cin >> t >> m;
    init(504);
    while(t--) solve();
    return 0;
}

// Written by xiezheyuan