【luogu P5320】勘破神机（数学）（数列特征方程）（第一类斯特林数）

勘破神机

题目链接：luogu P5320

题目大意

给你 $l,r,k$（其中 $l,r$ 很大，$k\leq 501$），求：
$\dfrac{1}{r-l+1}\sum\limits_{i=l}^rC_{f_i}^k$
$\dfrac{1}{r-l+1}\sum\limits_{i=l}^rC_{g_i}^k$
$f_n$ 为一个 $2*n$ 的网格往里面放 $1*2$ 的块填满的方案，$g_n$ 则是一个 $3*n$ 的网格。

思路

首先 $2$ 的答案我们可以看出答案其实是贴合与斐波那契数列 $F_i$。
具体看一看会发现是 $f_i=F_{i+1}$，所以记得 $l,r$ 加一。
然后斐波那契嘛，$f_i=f_{i-1}+f_{i-2}$，然后可以表示成那个带根号的通式嘛。

然后我们通过特征方程这个东西，我们知道如果有一个二阶递推式，我们可以把它表示成一个通项公式。
于是开始了下面的步骤：

2

$f_i=f_{i-1}+f_{i-2}$
$f[n]=\dfrac{\sqrt{5}}{5}[(\dfrac{1+\sqrt{5}}{2})^n-(\dfrac{1-\sqrt{5}}{2})^n]$
$A=\dfrac{\sqrt{5}}{5},\alpha=\dfrac{1+\sqrt{5}}{2},B=-\dfrac{\sqrt{5}}{5},\beta=\dfrac{1-\sqrt{5}}{2}$

3

思考要如何推出 $f$ 的值，首先把 $f_0,f_1$ 求出，因为可以作为边界，而且也是特征方程需要的。

$f_0=1,f_1=3$
$f_i$ 长度为 $2i$ 的答案（所以记得两个边界要改，但是注意 $r-l+1$ 要用之前的值）

于是进行思考，发现可以新放进来一个 $2*3$ 的，然后有 $3$ 中方法。
然后发现只有这样不够，因为可能会出现这样的：

而且呢，中间的粉色可能不止两个：

然后粉色可以是在 $1,2$ 行或者 $2,3$ 行，所以我们可以列出这样的式子：

$f_{i}=3f_{i-1}+2(f_{i-2}+f_{i-3}+...)$
$f_i=3f_{i-1}+2(\sum\limits_{j=0}^{i-2}f_j)$
$f_{i-1}=3f_{i-2}+2(\sum\limits_{j=0}^{i-3}f_j)$
$f_i-f_{i-1}=3f_{i-1}-f_{i-2}$
$f_i=4f_{i-1}-f_{i-2}$

搞出了二阶递推式，就直接开始求通项公式吧！

$f_i-4f_{i-1}+f_{i-2}=0$
$x^2-4x+1=0$
$x_1=\dfrac{4+2\sqrt{3}}{2},x_2=\dfrac{4-2\sqrt{3}}{2}$
$x_1=2+\sqrt{3},x_2=2-\sqrt{3}$
带入 $f_0=1,f_1=3$
$Ax_1^n+Bx_2^n=f_n$

$A+B=1$
$A(2+\sqrt{3})+B(2-\sqrt{3})=3$

$A(2+\sqrt3)+(1-A)(2-\sqrt3)=3$
$A(2+\sqrt3)+2-\sqrt{3}-A(2-\sqrt3)=3$
$2\sqrt3A=1+\sqrt3$
$A=\dfrac{1+\sqrt3}{2\sqrt3}=\dfrac{\sqrt3+3}{6}$
$A=\dfrac{3+\sqrt3}{6},B=\dfrac{3-\sqrt3}{6}$
$f[n]=\dfrac{3+\sqrt3}{6}(\dfrac{4+2\sqrt{3}}{2})^n+\dfrac{3-\sqrt3}{6}(\dfrac{4-2\sqrt{3}}{2})^n$
$A=\dfrac{3+\sqrt3}{6},\alpha=\dfrac{4+2\sqrt{3}}{2},B=\dfrac{3-\sqrt3}{6},\beta=\dfrac{4-2\sqrt{3}}{2}$
$A=\dfrac{3+\sqrt3}{6},\alpha=2+\sqrt{3},B=\dfrac{3-\sqrt3}{6},\beta=2-\sqrt{3}$

---

于是我们似乎把 $f_i$ 表示成了一个很好的形式，考虑开搞。

$f_i=A\alpha^i+B\beta^i$
$\sum\limits_{i=0}^n\binom{f_i}{k}$
$\sum\limits_{i=0}^n\dfrac{f_i!}{k!(f_i-k)!}$
$\sum\limits_{i=0}^n\dfrac{f_i^{\underline k}}{k!}$
$\dfrac{1}{k!}\sum\limits_{i=0}^n\sum\limits_{j=1}^k(-1)^{k-j}s(k,j)f_i^j$（斯特林数与下降幂的转化）
（补充：$x^{\overline k}=\sum\limits_{j=1}^ks(k,j)x^j$）
$\dfrac{1}{k!}\sum\limits_{i=1}^k(-1)^{k-i}s(k,i)\sum\limits_{j=0}^nf_j^i$

