什么是矩阵
矩阵，在数学上，矩阵是指纵横排列的二维数据表格，最早来自于方程组的系数及常数所构成的方阵。矩阵是高等代数学中的常见工具，也常见于统计分析等应用数学学科中。

长这个样子：



矢量也可以转为矩阵，可以看成nX1的行矩阵，或1Xn的矩阵。

矩阵列的运行比较复杂，下面就来一一探讨。
矩阵和标量的乘法
直接标量与各个分量相乘即可，不多废话了…同时kM=Mk即，谁在哪边都一样。
矩阵与矩阵的乘法
它会得到一个新的矩阵，而且维度与这两个矩阵有关系。

如A为4X3矩阵，B为3X6矩阵那么 AB维度就是4X6。

左矩阵的列数必须与右矩阵的行数想同，否则不能相乘。

矩阵不满足交换律：AB！=BA

满足结合律：(AB)C=A(BC) 甚至可以扩展至 ABCDE=((A(BC))D)E=(AB)(CD)E
方阵
方块矩阵，即行列数相同的矩阵。有一些运算和性质是只有方阵有具有，如对角元素

对角矩阵：



单位矩阵(I)：

单位矩阵乘完还等于原本的矩阵,设I为转：

MI=IM=M


转置矩阵(Mt)
对原矩阵的一种运算，即行变列，列变行。可以记作Mt



性制一：转两次就转回来了:

(Mt)t=M

性制二：矩阵串接转置，等于反射串接各矩阵

(AB)t=BtAt
逆矩阵(M-1)
这应该是这里最复杂的一种操作了。不是所有矩阵都有逆矩阵，它必须是一个方阵。

给定M-1来表示。最重要的特性就是M和M-1相乘会得到一个单位矩阵。也就是说：

MM-1=M-1M=I
并非所有有对应的逆矩阵，如果一个矩阵有对应的逆矩阵则这个矩阵称为是可逆的，否则称为不可逆的。

如果一个矩阵行列式不为0，那么它就是可逆的。
性质一：逆矩阵的逆矩阵就是它本身

（M-1）-1=M
性质二：单位矩阵的逆矩阵就是它本身

I-1=I
性制三：转置矩阵的逆矩阵是逆矩阵的转置

(Mt)-1=(M1)t
性质四：矩阵串接相乘后的逆矩阵等于反向串接各个矩阵的逆矩阵

(ABCD)-1=D-1C-1B-1A-1
性质五：允许我们还原这个变换

M-1(Mv)=(M-1M)v=Iv=v
正交矩阵
方正M和它的转置矩阵乘积为单位矩阵的话，它就是一个正交矩阵，即：

MMt=MtM=I
正交矩阵的逆矩阵和转置矩阵是一样的

Mt=M-1
三维变换中我们经常会需要作用逆矩阵来求解反射的变换。而逆矩阵的求解往往计算量很大，但转置矩阵就非常容易。


矩阵的每一行，即c1、c2、c3的是单位矢量，由于其相互垂直只有与自己点乘才能得到1，其他为0.
矩阵与矢量相乘
我们需要把矢量先转成行矩阵或是列矩阵，但要满足矩阵相乘的条件。通常我们使用右乘。

                    矩阵乘法的常数优化
                      philipsweng


虽然说作为键盘科学家，我们更应该关心程序的时间复杂度。但是一个写的不好的程序可能在实际运行会跟时间复杂度更差的程序差不了多少。我们我们也应该注意程序的常数优化。 

对于矩阵乘法来讲。我们实际上可以比较这两种打法。




实际上第二种打法在绝大多数情况下效率约为第一种的两倍。 

这是为什么？？？ 

注意到第二个程序的k放在了第二层。 

C[i][j] = a[i][k] * b[k][j] 

那么我们枚举的顺序就使得a[i]数据顺序枚举以及b[k]顺序枚举。 

那么我们数组的寻址速度会增加许多。 

然后对于取模操作，我们并不需要每次都取模，因为取模太缓慢。我们可以界定一个范围再取模（慎用。容易溢出）

今天我们不讲雷达，我们来讲讲协方差矩阵。为什么呢？因为如果我直接的去说超分辨估计的东西，放在一篇文章里面太长了，我自己写着都累，别人看着也累。同时协方差矩阵这个东西又比较基础，不仅雷达要用，其他各个工程领域也喜欢这个数学工具。所以今天我们就来把基础打好，学一学这个神奇又好玩的数学工具——协方差矩阵
一维数据的分析方法——方差
方差是在概率论和统计方差衡量随机变量或一组数据时离散程度的度量，样本方差的无偏估计公式为：
sigma = 1;
x = sigma .* randn(100000,1);
figure(1)
ksdensity(x);grid on;title('概率分布 sigma = 1')
方差能很好的解释了一维数据的分布特性，但是当我们开始分析二维数据的时候，完球了！
新问题——二维到N维数据
我们来构造一组二维数据：
x = 4.*randn(1,1000);
y = 1.*randn(1,1000);

