R语言中矩阵运算

最新推荐文章于 2023-12-09 16:14:13 发布

最古琴

最新推荐文章于 2023-12-09 16:14:13 发布

阅读量2.4k

点赞数

wwhl.wang

本文链接：https://blog.csdn.net/gdzjsubaoya/article/details/99900315

版权

撸了今年阿里、头条和美团的面试，我有一个重要发现.......>>>

推荐一款免费的自购省钱，分享赚钱的平台——赚赚熊，感兴趣的，特别是家里有婆娘专职带娃的，请看文末。

目录：
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

注意：生成矩阵是按列排列的。

1_2用array ( )函数构造多维数组
      用法为：array(data=NA,dim=length(data),dimnames=NULL)
            参数描述：data：是一个向量数据。
                      dim：是数组各维的长度，缺省时为原向量的长度。
                      dimname：是数组维的名字，缺省时为空。
例子：
> x=array(1:20,dim=c(4,5))
> x
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20

1_3用matrix()函数构造矩阵
    函数matrix)是构造矩阵(二维数组)的函数，其构造形式为
    matrix(data=NA，nrow=1，ncol=1，byrow=FALSE，dimnames=NULL)
    其中data是一个向量数据，nro、是矩阵的行数，ncol是矩阵的列数.当byrow=TRUE时，生成矩阵的数据按行放置，缺省时相当于byrow=FALSE，数据按列放置.dimname。是数组维的名字，缺省时为空.
    如构造一个3x5阶的矩阵
> A=matrix(1:15,nrow=3,byrow=TRUE)
> A
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15

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=2*1: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()是更为强大的外积运算函数，outer(x,y)计算向量二与y的外积，它等价于x %o%y
函数。outer()的一般调用格式为
      outer(x，y，fun=”*”)
     其中x, y矩阵(或向量)，fun是作外积运算函数，缺省值为乘法运算。函数outer()在绘制三维曲面时非常有用，它可生成一个x和y的网格。

3_5 矩阵的乘法
设A和B为两个矩阵，通常意义下的矩阵乘法是通过A%*%B来完成，crossprod(A,B)表示的是
t(A)%*%B，而tcrossprod(A,B)表示的是A%*%t(B)。最后我们通过运算知道x%*%A%*%x为二次型。

例子：
> A=array(1:9,dim=(c(3,3)))
> B=array(9:1,dim=(c(3,3)))
> A%*%B;
     [,1] [,2] [,3]
[1,]   90   54   18
[2,]  114   69   24
[3,]  138   84   30
> crossprod(A,B)==t(A)%*%B;
     [,1] [,2] [,3]
[1,] TRUE TRUE TRUE
[2,] TRUE TRUE TRUE
[3,] TRUE TRUE TRUE
> tcrossprod(A,B)==A%*%t(B);
     [,1] [,2] [,3]
[1,] TRUE TRUE TRUE
[2,] TRUE TRUE TRUE
[3,] TRUE TRUE TRUE

3_6 生成对角阵和矩阵取对角运算
    函数diag()依赖于它的变量，当v是一个向量时，diag(v)表示以v的元素为对角线元素的对角阵.当M是一个矩阵时，则diag(M)表示的是取M对角线上的元素的向量.如
> v=c(1,4,5);
> diag(v);
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    4    0
[3,]    0    0    5
> M=array(1:9,dim=c(3,3));
> diag(M);
[1] 1 5 9

3_7 解线性方程组和求矩阵的逆矩阵
若求解线性方程组Ax=b，其命令形式为solve(A,b)，求矩阵A的逆，其命令形式为solve(A).设矩阵A=t(array(c(1:8,10),dim=c(3,3))),b<-c(1,1,1),则解方程组Ax=b的解x和求矩阵A的逆矩阵的命令如下:

> 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

3_8 求矩阵的特征值与特征向量
函数eigen(Sm)是求对称矩阵Sm的特征值与特征向量，其命令形式为:ev=eigen(Sm),则ev存放着对称矩阵Sm特征值和特征向量，是由列表形式给出的，其中ev$values是Sm的特征值构成的向量，ev$vectors是Sm的特征向量构成的矩阵.如
> Sm=crossprod(A,A);
> ev=eigen(Sm);
> ev;
$values
[1] 303.19533618 0.76590739 0.03875643

