0x33同 余

定义

若整数 $a$ 和整数 $b$ 除以正整数 $m$ 的余数相等，则称 $a, b$ 模 $m$ 同余，记为 $a \equiv b \pmod{m}$ 。

同余类与剩余系

对于 $\forall a\in [0,m - 1]$ ，集合 $\{a + km\} (k\in \mathbb{Z})$ 的所有数模 $m$ 同余，余数都是 $a$ 。该集合称为一个模 $m$ 的同余类，简记为 $\overline{a}$ 。

模 $m$ 的同余类共有 $m - 1$ 个，分别为 $\overline{0},\overline{1},\overline{2},\dots ,\overline{m - 1}$ 。它们构成 $m$ 的完全剩余系。

$1 \sim m$ 中与 $m$ 互质的数代表的同余类共有 $\varphi(m)$ 个，它们构成 $m$ 的简化剩余系。例如，模 10 的简化剩余系为 $\{\overline{1}, \overline{3}, \overline{7}, \overline{9}\}$ 。

简化剩余系关于模 $m$ 乘法封闭。这是因为若 $a, b (1 \leq a, b \leq m)$ 与 $m$ 互质，则 $a * b$ 也不可能与 $m$ 含有相同的质因子，即 $a * b$ 也与 $m$ 互质。再由余数的定义即

可得到 $a * b \bmod m$ 也与 $m$ 互质，即 $a * b \bmod m$ 也属于 $m$ 的简化剩余系。

费马小定理

若 $p$ 是质数，则对于任意整数 $a$ ，有 $a^p \equiv a \pmod{p}$ 。

欧拉定理

若正整数 $a, n$ 互质，则 $a^{\varphi(n)} \equiv 1 \pmod{n}$ ，其中 $\varphi(n)$ 为欧拉函数。

证明：

设 $n$ 的简化剩余系为 $\{\overline{a_1},\overline{a_2},\dots ,\overline{a_{\varphi(n)}}\}$ 。对 $\forall a_i,a_j$ ，若 $a*a_{i}\equiv a*a_{j}(\mathrm{mod} n)$ 则 $a*(a_i - a_j)\equiv 0$ 。因为 $a,n$ 互质，所以 $a_{i} - a_{j}\equiv 0$ ，即 $a_{i}\equiv a_{j}$ 。故当 $a_{i}\neq a_{j}$ 时， $aa_{i},aa_{j}$ 也代表不同的同余类。

又因为简化剩余系关于模 $n$ 乘法封闭，故 $\overline{aa_i}$ 也在简化剩余系集合中。因此，集合 $\{\overline{a_1},\overline{a_2},\dots ,\overline{a_{\varphi(n)}}\}$ 与集合 $\{\overline{aa_1},\overline{aa_2},\dots ,\overline{aa_{\varphi(n)}}\}$ 都能表示 $n$ 的简化剩余系。综上所述：

因此 $a^{\varphi (n)}\equiv 1$ (mod $n$ )。

当 $p$ 是质数时， $\varphi(p) = p - 1$ ，并且只有 $p$ 的倍数与 $p$ 不互质。所以，只要 $a$ 不是 $p$ 的倍数，就有 $a^{p - 1} \equiv 1 \pmod{p}$ ，两边同乘 $p$ 就是费马小定理。若 $a$ 是 $p$ 的倍数，费马小定理显然成立。

证毕。

欧拉定理的推论

若正整数 $a, n$ 互质，则对于任意正整数 $b$ ，有 $a^b \equiv a^b \mod \varphi(n) \pmod{n}$ 。

证明：

设 $b = q*\varphi (n) + r$ ，其中 $0\leq r < \varphi (n)$ ，即 $r = b\bmod \varphi (n)$ 。于是：

证毕。

许多计数类的题目要求我们把答案对一个质数 $p$ 取模后输出。面对 $a + b, a - b, a * b$ 这样的算式，可以在计算前先把 $a, b$ 对 $p$ 取模。面对乘方算式，根据欧拉定理的推论，可以先把底数对 $p$ 取模、指数对 $\varphi(p)$ 取模，再计算乘方。

【例题】The Luckiest Number POJ3696

给定一个正整数 $L$ ， $L\leq 2*10^{9}$

问至少多少个8连在一起组成的正整数是 $L$ 的倍数？