D = [x;y];

plot(D(1,:),D(2,:),'.');
grid on;
axis equal;
我们构造的这组数据，一个是沿着X轴标准差为4和以及沿着Y轴标准差为1的数据，画图看到是这样的，但是这样子其实过于简单，大多数真实情况是比较复杂的，这种根本解决不了问题，为了增加点难度，我们引入了多元正态分布与线性变换的概念。
更真实的模拟——多元正态分布与线性变换
我们再搞点事情，用旋转矩阵把数据给他搅一下，实际上是我们在对原始数据做了一个变换 
 ，矩阵 
 被称为
变换矩阵(transformation matrix),为了将图中的点经过线性变换得到我们想要的，我们需要构造这个矩阵
例如有尺度矩阵
或者旋转矩阵
其中 
 是旋转的度数，那么在这个例子中我们用旋转矩阵来折腾数据，随手写个旋转矩阵函数
function out =rotate_2D(thetad) 
out =[cosd(thetad)-sind(thetad);sind(thetad)cosd(thetad)];
end
对原始数据进行旋转
N = 1000;
x = 4.*randn(1,N);
y = 1.*randn(1,N);

D = [x;y];

new_D = rotate_2d(40)*D;

plot(new_D(1,:),new_D(2,:),'.');grid on;axis equal
或者可以这么写，用矩阵的处理方法
N = 1000;
x = 1.*randn(1,N);
y = 1.*randn(1,N);

D = [x;y];
T = rotate_2d(40)*[4 0;0 1];
new_D = T*D;

plot(new_D(1,:),new_D(2,:),'.');grid on;axis equal
上面两种写法是等价的，都可以，不过为了理解方便，我用尺度矩阵和旋转矩阵继续往下整（挖坑了）。那么这个4和1就是标准差的缩放了，因此变换矩阵可以表示为
如果式中 
 ，我们会看到下图：
这下难度够了。二维数据，而且并没有像之前那样子是按照X和Y轴分布的，比较接近真实情况，构造数据的问题解决了，那么下来怎么分析呢？
新工具——协方差/协方差矩阵
看到这么一种数据，如果按照我们如果按照一维数据的分析思路方法，计算在X方向上的方差 
 和Y方向上的方差 
 ，肯定是会出问题。因为图形明显的告诉你，单纯的从水平方向和垂直方向去分析并不能对这个斜着分布的关系做出合理的解释，并且x值的增加同时y值也在增加，沿着这
两个轴分解产生了正相关。怎么办？
这时候我们就要扩展方差的概念了，为了捕捉这个特征，一帮闲的没事的人们就开始整活——协方差出现了：
对于二维的数据，我们可以得到的有这四个东西
我们用矩阵来表示，这个矩阵称为协方差矩阵：
由于X与Y的相关性和Y与X的相关性是一样的，因此 
协方差矩阵始终是一个对称阵，其对角线上的是方差，非对角线上的是协方差。二维正态分布的数据由它的均值和 2x2 的协方差矩阵就可以表示，我们将刚才的数据协方差矩阵打印出来：
Rxx = new_D *new_D './(N-1)
随着数据形态的变化，协方差矩阵的值也跟着在变，大家可以在代码里面修改一些参数试试，我就不上图了，那么问题又来了，协方差矩阵有了，我们要怎么分析，才能较好的还原数据的原始特性呢？
协方差矩阵的特征分解
从前面的部分可以知道，协方差矩阵定义了数据的离散程度(方差)和空间分布走向(协方差)。因此，我们可以认为，指向最大离散度方向的向量和其值(该值等于这个方向上的方差)，表示协方差矩阵。
假设前面的数据是 
 ，实际上我们刚才的操作其实是对数据 
 进行线性变换，可以得到：
其中 
 是含有旋转矩阵 
 和缩放矩阵 
 的变换矩阵：
矩阵 
 和 
 的定义我们再贴出来复习一下：
接下来我们就要讨论协方差矩阵 
 和线性变换 
 之间的关系。
我们来回忆一下刚才的数据构造流程，先从一组无缩放的数据开始，也就是标准差为1的无旋转，服从正态分布的数据，我们可以称之为"白数据”。
白数据的协方差矩阵，等于单位阵，因此方差和标准差为1，协方差为0：
现在，对 
 数据的标准差缩放到4：
