1.2 矩阵快速幂

1.2.1 矩阵快速幂的概念

矩阵A和矩阵B相乘的前提条件是矩阵A的列数和矩阵B的行数相等，相乘的结果矩阵的行数等于矩阵A的行数，列数等于矩阵B的列数；即，设A为m×p的矩阵，B为p×n的矩阵，则m×n的矩阵C为矩阵A与B的乘积，记作C=A×B，其中矩阵C中的第i行第j列元素。

所以，只有方阵（行数和列数相等）才能进行矩阵的幂运算。

根据线性代数的知识，n×n的方阵A和方阵B相乘，产生n×n的方阵C，方阵A、B和C分别用二维数组a、b和c表示，将二维数组c初始化为0，方阵相乘的程序段如下：

对于n×n的方阵A，给出正整数M，求方阵A的M次幂AM。结合快速幂算法和矩阵相乘算法，程序段如下，其中，函数mul(A,B)为方阵A和B相乘，res为当前的位权矩阵，初始时，res=A。

求Fibonacci数取模是矩阵快速幂的一个应用。Fibonacci数列递推公式为：f(0)=0，f(1)=1，f(2)=1，f(n)=f(n-1)+f(n-2)，其中n≥3。则对于n≥3，Fibonacci数列递推公式可以通过矩阵运算得到：。因此，求Fibonacci数列的第n项，可以通过对方阵通过矩阵快速幂求解n-2次获得。

或者，对于n≥3，，而=，所以，，即，对方阵通过矩阵快速幂求解n次，产生方阵的第一行第二列的项或者第二行第一列的项，就是f(n)。

矩阵快速幂还可以求解形如f(n)=a×f(n-1)+b×f(n-2)+c的递推式，其中a、b、c为常数，初始值f(1)和f(2)已知，则有。

矩阵快速幂也可以求解形如f(n)=cn-f(n-1)的递推式，其中c为常数，初始值f(1)已知，则有。

【1.2.1.1 Fibonacci】

在Fibonacci整数序列中，F0=0，F1=1，且对n≥2，Fn=Fn-1+Fn-2。例如，Fibonacci序列的前十项为0, 1, 1, 2, 3, 5, 8, 13, 21, 34,…

Fibonacci序列的另一个公式是：

给出整数n，请计算Fn的最后4位。

输入

输入包含多个测试用例。每个测试用例一行，给出n(其中0≤n≤1000000000)。输入用包含数字-1的一行结束。

输出

对于每个测试用例，输出Fn的最后4位。如果Fn的最后4位都是0，则输出“0”；否则，忽略前导零（即，输出Fn mod 10000）。

提示

本题和矩阵相乘关联，两个2×2方阵相乘：

任何2×2方阵的0次幂是单位矩阵：。

试题来源：Stanford Local 2006

在线测试：POJ 3070

试题解析

简述题意：给出n，输出Fn mod 10000的结果。

显然，由Fibonacci序列公式可知，本题就是要求通过矩阵快速幂进行运算，矩阵的右上角的数字就是Fn。

参考程序

【1.2.1.2 Matrix Power Series】

给出一个n×n的方阵A和一个正整数k，请您计算S=A+A2+A3+…+Ak。

输入

输入只包含一个测试用例。输入的第一行给出3个正整数n(n≤30)、k(k≤109)和m(m<104)。接下来的n行每行包含n个小于32768的非负整数，这是按行排列的矩阵A的元素。

输出

以矩阵A的方式给出矩阵S，矩阵S中的元素除以m取余。

试题来源：POJ Monthly--2007.06.03

在线测试：POJ 3233

试题解析

本题采用矩阵快速幂和二分法求解。首先，可以通过矩阵快速幂求Ak；其次，可以对幂次k进行二分，每次将规模减半，分k为奇偶两种情况：

●当k为偶数时，S(k)=(1+Ak/2)(A+A2+A3+…+Ak/2)=(1+Ak/2)×S(k/2)；

●当k为奇数时，S(k)=A+(A+A(k+1)/2)(A+A2+A3+…+Ak/2)=A+(A+A(k+1)/2)×S(k/2)。

例如，S(6)=A+A2+A3+…+A6=(1+A3)(A+A2+A3)=(1+A3)×S(3)，S(7)=A+A2+A3+…+A7=A+(A+A4)(A+A2+A3)=A+(A+A4)×S(3)。

参考程序

1.2.2 矩阵快速幂的应用

本节给出快速矩阵幂和其他算法相结合解题的实验。题目1.2.2.1是将快速矩阵幂和加法原理相结合进行解题的实验，题目1.2.2.2是将快速矩阵幂和调整的Floyd-Warshall算法相结合进行解题的实验，而题目1.2.2.3则采用统计分析法，由局部解推导出矩阵的数学模型。