$x$ 个8连在一起组成的正整数可写作 $8(10^{x} - 1) / 9$ 。题目就是让我们求出一个最小的 $x$ ，满足 $L\mid 8(10^{x} - 1) / 9$ 。设 $d = \gcd (L,8)$ 。

引理：

若正整数 $a, n$ 互质，则满足 $a^x \equiv 1 \pmod{n}$ 的最小正整数 $x_0$ 是 $\varphi(n)$ 的约数。

证明：

反证法。假设满足 $a^x \equiv 1 \pmod{n}$ 的最小正整数 $x_0$ 不能整除 $\varphi(n)$ 。

设 $\varphi(n) = qx_0 + r$ ( $0 < r < x_0$ )。因为 $a^{x_0} \equiv 1 \pmod{n}$ ，所以 $a^{qx_0} \equiv 1 \pmod{n}$ 。根据欧拉定理，有 $a^{\varphi(n)} \equiv 1 \pmod{n}$ ，所以 $a^r \equiv 1 \pmod{n}$ 。这与 $x_0$ 最小矛盾。故假设不成立，原命题成立。

证毕。

根据以上引理，我们只需求出欧拉函数 $\varphi(9L / d)$ ，枚举它的所有约数，用快速幂逐一检查是否满足条件即可。时间复杂度为 $O\left(\sqrt{L} \log L\right)$ 。

扩展欧几里得算法

Bézout定理

对于任意整数 $a, b$ ，存在一对整数 $x, y$ ，满足 $ax + by = \gcd(a, b)$ 。

证明：

在欧几里得算法的最后一步，即 $b = 0$ 时，显然有一对整数 $x = 1, y = 0$ ，使得 $a * 1 + 0 * 0 = \gcd(a, 0)$ 。

若 $b > 0$ ，则 $\gcd (a,b) = \gcd (b,a\bmod b)$ 。假设存在一对整数 $x,y$ ，满足 $b*x+(a\bmod b)*y = \gcd (b,a\bmod b)$ ，因为 $bx + (a\bmod b)y = bx + (a - b\lfloor a / b\rfloor)y = ay - b(x - \lfloor a / b\rfloor y)$ ，所以令 $x' = y, y' = x - \lfloor a / b\rfloor y$ ，就得到了 $ax' + by' = \gcd (a,b)$ 。

对欧几里得算法的递归过程应用数学归纳法，可知Bézout定理成立。

证毕。

Bézout 定理是按欧几里得算法的思路证明的，且上述证明同时给出了整数 $x$ 和 $y$ 的计算方法。这种计算方法被称为扩展欧几里得算法。

int exgcd(int a, int b, int &x, int &y) { if (b == 0) { x = 1, y = 0; return a; } }

int d = exgcd(b, a%b, x, y);
int z = x; x = y; y = z - y * (a / b);
return d;
}

定义变量 $d, x_0, y_0$ ，调用 $\mathrm{d} = \mathrm{exgcd}(a, b, x_0, y_0)$ 。注意在上述代码中， $x_0, y_0$ 是以引用的方式传递的。上述程序求出方程 $ax + by = \gcd(a, b)$ 的一组特解 $x_0, y_0$ ，并返回 $a, b$ 的最大公约数 $d$ 。

对于更为一般的方程 $ax + by = c$ ，它有解当且仅当 $d|c$ 。我们可以先求出 $ax + by = d$ 的一组特解 $x_0, y_0$ ，然后令 $x_0, y_0$ 同时乘上 $c / d$ ，就得到了 $ax + by = c$ 的一组特解 $(c / d)x_0, (c / d)y_0$ 。

事实上，方程 $ax + by = c$ 的通解可以表示为：

其中 $k$ 取遍整数集合， $d = \gcd(a, b)$ ， $x_0, y_0$ 是 $ax + by = \gcd(a, b)$ 的一组特解。乘法逆元

若整数 $b, m$ 互质，并且 $b \mid a$ ，则存在一个整数 $x$ ，使得 $a / b \equiv a * x \pmod{m}$ 。称 $x$ 为 $b$ 的模 $m$ 乘法逆元，记为 $b^{-1} \pmod{m}$ 。

因为 $a / b\equiv a*b^{-1}\equiv a / b*b*b^{-1}(\mathrm{mod} m)$ ，所以 $b*b^{-1}\equiv 1(\mathrm{mod}m)$ 。