变换后的数据 
 如图所示：
 的协方差矩阵 
 为：
故 
 的协方差矩阵 
 ，与作用在原始数据上的线性变换矩阵 
 的关系为 
 ,并且：
但是这个式子是否适用于旋转变换呢？我们将协方差矩阵分解为旋转和缩放矩阵的乘积。
协方差矩阵可以用特征向量和特征值表示，在二维情况下，有2个特征向量和2个特征值。通过矩阵来定义两个等式的关系：
式中 
 是由 
 的各特征向量构成的矩阵， 
 为相应特征向量对应的非0特征值组成的对角矩阵，上式可以写为我们常见的形式：
这个式子我们称为协方差矩阵的特征分解。由于特征向量表示了最大方差的方向，特征值表示了该方向上方差的值，也就是说，
 实际上就是旋转矩阵，而 
 就是缩放矩阵，
协方差矩阵可以进一步表示为:
式中 
 为旋转矩阵， 
 为缩放矩阵。因为 
 为对角缩放矩阵，所以 
 ，并且由于 
 为正交矩阵，所以 
。那么协方差矩阵最后的形式可以写为：
也就是说，协方差矩阵特征分解后的特征向量，总是指向数据最大方差的方向，并规定了其走向。而且由于旋转矩阵的正交性，次最大特征值对应的特征向量，总是与最大特征值的特征向量正交。
我们把代码补完，验证一下：
N = 1000;
x = 1.*randn(1,N);
y = 1.*randn(1,N);

D = [x;y];
T = rotate_2d(40)*[4 0;0 1];
new_D = T*D;


Rxx = new_D*new_D'./(N-1);
[V,D] = eig(Rxx)

figure(1)
plot(new_D(1,:),new_D(2,:),'.');grid on;axis equal;hold on;
quiver(0,0,V(1,1).*sqrt(D(1,1)),V(2,1).*sqrt(D(1,1)),'LineWidth',3)
quiver(0,0,V(1,2)*sqrt(D(2,2)),V(2,2)*sqrt(D(2,2)),'LineWidth',3)

legend('数据分布','特征向量1','特征向量2')
hold off
我们看到特征向量其实就是原来的旋转矩阵，而特征值的结果也符合原来数据的方差，即1和16，即为原来的标准差1和4。
那么更高维度的数据是否也符合要求呢？
我们构造一个三维数据试一试：
旋转函数：
function yaw = rotate_yd(theta)
yaw = [cosd(theta) 0 sind(theta);0 1 0;-sind(theta) 0 cosd(theta)];
end
代码部分：
clear;clc;format longG

%--------------------------------------------------------------------------
%   一组标准的概率分布，构造一个二维数据
%--------------------------------------------------------------------------
N = 2000;
x = 5.*randn(1,N);
y = 1.*randn(1,N);
z = 3.*randn(1,N);

A = [x;y;z];

A = rotate_zd(0)*rotate_yd(20)*rotate_xd(0)*A;


Rxx = A*A'./(N-1);
[V,D] = eig(Rxx);

figure(1)
plot3(A(1,:),A(2,:),A(3,:),'.');grid on;axis equal;hold on;

quiver3(0,0,0,V(1,1).*sqrt(D(1,1)).*5,V(2,1).*sqrt(D(1,1)).*5,V(3,1).*sqrt(D(1,1)).*5,'LineWidth',2)
quiver3(0,0,0,V(1,2)*sqrt(D(2,2)).*5,V(2,2)*sqrt(D(2,2)).*5,V(3,2)*sqrt(D(2,2)).*5,'LineWidth',2)
quiver3(0,0,0,V(1,3)*sqrt(D(3,3)).*5,V(2,3)*sqrt(D(3,3)).*5,V(3,3)*sqrt(D(3,3)).*5,'LineWidth',2)
hold off
适当的缩放了一下特征向量的长度，我乘了个5进去，为了画图好看，我们看看图形：
总结
数据的协方差矩阵，直接与数据的线性变换有关，该线性变换完全由特征向量和特征值定义，特征向量表示了矩阵的旋转，特征值对应各个维度上的缩放因子的平方。
参考：
Xinyu Chen：如何直观地理解「协方差矩阵」？​zhuanlan.zhihu.com
A geometric interpretation of the covariance matrix​www.visiondummy.com
翻译：协方差矩阵的几何解释​njuferret.github.io
作为2020年的最后一天，以这篇文章完满结束收尾。我是从今年才开始决定写文章，想把自己的一些学习心得更新上来，但是自己又比较拖延，速度比较慢，希望大家能够见谅，写一篇文章下来需要好多天，每篇都保证较好的图，文字和公式说明，所以确实是一件花心思的事情。
那么希望在2021年里，能够继续保持下去，持续的把自己在工作中学习中的一些理解分享给大家，能够帮助到你们就是最大的快乐。
最后，有不足之处欢迎指正，评论区不能及时回复，望见谅，转载请注明知乎来源，谢谢~

