线性代数与理性打表
1. 前言和前置技能
1.1 一些扯淡
有这样一个神仙，叫刘承奥。
他出了这样一个题，叫某少女附中的体育课。
这道题的题解的画风是这个样子的。
后来我就走上了一条不归路。
以下大部分内容都没有证明。
即得易见平凡，仿照上例显然，留作习题答案略，读者自证不难。
反之亦然同理，推论自然成立，略去过程Q.E.D.，由上可知证毕。
1.2 前置技能
- 多项式全家桶
2. 抽象代数
2.1 代数系统
我为什么要专门介绍这个和 OI 和高考基本上没什么关系的东西呢？
主要是为了能听懂 LCA 这种人说话。
以及抽象代数是数学的另一大支柱级学科，其地位甚至比微积分还要重要，大部分数学概念都是用抽象代数的知识来定义的，如果有人张口闭口就是什么环啊域啊之类的东西，咱们起码得能听懂是不是？而且我觉得最近不少 OIER 喜欢用抽象代数里的知识来给某个东西下一个严谨的定义，尤其是 LCA ，起码域这个词已经出现在 vfk 的论文里了。
而且我猜大家逃不过学他的命运，还是看看吧。
并没有什么奇奇怪怪的式子，都是些定义。证明啥的我也不会证，明白都是个什么东西就行。
抽象代数是研究某种“结构”的学科，正如我们每天研究数据结构一样。
数据结构研究的是数据与数据之间的“关系“：某个节点上的信息，它的”父亲“的信息，它的”儿子“的信息，他们之间的所有关系，我们称之为一个”结构“。
而抽象代数是研究数与数之间的“关系”的学科。换句话说说，它研究的是我们的”代数系统”。
我们在日常做数学计算的时候，有没有想过，为什么我们能做加减乘除？我们能不能自己定义一套用来计算的东西，包括我们计算的对象，以及运算法则呢？
难道不是日常吗？
typedef std::vector<int> Poly;
inline Poly operator + (const Poly&, const Poly&);
inline Poly operator - (const Poly&, const Poly&);
inline Poly operator * (const Poly&, const Poly&);
inline Poly operator / (const Poly&, const Poly&);
如果我们有一个用来运算的集合 (A) ，和一个 (n) 元函数 (f) ：
表示 (f) 以 (n) 个 (A) 集合里的元素作为参数，返回一个 (A) 集合里的数，我们就有了一个“代数系统”。我们不用去研究 (A) 里的元素如何如何，我们只对“这个系统有什么性质”以及“我能在这个系统里做什么”感兴趣。
一般来说，这个函数 (f) 是二元函数，比如加减乘除，我们用 (circ) 代替，然后用 (langle A,circ angle)来表示这个代数系统。
以下所有的代数系统都得满足封闭性，你不能有 undefined behavior。
2.2 运算律
抽象代数要从娃娃抓起。
下面，我们来回忆一下小学一年级的知识：运算律。
没错，多项式是七年级上册学的，抽象代数是从一年级上册学的。
回忆一下我们一年级都学了啥：
假设我们有种二元运算，不妨就是加法 (+) 吧！
它满足交换律：
结合律：
结合律让我们的结果只和序列有关而与运算顺序无关，交换律让我们的结果只和运算的变量有关。
如果我们还有一种二元运算，不妨就是乘法 ( imes) 吧！
它和加法一起满足分配律：
这就是我们日常所使用的代数系统的基石法则。将它放在一年级讲，也不无道理。
但是如果我们的系统不满足以上的性质呢？就像矩阵乘法不满足交换律一样。
如果一个代数系统只满足结合律，我们称之为半群。
但是代数系统还有别的性质。
2.3 单位元和逆元
如果一个代数系统存在这样一个元素 (e) 满足对任意元素 (ex= x) 我们称之为单位元。
如果一个代数系统存在这样一个元素 (x^{-1}) 使得 (x^{-1}x=e) 我们称 $x^{-1} $ 为 (x) 的逆元。
存在单位元的代数系统我们称之为幺半群。
满足结合律，存在单位元和逆元的代数系统我们称之为群。
满足交换律的群我们称之为交换群，也叫阿贝尔群，以纪念某位神仙。
想想我们之前的置换群，似乎也满足这样三个性质。
特别的，一个元素 (a) 自己与自己进行多次运算生成了一个群，那么这个群就是循环群，(a) 叫做这个群的生成元。有没有想起来原根？原根就是模 (p) 剩余系下的生成元，而且这个剩余系还是交换群。
置换群也是一个单独的类别，还有对称群，在此不表，感兴趣的同学可参阅普通《高中课程标准实验教科书 数学 选修三系列》。(然而并不是高考内容，大概连书都见不到吧。)
2.4 环和域
如果现在我们的群有两种运算 (+) 和 ( imes) ，而且 (+) 满足群的所有性质， (+) 和 ( imes) 满足分配率，我们就称之为环。
注意这个 ( imes) 不一定满足群的性质，它只需要是个半群。
更一般的，如果 (+) 和 ( imes) 都满足群所有的性质，而且还满足交换律，我们就称之为域。
我们日常与之打交道的东西除了整数集之外都是域 ，有理数域 (mathbb Q) ，实数域 (mathbb R) ，复数域 (mathbb C) ，还有模 (p) 剩余系 (mathbb F_p)。域的限制最少，所以我们一般定义的操作与数都是在域下定义的。
当然域还可以扩张，例如可以这样来扩张有理数域 ：
其中 (w) 是某个平方根，可能是无理数，也可能是虚数。这样就能让某个方程在这个数域内有解。
这就是二次域啦！有个求解二次剩余的算法( Cipolla )就是以此为基础的。感兴趣的同学可以参考 WC2019 朱老大讲的《简单数论》。
群论和抽象代数是一门十分博大精深的学科，这里只是解释一点点名词，更加详细的就要等到大学里去学了。希望没有带偏大家
大概知道上面这些就能听懂 LCA 这样的人说话了吧。
3. 线性代数
只讲行列式。因为我不会别的。
有时间的可以去 b 站看一看 3Bule1Brown 的视频。
3.1 行列式
矩阵的相关运算大家已经非常熟悉了。然而线性代数中还有一类重要的东西——行列式。
行列式是形如下面这个样子的东西（注意是方的）
它是一个式子，或者说并不是像矩阵那样是个存数的容器，它是一个表达式的另一种写法，是一个值。
如果我们有一个矩阵 (A) ，那么矩阵 $ A $ 的行列式可以记为 (det A) 或者 (|A|)，表示将 (A) 里的数排成上面那样的行列式。
道理我都懂，可是行列式到底应该怎么算？
其实行列式有许多定义的方法，这里介绍两种。他们本质上是一样的。
如果我们有一个排列 ((j_1j_2cdots j_n)) ，那么设函数 ( au(j_1j_2cdots j_n)) ，当 ((j_1j_2cdots j_n)) 的逆序对数为偶数时，( au(j_1j_2cdots j_n) = 0) ，否则 ( au(j_1j_2cdots j_n) = 1) 。
那么行列式的值就定义为：
显然，展开式中各项系数正负参半。
但是我们还有一个更容易理解的定义。
行列式最早是用来解多元一次方程组的。举个最简单的例子，二元一次方程组。我们虽然能很轻易的用高斯消元解出答案，但是如果这是在考场上，你在联立两个直线方程，然后用高消解？或者你在计算几何里高消出来直线交点？一个是慢，另一个是炸精。
为啥不用个 (O(1)) 的式子呢？
现在有方程
我们可以很容易地用代入消元法解出来
人（神仙）们从上式中发现，如果记
那么解就能表示成
我们把 (D) 称为二阶行列式。
而对于 (9) 个元素 (a_{ij}(i, j=1,2,3)) ，定义三阶行列式
我们会发现它其实是三条主对角线乘积之和减去三条副对角线乘积之和。
主对角线就是左上到右下的三个元素，不够三个的从另一边继续往下数就好了。副对角线就是右上到左下。
特别的，如果我们将方程写成这个样子：
然后将系数提出来作为一个行列式
那么 (x_j) 就等于将 (|A|) 的第 (j) 列替换成 ((b_1dots b_n)) 后的行列式的值除以 (|A|) 。
这就是克拉默(Cramer)法则。读者自证不难。
于是我们就得到了 (O(1)) 解方程的式子啦！个鬼！
首先上面三阶行列式的算法在 (n > 3) 的时候是假的！
按照最开始的定义，最后出来的式子得有 (n!) 项，然而那样算只有 (2n) 项，肯定是假的呀！
然后如果用行列式的定义来暴力算的话，解一项的花费是 (n! nlog n) ，因为还得算逆序对数。于是我们就有了一 (O(n!n^2log n)) 的优秀算法。刺不刺激？
那我们要这玩意有啥用啊！
先别急，我们先把那个 (log n) 去了再说。
我买的那本 《线性代数》上实际上是这样递归定义行列式的：
已知行列式 (D) ，设 (A_{ij}) 为 (D) 的代数余子式，表示将 (D) 去掉第 (i) 行和第 (j) 列后剩下的行列式乘上 ((-1)^{i+j}) 的值。
定义
这样定义的话，很多性质的证明就能用数学归纳法，而且可以简化某些计算，比最开始的定义更加实用。
然而直接这样暴力算的话，依然是 (O(n!)) 的。能不能继续优化呢？
当然可以。不过我们得先知道一些行列式的性质。
以下性质，读者自证不难。
性质 1 行列式和行与列互换，其值不变。即
性质 2 行列式对任意一行进行如下展开，其值不变。
性质 3 （线性性质）有以下两条
推论 1 某行全为 (0) 的行列式值为 (0) 。
性质 4 行列式中，两行对应元素全相等，其值为 (0) 。
推论 2 行列式中，两行对应元素成比例，其值为 (0) 。
性质 5 行列式中，把某行的元素分别乘常数 (k) 加到另一行的对应元素上，其值不变。
性质 6 行列式中，对换两行，行列式的值反号。
其中性质 5， 6 是最重要的性质。
特别的，一个下三角行列式
的值为 (prodlimits_{i=1}^na_{ii}) 。直接按定义展开即可。上三角同理。
这样，根据性质 5， 6 ，我们可以用高斯消元将行列式消至上（下）三角，即可求出来行列式的值。我们就有了一个 (O(n^3)) 计算行列式的值的算法。
当然行列式的性质还有许许多多，我们现在还用不到。而计算行列式的值，本来也是用不到的，只是接下来的某个算法需要用到而已。
4. 矩阵树定理
用来计算有标号无向图生成树个数的利器。
是由基尔霍夫(Kirchhoff)证明的，我原来还奇怪为什么完全搜不到这位基尔霍夫的资料，而只能搜到一个物理学家。后来才知道，他们就是一个人。物理学家基尔霍夫在研究电路的时候，顺便搞出来了这个东西。而基尔霍夫定律是研究电路的基本定律之一，大概高考会讲吧。大家日后见到它的时候应该会有熟悉感吧，因为其中一个电流定律 KCL 完全就是网络流的流量守恒，我之前还用这玩意做过物理题。（显然是不怎么优的解法）。
内容：
对于一个无向图 (G) ，我们有它的邻接矩阵 (A) 和度数矩阵 (D) ,邻接矩阵就是若 (i, j) 有边则 (a_{ij} = 1)，度数矩阵就是 (d_{ii}) 等于 (i) 的度数。
然后我们定义 (G) 的基尔霍夫矩阵 (C = D - A) 。(G) 的生成树个数即为 (C) 的任意一个 ((n-1)) 阶主子式的行列式的值。
矩阵的 ((n-1)) 阶主子式与行列式余子式类似，就是去掉这个矩阵的某一行和某一列剩下的部分。
注意主子式和余子式并不是差不多的东西。只是((n-1)) 阶主子式与余子式长得差不多而已。
顺带提一下有向图的矩阵树定理。
如果要求以 (i) 为根的外向生成树的数量，就是出度矩阵 - 邻接矩阵，然后去掉第 (i) 行 (i) 列求行列式。
内向就入度矩阵 - 邻接矩阵。
正确性？读者自证不难
然后我们就有一个 (O(n^3)) 求图的生成树个数的算法啦！
可以直接拿到生成树计数的 (60) 分是不是很棒棒！
实际上那个题已经给了矩阵树定理和行列式的计算方法了，不过是那个 (O(n!n^2log n)) 的方法。
怎么和某 THUWC 题有异曲同工之妙。
然而好戏才刚刚开始。
5. 常系数齐次线性递推式第n项的快速计算
没有什么朗朗上口的算法名。
如果你的 DP 能写成 (k) 阶矩阵转移的形式，那就能 (O(k^3log n)) 求解了对不对？
看起来很棒棒，然而当转移矩阵一大就还不如暴力。
于是我们就需要这个黑科技！
它能将复杂度降到 (O(k^2log n)) 乃至 (O(klog k log n)) 。
于是有些毒瘤出题人就是要用这玩意最后恶心你一下。
比如NOI2017 泳池。
好在它的编程复杂度并没有多高，起码 (O(k^2log n)) 和矩阵快速幂一样好写，还不用写什么矩阵结构体，何乐而不为呢？
首先它是用来干什么的呢？
顾名思义，它是用来快速计算常系数齐次线性递推式第n项的。（逃
啥是递推式？就是
啥是常系数？递推式中 (a_i) 与 (n) 无关。
啥是线性？递推式中没有乱七八糟的二次三次项。
啥是齐次？递推式没有常数项。
5.1 暴力
这玩意很显然能拿矩阵快速幂做吧。我们构造这样一个转移矩阵 (A)
然后愉快的 (O(k^3log n)) 啦！
最终出题人将数据范围定为 (k leq 3e4) 。
等会这 (O(k^2log n)) 也跑不过去吧！
没事反正我们还有 (O(klog k log n)) 。
先不管那个，我们先来优化一下。
5.2 线代魔法
不妨设我们最开始的向量为 (v) 。
首先我们并不关心这个矩阵的 (n) 次幂是什么样子的。然后我们把它求了出来。
然后我们也不关心这个向量转移 (n) 次之后整个是什么样子的。我们只想要第一项。然后我们整个求了出来。
于是我们浪费了许多资源在无关紧要的地方。
优化的基本策略不就是不求我们不想要的东西嘛？
假设我们构造了一个奥妙重重的序列 (c) 使得
然后我们就可以愉快的 (O(k^4)) 的计算 (A^n) 啦！
嗯好像浪费了更多的资源。
冷静一下，一定是哪里出了问题。
我们不是只需要求 (A^n v) 吗？
根据矩阵乘法的分配率，我们得到
我们好像不需要整个向量，我们只关心第一个值。于是
然后，((A^iv)_0) 是啥啊？转移了 (i) 次的 (v) 的第 (0) 项。。。不就是 (v_i) 嘛？
于是
(O(k)) 。完美。
5.3 多项式科技
那么问题来了，怎么构造 (c) 呢？
回顾一下这个式子
左边，(n) 次。右边，(k-1) 次。
喵喵喵？
这能相等？
次数都不一样诶！
一定是发生了什么奥妙重重的事情。
右边原来肯定是一个 (n) 次的多项式，只是因为某些不明的原因，高次项系数全都成 (0) 了。
于是我们先复原一下这个式子，令 (c) 的多项式为 (C(x)) ，然后将它写成
其中 (Q(A)) 和 (G(A)) 是两个玄妙的多项式，满足 (Q(A)G(A) = 0) 。
为啥这样构造呢？
巧了，我们恰好能找出来这样一个次数为 (k) 的多项式使得 (G(A) = 0) 。怎么找待会儿讲。
然后就是一个多项式取模了。
然后我们构造一个多项式 (F(x) = x) ，快速幂加多项式取模，就可以愉快的算出来 (C(x)) 啦！
总时间复杂度 (O(klog k log n)) 。然而常数爆炸。你都全家桶了还在意啥常数
5.4 Cayley-Hamilton 定理
5.4.1 特征多项式
设 (A) 是 (n) 阶矩阵，若数 (lambda) 和非零列向量 (x) 使关系式
成立，那么这样的数 (lambda) 称为矩阵 $ A $ 的特征值，非零向量 (x) 称为 (A) 的对应于 (lambda) 的特征向量。
然后我们给它移一下项，就变成了
其中 (I) 为单位矩阵。注意 (lambda) 是个标量，得乘一个单位矩阵才能和 (A) 相加减。
存在 (x ot= vec 0) 的充要条件是 (det (lambda I - A) = 0) 。读者自证不难。
(f(lambda) = det (lambda I - A)) 又被称为 (A) 的特征多项式，其零点即为特征值。
5.4.2 Cayley-Hamilton 定理
若 (f(lambda)) 是 (A) 的特征多项式，则 (f(A) = 0) 。
然而证明似乎不能直接将 (A) 代进去。读者自证不难。
于是我们不就构造出一个 (G(x)) 了嘛！
桥豆麻袋，这行列式怎么算啊！
直接暴力算是 (O(k^3)) 的。显然不行。
这个时候我们需要用到数学题的思想：用手玩来降低时间复杂度。
我们先看看这个行列式长什么样子。
然后我们给第一行按照定义展开，就得到了一些近似于下三角的东西。
稍微手玩一下，就能得到
然后我们就真的只需要 std::reverse 一下递推式系数就能得到 (G(x)) 啦！
于是问题圆满解决。
5.4.3 实现
直接快速幂加多项式取模就行。
然而大部分时间没必要用多项式，(k^2) 暴力就够了。
可是这暴力怎么写呢？
现在我们将两个多项式 (f) 和 (g) 卷了起来，应该是一个 (2k-1) 次的多项式 (h) 。我们要对它进行取模，使其次数降为 (k-1) 。注意到 (G(x)) 的最高次项的系数为 (1) ，我们只需要从 (2k-1) 循环到 (k) ，然后将 (G(x)) 的最高项与之对齐，乘上当前项的系数，一减即可。而且由于 (G(x)) 中除了最高次项之外其他项的系数的符号都是负，我们甚至不需要预处理 (G(x)) ，直接用给的转移系数做就行。
这个过程就是在模拟手动取模的过程，不理解的可以手玩一下。记得处理 (k = 1) 的情况。
6. Berlekamp-Massey 算法
有的时候我们并不能很轻易地推出来递推系数。就是打表打不出来的那种。
那我们能不能通过某个算法让计算机给我们找规律呢？
还真有。这就是 Berlekamp-Massey 算法。
不过这个算法局限性也是巨大的，应用范围非常小。
好就好在十分好写，见题直接先来一个素质二连看看能不能打出表来。
Berlekamp-Massey 算法，简称 BM 算法，是用来求解一个数列的最短线性递推式的算法。
BM 算法可以在 (O(k^2)) 的时间内求解一个长度为 (k) 的数列的最短线性递推式。
这个算法毛爷爷的2017论文讲了，不过那篇论文基本上不能看吧。。。
实际上 BM 算法的过程还是很简单的。
对于数列 ({a_1,a_2,dots,a_n}) ，我们称数列 ({r1, r2,dots,r_k}) 为其线性递推式当且仅当
对于任意 (k+1 leq i leq n) 成立。
若数列 ({a_1,a_2,dots,a_n}) 的线性递推式 ({r1, r2,dots,r_k}) 还满足 (k) 是所有该数列的线性递推式中最小的，则称 ({r1, r2,dots,r_k}) 为 ({a_1,a_2,dots,a_n}) 的最短线性递推式。
我们定义一个递推式代入一个数列的第 (i) 项的结果为
我们现在想要求递推式。怎么求呢？我们开始随便口胡一个递推式，就 ({}) 好了。
我们一个个地将数列里的数代入递推式，直到 (Delta_i ot =0) 。说明我们找到了数列中第一个非零的值。
这个时候我们能做些什么呢？我们好像什么都做不了。那就给递推式填上 (i) 个 (0) 吧！
然后我们继续代入，直到下一个 (Delta_i ot =0) 的位置。这个时候我们就不能草率地全填 (0) 了。毕竟最后得到一个和原数组等长的递推式并不是我们想要的结果。
于是我们得考虑修复这个递推式。
我们现在有递推式 (g) 对于 (1, 2, dots,i-1) 都成立，然后挂在了 第 (i) 项。如果我们能凭空变出一个递推式 (f) 使得它代进前 (i-1) 项都是 (0)，代进第 (i) 项是 (1) ，那么我们直接将 (g) 减去 (Delta_i f) 不就行了？
可是从哪儿去找这么一个奥妙重重的 (f) 呢？
根据字符串那一套理论，失败即是经验。按照以上的过程，我们一定不是第一次遇到这种情况。之前肯定也有一个递推式 (f) ，它满足代入前 (j-1) 项的结果都为 (0) ，而代入第 (j) 项的结果为 $ Delta_j ot = 0 $ 。我们将 (f) 和当前递推式 (g) 对齐，并让 (g) 减去 (f) ，我们就可以成功修补 (f) 啦！
具体的，设 (f) 第一次代入出了偏差的位置为 (j) ，那么我们就有
设 (mul = frac {Delta_i} {Delta_j}) ，我们将 (f) 的前面填入 (i - j - 1) 个 (0) ，然后填入 (-mul) ，和 (f) 的各项乘 (mul) 的结果。于是它代入前 (i-1) 项的结果都是 (0) ，而代入第 (i) 项的结果就是 (Delta_i) 。然后一减就好啦！
而我们有多个 (f) ，选哪一个呢？
选最近的是假的，这样跑不出来最短递推式，虽然并没有什么区别的样子。
我们可以维护一下当前最优的一个递推式，使得它对齐之后的长度最少，随时更新一下就行。
跑完一遍 BM 出来的就是递推式，搭配常系数齐次线性递推式第n项的快速计算我们就可以闷声 (O(k^2log n)) 啦！是不是很愉快？
由于每一项都可能会重新计算递推式，于是复杂度就是 (O(k^2)) 的，(k) 是你扔进去的数的数量。不过肯定远远跑不满就是了。
7. 实现
虽然一套板子下来并没有多少，但是如同 Dinic 一样，这个代码细节还是很多的，并不是太好背。
虽然比 Dinic还要黑盒一些，但是一不留神写假了，就会跑出来各种奇奇怪怪的结果。有的时候是与答案错开几项的数，有的时候是完全不对的数。
我们先来讲一下如何素质三连碾这个题。生成树计数
当时这个题出出来的时候，人类科技并不是很发达。于是出题人像模像样地给了矩阵树定理，企图诱导选手们想歪。会点线代的人大概就打个 (60) 分跑了吧。结果正解是个状压 DP 。
然而 (12) 年过去了，人类的科技有了突飞猛进的发展。
直接矩阵树定理跑出来 (N leq 100) 的点，然后扔 BM 里跑出一个递推式，再用常系数齐次线性递推直接出解。其实也用不到常系数齐次线性递推，矩阵快速幂能跑出来，递推式只有 (46) 阶。不过常系数齐次线性递推真的比矩阵快速幂好写。
粘个代码理解一下细节。
#include <cstdio>
#include <cstring>
#include <vector>
#include <cmath>
const int MOD = 65521;
inline int Qpow(int a, int p) {
int ret = 1;
while (p) {
if (p & 1) ret = 1ll * ret * a % MOD;
a = 1ll * a * a % MOD; p >>= 1;
}
return ret;
}
const int MAXN = 210;
namespace Linear_Seq {
std::vector<int> BM(const std::vector<int>& f) {
std::vector<int> best, last;
int lf = 0, ld = 0;
for (size_t i = 1; i < f.size(); ++i) {
int dt = MOD - f[i];
for (size_t j = 0; j < last.size(); ++j)
dt = (dt + 1ll*f[i-j-1]*last[j]%MOD) % MOD;
if (!dt) continue;
if (last.empty()) {
last.resize(i), lf = i, ld = dt;
continue;
}
int k = 1ll*(MOD-dt)*Qpow(ld, MOD-2) % MOD;
std::vector<int> cur(i - lf - 1);
cur.push_back(MOD-k);
for (size_t j = 0; j < best.size(); ++j)
cur.push_back(1ll*best[j]*k%MOD);
if (cur.size() < last.size()) cur.resize(last.size());
for (size_t j = 0; j < last.size(); ++j)
(cur[j] += last[j]) %= MOD;
if (cur.size() <= i-lf+best.size())
best = last, lf = i, ld = dt;
last = cur;
}
return last;
}
void mul(int* a, int* b, const std::vector<int>& g) {
static int tmp[MAXN];
int n = g.size();
for (int i = 0; i < (n << 1); ++i) tmp[i] = 0;
for (int i = 0; i < n; ++i) if (a[i])
for (int j = 0; j < n; ++j)
tmp[i+j] = (tmp[i+j]+1ll*a[i]*b[j]%MOD)%MOD;
for (int i = 2*n-1; i >= n; --i) if (tmp[i])
for (int j = n-1; ~j; --j)
tmp[i-j] = (tmp[i-j]+1ll*tmp[i]*g[j]%MOD)%MOD;
for (int i = 0; i < n; ++i) a[i] = tmp[i];
}
int a[MAXN], b[MAXN];
int Calc(const std::vector<int>& f, long long n) {
if (1u*n < f.size()) return f[n];
std::vector<int> v = BM(f);
v.insert(v.begin(), 0);
a[0] = 1, b[1] = 1;
while (n) {
if (n & 1) mul(a, b, v);
mul(b, b, v); n >>= 1;
}
int ans = 0;
for (size_t i = 0; i < v.size(); ++i)
ans = (ans + 1ll*a[i]*f[i]%MOD) % MOD;
return ans;
}
}
int G[MAXN][MAXN];
int Gauss_Det(int n) {
int cnt = 0;
for (int i = 1; i <= n; ++i) {
int p = i;
for (int j = i; j <= n; ++j) if (G[j][i]) { p = j; break; }
if (p != i) {
++cnt;
for (int j = 1; j <= n; ++j) std::swap(G[i][j], G[p][j]);
}
int inv = Qpow(G[i][i], MOD-2);
for (int j = i+1; j <= n; ++j) {
int rate = 1ll * G[j][i] * inv % MOD;
for (int k = 1; k <= n; ++k)
G[j][k] = (G[j][k]+MOD-1ll*rate*G[i][k]%MOD) % MOD;
}
}
int ans = cnt ? MOD-1 : 1;
for (int i = 1; i <= n; ++i) ans = 1ll * ans * G[i][i] % MOD;
return ans;
}
int K;
long long N;
std::vector<int> f(1);
int main() {
scanf("%d%lld", &K, &N);
for (int n = 1; n <= 100; ++n) {
memset(G, 0, sizeof G);
for (int i = 1; i <= n; ++i)
for (int j = std::max(1, i-K); j <= n && j <= i + K; ++j)
if (!G[i][j])
G[i][i]++, G[j][j]++, G[i][j]--, G[j][i]--;
for (int i = 1; i <= n; ++i)
for (int j = 1; j <= n; ++j) {
if (G[i][j] < 0) G[i][j] += MOD;
}
f.push_back(Gauss_Det(n-1));
}
printf("%d
", Linear_Seq::Calc(f, N));
}
于是我们就可以理性打表了！
8. 尾声
然而 2019 年集训队论文，神仙 zzq 搞了一些更牛逼的递推法。感兴趣的同学可以去了解一下。