【1.2.2.1 Blocks】

Panda接到了一项任务，要给一排积木涂上颜色。Panda是一个聪明的孩子，他善于思考数学问题。假设一排积木有N块，每块积木可以涂成红色、蓝色、绿色或黄色。出于某些特殊的原因，Panda希望红色的积木和绿色的积木数量都是偶数。在这种情况下，Panda想知道有多少种不同的方法来给这些积木涂上颜色。

输入

输入的第一行给出整数T(1≤T≤100)，表示测试用例的数目。接下来的T行每行给出一个整数N(1≤N≤109)，表示积木的数目。

输出

对于每个测试用例，在一行输出给积木涂上颜色的方法数。因为答案可能很大，所以答案要对10007取余。

试题来源：PKU Campus 2009 (POJ Monthly Contest-2009.05.17), Simon

在线测试：POJ 3734

试题解析

本题给出要涂上颜色的N块积木，这N块积木排成一排，每块积木可以涂成红色、蓝色、绿色或黄色，问涂色结束后，红色的积木和绿色的积木数量都是偶数的涂色方法有多少种（结果对10007取余）。

分析Panda涂色到第i块积木时的情况。设Panda涂色到第i块积木时，红色的积木和绿色的积木数量都是偶数的方法数为ai，红色的积木和绿色的积木数量一奇一偶的方法数为bi，红色的积木和绿色的积木数量都是奇数的方法数为ci，则可以推导出如下公式：

●ai+1=2ai+bi，即：在红色的积木和绿色的积木数量都是偶数的情况下，接下来的一块积木不涂红色或绿色，涂蓝色或黄色，根据加法原理，方法数为2ai；在红色的积木和绿色的积木数量一奇一偶的情况下，涂颜色数为奇数的积木，方法数为bi；而在红色的积木和绿色的积木数量都是奇数的情况下，仅涂一块积木，无法使红色的积木和绿色的积木数量都是偶数。所以，根据加法原理，ai+1=2ai+bi。

●bi+1=2ai+2bi+2ci，即：在红色的积木和绿色的积木数量都是偶数的情况下，接下来的一块积木涂红色或绿色，方法数为2ai；在红色的积木和绿色的积木数量一奇一偶的情况下，接下来的一块积木涂蓝色或黄色，方法数为2bi；在红色的积木和绿色的积木数量都是奇数的情况下，接下来的一块积木涂红色或绿色，方法数为2ci。根据加法原理，bi+1=2ai+2bi+2ci。

●ci+1=bi+2ci，即：在红色的积木和绿色的积木数量都是偶数的情况下，仅涂一块积木，无法使红色的积木和绿色的积木数量都是奇数；在红色的积木和绿色的积木数量一奇一偶的情况下，涂颜色数为偶数的积木，方法数为bi；在红色的积木和绿色的积木数量都是奇数的情况下，接下来的一块积木涂蓝色或黄色，方法数为2ci。根据加法原理，ci+1=bi+2ci。

上述公式用矩阵表示如下：

所以，本题可以采用矩阵的快速幂运算求解。

参考程序

【1.2.2.2 Cow Relays】

N(2≤N≤1000000)头奶牛决定在牧场使用T(2≤T≤100)条奶牛赛道进行接力赛，以进行它们的体能训练计划。

每条赛道连接着两个不同的交点I1i(1≤I1i≤1000)和I2i(1≤I2i≤1000)，每个交点是至少两条赛道的节点。奶牛们知道每条赛道的长度lengthi(1≤lengthi≤1000)、赛道连接的两个交点，两个交点直接连接着两条不同的赛道的情况是不存在的。这些赛道构成一种数学上称为图的结构。

为了进行接力赛，N头奶牛要站在不同的交点上，在有些交点上可能有多头奶牛。奶牛要站在适当的位置，这样就可以一头奶牛接一头奶牛地把接力棒传递下去，最后到达终点。

请编写一个程序来确定奶牛所站的适当的位置，要找到连接起始交点（S）和结束交点（E）的最短路径，并正好通过N条奶牛赛道。

输入

第1行：由四个空格分隔的整数N、T、S和E。

第2～T+1行：第i+1行为用三个空格分隔的整数lengthi、I1i和I2i，描述第i条赛道。

输出

输出1行，给出一个整数，从交点S到交点E，且正好通过N条奶牛赛道的最短距离。

试题来源：USACO 2007 November Gold

在线测试：POJ 3613

试题解析

