矩阵乘法

2019-04-12 15:47发布

矩阵乘法是个很高端的东西。
注意事项: 以下的讲述下标都从0开始

矩阵是什么?

不解释

矩阵加减法?

矩阵加矩阵,就将对应的加上。
矩阵加常数,就将每一个元素都加上这个常数。
减法同理。

矩阵乘法

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

一个mn的的A矩阵,和一个np的B矩阵相乘,将得到一个mp的矩阵C
C(i,j)=k=1nA(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 operator*(const Matrix& x,const Matrix& y) { Matrix ret; int i,j,k; for (i=0;ifor (j=0;jfor (k=0;kreturn 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=[f(n2)f(n1)]
现在我们要通过A推出B
B=[f(n1)f(n)]=[f(n1)f(n2)+f(n1)]
我们要求出一个2*2的 转移矩阵 T,满足
AT=B

[f(n2)f(n1)]T=[f(n1)f(n2)+f(n1)]
我们知道,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=[0111]
显然这个式子是正确的。
[f(n2)f(n1)][0111]=[f(n1)f(n2)+f(n1)]
求出这个矩阵后,就可以愉快地快速幂了。
时间复杂度O(lgN)

代码

#include #include #include 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 operator*(const Matrix& x,const Matrix& y) { Matrix ret; int i,j,k; for (i=0;ifor (j=0;jfor (k=0;kreturn 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 ",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)θx=(xa)cosθ(yb)sinθ+ay=(xa)si
登录 后发表评论
0条评论
还没有人评论过~