矩阵乘法是个很高端的东西。 

注意事项： 以下的讲述下标都从0开始



矩阵是什么
矩阵加减法
矩阵乘法

矩阵乘法的运算定律
模板
矩阵乘法的应用


例题一 斐波那契数列

描述

题目背景
题目描述
输入格式
输出格式
输入样例1
输出样例1
输入样例2
输出样例2
说明


解法
代码


例题二 NOIP2013模拟联考14图形变换


Description
Input
Output
Sample Input
Sample Output
Data Constraint
Hint


解法
代码


例题三 SDOI2009HH去散步

描述

题目描述
输入格式
输出格式
输入样例1
输出样例1
说明


解法
代码


总结
补充







矩阵是什么？

不解释



矩阵加减法？

矩阵加矩阵，就将对应的加上。 

矩阵加常数，就将每一个元素都加上这个常数。 

减法同理。



矩阵乘法

一张好图，在这里发现的，方便理解矩阵乘法。 

 

一个$m*n$的的A矩阵，和一个$n*p$的B矩阵相乘，将得到一个$m*p$的矩阵C 

 $$C(i,j)=\sum_{k=1}^n A(i,k)*B(k,j)$$  

可以简单地理解为，A中i行的元素，与B中j列的元素，对应相乘得到的和。




矩阵乘法的运算定律

（1） 不满足交换律 

（2） 满足结合律：(AB)C=A(BC) 

（3） 满足分配律：(A+B)C=AC+BC A(B+C)=AB+AC



模板



template <typename T,int M,int N>
    struct Matrix
    {
        T mat[M][N];
        inline Matrix(){memset(mat,0,sizeof mat);}
        inline T* operator[](int i){return mat[i];}
    };
template <typename T,int M,int N,int P>
    inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y)
        {
            Matrix<T,M,P> ret;
            int i,j,k;
            for (i=0;i<M;++i)
                for (j=0;j<P;++j)
                    for (k=0;k<N;++k)
                        ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j];
            return ret;
        }
template <typename T,typename T_>
    inline T pow(T x,T_ y)
    {
        T ret=x;
        --y;//这两句话可以忽略初值问题，但y<=0时就悲剧了
        while (y)
        {
            if (y&1)
                ret=ret*x;
            y>>=1;
            x=x*x;  
        }
        return ret;
    }



矩阵乘法的应用

矩阵乘法可以优化时间。 

将某些操作转化为矩阵乘法的形式，就可以用快速幂减少时间。因为矩阵乘法满足结合律。





例题一 斐波那契数列



描述


  
题目背景
  
  
大家都知道，斐波那契数列是满足如下性质的一个数列： 

  • f(1) = 1 

  • f(2) = 1 

  • f(n) = f(n-1) + f(n-2) (n ≥ 2 且 n 为整数)
  
  
题目描述
  
  
请你求出 f(n) mod 1000000007 的值。
  
  
输入格式：
  
  
·第 1 行：一个整数 n
  
  
输出格式：
  
  
第 1 行： f(n) mod 1000000007 的值
  
  
输入样例#1：
  
  
5
  
  
输出样例#1：
  
  
5
  
  
输入样例#2：
  
  
10
  
  
输出样例#2：
  
  
55
  
  
说明
  
  
对于 60% 的数据： n ≤ 92 

  对于 100% 的数据： n在long long(INT64)范围内。




解法

这是一道裸斐波拉契数列。传统递推是O(N)的，显然会炸。 

斐波拉契数列有个通项，但我们在这里用矩阵乘法解决。 

我们知道f(n)是由f(n-1)和f(n-2)推过来的，不妨设一个1*2的矩阵 

 $$A=\left[\begin{matrix}
f(n-2) & f(n-1)
\end{matrix}\right]$$  

现在我们要通过A推出B 

 $$B=\left[\begin{matrix}
f(n-1) & f(n)
\end{matrix}\right]=
\left[\begin{matrix}
f(n-1) & f(n-2)+f(n-1)
\end{matrix}\right]$$  

我们要求出一个2*2的 转移矩阵 T，满足 

 $$A*T=B$$  

即 

 $$\left[\begin{matrix}
f(n-2) & f(n-1)
\end{matrix}\right]*T=\left[\begin{matrix}
f(n-1) & f(n-2)+f(n-1)
\end{matrix}\right]$$  