$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阶方阵A，x为标量，v是非零的N维列向量，且满足Ax=xv ，则称x为矩阵A的特征值，v 是相对应于x 的特征向量。特征值的全体成为A的谱。

(2).在r中的实现：在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).特征值分解的性质：我们知道当所求的的特征向量构成的矩阵可逆时会满足solve(vectors)%*%A%*%vectors=diag(values),下面进行验证。

> 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 矩阵的奇异值分解
函数svd(A)是对矩阵A作奇异值分解，即A =U%*%D%*%t(V)，其中U, V是正交阵，D为对角阵，也就是矩阵A的奇异值.svd(A)的返回值也是列表，svd(A)$d表示矩阵A的奇异值，即矩阵D的对角线上的元素.svd(A)$u对应的是正交阵U, svd(A) $v对应的是正交阵V.例如，
> A<-t(array(c(1:8,10),dim=c(3,3)))
> SVD=svd(A);
> SVD;
$d
[1] 17.4125052 0.8751614 0.1968665

$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为m*n矩阵，如果存在m*m酉矩阵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

4_4 Schur分解
引言：
    从特征值的分解中可以看出，特征值的分解是有条件的，如果特征向量不是线性无关的，那么对于一个矩阵来说便不能采用特征值分解的方法对矩阵进行分解。例如对于矩阵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)));
> A;
     [,1] [,2] [,3]
[1,]    6   12   19
[2,]   -9  -20  -33
[3,]    4    9   15
> det(A);
[1] -1
> W=eigen(A);
> W;
$values
[1]  1+0i  1-0i -1+0i

$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=USU(H)，其中U是标准的正交矩阵（即满足UU(H)=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=(x1,...xn) 都有
t(X)AX>0，就称A正定(Positive Definite)。正定矩阵在相合变换下可化为标准型，即单位矩阵。

Cholesky分解：
对任意的正定矩阵A，存在上三角矩阵R，使A=t(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

==============以下转自百度百科=================================

===================以下转自百度百科=================================

==============以下转自百度百科=================================

http://baike.baidu.com/link?url=18IT0KMVS7e94kG9o_YNUr_SAmsCFFkKRpsQtdjcG92nuGLI7qrRO-K2XQnwFb2o0TNWGlgjGu-p8wOeLCO3TWKcYSSqjndZozxTUk6F2LFkRDD9cqYrzrYWEEFTlaGn

矩阵乘法

矩阵相乘最重要的方法是一般矩阵乘积。它只有在第一个矩阵的列数（column）和第二个矩阵的行数（row）相同时才有意义[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=BTAT．
矩阵乘法一般不满足交换律[3] 。

Hadamard乘积

矩阵

与

矩阵

的Hadamard积记为

。其元素定义为两个矩阵对应元素的乘积

的m×n矩阵[2] 。例如，

Kronecker乘积

Kronecker积是两个任意大小的矩阵间的运算，表示为

。克罗内克积也成为直积或张量积[4] .以德国数学家利奥波德·克罗内克命名。计算过程如下例所示：

实现

C++代码

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);

}

实际应用

数据统计

某公司有四个工厂，分布在不同地区，同时三种产品，产量（单位；t），试用矩阵统计这些数据。

工厂\产品	P1	P2	P3
甲	5	2	4
乙	3	8	2
丙	6	0	4
丁	0	1	6

可用下述矩阵描述

，其中四行分别表示甲乙丙丁四个工厂的生产情况，三列分布表示三种产品P1，P2，P3的产量。

再设矩阵

，其中第一列表示三种产品的单件利润，第二列表示三种产品的单件体积。

矩阵C的第一列数据分别表示四个工厂的利润，第二列分别表示四个工厂产品需要的存储空间。

VOJ1067

我们可以用上面的方法二分求出任何一个线性递推式的第n项，其对应矩阵的构造方法为：在右上角的(n-1)*(n-1)的小矩阵中的主对角线上填1，矩阵第n行填对应的系数，其它地方都填0。例如，我们可以用下面的矩阵乘法来二分计算f(n) = 4f(n-1) - 3f(n-2) + 2f(n-4)的第k项：

利用矩阵乘法求解线性递推关系的题目我能编出一卡车来。这里给出的例题是系数全为1的情况。