如果 $m$ 是质数（此时我们用符号 $p$ 代替 $m$ ）并且 $b < p$ ，根据费马小定理， $b^{p - 1} \equiv 1 \pmod{p}$ ，即 $b * b^{p - 2} \equiv 1 \pmod{p}$ 。因此，当模数 $p$ 为质数时， $b^{p - 2}$ 即为 $b$ 的乘法逆元。

如果只是保证 $b, m$ 互质，那么乘法逆元可通过求解同余方程 $b * x \equiv 1 \pmod{m}$ 得到。下个部分我们就来介绍线性同余方程及其求解方法。

有了乘法逆元，我们在计数类问题中即使遇到 $a / b$ 这样的除法算式，也可以先把 $a, b$ 各自对模数 $p$ 取模，再计算 $a * b^{-1} \bmod p$ 作为最终的结果。当然，前提是必须保证 $b, p$ 互质（当 $p$ 是质数时，等价于 $b$ 不是 $p$ 的倍数）。

到此为止，我们在模 $p$ 运算下对加、减、乘、除、乘方运算都已经有了适当的处理方法，请读者回顾并加以思考。

【例题】Sumdiv POJ1845

求 $A^B$ 的所有约数之和 mod 9901 $(1 \leq A, B \leq 5 * 10^7)$ 。

在 $0 \times 03$ 节中，我们已经用递归和分治算法解决了取模乘法下等比数列的求和问题。在本题中我们将直接借助乘法逆元计算等比数列的求和公式。

把 $A$ 分解质因数，表示为 $p_1^{c_1} * p_2^{c_2} * \dots * p_n^{c_n}$ ，由算术基本定理的“约数和”推论， $A^B$ 的所有约数之和为： $\left(1 + p_1 + p_1^2 + \dots + p_1^{B*c_1}\right) * \left(1 + p_2 + p_2^2 + \dots + p_2^{B*c_2}\right) * \dots * \left(1 + p_n + p_n^2 + \dots + p_n^{B*c_n}\right)$ 。

上式的每一项都是一个等比数列。以第一项为例， $1 + p_1 + p_1^2 + \dots + p_1^{B*c_1} = (p_1^{B*c_1 + 1} - 1) / (p_1 - 1)$ 。可以先用快速幂计算分子 $\left(p_1^{B*c_1 + 1} - 1\right) \mod 9901$ 和分母 $(p_1 - 1) \mod 9901$ 。因为9901是质数，只要 $p_1 - 1$ 不是9901的倍数，就只需计算 $p_1 - 1$ 的乘法逆元 $inv$ ，用乘 $inv$ 代替除以 $(p_1 - 1)$ ，直接计算出等比数列求和公式的结果。

特别地，若 $p_1 - 1$ 是9901的倍数，此时乘法逆元不存在，但是 $p_1 \mod 9901 = 1$ ，所以 $1 + p_1 + p_1^2 + \dots + p_1^{B*c_1} \equiv 1 + 1 + 1^2 + \dots + 1^{B*c_1} \equiv B*c_1 + 1 \pmod{9901}$ 。