我们知道，B[0][0]=A[0][1] B[0][1]=A[0][0]+A[0][1] 

对于T[0][0]，A[0][0]对B[0][0]没有贡献，所以T[0][0]=0 

对于T[1][0]，A[0][1]对B[0][0]有贡献，B[0][0]=A[0][1]*1，所以T[1][0]=1 

对于T[0][1]，A[0][0]对B[0][1]有贡献，B[0][1]=A[0][0]*1+A[0][1]*1，所以T[0][1]=1 

对于T[1][1]，A[1][1]对B[0][1]有贡献，B[1][1]=A[0][0]*1+A[0][1]*1，所以T[1][1]=1 

综上所述， 

 $$T=\left[\begin{matrix}
0 & 1 \\
1 & 1
\end{matrix}\right]$$  

显然这个式子是正确的。 

 $$\left[\begin{matrix}
f(n-2) & f(n-1)
\end{matrix}\right]*\left[\begin{matrix}
0 & 1 \\
1 & 1
\end{matrix}\right]=\left[\begin{matrix}
f(n-1) & f(n-2)+f(n-1)
\end{matrix}\right]$$  

求出这个矩阵后，就可以愉快地快速幂了。 

时间复杂度O(lgN)




代码



#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define mod 1000000007
template <typename T,int M,int N>
    struct Matrix
    {
        T mat[M][N];
        inline Matrix(){memset(mat,0,sizeof mat);}
        inline T* operator[](int i){return mat[i];}
    };
template <typename T,int M,int N,int P>
    inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y)
        {
            Matrix<T,M,P> ret;
            int i,j,k;
            for (i=0;i<M;++i)
                for (j=0;j<P;++j)
                    for (k=0;k<N;++k)
                        (ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j])%=mod;
            return ret;
        }
template <typename T,typename T_>
    inline T pow(T x,T_ y)
    {
        T ret=x;
        --y;
        while (y)
        {
            if (y&1)
                ret=ret*x;
            y>>=1;
            x=x*x;  
        }
        return ret;
    }
int main()
{
    long long n;
    scanf("%lld",&n);
    if (n==1)
        putchar('1');
    else
    {
        Matrix<long long,2,2> tmp;
        tmp[0][1]=tmp[1][0]=tmp[1][1]=1;
        tmp=pow(tmp,n-1);
        Matrix<long long,1,2> a;
        a[0][1]=1;//a=[f(0) f(1)]
        a=a*tmp;//[f(0) f(1)]*(T^(n-1))=[f(n-1) f(n)]
        printf("%lld\n",a[0][1]);
    }
}



例题二 【NOIP2013模拟联考14】图形变换


  
Description
  
  
翔翔最近接到一个任务，要把一个图形做大量的变换操作，翔翔实在是操作得手软，决定写个程序来执行变换操作。 

  翔翔目前接到的任务是，对一个由n个点组成的图形连续作平移、缩放、旋转变换。相关操作定义如下： 

  Trans(dx,dy) 表示平移图形，即把图形上所有的点的横纵坐标分别加上dx和dy； 

  Scale(sx,sy) 表示缩放图形，即把图形上所有点的横纵坐标分别乘以sx和sy； 

  Rotate(θ,x0,y0) 表示旋转图形，即把图形上所有点的坐标绕(x0,y0)顺时针旋转θ角度 

  由于某些操作会重复运行多次，翔翔还定义了循环指令： 

  Loop(m) 

  …  

  End  

  表示把Loop和对应End之间的操作循环执行m次，循环可以嵌套。
  
  
Input
  
  
第一行一个整数n(n<=100)表示图形由n个点组成； 

  接下来n行，每行空格隔开两个实数xi，yi表示点的坐标； 

  接下来一直到文件结束，每行一条操作指令。保证指令格式合法，无多余空格。
  
  
Output
  
  
输出有n行，每行两个空格隔开实数xi，yi表示对应输入的点变换后的坐标。 

  本题采用Special Judge判断，只要你输出的数值与标准答案误差不能超过1即可。
  
  
Sample Input
  
  
3 

  0.5 0 

  2.5 2 

  -4.5 1 

  Trans(1.5,-1) 

  Loop(2) 

  Trans(1,1) 

  Loop(2) 

  Rotate(90,0,0) 

  End 

  Scale(2,3) 

  End
  
  
Sample Output
  
  
10.0000 -3.0000 

  18.0000 15.0000 

  -10.0000 6.0000
  
  
Data Constraint
  
  
保证操作中坐标值不会超过double范围，输出不会超过int范围； 

  指令总共不超过1000行； 

  对于所有的数据,所有循环指令中m<=1000000； 

  对于60%的数据,所有循环指令中m<=1000； 

  对于30%的数据不含嵌套循环。
  
  