给定一个有向图，问从A点恰好走k步（允许重复经过边）到达B点的方案数mod p的值

把给定的图转为邻接矩阵，即A(i,j)=1当且仅当存在一条边i->j。令C=A*A，那么C(i,j)=ΣA(i,k)*A(k,j)，实际上就等于从点i到点j恰好经过2条边的路径数（枚举k为中转点）。类似地，C*A的第i行第j列就表示从i到j经过3条边的路径数。同理，如果要求经过k步的路径数，我们只需要二分求出A^k即可。

#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为例进行讲解。假设我们把这个矩形横着放在电脑屏幕上，从右往左一列一列地进行填充。其中前n-2列已经填满了，第n-1列参差不齐。现在我们要做的事情是把第n-1列也填满，将状态转移到第n列上去。由于第n-1列的状态不一样（有8种不同的状态），因此我们需要分情况进行讨论。在图中，我把转移前8种不同的状态放在左边，转移后8种不同的状态放在右边，左边的某种状态可以转移到右边的某种状态就在它们之间连一根线。注意为了保证方案不重复，状态转移时我们不允许在第n-1列竖着放一个多米诺骨牌（例如左边第2种状态不能转移到右边第4种状态），否则这将与另一种转移前的状态重复。把这8种状态的转移关系画成一个有向图，那么问题就变成了这样：从状态111出发，恰好经过n步回到这个状态有多少种方案。比如，n=2时有3种方案，111->011->111、111->110->111和111->000->111，这与用多米诺骨牌覆盖3x2矩形的方案一一对应。这样这个题目就转化为了我们前面的例题8。

后面我写了一份此题的源代码。你可以再次看到位运算的相关应用。

经典题目10

POJ2778

题目大意是，检测所有可能的n位DNA串有多少个DNA串中不含有指定的病毒片段。合法的DNA只能由ACTG四个字符构成。题目将给出10个以内的病毒片段，每个片段长度不超过10。数据规模n<=2 000 000 000。

下面的讲解中我们以ATC,AAA,GGC,CT这四个病毒片段为例，说明怎样像上面的题一样通过构图将问题转化为例题8。我们找出所有病毒片段的前缀，把n位DNA分为以下7类：以AT结尾、以AA结尾、以GG结尾、以?A结尾、以?G结尾、以?C结尾和以??结尾。其中问号表示“其它情况”，它可以是任一字母，只要这个字母不会让它所在的串成为某个病毒的前缀。显然，这些分类是全集的一个划分（交集为空，并集为全集）。现在，假如我们已经知道了长度为n-1的各类DNA中符合要求的DNA个数，我们需要求出长度为n时各类DNA的个数。我们可以根据各类型间的转移构造一个边上带权的有向图。例如，从AT不能转移到AA，从AT转移到??有4种方法（后面加任一字母），从?A转移到AA有1种方案（后面加个A），从?A转移到??有2种方案（后面加G或C），从GG到??有2种方案（后面加C将构成病毒片段，不合法，只能加A和T）等等。这个图的构造过程类似于用有限状态自动机做串匹配。然后，我们就把这个图转化成矩阵，让这个矩阵自乘n次即可。最后输出的是从??状态到所有其它状态的路径数总和。

题目中的数据规模保证前缀数不超过100，一次矩阵乘法是三方的，一共要乘log(n)次。因此这题总的复杂度是100^3 * log(n)，AC了。

最后给出第9题的代码供大家参考（今天写的，熟悉了一下C++的类和运算符重载）。为了避免大家看代码看着看着就忘了，我把这句话放在前面来说：

Matrix67原创，转贴请注明出处。[5]

#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次(k是所给矩阵的行数)所乘的矩阵,那么我们就可以使用快速幂了,对于余数就暴力算就可以啦.

#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

http://baike.baidu.com/link?url=18IT0KMVS7e94kG9o_YNUr_SAmsCFFkKRpsQtdjcG92nuGLI7qrRO-K2XQnwFb2o0TNWGlgjGu-p8wOeLCO3TWKcYSSqjndZozxTUk6F2LFkRDD9cqYrzrYWEEFTlaGn

----------------------下方为免费的广告---------------------

1.识别下图中的二维码，跳到步骤2的界面；