int a, b, m, ans = 1, mod = 9901;  
int p[20], c[20];  
void divide(int n) {  
    m = 0;  
    for (int i = 2; i * i <= n; i++) {  
        if (n % i == 0) {  
            p[++m] = i, c[m] = 0;  
        }  
    }  
    if (n > 1)  
        p[++m] = n, c[m] = 1;  
}  
int power(int a, long long b) {  
    int c = 1;  
    for (; b; b >= 1) {  
        if (b & 1) c = (long long)c * a % mod;  
        a = (long long)a * a % mod;  
    }  
    return c;  
}  
int main() {  
    cin >> a >> b;

divide(a); //分解质因数  
for (int i = 1; i <= m; i++) {  
    if ((p[i] - 1) % mod == 0) { //没有逆元时，特判  
        ans = ((long long)b * c[i] + 1) % mod * ans % mod;  
        continue;  
    }  
    int x = power(p[i], (long long)b * c[i] + 1); //分子  
    x = (x - 1 + mod) % mod;  
    int y = p[i] - 1; //分母  
    y = power(y, mod - 2); //根据费马小定理求乘法逆元  
    ans = (long long)ans * x % mod * y % mod;  
}  
cout << ans << endl;

$\spadesuit$ 线性同余方程

给定整数 $a, b, m$ ，求一个整数 $x$ 满足 $a * x \equiv b \pmod{m}$ ，或者给出无解。因为未知数的指数为1，所以我们称之为一次同余方程，也称线性同余方程。

$a * x \equiv b \pmod{m}$ 等价于 $a * x - b$ 是 $m$ 的倍数，不妨设为 $-y$ 倍。于是，该方程可以改写为 $a * x + m * y = b$ 。

根据Bézout定理及其证明过程，线性同余方程有解当且仅当 $\gcd (a,m)|b$ 。

在有解时，先用欧几里得算法求出一组整数 $x_0, y_0$ ，满足 $a * x_0 + m * y_0 = \gcd(a, m)$ 。然后 $x = x_0 * b / \gcd(a, m)$ 就是原线性同余方程的一个解。

方程的通解则是所有模 $m / \gcd (a,m)$ 与 $x$ 同余的整数。

【例题】同余方程 NOIP2012/CODEVS1200

求关于 $x$ 的同余方程 $a*x \equiv 1 \pmod{b}$ 的最小正整数解。输入数据保证一定有解。

由上述线性同余方程的知识可得， $a*x\equiv 1(\mathrm{mod}b)$ 有解当且仅当 $\gcd (a,b) = 1$ 。方程可改写为 $a*x + b*y = 1$ ，用欧几里得算法求出一组特解 $x_0,y_0$ ，则 $x_0$ 就是原方程的一个解，通解为所有模 $b$ 与 $x_0$ 同余的整数。通过取模操作把解的范围移动到 $1\sim b$ 之间，就得到了最小正整数解。

long long a, b, x, y;
long long exgcd(long long a, long long b,

long long &x, long long &y) { if (!b) { x=1, y=0; return a; } long long d = exgcd(b, a%b, x, y); long long z = x; x = y; y = z - y*(a/b); return d; } int main() { cin >> a >> b; exgcd(a, b, x, y); cout << (x % b + b) % b << endl; }

中国剩余定理

设 $m_{1}, m_{2}, \dots, m_{n}$ 是两两互质的整数， $m = \prod_{i=1}^{n} m_{i}$ ， $M_{i} = \frac{m}{m_{i}}$ ， $t_{i}$ 是线性同余方程 $M_{i} t_{i} \equiv 1 \pmod{m_{i}}$ 的一个解。对于任意的 $n$ 个整数 $a_{1}, a_{2}, \dots, a_{n}$ ，方程组

有整数解，解为 $x = \sum_{i=1}^{n} a_i M_i t_i$ 。

证明：

因为 $M_{i} = m / m_{i}$ 是除 $m_{i}$ 之外所有模数的倍数，所以 $\forall k\neq i,a_{i}M_{i}t_{i}\equiv$ 0 (mod $m_{k})$ 。又因为 $a_{i}M_{i}t_{i}\equiv a_{i}$ (mod $m_{i})$ ，所以代入 $x = \sum_{i = 1}^{n}a_{i}M_{i}t_{i}$ ，原方程组成立。

证毕。

中国剩余定理给出了模数两两互质的线性同余方程组的一个特解。方程组的通解可以表示为 $x + km (k \in \mathbb{Z})$ 。有些题目要求我们求出最小的非负整数解，只需把 $x$ 对 $m$ 取模，并让 $x$ 落在 $0 \sim m - 1$ 的范围内即可。

另外，即使模数不满足两两互质，我们也有方法判断线性同余方程组是否有解，并求出方程组的解。

【例题】Strange Way to Express Integers POJ2891

给定 $2n$ 个正整数 $a_1, a_2, \dots, a_n$ 和 $m_1, m_2, \dots, m_n$ ，求一个最小的正整数 $x$ ，满足 $\forall i \in [1, n]$ ， $x \equiv a_i \pmod{m_i}$ ，或者给出无解。

本题中 $m_{i}$ 不一定两两互质，中国剩余定理不再适用。可以考虑使用数学归纳法，

假设已经求出了前 $k - 1$ 个方程构成的方程组的一个解 $x$ 。记 $m = \prod_{i=1}^{k-1} m_i$ ，则 $x + i * m (i \in \mathbb{Z})$ 是前 $k - 1$ 个方程的通解。

考虑第 $k$ 个方程，求出一个整数 $t$ ，使得 $x + t * m \equiv a_k \pmod{m_k}$ 。该方程等价于 $m * t \equiv a_k - x \pmod{m_k}$ ，其中 $t$ 是未知量。这就是一个线性同余方程，可以用扩展欧几里得算法判断是否有解，并求出它的解。若有解，则 $x' = x + t * m$ 就是前 $k$ 个方程构成的方程组的一个解。

综上所述，我们使用 $n$ 次扩展欧几里得算法，就求出了整个方程组的解。

$\spadesuit$ 高次同余方程

关于高次同余方程，有 $a^x \equiv b (\bmod p)$ 和 $x^a \equiv b (\bmod p)$ 两类问题。不过后者超出了我们的探讨范围，学有余力的读者可以查阅“原根”“阶”“指标”的相关资料。我们重点来解决前者。

问题：

给定整数 $a, b, p$ ，其中 $a, p$ 互质，求一个非负整数 $x$ ，使得 $a^x \equiv b (\bmod p)$ 。

Baby Step, Giant Step 算法

因为 $a, p$ 互质，所以可以在模 $p$ 意义下执行关于 $a$ 的乘、除法运算。

设 $x = i * t - j$ ，其中 $t = \left\lceil \sqrt{p} \right\rceil, 0 \leq j \leq t - 1$ ，则方程变为 $a^{i * t - j} \equiv b (\bmod p)$ 。即 $(a^t)^i \equiv b * a^j (\bmod p)$ 。

对于所有的 $j \in [0, t - 1]$ ，把 $b * a^j \mod p$ 插入一个 Hash 表。

枚举 $i$ 的所有可能取值，即 $i \in [0, t]$ ，计算出 $(a^t)^i \bmod p$ ，在 Hash 表中查找是否存在对应的 $j$ ，更新答案即可。时间复杂度为 $O(\sqrt{p})$ 。

下面的参考程序实现了 BabyStep, GiantStep 算法，计算同余方程 $a^x \equiv b (\bmod p)$ 的最小非负整数解，无解时返回 -1。

int baby_step_giant_step(int a, int b, int p) {
    map<int, int> hash;
    hash.clear();
    b % = p;
    int t = (int)sqrt(p) + 1;
    for (int j = 0; j < t; j++) {
        int val = (long long)b * power(a, j, p) % p; // b*a^j
        hash[val] = j;
    }
    a = power(a, t, p); // a^t

if  $(a == 0)$  return  $b == 0$  ? 1 : -1;  
for (int i = 0; i <= t; i++) {  
    int val = power(a, i, p); // (a^t)^i  
    int j = hash.find(val) == hash.end() ? -1 : hash[val];  
    if (j >= 0 && i * t - j >= 0) return i * t - j;  
}  
return -1;

$0 \times 34$ 矩阵乘法

一个 $n*m$ 的矩阵可看作一个 $n*m$ 的二维数组。矩阵的加法和减法就是把两个矩阵对应位置上的数相加减，即 $C = A + B \Leftrightarrow \forall i \in [1,n], \forall j \in [1,m], C_{i,j} = A_{i,j} \pm B_{i,j}$ ，其中 $A, B, C$ 都是 $n*m$ 矩阵。

矩阵乘法的定义略微复杂一些。设 $A$ 是 $n*m$ 矩阵， $B$ 是 $m*p$ 矩阵，则 $C = A*B$ 是 $n*p$ 矩阵，并且 $\forall i \in [1,n], \forall j \in [1,p]$ ：

也就是说，参与矩阵乘法运算的第一个矩阵的列数必须等于第二个矩阵的行数，所得到的结果矩阵的行数、列数分别等于第一个矩阵的行数、第二个矩阵的列数。结果矩阵 $C$ 第 $i$ 行第 $j$ 列的数，是由 $A$ 的第 $i$ 行的 $m$ 个数与 $B$ 的第 $j$ 列的 $m$ 个数分别相乘再相加得到的。

矩阵乘法满足结合律，即 $(A*B)*C = A*(B*C)$ ，满足分配律，即 $(A + B)*C = A*C + B*C$ ，不满足交换律，即 $A*B$ 不一定等于 $B*A$ 。

考虑矩阵乘法中一种特殊的情形， $F$ 是 $1*n$ 矩阵， $A$ 是 $n*n$ 矩阵，则 $F' = F*A$ 也是 $1*n$ 矩阵。 $F$ 和 $F'$ 可看作一维数组，省略它们的行下标 1。按照矩阵乘法的定义， $\forall j \in [1,n]$ ，有 $F'_j = \sum_{k=1}^{n} F_k * A_{k,j}$ 。它的含义是，经过一次与 $A$ 的矩阵乘法运算后， $F$ 数组中的第 $k$ 个值会以 $A_{k,j}$ 为系数累加到 $F'$ 数组的第 $j$ 个值上，等价于在一个向量的 $k,j$ 两个状态之间发生了递推。

【例题】Fibonacci POJ3070

在斐波那契数列中， $Fib_{0} = 0, FIB_{1} = 1, Fib_{n} = Fib_{n - 1} + Fib_{n - 2}(n > 1)$ 。现给定整数 $n,m$ $(0\leq n\leq 2*10^{9},m = 10000)$ ，求 $Fib_{n}\mod m_{\circ}$

直接循环递推计算斐波那契数列，时间复杂度显然为 $O(n)$ 。不过， $Fib_{n}$ 只与

$F i b_{n - 1}$ 和 $F i b_{n - 2}$ 有关，我们在递推时只需要保存最近的两个斐波那契数，即可得到下一个斐波那契数。

设 $F(n)$ 表示一个 $1*2$ 的矩阵， $F(n) = [Fib_{n} \quad Fib_{n+1}]$ 。

我们希望根据 $F(n - 1) = [Fib_{n - 1} \quad Fib_n]$ 计算出 $F(n)$ 。也就是说，要把 $F(n - 1)$ 第2列上的数作为 $F(n)$ 第1列上的数，并把 $F(n - 1)$ 第1、2列上的数都累加到 $F(n)$ 的第2列上。因此，我们令矩阵 $A$ 第2行第1列、第1行第2列、第2行第2列的数都是1，第1行第1列的数是0。

请读者根据矩阵乘法的定义对上式进行验证。

上式是一个矩阵乘法递推式，初值为 $F(0) = [0 \quad 1]$ ，目标为 $F(n) = F(0) * A^n$ 。因为矩阵乘法满足结合律，所以我们可以用快速幂计算 $F(0) * A^n$ ，得到的矩阵第1列上的数就是 $Fib_n$ 。算法的时间复杂度为 $O(2^3 \log n)$ ，其中 $2^3$ 是矩阵乘法所消耗的时间。

题目还要求我们对 $m$ 取模，可以直接在执行矩阵乘法时，对各个整数的加法、乘法运算加入取模操作。详细过程请参阅下面的参考程序：

const int mod = 10000;  
int n;  
void mul(int f[2], int a[2][2]) {  
    int c[2];  
    memset(c, 0, sizeof(c));  
    for (int j = 0; j < 2; j++)  
        for (int k = 0; k < 2; k++)  
            c[j] = (c[j] + (long long)f[k] * a[k][j]) % mod;  
        memcpy(f, c, sizeof(c));  
}  
void mulself(int a[2][2]) {  
    int c[2][2];  
    memset(c, 0, sizeof(c));  
    for (int i = 0; i < 2; i++)  
        for (int j = 0; j < 2; j++)  
            for (int k = 0; k < 2; k++)  
                c[i][j] = (c[i][j] + (long long)a[i][k]*a[k][j])%mod;  
       memcpy(a, c, sizeof(c));  
}  
int main() {

while (cin >> n && n != -1) {  
    int f[2] = {0,1};  
    int a[2][2] = {{0,1}, {1,1}};  
    for (; n; n >= 1) {  
        if (n & 1) mul(f, a);  
            mulself(a);  
    }  
    cout << f[0] << endl;  
}

这道题使我们对“矩阵乘法加速递推”有了初步的认识。一般来说，如果一类问题具有以下特点：

1. 可以抽象出一个长度为 $n$ 的一维向量，该向量在每个单位时间发生一次变化。
   2.变化的形式是一个线性递推（只有若干“加法”或“乘一个系数”的运算）。

2. 该递推式在每个时间可能作用于不同的数据上，但本身保持不变。

3. 向量变化时间（即递推的轮数）很长，但向量长度 $n$ 不大。

那么可以考虑采用矩阵乘法进行优化。我们把这个长度为 $n$ 的一维向量称为“状态矩阵”，把用于与“状态矩阵”相乘的固定不变的矩阵 $A$ 称为“转移矩阵”。若状态矩阵中的第 $x$ 个数对下一单位时间状态矩阵中的第 $y$ 个数产生影响，则把转移矩阵的第 $x$ 行第 $y$ 列赋值为适当的系数。

矩阵乘法加速递推的关键在于定义出“状态矩阵”，并根据递推式构造出正确的“转移矩阵”，之后就可以利用快速幂和矩阵乘法完成程序实现了。时间复杂度为 $O(n^{3}\log T)$ ，其中 $T$ 为递推总轮数。

【例题】石头游戏 BZOJ2973

石头游戏在一个 $n$ 行 $m$ 列 $(1 \leq n, m \leq 8)$ 的网格上进行，每个格子对应一种操作序列，操作序列至多有10种，分别用 $0 \sim 9$ 这10个数字指明。

操作序列是一个长度不超过6且循环执行、每秒执行一个字符的字符串。每秒钟，所有格子同时执行各自操作序列里的下一个字符。序列中的每个字符是以下格式之一：

1. 数字 $0 \sim 9$ ：表示拿 $0 \sim 9$ 个石头到该格子。

2. NWSE：表示把这个格子内所有的石头推到相邻的格子，N表示上方，W表示左方，S表示下方，E表示右方。
   3.D：表示拿走这个格子的所有石头。

给定每种操作序列对应的字符串，以及网格中每个格子对应的操作序列，求石头游戏进行了 $t$ 秒之后，石头最多的格子里有多少个石头。在游戏开始时，网格是空的。

设 $\forall i\in [1,n],j\in [1,m]$ ， $\mathrm{num}(i,j) = (i - 1)*m + j$ 。

把网格看作长度为 $n * m$ 的一维向量，定义1行 $n * m + 1$ 列的“状态矩阵” $F$ 下标为 $0 \sim n * m$ ，其中 $F[\mathrm{num}(i,j)]$ 记录格子 $(i,j)$ 中石头的个数。

特别地，令 $F[0]$ 始终等于1。

$F$ 会随着时间的增长而不断变化。设 $F_{k}$ 表示 $k$ 秒之后的“状态矩阵”。在游戏开始时，根据题意和 $F$ 的定义，有 $F_{0} = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \end{bmatrix}$ 。

注意到操作序列的长度不超过6，而 $1\sim 6$ 的最小公倍数是60，所以每经过60秒，所有操作序列都会重新处于最开始的字符处。我们可以统计出第 $k$ 秒 $(1\leq k\leq 60)$ 各个格子执行了什么字符，第 $k + 60$ 秒执行的字符与第 $k$ 秒一定是相同的。

对于 $1\sim 60$ 之间的每个 $k$ ，各个格子在第 $k$ 秒执行的操作字符可以构成一个“转移矩阵” $A_{k}$ ，矩阵行、列下标都是 $0\sim n*m$ 。构造方法如下：

1. 若网格 $(i,j)$ 第 $k$ 秒的操作字符为“N”，且 $i > 1$ ，则令 $A_{k}[\mathrm{num}(i,j),\mathrm{num}(i - 1,j)] = 1$ ，表示把石头推到上边的格子。字符“W”“S”“E”类似。

2. 若网格 $(i,j)$ 第 $k$ 秒的操作字符是一个数字 $x$ ，则令 $A_{k}[0, \mathrm{num}(i,j)] = x$ 。

3. 令 $A_{k}[0,0] = 1$ 。

4. $A_{k}$ 的其他部分均赋值为0。

上述构造实际上把“状态矩阵”的第0列作为“石头来源”。 $A_{k}[0,0] = 1$ 保证了 $F[0]$ 始终等于1。 $A_{k}[0,\mathrm{num}(i,j)] = x$ 相当于从“石头来源”获取 $x$ 个石头，放到格子 $(i,j)$ 中。使用矩阵乘法加速递推，遇到常数项时，经常需要在“状态矩阵”中添加一个额外的位置，始终存储常数1，并乘上“转移矩阵”中适当的系数，累加到“状态矩阵”的其他位置。

设 $A = \prod_{i=1}^{60} A_i$ , $t = q * 60 + r$ ( $0 \leq r < 60$ )。根据上面的讨论：

使用矩阵乘法和快速幂计算该式， $F_{t}$ 第 $1\sim n*m$ 列中的最大值即为所求。