Hint
  
  
【友情提醒】 

  pi的值最好用系统的值。C++的定义为：#define Pi M_PI 

  Pascal就是直接为：pi 

  不要自己定义避免因为pi带来的误差。




解法

旋转公式： 

 $$绕(a,b)逆时针转\theta
\\
x'=(x-a)*cos\theta-(y-b)*sin\theta+a
\\
y'=(x-a)*sin\theta+(y-b)*cos\theta+b$$  

暴力铁定超时。 

先考虑使用2*2的转移矩阵T，但是只能对付乘法，加法和旋转就要用矩阵加法了。矩阵乘法和矩阵加法混在一起很麻烦（矩阵套矩阵？）。 

我们可以新增一个常数项 

 $$A=\left[\begin{matrix}
x & y & 1
\end{matrix}\right]$$  

现在要求出T1、T2、T3满足 

 $$\left[\begin{matrix}
x & y & 1
\end{matrix}\right]*T1=\left[\begin{matrix}
x+a & y+b & 1
\end{matrix}\right]
\\
\left[\begin{matrix}
x & y & 1
\end{matrix}\right]*T2=\left[\begin{matrix}
x*a & y*b & 1
\end{matrix}\right]
\\
\left[\begin{matrix}
x & y & 1
\end{matrix}\right]*T3=\left[\begin{matrix}
(x-a)*cos\theta-(y-b)*sin\theta+a
&
(x-a)*sin\theta+(y-b)*cos\theta+b
&
1
\end{matrix}\right]$$  

这样加、乘就特别容易，将旋转公式用乘法分配律拆开，也可以得出转移矩阵 

 $$T1=\left[\begin{matrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
a & b & 1
\end{matrix}\right]
\\
T2=\left[\begin{matrix}
a & 0 & 0 \\
0 & b & 0 \\
0 & 0 & 1
\end{matrix}\right]
\\
T3=\left[\begin{matrix}
\cos\theta & \sin\theta & 0 \\
-\sin\theta & \cos\theta & 0 \\
a-a*\cos\theta+b*\sin\theta & b-a*\sin\theta-b*\cos\theta & 1 \\
\end{matrix}\right]$$  

可以代进去看看，是成立的。 

于是我们可以将所有操作乘起来（有循环时递归，算出这个循环里一次的乘积后快速幂），最后在将每个坐标乘这个积，就是答案。




代码

英语不好，将theta打成xida，懒得改，注意一下就好了。



#include <cstdio>
#include <cstring>
#include <cmath>
#include <cassert>
#include <algorithm>
using namespace std;
#define I_O(x) freopen(""#x".in","r",stdin);freopen(""#x".out","w",stdout)
#define Pi 3.14159265358979323846
template <typename T,int M,int N>
    struct Matrix
    {
        T mat[M][N];
        inline Matrix(){memset(mat,0,sizeof mat);}
        inline T* operator[](int i){return mat[i];}
    };
template <typename T,int M,int N,int P>
    inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y)
        {
            Matrix<T,M,P> ret;
            int i,j,k;
            for (i=0;i<M;++i)
                for (j=0;j<P;++j)
                    for (k=0;k<N;++k)
                        ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j];
            return ret;
        }
template <typename T,typename T_>
    inline T pow(T x,T_ y)
    {
        T ret=x;
        --y;
        while (y)
        {
            if (y&1)
                ret=ret*x;
            y>>=1;
            x=x*x;  
        }
        return ret;
    }