设 $g_{n,j}=\sum\limits_{i=0}^nf_i^j$
$f_i=A\alpha^i+B\beta^i$
$g_{n,j}=\sum\limits_{i=0}^n(A\alpha^i+B\beta^i)^j$
$g_{n,j}=\sum\limits_{k=0}^j\binom{j}{k}\sum\limits_{i=0}^n(A\alpha^i)^k(B\beta^i)^{j-k}$
$g_{n,j}=\sum\limits_{k=0}^j\binom{j}{k}A^kB^{j-k}\sum\limits_{i=0}^n\alpha^{ik}\beta^{i(j-k)}$
$g_{n,j}=\sum\limits_{k=0}^j\binom{j}{k}A^kB^{j-k}\dfrac{1-\alpha^{(n+1)k}\beta^{(n+1)(j-k)}}{1-\alpha^k\beta^{j-k}}$（等比数列求和）

$\dfrac{1}{k!}\sum\limits_{i=1}^k(-1)^{k-i}s(k,i)\sum\limits_{j=0}^nf_j^i$
$\dfrac{1}{k!}\sum\limits_{i=1}^k(-1)^{k-i}s(k,i)g_{n,i}$
$\dfrac{1}{k!}\sum\limits_{i=1}^k(-1)^{k-i}s(k,i)(\sum\limits_{p=0}^i\binom{i}{p}A^pB^{i-p}\dfrac{1-(\alpha^{p}\beta^{i-p})^{n+1}}{1-\alpha^p\beta^{i-p}})$

---

小提示：
记得特判等比数列公比为 $1$。

代码

#include<cstdio>
#define ll long long
#define mo 998244353

using namespace std;

const int K = 550;
int gen;
int T, op, k;
ll C[K][K], S[K][K], l, r;

ll ksm(ll x, ll y) {
	ll re = 1;
	while (y) {
		if (y & 1) re = re * x % mo;
		x = x * x % mo; y >>= 1;
	}
	return re;
}

struct complex {
	ll x, y;
	
	complex(ll a = 0, ll b = 0) {
		x = a; y = b;
	}
}A, B, alpha, beta;//x+sqrt{y}
complex operator +(complex x, complex y) {
	return (complex){(x.x + y.x) % mo, (x.y + y.y) % mo};
}
complex operator -(complex x) {return (complex){(mo - x.x) % mo, (mo - x.y) % mo};}
complex operator -(complex x, complex y) {return x + (-y);}
complex operator *(complex x, complex y) {
	return (complex){(x.x * y.x + gen * x.y * y.y) % mo, (x.x * y.y + x.y * y.x) % mo};
}
complex inv(complex x) {
	ll fm = ksm((((x.x * x.x - gen * x.y * x.y) % mo + mo) % mo) % mo, mo - 2);
	return (complex){x.x * fm % mo, (mo - x.y) * fm % mo};
}

complex ksm(complex x, ll y) {
	complex re = complex(1, 0);
	while (y) {
		if (y & 1) re = re * x;
		x = x * x; y >>= 1;
	}
	return re;
}

void Init() {
	C[0][0] = 1;
	for (int i = 1; i < K; i++) {
		C[i][0] = 1; for (int j = 1; j < K; j++) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mo;
	}
	S[0][0] = 1;
	for (int i = 1; i < K; i++) {
		for (int j = 1; j < K; j++)
			S[i][j] = (S[i - 1][j - 1] + (i - 1) * S[i - 1][j] % mo) % mo;
	}
	
	gen = (op == 2) ? 5 : 3;
	if (op == 2) {
		A = (complex){0, ksm(5, mo - 2)}; alpha = (complex){ksm(2, mo - 2), ksm(2, mo - 2)};
		B = (complex){0, (mo - ksm(5, mo - 2)) % mo}; beta = (complex){ksm(2, mo - 2), (mo - ksm(2, mo - 2)) % mo};
	}
	if (op == 3) {
		A = (complex){ksm(2, mo - 2), ksm(6, mo - 2)}; alpha = (complex){2, 1};
		B = (complex){ksm(2, mo - 2), (mo - ksm(6, mo - 2)) % mo}; beta = (complex){2, mo - 1};
	}
}

ll slove(ll n) {
	ll ans = 0, jc = 1; ll di = ((k & 1) ? 1 : mo - 1);
	for (int i = 1; i <= k; i++, di = mo - di) {
		complex sum; jc = jc * i % mo;
		for (int p = 0; p <= i; p++) {
			complex alphabeta = ksm(alpha, p) * ksm(beta, i - p);
			//你是一个一个一个特判啊啊啊！！！
			if (alphabeta.x == 1 && alphabeta.y == 0) {//特判公比为1的情况
				sum = sum + ((complex){C[i][p], 0} * ksm(A, p) * ksm(B, i - p) * (complex){(n + 1) % mo, 0});
				continue;
			}
			sum = sum + ((complex){C[i][p], 0} * ksm(A, p) * ksm(B, i - p) * ((complex){1, 0} - ksm(alphabeta, n + 1)) * inv((complex){1, 0} - alphabeta));
		}
		if (sum.y) {printf("f\n"); return -1;}
		(ans += di * S[k][i] % mo * sum.x % mo) %= mo;
	}
	return ans * ksm(jc, mo - 2) % mo;
}

int main() {
	scanf("%d %d", &T, &op);
	Init();
	while (T--) {
		scanf("%lld %lld %d", &l, &r, &k); ll sz = (r - l + 1) % mo;
		if (op == 2) l++, r++;
			else l = (l + 1) / 2, r = r / 2;
		printf("%lld\n", (slove(r) - slove(l - 1) + mo) % mo * ksm(sz, mo - 2) % mo);
	}
	
	return 0;
} 

__EOF__

本文作者：あおいSakura
本文链接：https://www.cnblogs.com/Sakura-TJH/p/luogu_P5320.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！