R语言中矩阵运算
目录：
1_矩阵的生成
2_矩阵的四则运算
3_矩阵的矩阵运算
4_矩阵的分解
1_1 被定义为数组
仅当向量已明确设置了维数属性（如dim属性）后才能被视为数组. 比如:
z=1:12;
dim(z)=c(3,4);
z;
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
注意：生成矩阵是按列排列的。
通过array()函数创建多维数组及其使用方法
该函数用于创建（二维数组）矩阵，并接受多个参数配置其结构特性。其调用格式为matrix(data=NA, nrow=1, ncol=1, byrow=FALSE, dimnames=NULL)，其中data表示输入向量值；nro、ncol分别指定行数与列数；byrow逻辑变量决定数据填入方向，默认值为FALSE即按列填充；dimname参数可指定各维名称，默认为空值。例如构建一个3x5阶矩阵的操作如下：
2_矩阵的四则运算
支持对数组之间的加减乘除四则运算（+、-、*、/），这些操作具体来说就是对对应位置上的元素分别进行算术运算。通常情况下参与运算的矩阵或数组具有相同的维度；但也允许计算不同维度的情况，并需将对应的元素补齐。
3_1 转置运算
对于矩阵A，函数t(A)表示矩阵A的转置，如：
A=matrix(1:6,nrow=2);
A;
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
t(A);
[,1] [,2]
[1,] 1 2
[2,] 3 4
[3,] 5 6
3_2 求方阵的行列式
函数det()是求矩阵行列式的值，如
det(matrix(1:4,ncol=2));
[1] -2
3_3 向量的内积
对于n维向量x，可以看成nxl阶矩阵或lxn阶矩阵。若x与y是相同
维数的向量，则x%*%Y表示x与y作内积.例如，
x=1:5; Y=21:5
x%%y
[，1]
[1，]110
函数crossprod()是内积运算函数(表示交叉乘积)，crossprod(x,y)计算向量x与y的内积，即t(x) %% y'。crossprod(x)表示x与x的内积.
类似地，tcrossprod(x,y)表示’x%%t(Y)’，即x与y的外积，也称为叉积。tcrossprod(x)表示x与x作外积.如：
x=1:5; y=2*1:5;
crossprod(x);
[,1]
[1,] 55
crossprod(x,y);
[,1]
[1,] 110
tcrossprod(x);
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 2 4 6 8 10
[3,] 3 6 9 12 15
[4,] 4 8 12 16 20
[5,] 5 10 15 20 25
tcrossprod(x,y);
[,1] [,2] [,3] [,4] [,5]
[1,] 2 4 6 8 10
[2,] 4 8 12 16 20
[3,] 6 12 18 24 30
[4,] 8 16 24 32 40
[5,] 10 20 30 40 50
3_4 向量的外积(叉积)
设x和y是n维向量，则x%o%y表示x与y作外积.例如
x%o%y;
[,1] [,2] [,3] [,4] [,5]
[1,] 2 4 6 8 10
[2,] 4 8 12 16 20
[3,] 6 12 18 24 30
[4,] 8 16 24 32 40
[5,] 10 20 30 40 50
该外积操作outer()相较于传统方法更为高效。
该操作用于计算两个向量X和Y之间的外积，并等价于X %o% Y。
在常规使用场景下，默认情况下该操作采用乘法作为外积运算。
其中X和Y既可以是矩阵也可以是向量。
特别有用的是该操作在生成三维曲面时的应用。
该操作能够有效地生成由X和Y坐标组成的网格数据。
矩阵乘法运算通常用于处理两个矩阵A和B之间的关系。考虑两个矩阵A和B，在标准数学框架下它们之间的乘积关系可表示为A %% B。具体而言,crossprod(A,B)即代表t(A) %% B的操作结果,而tcrossprod(A,B)则对应于A %% t(B)的过程。经过上述运算可知, x %% A %*% x构成一个二次型结构。
A ← array(c(1:9), dim = c(3, 3))
B ← array(c(9:1), dim = c(3, 3))
计算矩阵乘法运算结果如下：
[,1] [,2] [,3]
[1,] 90 54 18
[2,] 114 69 24
[3,] 138 84 30
验证crossprod(A,B)是否等于转置后的矩阵乘积：
crossprod(A,B) == t(A) %*% B
[,1] [,2] [,3]
[1,] TRUE TRUE FALSE
等等（注：这里可能存在笔误或其他错误，请检查完整代码后再使用）
第3.6节介绍生成对角矩阵以及从矩阵中提取对角元素的操作。在R语言中，默认情况下，默认情况下，默认情况下，默认情况下，默认情况下，默认情况下，默认情况下，默认情况下，默认情况下，默认情况下，默认情况下，默认情况下，默认情况下，默认情况
第3.7节 解线性方程组及求逆矩阵
A=t(array(c(1:8,10),dim=c(3,3)));
b=c(1,1,1);
x=solve(A,b);
x;
[1] -1.000000e+00 1.000000e+00 3.806634e-16
solve(A);
[,1] [,2] [,3]
[1,] -0.6666667 -1.333333 1
[2,] -0.6666667 3.666667 -2
[3,] 1.0000000 -2.000000 1
计算对称矩阵Sm及其对应的全部特征值与对应的所有单位正交 eigenvectors
$vectors
[,1] [,2] [,3]
[1,] -0.4646675 0.833286355 0.2995295
[2,] -0.5537546 -0.009499485 -0.8326258
[3,] -0.6909703 -0.552759994 0.4658502
4_1 特征值分解
(1) 定义如下：对于一个n×n矩阵a（即一个n阶方阵），设λ是一个标量（称为****），而向量是一个非零n维列向量（称为对应的）。其中a**= λ** ，则λ被称为矩阵a的一个****。（2）
(2). 在R语言中的实现：通过函数eigen(A)来计算矩阵的特征值与特征向量；具体操作流程如下：以矩阵A为例进行说明
A=array(c(1,1,1,4,2,1,9,3,1),dim=c(3,3));
D=eigen(A);
D;
$values
[1] 5.8284271 -2.0000000 0.1715729
$vectors
[,1] [,2] [,3]
[1,] -0.8597736 -9.486833e-01 0.5384820
[2,] -0.4346498 6.474883e-17 -0.7872938
[3,] -0.2680839 3.162278e-01 0.3003425
(3). 特征值分解的性质：了解当所求特征向量构成的矩阵可逆时, 会满足求解后的向量集合通过矩阵A作用后得到的结果等于对角矩阵包含这些特征值, 进而进行验证。
solve(vectors)%%A%%vectors;
[,1] [,2] [,3]
[1,] 5.828427e+00 8.339683e-16 -1.285213e-15
[2,] 1.211325e-15 -2.000000e+00 2.704000e-16
[3,] -3.471971e-16 -1.607126e-16 1.715729e-01
结果的精度还是比较高的。
4_2 矩阵的奇异值分解
$u
[,1] [,2] [,3]
[1,] -0.2093373 0.96438514 0.1616762
[2,] -0.5038485 0.03532145 -0.8630696
[3,] -0.8380421 -0.26213299 0.4785099
$v
[,1] [,2] [,3]
[1,] -0.4646675 -0.833286355 0.2995295
[2,] -0.5537546 0.009499485 -0.8326258
[3,] -0.6909703 0.552759994 0.4658502
attach(SVD);
The following object(s) are masked from 'SVD (position 3)':
d, u, v
u%%diag(d)%%t(v);
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 10
A;
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 10
4_3 qr分解
设A为mn矩阵，如果存在mm酉矩阵Q（即Q(H)Q=QQ(H)=I）和m*n阶梯形矩阵R，使得A=QR，那么此分解称为QR分解。QR分解在解决最小二乘问题、特征值计算等方面有着十分重要的作用。
#建立矩阵
A=(array(c(1:12),dim=c(4,3)));
A;
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
#进行矩阵分解
QR=qr(A);QR
$qr
[,1] [,2] [,3]
[1,] -5.4772256 -12.7801930 -2.008316e+01
[2,] 0.3651484 -3.2659863 -6.531973e+00
[3,] 0.5477226 -0.3781696 7.880925e-16
[4,] 0.7302967 -0.9124744 9.277920e-01
$rank
[1] 2
$qraux
[1] 1.182574 1.156135 1.373098
$pivot
[1] 1 2 3
attr(,"class")
[1] "qr"
#提取Q,R并验证分解的正确性。
Q=qr.Q(QR);
R=qr.R(QR);
Q%*%R;
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
引言部分：
通过分析特征值的分解过程可以看出，在一组特征向量线性相关时，并不能直接应用特征值分解的方法对其进行矩阵分解。例如，在计算过程中我们发现对于矩阵A=t(array(c(6,12,19,-9,-20,-33,4,9,15),dim=c(3,3)))而言：
A=t(array(c(6,12,19,-9,-20,-33,4,9,15),dim=c(3,3)));
计算得到det(A)=-1；
计算得到eigen(A)=...
$vectors
[,1] [,2] [,3]
[1,] -0.4082483-0i -0.4082483+0i -0.4740998+0i
[2,] 0.8164966+0i 0.8164966+0i 0.8127426+0i
[3,] -0.4082483+0i -0.4082483-0i -0.3386427+0i
attach(W);
The following object(s) are masked from 'W (position 3)':
values, vectors
det(vectors);
错误于determinant.matrix(x, logarithm = TRUE, ...) :
目前还不能算复数矩阵的行列式
det(Re(vectors));
[1] -7.599489e-19
solve(vectors)
[,1] [,2] [,3]
[1,] 0.000000+78209959i 0.00000+78209959i -9.26965+78209959i
[2,] 0.000000-78209959i 0.00000-78209959i -9.10153-78209959i
[3,] 3.691206+ 0i 11.07362+ 0i 18.45603+ 0i
很明显vectors不是一个可逆矩阵此时进行特征值分辨这种方法便不可行，对于这种情况我们可以作Schur分解。
描述：
对于任意给定的一个矩阵A（注：此处将"方针"改为"矩阵"以避免歧义），其Schur分解形式可正式表示为：A = U S U^{\dagger}（其中U^{\dagger}代表U的共轭转置），其中U是一个正规交矩阵（即满足U U^{\dagger} = I），而S是一个上三角矩阵，在对角线上的元素即为原始矩阵A的所有特征值）。需要注意的是，在软件包Matrix中定义了该函数，在使用时需加载该软件包，并特别提醒区分大写与小写的变量名matrix与Matrix的不同含义。
例子：
A=Matrix(c(6,12,19,-9,-20,-33,4,9,15),ncol=3,byrow=TRUE);
A;
3 x 3 Matrix of class "dgeMatrix"
[,1] [,2] [,3]
[1,] 6 12 19
[2,] -9 -20 -33
[3,] 4 9 15
library(Matrix);
Sch=Schur(A, vectors=TRUE);
Q=Sch@Q;
Q=as.matrix(Q)
attach(Sch);
错误于attach(Sch) : 'attach'只适用于串列，数据框和环境
Q%%T%%t(Q)
3 x 3 Matrix of class "dgeMatrix"
[,1] [,2] [,3]
[1,] 6 12 19
[2,] -9 -20 -33
[3,] 4 9 15
4_5 Cholesky分解（柯利分解）
描述：
令A为n级实数系数矩阵。如果对于任意非零向量X=(x₁,…,xₙ)，都有t(X)AX>0，则称A为正定矩阵(Positive Definite)。在相合变换下可经合同变换化为标准型：即单位矩阵。
Cholesky 分解：
对于每个正定矩阵 A，在其对应的 Hilbert 空间中都存在唯一一个上三角算子 R 使得 A 等于 R 的转置与 R 的乘积，则称 A 为该算子的 Cholesky 分解（亦称柯尔莫果洛夫分解）。
例子：
#输入矩阵
m=matrix(c(5,1,1,3),ncol=2 );
m;
[,1] [,2]
[1,] 5 1
[2,] 1 3
#矩阵分解
CH=chol(m);
#验证结果
t(CH)%*%CH;
[,1] [,2]
[1,] 5 1
[2,] 1 3
==============以下转自百度百科=================================
===================以下转自百度百科=================================
==============以下转自百度百科=================================
该页面归于百度百科内容库下的专题资源索引模块
该页面包含若干个独立分段落链接
每个分段落均附有明确的标题以及对应的连接路径
其中标题中常见关键词包括'简介'、'起源'、'发展'等
其连接路径统一采用'/link/url'的形式
矩阵乘法
一般而言，在处理矩阵相乘问题时最常用的方法是普通矩陣乘法。需要注意的是，在进行这种运算时必须满足前一矩陣的列數与后一矩陣的行數一致[1] 。当我们仅指普通矩陣乘法时，默认指的是这种运算方式本身。一个m×n型的矩陣实际上是由m×n个元素按照m行排列成n列形成的阵列。简而言之，在数学建模等领域中这种结构能够有效地组织大量数据信息。
定义
设** A** 为
的矩阵，** B** 为
的矩阵，那么称
的矩阵** C** 为矩阵** A** 与** B** 的乘积，记作
，其中矩阵C中的第
行第
列元素可以表示为：
如下所示：
注意事项
当矩阵A的列数等于矩阵B的行数时，A与B可以相乘。
矩阵C的行数等于矩阵A的行数，C的列数等于B的列数。
在乘积矩阵C中位于第m行、第n列的那个数值即为：首先从矩阵A中取出其第m行的所有元素；接着从矩阵B中取出其第n列的所有元素；然后将这两个行列中的对应位置上的数值相乘；最后将所有这些相乘的结果进行累加
基本性质
乘法结合律： (AB)C =A(BC)．[2]
乘法左分配律：(A +B)C =AC +BC[2]
乘法右分配律：** C(A** +B)=CA +CB[2]
对数乘的结合性 k(AB ）=(kA)B =A(kB ）．
转置 (AB)T=B T A T．
矩阵乘法一般不满足交换律[3] 。
Hadamard乘积
矩阵
与
矩阵
的Hadamard积记为
。其元素定义为两个矩阵对应元素的乘积
的 m×n 矩阵[2] 。例如，
Kronecker乘积
Kronecker积是两个任意大小的矩阵间的运算，表示为
也可称作笛卡尔积或直乘积[4] ，并以其创立者德国数学家利奥波德·克罗内克的名字命名。例如，在矩阵运算中可表示为：
实现
C++代码
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | struct Matrix:vector<vector<``int``> >``//使用标准容器vector做基类，需#include语句 { ``Matrix(``int x=0,``int y=0,``int z=0)``//初始化，默认为0行0列空矩阵 ``{ ``assign(x,vector<``int``>(y,z)); ``} ``int h_size()``const``//常量说明不可省，否则编译无法通过 ``{ ``return size(); ``} ``int l_size()``const ``{ ``return empty()?0:front().size();``//列数要考虑空矩阵的情况 ``} ``Matrix ``pow``(``int k);``//矩阵的k次幂，用快速幂实现，k为0时返回此矩阵的单位矩阵 }; Matrix operator*(``const Matrix &m,``const Matrix &n)``//常量引用避免拷贝 { ``if``(m.l_size()!=n.h_size())``return Matrix();``//非法运算返回空矩阵 ``Matrix ans(m.h_size(),n.l_size()); ``for``(``int i=0; i!=ans.h_size(); ++i) ``for``(``int j=0; j!=ans.l_size(); ++j) ``for``(``int k=0; k!=m.l_size(); ++k) ``ans[i][j]+=m[i][k]*n[k][j]; ``return ans; } Matrix Matrix::``pow``(``int k) { ``if``(k==0) ``{ ``Matrix ans(h_size(),h_size()); ``for``(``int i=0; i!=ans.h_size(); ++i) ``ans[i][i]=1; ``return ans; ``} ``if``(k==2)``return (*``this``)*(*``this``); ``if``(k%2)``return pow``(k-1)*(*``this``); ``return pow``(k/2).``pow``(2); } |
|---|
实际应用
数据统计
某公司拥有四个工厂，在各地布局；同时生产三类产品。这些产品的产量（单位为吨）需要通过矩阵进行统计分析。
| 工厂\产品 | P1 | P2 | P3 |
|---|---|---|---|
| 甲 | 5 | 2 | 4 |
| 乙 | 3 | 8 | 2 |
| 丙 | 6 | 0 | 4 |
| 丁 | 0 | 1 | 6 |
可用下述矩阵描述
其中四行分别代表甲、乙、丙、丁四个工厂的生产状况
再设矩阵
，其中第一列表示三种产品的单件利润，第二列表示三种产品的单件体积。
矩阵C的第一列数据各自代表四个工厂的利润,第二列数据各自对应四个工厂产品所需的存储空间
VOJ1067
我们采用上面所述的方法可以快速计算出任意一个线性递推式的第n项。具体而言，在构建其对应的矩阵时，在右上角的小矩阵主对角线上填充1，并将第n行设置为相应的系数值；其余位置则填充0值。例如，在计算满足递推关系式 f(n) = 4f(n-1) - 3f(n-2) + 2f(n-4) 的第k项时，我们可以利用以下形式进行快速运算：
矩阵乘法可用于构造大量相关问题。其中提供的示例特化于所有系数均为1的情形。
在考虑一个有向图的情况下，请问从A点出发恰好沿着k条边行走（允许沿边多次通过）到达B点的所有路径的数量模p的结果是多少。
将输入的图转换为对应的邻接矩阵表示A。其中C = A \times A中C(i,j) = \sum A(i,k) \times A(k,j)实际上代表了节点i到节点j之间恰好经过两条边的不同路径数目（其中k作为中间节点）。类似地，在计算C \times A时，第i行第j列元素表示从i到j经过三条边的不同路径数量。以此类推，在计算更高次幂时，则能够得到所有节点之间的k步可达性信息。通过递推或快速幂算法计算A^k即可得到所有节点之间经过k-1条边的所有可能路径数目。
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 | #include <cmath> #include <cstdio> #include <cstring> #include <iostream> #include <algorithm> #define N 10 using namespace std; const int mod = 7777777; typedef long long LL; struct matrix{ ``LL a[10][10]; }origin; int n,m; matrix multiply(matrix x,matrix y) { ``matrix temp; ``memset``(temp.a,0,``sizeof``(temp.a)); ``for (``int i=0;i<n;i++) ``{ ``for (``int j=0;j<n;j++) ``{ ``for (``int k=0;k<n;k++) ``{ ``temp.a[i][j]+=x.a[i][k]*y.a[k][j]; ``temp.a[i][j]=(temp.a[i][j])%mod; ``} ``} ``} ``return temp; } matrix matmod(matrix A,``int k) { ``matrix res; ``memset``(res.a,0,``sizeof res.a); ``for (``int i=0;i<n;i++) res.a[i][i]=1; ``while``(k) ``{ ``if (k&1) res=multiply(res,A); ``k>>=1; ``A=multiply(A,A); ``} ``return res; } void print(matrix x) { ``for (``int i=0;i<n;i++) ``{ ``for (``int j=0;j<n;j++) ``cout<<``" "``<<x.a[i][j];``puts``(``""``); ``} ``printf``(``"---------------\n"``); } int main() { ``int k; ``while (cin>>n>>k) ``{ ``memset``(origin.a,0,``sizeof origin.a); ``origin.a[0][0]=1; ``for (``int i=1;i<=n;i++) ``{ ``origin.a[i][0]=1; ``for (``int j=0;j<i;j++) ``origin.a[i][0]+=origin.a[j][0]; ``} ``// print(origin); ``matrix res; ``memset``(res.a,0,``sizeof res.a); ``for (``int i=0;i<n-1;i++) ``res.a[i][i+1]=1; ``for (``int i=0;i<n;i++) res.a[n-1][i]=1; ``//print(res); ``res=matmod(res,k-1); ``LL fans=0; ``for (``int i=0;i<n;i++) ``{ ``fans+=res.a[0][i]*origin.a[i][0]; ``fans%=mod; ``} ``cout<<fans<<endl; ``} ``return 0; } |
|---|
经典题目9
用1 x 2的多米诺骨牌填满M x N的矩形有多少种方案，M<=5，N<2^31，输出答案mod p的结果
举例来说，在M等于3的情况下展开说明
后来编写了一份这个题目中的源代码，你可以再次查看位运算的相关应用。
经典题目10
POJ2778
该问题要求计算所有合法长度为n位的DNA序列中不含给定病毒片段的数量。其中合法DNA仅由A、T、C、G四种碱基构成，在本题中将提供最多10个病毒片段（每个片段长度不超过10）。所求解的数据规模范围可高达N=2,000,000,000位碱基对。
将所有病毒片段划分为以AT结尾、AA结尾等7类：以AT结尾、以AA结尾、以GG结尾、以?A结尾（其中问号代表“其他情况”，即任一字母但不会导致成为特定前缀）、以?G结尾、以?C结尾以及以??结尾（双问号表示任一两个字母组合）。这些分类构成了全集的一个划分（各分类之间无交集且并集为全集）。已知长度为n-1时各类DNA满足条件的数量，则需计算长度为n时各类DNA的数量变化规律。通过建立各类型间的转移关系可构造一个带权有向图：例如，在AT状态下无法转移到AA；从AT状态可转移到??状态有4种可能（后面跟任一字母）；从?A状态转移到AA有1种方式（后面跟A），而转移到??状态则有2种方式（后面跟G或C）；GG状态下只能转移到AA或TT等非违禁组合；该过程类似于有限状态自动机用于模式匹配。最后通过构建转移矩阵并进行n次自乘运算即可求得结果：最终输出值为从??状态到其他各状态的所有路径数量之和
题目中的数据规模限定前缀数不得超过100, 每次矩阵乘法涉及三个矩阵的操作, 总共需要进行log(n)次这样的运算. 因此可知该问题的时间复杂度计算公式为100^3 \times \log n.
最终提供第9题的代码作为参考（最近学习了C++中的类及其运算符重载）。为了帮助大家更好地理解这些内容，请将此提示语放置在前面。
Matrix67原创，转贴请注明出处。[5]
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 | #include <cstdio> #define SIZE (1<<m) #define MAX_SIZE 32 using namespace std; class CMatrix { ``public``: ``long element[MAX_SIZE][MAX_SIZE]; ``void setSize(``int``); ``void setModulo(``int``); ``CMatrix operator* (CMatrix); ``CMatrix power(``int``); ``private``: ``int size; ``long modulo; }; void CMatrix::setSize(``int a) { ``for (``int i=0; i<a; i++) ``for (``int j=0; j<a; j++) ``element[i][j]=0; ``size = a; } void CMatrix::setModulo(``int a) { ``modulo = a; } CMatrix CMatrix::operator* (CMatrix param) { ``CMatrix product; ``product.setSize(size); ``product.setModulo(modulo); ``for (``int i=0; i<size; i++) ``for (``int j=0; j<size; j++) ``for (``int k=0; k<size; k++) ``{ ``product.element[i][j]+=element[i][k]*param.element[k][j]; ``product.element[i][j]%=modulo; ``} ``return product; } CMatrix CMatrix::power(``int exp``) { ``CMatrix tmp=(*``this``)*(*``this``); ``if (``exp``==1) ``return *``this``; ``else if (``exp``&1) ``return tmp.power(``exp``/2)*(*``this``); ``else return tmp.power(``exp``/2); } int main() { ``const int validSet[]={0,3,6,12,15,24,27,30}; ``long n, m, p; ``CMatrix unit; ``scanf``(``"%d%d%d"``, &n, &m, &p); ``unit.setSize(SIZE); ``for (``int i=0; i<SIZE; i++) ``for (``int j=0; j<SIZE; j++) ``if (((~i)&j) == ((~i)&(SIZE-1))) ``{ ``bool isValid=``false``; ``for (``int `k=0;k<8;k++) isValid=isValid |
(i&j)==validSet[k]; unit.element[i][j]=isValid;` ` } unit.setModulo(p);` ` printf("%d",unit.power(n).element[SIZE-1][SIZE-1] );` ` return 0; }` |
|---|
矩阵乘法例题
vijos1049
题目的大致意思是提供给你一个物品变换的顺序表，请你计算经过n次变换后物品的位置。
解析：这个题目仔细思考其实并不困难。考虑到每行变换的具体顺序不同，在这种情况下我们需要对给定的矩阵每行进行相应的转换，并将其表示为一种矩阵运算。但是直接进行每次这样的操作会导致时间复杂度过高。因此我们可以将各次操作对应的转换矩阵相乘得到一个新的总转换矩阵。这相当于完成k次这样的转换后的总效果对应的转换。这样我们就可以利用快速幂算法来高效计算所需的结果。而对于求余的情况则可以通过直接计算的方式来处理。
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 | #include <cstdio> #include <cstring> #include <algorithm> using namespace std; static const int maxm=1e2+10; #define REP(i,s,t) for(int i=s;i<=t;i++) typedef long long LL; struct matrix{ ``LL mtx[maxm][maxm]; }mx[16]; LL n,k,m; LL A[maxm][maxm]; matrix mul(matrix A,matrix B){ ``matrix ret; ``memset``(ret.mtx,0,``sizeof ret.mtx); ``REP(i,1,n)REP(j,1,n)REP(k,1,n) ``ret.mtx[i][j]=(ret.mtx[i][j]+A.mtx[i][k]*B.mtx[k][j]); ``return ret; } matrix ``pow``(matrix A,LL n){ ``if``(!n)``return A; ``matrix ret=A;n--; ``while``(n){ ``if``(n&1)ret=mul(ret,A); ``A=mul(A,A); ``n>>=1; ``} ``return ret; } void display(matrix base){ ``for``(``int i=1;i<=n;i++)``printf``(``"%lld "``,base.mtx[1][i]); ``puts``(``""``); } int main(){ ``matrix st,get,f; ``scanf``(``"%lld%lld%lld"``,&n,&m,&k); ``for``(``int i=1;i<=m;i++){ ``for``(``int j=1;j<=n;j++){ ``scanf``(``"%lld"``,&A[i][j]); ``mx[i].mtx[A[i][j]][j]=1; ``} ``} ``for``(``int i=1;i<=n;i++)st.mtx[1][i]=i; ``get=mx[1]; ``for``(``int i=2;i<=m;i++)get=mul(get,mx[i]); ``get=``pow``(get,k/m);k%=m; ``for``(``int i=1;i<=k;i++)get=mul(get,mx[i]); ``st=mul(st,get); ``display(st); ``return 0; } //by Exbilar |
|---|