int n;
Matrix<double,1,3> d[100];
char cz[101];
Matrix<double,3,3> dfs(int);
int main()
{
    I_O(transform);
    scanf("%d",&n);
    int i;
    double x,y;
    for (i=0;i<n;++i)
    {
        scanf("%lf%lf",&d[i][0][0],&d[i][0][1]);
        d[i][0][2]=1;
    }
    Matrix<double,3,3> a,tmp1,tmp2,tmp3;//分三个变量是为避免多次赋值
    a[0][0]=a[1][1]=a[2][2]=1;//a的初值
    tmp1[0][0]=tmp1[1][1]=tmp1[2][2]=1;
    tmp2[2][2]=1;
    tmp3[2][2]=1;
    double xida,tmp_sin,tmp_cos;
    int t;
    while (scanf("%s",cz)==1)
    {
        switch (*cz)
        {
            case 'T':
//              1 0 0
//              0 1 0
//              a b 1
                sscanf(cz,"Trans(%lf,%lf",&tmp1[2][0],&tmp1[2][1]);
                a=a*tmp1;
            break;
            case 'S':
//              a 0 0
//              0 b 0
//              0 0 1
                sscanf(cz,"Scale(%lf,%lf",&tmp2[0][0],&tmp2[1][1]);
                a=a*tmp2;
            break;
            case 'R':
                sscanf(cz,"Rotate(%lf,%lf,%lf",&xida,&x,&y);
                xida=-xida/180*Pi;//顺时针角度转成逆时针弧度
//              cos           sin           0
//              -sin          cos           0
//              a-a*cos+b*sin b-a*sin-b*cos 1
                tmp_sin=sin(xida);
                tmp_cos=cos(xida);
                tmp3[0][0]=tmp_cos; tmp3[0][1]=tmp_sin; tmp3[0][2]=0;
                tmp3[1][0]=-tmp_sin; tmp3[1][1]=tmp_cos; tmp3[1][2]=0;
                tmp3[2][0]=x-x*tmp_cos+y*tmp_sin; tmp3[2][1]=y-x*tmp_sin-y*tmp_cos;
                a=a*tmp3;
            break;
            default:
                sscanf(cz,"Loop(%d",&t);
                a=a*dfs(t);
        }
    }
    Matrix<double,1,3> tmp;
    for (i=0;i<n;++i)
    {
        tmp=d[i]*a;
        printf("%.4lf %.4lf\n",tmp[0][0],tmp[0][1]);
    }
}
Matrix<double,3,3> dfs(int t)
{
    double x,y;
    Matrix<double,3,3> a,tmp1,tmp2,tmp3;
    a[0][0]=a[1][1]=a[2][2]=1;
    tmp1[0][0]=tmp1[1][1]=tmp1[2][2]=1;
    tmp2[2][2]=1;
    tmp3[2][2]=1;
    int times;
    double xida,tmp_sin,tmp_cos;
    for (scanf("%s",cz);*cz!='E';scanf("%s",cz))
    {
        switch (*cz)
        {
            case 'T':
                sscanf(cz,"Trans(%lf,%lf",&tmp1[2][0],&tmp1[2][1]);
//              1 0 0
//              0 1 0
//              a b 1
                a=a*tmp1;
            break;
            case 'S':
                sscanf(cz,"Scale(%lf,%lf",&tmp2[0][0],&tmp2[1][1]);
//              a 0 0
//              0 b 0
//              0 0 1
                a=a*tmp2;
            break;
            case 'R':
                sscanf(cz,"Rotate(%lf,%lf,%lf",&xida,&x,&y);
                xida=-xida/180*Pi;
//              cos           sin           0
//              -sin          cos           0
//              a-a*cos+b*sin b-a*sin-b*cos 1
                tmp_sin=sin(xida);
                tmp_cos=cos(xida);
                tmp3[0][0]=tmp_cos; tmp3[0][1]=tmp_sin; tmp3[0][2]=0;
                tmp3[1][0]=-tmp_sin; tmp3[1][1]=tmp_cos; tmp3[1][2]=0;
                tmp3[2][0]=x-x*tmp_cos+y*tmp_sin; tmp3[2][1]=y-x*tmp_sin-y*tmp_cos;
                a=a*tmp3;
            break;
            default:
                sscanf(cz,"Loop(%d",&times);
                a=a*dfs(times);
        }
    }
    return pow(a,t);
}



例题三 【SDOI2009】HH去散步



描述


  
题目描述
  
  
HH有个一成不变的习惯，喜欢饭后百步走。所谓百步走，就是散步，就是在一定的时间 内，走过一定的距离。 但是同时HH又是个喜欢变化的人，所以他不会立刻沿着刚刚走来的路走回。 又因为HH是个喜欢变化的人，所以他每天走过的路径都不完全一样，他想知道他究竟有多 少种散步的方法。 

  现在给你学校的地图（假设每条路的长度都是一样的都是1），问长度为t，从给定地 点A走到给定地点B共有多少条符合条件的路径
  
  
输入格式：
  
  
第一行：五个整数N，M，t，A，B。其中N表示学校里的路口的个数，M表示学校里的 路的条数，t表示HH想要散步的距离，A表示散步的出发点，而B则表示散步的终点。 

  接下来M行，每行一组Ai，Bi，表示从路口Ai到路口Bi有一条路。数据保证Ai != Bi，但 不保证任意两个路口之间至多只有一条路相连接。 路口编号从0到N − 1。 同一行内所有数据均由一个空格隔开，行首行尾没有多余空格。没有多余空行。 答案模45989。
  
  
输出格式：
  
  
一行，表示答案。
  
  
输入样例#1：
  
  
4 5 3 0 0 

  0 1 

  0 2 

  0 3 

  2 1 

  3 2
  
  