本题给出一个T条边的带权无向图G，求从起点S到终点E恰好经过N条边的最短路径，其中，2≤T≤100，2≤N≤1000000。所以，从起点S到终点E恰好经过的路径，边可以重复。

对于图的邻接矩阵A，Ak[i][j]表示点i到点j恰好经过k条边的路径数。基于此，设带权无向图G表示为邻接矩阵A，图G的节点数为n，并设从节点i出发到节点j的长度为k的最短路径为A(k)[i][j]，则有。对Floyd-Warshall算法进行调整如下：

这里要注意，Floyd-Warshall算法和调整的Floyd-Warshall算法有以下两处不同。

●在Floyd-Warshall算法中，中间节点变量j1作为for三层循环的最外层的循环变量；而在调整的Floyd-Warshall算法中，中间节点变量j1作为for三层循环的最内层的循环变量。

●Floyd-Warshall算法用于计算图中所有节点对之间的最短路径；对邻接矩阵A，一次调整的Floyd-Warshall算法运算，可以视为计算A2，求出的是每对节点间经过两条边的最短路径，所以，要计算每对节点间经过N条边的最短路径，就要计算AN。

对于本题，因为N的数据范围较大，所以AN的计算采用矩阵快速幂运算完成。此外，把INF设为0x3fffffff，即(1<<30)-1。

参考程序

【1.2.2.3 Training little cats】

Facer的宠物猫刚生下一窝小猫。考虑到这些可爱的猫的健康状况，Facer决定让这些猫做些运动。Facer为他的猫精心设计了一套动作。他现在要求你监督猫完成练习。Facer对猫设计的练习包含以下三种不同的动作。

●gi：给第i只猫一颗花生。

●ei：让第i只猫吃掉它所拥有的花生。

●sij：让第i只猫和第j只猫交换它们的花生。

所有的猫都要执行一个动作的序列，而且必须重复这个序列m次！可怜的猫！也只有Facer能想出这种馊主意。

请确定最终每只猫拥有的花生的数量，并输出这些数量，以便保存这些值。

输入

输入由多个测试用例组成，以三个零（0 0 0）结束。对于每个测试用例，首先给出三个整数n、m和k，其中n是猫的数目，k是一个动作序列的长度。接下来的k行给出动作序列。其中m≤1000000000、n≤100、k≤100。

输出

对于每个测试用例，在一行中输出n个数字，表示最终每只猫拥有的花生的数量。

试题来源：PKU Campus 2009 (POJ Monthly Contest-2009.05.17), Facer

在线测试：POJ 3735

试题解析

本题给出n只猫，开始时每只猫有0颗花生，猫能做三种动作。给出一个由k个动作组成的动作序列，将一个动作序列做m次后，问：最终每只猫有多少颗花生？

对于本题，由于要重复一个动作序列m次，而m≤1000000000，因此不能采用模拟的方法解题。

本题采用矩阵表示动作。初始矩阵A为一个n+1元组，编号为0～n，0号元素固定为1，1～n号元素分别为对应的猫所拥有的花生数。以3只猫为例，初始矩阵A=(1 0 0 0)，再构造一个(n+1)×(n+1)的单位矩阵M，在单位矩阵M上表示动作。

gi：给第i只猫一颗花生，在M上使M[0][i]变为1。例如g 2，则：

，A×M=(1 0 1 0)，也就是第2只猫有一颗花生，其他猫没有花生。

ei：让第i只猫吃掉它所拥有的花生，在M上使M[i][i]变为0。例如，第1、第2、第3只猫分别有2、3、4颗花生，则当前矩阵A=(1 2 3 4)，当前动作为e 3，则：

，A×M=(1 2 3 0)，也就是第1、第2只猫分别有2颗和3颗花生，而第4只猫没有花生。

sij：让第i只猫和第j只猫交换它们的花生，在M上使第i列与第j互换。例如，第1、第2、第3只猫分别有2、3、4颗花生，则当前矩阵A=(1234)，当前动作为s 1 3，则

，A×M=(1 4 3 2)，也就是第1、第2、第3只猫分别有4、3、2颗花生，第1只猫和第3只猫交换了它们的花生。

对于k个动作，按动作顺序，根据动作种类，对(n+1)×(n+1)的单位矩阵M处理如下。

●gi：给第i只猫一颗花生，M[0][i]++。

●ei：让第i只猫吃掉它所拥有的花生，将M上第i列的所有元素置为0。

●sij：让第i只猫和第j只猫交换它们的花生，在M上使第i列与第j互换。

这样就得到了表示一个动作序列的矩阵T。由于要重复一个动作序列m次，因此最终每只猫拥有的花生的数量为A×Tm。

由于稀疏矩阵，在做矩阵乘法时要进行优化。

参考程序