输出样例#1：
  
  
4
  
  
说明
  
  
对于30%的数据，N ≤ 4，M ≤ 10，t ≤ 10。 

  对于100%的数据，N ≤ 50，M ≤ 60，t ≤ 2^30，0 ≤ A,B




解法

先说一下大意（我等好久才弄懂）。 

给你一张图，要从A走到B，不能走上一次走的路，要走t步，问方案数。 

不能走上一次走的路，一次。 

举个例子： 

 

路径可以这样：0->1->2->3->1->0 

所以我们可以想到DP 

设f[i][j]表示第i步走了j这条边的方案数，这样在转移时，就可以方便地判断，避免走上一次走的路。 

转移：$f[i][bc]=\sum_{ab\neq bc} f[i-1][ab]$ 

但是，我们发现，每次的决策点都是一样的，很明显能直接用矩阵乘法来转移，快速幂即可。



代码



#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
using namespace std;
template <typename T,int M,int N>
    struct Matrix
    {
        T mat[M][N];
        inline Matrix(){memset(mat,0,sizeof mat);}
        inline T* operator[](int i){return mat[i];}
    };
template <typename T,int M,int N,int P>
    inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y)
        {
            Matrix<T,M,P> ret;
            int i,j,k;
            for (i=0;i<M;++i)
                for (j=0;j<P;++j)
                    for (k=0;k<N;++k)
                        ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j];
            for (i=0;i<M;++i)
                for (j=0;j<P;++j)
                    ret.mat[i][j]%=45989;//减少mod运算（mod运算很慢的）
            return ret;
        }
template <typename T,typename T_>
    inline T pow(T x,T_ y)
    {
        T ret=x;
        --y;
        while (y)
        {
            if (y&1)
                ret=ret*x;
            y>>=1;
            x=x*x;  
        }
        return ret;
    }
int n,m,t,u,v;
struct EDGE
{
    int from,to;
    EDGE* las;
    EDGE* rev;
} e[120];
EDGE* last[20];
Matrix<long long,1,120> f;
Matrix<long long,120,120> T;
int li[60];//用于记录连接终点的边
int nl;
int main()
{
    scanf("%d%d%d%d%d",&n,&m,&t,&u,&v);
    int i,j=-1,x,y;
    for (i=1;i<=m;++i)
    {
        scanf("%d%d",&x,&y);
        ++j;
        e[j]={x,y,last[x],e+j+1};
        last[x]=e+j;
        if (y==v)
            li[nl++]=j;
        ++j;
        e[j]={y,x,last[y],e+j-1};
        last[y]=e+j;
        if (x==v)
            li[nl++]=j;
    }
    EDGE* ei;
    for (ei=last[u];ei;ei=ei->las)
        f[0][int(ei-e)]=1;
    m<<=1;
    EDGE *ej,*end=e+m;
    for (ei=e;ei<end;++ei)
        for (ej=last[ei->to];ej;ej=ej->las)
            if (ei->rev!=ej)
                T[int(ei-e)][int(ej-e)]=1;//生成转移矩阵
    f=f*pow(T,t-1);
    int ans=0; 
    for (i=0;i<nl;++i)
        ans+=f[0][li[i]];
    printf("%d\n",ans%45989);
}



总结


遇到矩阵乘法的题，先将转移后的式子拆开，再帮它配转移矩阵。
如果用普通的方式配矩阵必须出现加法，那么可以尝试加常数项上去。
若有某些DP（还是递推？）的题目的决策点固定，那么可以在程序中生成转移矩阵，直接快速幂。




补充

矩阵乘法模板2.0



template <typename T,int M,int N>
    struct Matrix
    {
        int m,n;
        T mat[M][N];
        inline Matrix(){m=M;n=N;memset(mat,0,sizeof mat);}
        inline void SetUpSize(int _m,int _n){m=_m;n=_n;} 
        inline T* operator[](int i){return mat[i];}
    };
template <typename T,int M,int N,int P>
inline Matrix<T,M,P> operator*(const Matrix<T,M,N>& x,const Matrix<T,N,P>& y)
    {
        Matrix<T,M,P> ret;
        int i,j,k;
        for (i=0;i<x.m;++i)
            for (j=0;j<y.n;++j)
                for (k=0;k<x.n;++k)
                    ret.mat[i][j]+=x.mat[i][k]*y.mat[k][j];
        return ret;
    }
template <typename T,typename T_>
    inline T pow(T x,T_ y)
    {
        T ret=x;
        --y;//这两句话可以忽略初值问题，但y<=0时就悲剧了
        while (y)
        {
            if (y&1)
                ret=ret*x;
            y>>=1;
            x=x*x;  
        }
        return ret;
    }