一种混合负二进制编码的光学矩阵乘法系统

一种混合负二进制编码的光学矩阵乘法系统（完整版）实用资料(可以直接使用，可编辑完整版实用资料，欢迎下载）

第5卷第3期2007年6月一种混合负二进制编码的光学矩阵乘法系统（完整版）实用资料(可以直接使用，可编辑完整版实用资料，欢迎下载）光学与光电技术OPTICS&OPTOELECTRONICTECHNOLOGYVol.5,No.3June,2007收稿日期2006207226;收到修改稿日期2006210213文章编号:167223392(20070320066203一种混合负二进制编码的光学矩阵乘法系统张锐1,2李修建1杨建坤1胡文华1(1国防科学技术大学技术物理研究所,长沙410073;266222,摘要将数学中的编码方法与光学原理相结合,,实现矩阵2,分析了其工作性能。研究表明,,,对光计算;;光折变晶体O438A1引言在电子计算机中许多复杂的运算都可以分解为乘法和卷积运算,用光学的方法来高效地实现。光学矩阵乘法一直是光计算领域中的一个很重要的部分。早期的一些研究主要都是进行光学模拟式的矩阵乘法系统,这类系统大都精度低。例如,Heinz最早提出用相干空间滤波技术实现矩阵乘法[1],由于光输出时损失太大,它的技术很难实现大规模的矩阵乘法。后来由Schneider提出一种由球面柱面透镜组与光孔阵列组合成的非相干系统,当矩阵维数增加时,衍射现象限制了系统的精度[2]。为了提高其精度,引入了许多的编码方法,实现了数字式的矩阵乘法系统。例如,由D.saltis提出的二进制乘法转化技术[3],周常河教授提出二进制补码编码技术[4],刘立人教授提出运用混合负二进制编码实现向量矩阵乘法[5]等。本文跟踪目前国际上的研究热点———光折变晶体的研究[6],利用其光学性质,提出了一种运用数学编码技术实现混合负二进制可编码的矩阵乘法的光学系统。2混合负二进制编码光学矩阵乘法系统2.1编码技术在光学矩阵乘法研究领域中,编码问题一直是人们研究的热点。光强编码采用数字编码数制有:二进制、二进制补码、混合负二进制等。二进制编码是利用计算机算法的思路衍变过来的一种编码方法,利用这种方法来实现光学矩阵乘法运算,必须考虑进位和反馈回路等问题,不但结构复杂而且误码率高。在此基础上提出的用二进制补码编码技术实现矩阵乘法运算,无需进位。但是由于此类系统对透镜的尺寸要求很高,因此很难实现大规模的矩阵乘法。结合以上两种编码方法的优缺点,人们又提出了混合负二进制编码。任何一个十进制的数a都可以表示为a=∑Mn=-Nan(-2n(an是混合负二进制中的第n位数,N和M都是整数。这种编码算法无需进位,无需正负号,既能表示实数,也能表示复数。更重要的是,它便于光路的实现,是矩阵乘法系统编码的最佳选择。2.2光折变晶体的光学特性光折变晶体在光辐射下由光强的空间分布会引起材料折射率的相关变化,典型的应用就是相位共轭装置[7]。当三束相干的入射光同时照射在光折变晶体中某点时,会产生第四束相干光波,这也就是四波混频效应,其原理如图1如所示。其中k1和k2是同频率、同强度、方向相反的写入光,k3是与之同频但强度略小的读出光。当k3以特定的角度入射时,就会同时出现反射光k4和共轭光k5,最后经半透半反镜得到输出光k6。只要当k1、k2、k3其中一个不存在时,就不会有输出光产生。第3期张锐等:一种混合负二进制编码的光学矩阵乘法系统图1相位共轭器原理图Fig.1Schematicdiagramofphaseconjugatesystem2.3光学矩阵乘法系统在矩阵2矩阵乘法系统中实现两个矩阵A的相乘。上,,成相乘,。若A=12…a1qa21a22…a2q……am1am2…am×qB=b11b12…b1qb21b22…b2q……bn1bn2…bq×n则C=AB=c11c12…c1nc21c22…c2n……cm1cm2…cm×n其中,矩阵C中每个元素cik=ai1b1k+ai2b2k+…+aikbik(i=1,2,3,…,p;k=1,2,3,…,n。当A和B中的元素都用混合负二进制表示时,所得到的结果也必定是混合二进制的。现以A=521-,B=-10-74-为例,为了实现两个矩阵的乘法,首先将矩阵的信息分别加载到两个大小相等的矩阵平面上,如图2所示。再对矩阵A和B中元素进行混合负二进制编码。同时将编码后的信息按图3、4的形式重复排列,其中1代表发光,0和空格代表不发光。利用这样的编排方法,只要将两个矩阵面中的矩阵元素对应相乘即可,无需考虑进位问题。将光折变晶体相对平行地置于两个矩阵面正中间,且其中元素与矩阵面的元素一一对应,最后经透镜组传输到接收平面。系统原理如图5所示。a11a11a12a12a21a21a22a22b11b11b21b22b11b12b21b22图2矩阵A和矩阵BArrangementof图3矩阵A混合负二进制编码图Fig.3NegativebinaryofmatrixA图4矩阵B混合负二进制编码图Fig.4NegativebinaryofmatrixB图5混合负二进制编码矩阵乘法系统原理图Fig.5Schematicdiagramofopticalmatrixmultiplicationschemewithnegativebinaryencoding若在晶体某个位置同时有方向相反频率相同的写入光时,则当读出光以合适角度照射时,会出现共轭光。以a11b11+a12b21=c11为例,若两个矩阵面的相互对应的点同时为1时,则在这两点中间的光折变晶体面的对应位置上会有读出光的共轭光束产生。例如第2行第4、6列,第4行第2、4列,76光学与光电技术第5卷第7行第3、4列位置。同一列的8束共轭光经过柱面镜后会聚为一点,在光折变晶体C11位置的7列就会在接收面C11位置处得到7个点的光强分布信息,如图6所示。第2、3、6个点的光强基本相同,而第4个点的光强是它们的3倍,因此光强信息代表混合负二进制0113010,转化为十进制为-42,且符合乘法规则。同理,可以得到接收面其他各点的光强信息。所以此光学系统验证式(1成立,实现了矩阵2矩阵乘法。C=A×B=521--10-74十进制-42-39-428(1图6C11处的光强分布Fig.6ArrangementofopticalintensityonC113结论利用了光折变晶体的特性与混合负二进制编码的方法,提出了一种矩阵乘法系统。不但可以实现实数的编码而且可以实现复数的编码,无需光电转换即可实现高度并行的矩阵2矩阵乘法运算,是一种全光学光计算系统。由于矩阵面编码信息简单,便于扩展,可用于实现大规模矩阵乘法运算。但是仍有许多待解决的问题,最重要的就是工艺实现。要产生共轭光,读出光必须在合适的角度入射,入射角度和位置需做精确调整。另外,两矩阵。、“非”等逻辑运。参考文献[1]RAHeinz,JOArtman,SHLee.Matrixmultipli2cationbyopticalmethods[J].Appl.Opt.,1970(9:216122164.[2]WSchneider,WFink.Incoherentopticalmatrixmultiplication[J].Opt.Acta.,1975(22:8792883.[3]DPsaltis,DCasasent.Accuratenumericalomputa2tionbyopticalconvolutionphoto2opt[J].Instrum.Eng.,1980,(232:1512153.[4]ChangheZhou,LirenLiu,ZhijiangWang.Binary2enco2dedvector2matrixmultiplica2tionarchitecture[J].Opt.Let.,1992,(17:180021802.[5]LirenLiu,GuoqiangLi,YaozuYin.Opticalcomplexmatrix2vectormult2iplication[J].Opt.Let.,1994,(19:175921761.[6]AKDas,SMukhopaadhyay.All2opticaltime2domainmultiplexingde2multiplexingschemewithnonlinearma2terial[J].Opt.Eng.,2005,(44:2021203[7]ArchanKumarDas,SourangshuMukhopadhyay.Anall2opticalmatrixmultiplicationschemewithnon2line2armaterialbasedswithchingsystem[J].ChineseOpt.Let.,2005,(30:1722175.OpticalMatrixMultiplicationSystemWithNegativeBinaryEncodingZHANGRui1,2LIXiu2jian1YANGJian2kun1HUWen2hua1(1Tech2PhysicalResearchCenter,NationalUniversityofDefenceTechnology,Changsha410073,China;266222ArmyofBeijingMilitaryCommand,Beijing102202,ChinaAbstractInthispaperweapplytheencodinginmathematicstotheopticalprincipleandproposeanewopticalmatrixmulti2plicationsystemwithnegativebinaryencodingbyusingfour2wavemixingofphotorefractivecrystal.Itsprincipleisintro2ducedanditsoperatingperformanceisanalyzed.Theresearchshowsthatthematrixencodingofthesystemiseasyanditcanrealizethematrixmultiplicationcosmically.Anditisusefulfortheresearchonopticalcalculation.Keywordsopticalcalculation;opticalmatrixmultiplication;negativebinary;photorefractivecrystal86fibonacci斐波那契数列(矩阵复习）2021-11-0221:03描述Description著名的斐波那契数列

f[n]=1n=1,2

f[n-1]+f[n-2]n>2

现在求第n项，由于f[n]可能很大，你只需要输出mod32768的值即可。输入格式InputFormat仅有一行，为n,1<=n<=maxlongint输出格式OutputFormat仅有一个数字，即答案样例输入SampleInput3样例输出SampleOutput2时间限制TimeLimitation各个测试点1s注释Hint各个测试点1s===================================================================================var

i,j,k,m,n:longint;

a,c:array[1..1,1..2]oflongint;

b,cc:array[1..2,1..2]oflongint;

begin

a[1,1]:=1;

a[1,2]:=1;

b[1,1]:=0;

b[1,2]:=1;

b[2,1]:=1;

b[2,2]:=1;

readln(n);

dec(n,2);

whilen>0do

begin

ifodd(n)then

begin

fillchar(c,sizeof(c),0);

fork:=1to2do

fori:=1to1do

forj:=1to2do

c[i,j]:=(c[i,j]+a[i,k]*b[k,j])and32767;

a:=c;

end;

fillchar(cc,sizeof(cc),0);

fork:=1to2do

fori:=1to2do

forj:=1to2do

cc[i,j]:=(cc[i,j]+b[i,k]*b[k,j])and32767;

b:=cc;

n:=nshr1;

end;

writeln(a[1,2]);

end.矩阵乘法百科名片矩阵乘法是一种高效的算法可以把一些一维递推优化到log（n)，还可以求路径方案等，所以更是是一种应用性极强的算法。矩阵，是线性代数中的基本概念之一。一个m×n的矩阵就是m×n个数排成m行n列的一个数阵。由于它把许多数据紧凑的集中到了一起，所以有时候可以简便地表示一些复杂的模型。矩阵乘法看起来很奇怪，但实际上非常有用，应用也十分广泛。目录基本定义它是这样定义的，只有当矩阵A的列数与矩阵B的行数相等时A×B才有意义。一个m×n的矩阵a（m,n)左乘一个n×p的矩阵b（n,p)，会得到一个m×p的矩阵c（m,p)，满足矩阵乘法满足结合率，但不满足交换率一般的矩乘要结合快速幂才有效果``一般矩乘的代码：functionmul(a,b:Tmatrix):Tmatrix;vari,j,k:longint;c:Tmatrix;beginfillchar(c,sizeof(c),0);fork:=0tondofori:=0tomdoforj:=0topdobegininc(c[i,j],a[i,k]*b[k,j]);ifc[i,j]>rathenc[i,j]:=c[i,j]modra;end;mul:=c;end;编辑本段经典题目好像目前还没有这方面题目的总结。这几天连续看到四个问这类题目的人，今天在这里简单写一下。这里我们不介绍其它有关矩阵的知识，只介绍矩阵乘法和相关性质。不要以为数学中的矩阵也是黑色屏幕上不断变化的绿色字符。在数学中，一个矩阵说穿了就是一个二维数组。一个n行m列的矩阵可以乘以一个m行p列的矩阵，得到的结果是一个n行p列的矩阵，其中的第i行第j列位置上的数等于前一个矩阵第i行上的m个数与后一个矩阵第j列上的m个数对应相乘后所有m个乘积的和。比如，下面的算式表示一个2行2列的矩阵乘以2行3列的矩阵，其结果是一个2行3列的矩阵。其中，结果的那个4等于2*2+0*1：下面的算式则是一个1x3的矩阵乘以3x2的矩阵，得到一个1x2的矩阵：矩阵乘法的两个重要性质：一，矩阵乘法不满足交换律；二，矩阵乘法满足结合律。为什么矩阵乘法不满足交换律呢？废话，交换过来后两个矩阵有可能根本不能相乘。为什么它又满足结合律呢？仔细想想你会发现这也是废话。假设你有三个矩阵A、B、C，那么(AB)C和A(BC)的结果的第i行第j列上的数都等于所有A(ik)*B(kl)*C(lj)的和（枚举所有的k和l）。经典题目1给定n个点，m个操作，构造O(m+n)的算法输出m个操作后各点的位置。操作有平移、缩放、翻转和旋转这里的操作是对所有点同时进行的。其中翻转是以坐标轴为对称轴进行翻转（两种情况），旋转则以原点为中心。如果对每个点分别进行模拟，那么m个操作总共耗时O(mn)。利用矩阵乘法可以在O(m)的时间里把所有操作合并为一个矩阵，然后每个点与该矩阵相乘即可直接得出最终该点的位置，总共耗时O(m+n)。假设初始时某个点的坐标为x和y，下面5个矩阵可以分别对其进行平移、旋转、翻转和旋转操作。预先把所有m个操作所对应的矩阵全部乘起来，再乘以(x,y,1)，即可一步得出最终点的位置。经典题目2给定矩阵A，请快速计算出A^n（n个A相乘）的结果，输出的每个数都modp。由于矩阵乘法具有结合律，因此A^4=A*A*A*A=(A*A)*(A*A)=A^2*A^2。我们可以得到这样的结论：当n为偶数时，A^n=A^(n/2)*A^(n/2)；当n为奇数时，A^n=A^(n/2)*A^(n/2)*A（其中n/2取整）。这就告诉我们，计算A^n也可以使用二分快速求幂的方法。例如，为了算出A^25的值，我们只需要递归地计算出A^12、A^6、A^3的值即可。根据这里的一些结果，我们可以在计算过程中不断取模，避免高精度运算。经典题目3POJ3233(感谢rmq)题目大意：给定矩阵A，求A+A^2+A^3+...+A^k的结果（两个矩阵相加就是对应位置分别相加）。输出的数据modm。k<=10^9。这道题两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：A+A^2+A^3+A^4+A^5+A^6=(A+A^2+A^3)+A^3*(A+A^2+A^3)应用这个式子后，规模k减小了一半。我们二分求出A^3后再递归地计算A+A^2+A^3，即可得到原问题的答案。经典题目4VOJ1049题目大意：顺次给出m个置换，反复使用这m个置换对初始序列进行操作，问k次置换后的序列。m<=10,k<2^31。首先将这m个置换“合并”起来（算出这m个置换的乘积），然后接下来我们需要执行这个置换k/m次（取整，若有余数则剩下几步模拟即可）。注意任意一个置换都可以表示成矩阵的形式。例如，将1234置换为3124，相当于下面的矩阵乘法：置换k/m次就相当于在前面乘以k/m个这样的矩阵。我们可以二分计算出该矩阵的k/m次方，再乘以初始序列即可。做出来了别忙着高兴，得意之时就是你灭亡之日，别忘了最后可能还有几个置换需要模拟。经典题目5《算法艺术与信息学竞赛》207页（2.1代数方法和模型，[例题5]细菌，版次不同可能页码有偏差）大家自己去看看吧，书上讲得很详细。解题方法和上一题类似，都是用矩阵来表示操作，然后二分求最终状态。经典题目6给定n和p，求第n个Fibonacci数modp的值，n不超过2^31根据前面的一些思路，现在我们需要构造一个2x2的矩阵，使得它乘以(a,b)得到的结果是(b,a+b)。每多乘一次这个矩阵，这两个数就会多迭代一次。那么，我们把这个2x2的矩阵自乘n次，再乘以(0,1)就可以得到第n个Fibonacci数了。不用多想，这个2x2的矩阵很容易构造出来：经典题目7VOJ1067我们可以用上面的方法二分求出任何一个线性递推式的第n项，其对应矩阵的构造方法为：在右上角的(n-1)*(n-1)的小矩阵中的主对角线上填1，矩阵第n行填对应的系数，其它地方都填0。例如，我们可以用下面的矩阵乘法来二分计算f(n)=4f(n-1)-3f(n-2)+2f(n-4)的第k项：利用矩阵乘法求解线性递推关系的题目我能编出一卡车来。这里给出的例题是系数全为1的情况。经典题目8给定一个有向图，问从A点恰好走k步（允许重复经过边）到达B点的方案数modp的值把给定的图转为邻接矩阵，即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即可。经典题目9用1x2的多米诺骨牌填满MxN的矩形有多少种方案，M<=5，N<2^31，输出答案modp的结果我们以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。后面我写了一份此题的源代码。你可以再次看到位运算的相关应用。经典题目10POJ2778题目大意是，检测所有可能的n位DNA串有多少个DNA串中不含有指定的病毒片段。合法的DNA只能由ACTG四个字符构成。题目将给出10个以内的病毒片段，每个片段长度不超过10。数据规模n<=2000000000。下面的讲解中我们以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原创，转贴请注明出处。#include<cstdio>#defineSIZE(1<<m)#defineMAX_SIZE32usingnamespacestd;classCMatrix{public:longelement[MAX_SIZE][MAX_SIZE];voidsetSize(int);voidsetModulo(int);CMatrixoperator*(CMatrix);CMatrixpower(int);private:intsize;longmodulo;};voidCMatrix::setSize(inta){for(inti=0;i<a;i++)for(intj=0;j<a;j++)element[i][j]=0;size=a;}voidCMatrix::setModulo(inta){modulo=a;}CMatrixCMatrix::operator*(CMatrixparam){CMatrixproduct;product.setSize(size);product.setModulo(modulo);for(inti=0;i<size;i++)for(intj=0;j<size;j++)for(intk=0;k<size;k++){product.element[i][j]+=element[i][k]*param.element[k][j];product.element[i][j]%=modulo;}returnproduct;}CMatrixCMatrix::power(intexp){CMatrixtmp=(*this)*(*this);if(exp==1)return*this;elseif(exp&1)returntmp.power(exp/2)*(*this);elsereturntmp.power(exp/2);}intmain(){constintvalidSet[]={0,3,6,12,15,24,27,30};longn,m,p;CMatrixunit;scanf("%d%d%d",&n,&m,&p);unit.setSize(SIZE);for(inti=0;i<SIZE;i++)for(intj=0;j<SIZE;j++)if(((~i)&j)==((~i)&(SIZE-1))){boolisValid=false;for(intk=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]);return0;}理解矩阵乘法大多数人在高中，或者大学低年级，都上过一门课《线性代数》。这门课其实是教矩阵。刚学的时候，还蛮简单的，矩阵加法就是相同位置的数字加一下。矩阵减法也类似。矩阵乘以一个常数，就是所有位置都乘以这个数。但是，等到矩阵乘以矩阵的时候，一切就不一样了。这个结果是怎么算出来的？教科书告诉你，计算规则是，第一个矩阵第一行的每个数字（2和1），各自乘以第二个矩阵第一列对应位置的数字（1和1），然后将乘积相加（2x1+1x1），得到结果矩阵左上角的那个值3。也就是说，结果矩阵第m行与第n列交叉位置的那个值，等于第一个矩阵第m行与第二个矩阵第n列，对应位置的每个值的乘积之和。怎么会有这么奇怪的规则？我一直没理解这个规则的含义，导致《线性代数》这门课就没学懂。研究生时发现，线性代数是向量计算的基础，很多重要的数学模型都要用到向量计算，所以我做不了复杂模型。这一直让我有点伤心。前些日子，受到一篇文章的启发，我终于想通了，矩阵乘法到底是什么东西。关键就是一句话，矩阵的本质就是线性方程式，两者是一一对应关系。如果从线性方程式的角度，理解矩阵乘法就毫无难度。下面是一组线性方程式。矩阵的最初目的，只是为线性方程组提供一个简写形式。老实说，从上面这种写法，已经能看出矩阵乘法的规则了：系数矩阵第一行的2和1，各自与x和y的乘积之和，等于3。不过，这不算严格的证明，只是线性方程式转为矩阵的书写规则。下面才是严格的证明。有三组未知数x、y和t，其中x和y的关系如下。x和t的关系如下。有了这两组方程式，就可以求y和t的关系。从矩阵来看，很显然，只要把第二个矩阵代入第一个矩阵即可。从方程式来看，也可以把第二个方程组代入第一个方程组。上面的方程组可以整理成下面的形式。最后那个矩阵等式，与前面的矩阵等式一对照，就会得到下面的关系。矩阵乘法的计算规则，从而得到证明。阅读材料：两个具体矩阵合成的乘法公式的验证通过前面的学习，我们知道两个变换A和B的合成是一个新的变换，比如：两个旋转变换矩阵A=⎛0-1⎫⎛0-1⎫⎛1⎫⎪⎪B=α=和都将向量⎪10⎪2⎪⎪逆时针旋转90度，那么很自然10⎝⎭⎝⎭⎝⎭的想法就是它们的合成变换AB将向量α逆时针旋转180度，而将向量α逆时针旋转180度的变换我们可以用矩阵⎛-10⎫⎪来表示（为什么？）。于是，我们思考对于一般的两个⎪⎝0-1⎭矩阵变换A和B，它们的合成变换AB是否能用一个新的矩阵表示出来，接下来我们就用具体的几个例子来验证：两个一般矩阵变换合成的乘法公式。即对于任意两个矩阵变换⎛ab⎫⎛uv⎫⎪⎪A=和B=⎪，则它们的合成是一个新变换C，记C=AB，且满足⎪stcd⎝⎭⎝⎭⎛ab⎫⎛uv⎫⎛au+bsav+bt⎫C=cd⎪⎪st⎪⎪=cu+dscv+dt⎪⎪——————（*）⎝⎭⎝⎭⎝⎭范例1：两次旋转⎛cosθ我们取矩阵A=sinθ⎝-sinθ⎫⎛cosϕ⎪，B=⎪sinϕcosθ⎭⎝-sinϕ⎫⎪，它们分别表示把平面上cosϕ⎪⎭的任意向量α逆时针旋转θ度和逆时针旋转ϕ度的变换。利用A和B，对平面上的任意向量连续做两次变换，即先逆时针旋转θ度，再逆时针旋转ϕ度。实际上连续完成这两个变换，变换的效果可以用一个变换来表示，即使向量α逆时针旋转θ+ϕ度的变换，矩阵θ+ϕ)-sin(θ+ϕ)⎫⎛cos(⎪C=恰好表示了这一变换。让我们来看下面的几何解释：（其sin(θ+ϕ)cos(⎪θ+ϕ)⎭⎝中向量α就是向量OD）旋转θ+ϕ度到D'点。它们的合成是逆时针旋转θ+ϕ度的旋转变换下面我们再从代数上验证，通过前面的学习我们知道所谓两个矩阵变换相等就是指它们作用在任意平面向量上的结果一样，并且为了验证的简便，根据矩阵的线性性质，我们只需要验证它们在单位向量再用等式y⎪⎪=x0⎪⎪+y1⎪⎪即可说明任意向0⎪⎪和1⎪⎪上的作用，⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭量的正确性。⎛1⎫⎛0⎫⎛x⎫⎛1⎫⎛0⎫⎛⎛cosθ⎛1⎫⎛1⎫⎪⎪BA=B(A)=B0⎪0⎪sinθ⎝⎭⎝⎭⎝⎝-sinθ⎫⎛1⎫⎫⎛cosθ⎫⎪⎪⎪⎪=B⎪⎪⎪⎪cosθ⎭⎝0⎭⎭⎝sinθ⎭⎛cosϕ=sinϕ⎝-sinϕ⎫⎛cosθ⎫⎛cos(θ+ϕ)⎫⎪⎪=sinθ⎪sin(θ+ϕ)⎪⎪cosϕ⎪⎭⎝⎭⎝⎭θ+ϕ)-sin(θ+ϕ)⎫⎛1⎫⎛cos(θ+ϕ)⎫⎛1⎫⎛cos(⎛1⎫⎛1⎫⎪⎪⎪⎪⎪⎪C==BA=C，于是0⎪sin(θ+ϕ)cos(⎪⎪⎪⎪⎪，θ+ϕ)0sin(θ+ϕ)00⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭⎛0⎫⎛0⎫⎪BA=B(A1⎪1⎪⎪)=⎝⎭⎝⎭⎛⎛cosθB⎝⎝sinθ-sinθ⎫⎛0⎫⎫⎪1⎪⎪⎪⎪=cosθ⎪⎭⎝⎭⎭⎛-sinθ⎫Bcosθ⎪⎪⎝⎭⎛cosϕ=sinϕ⎝⎛0⎫⎝1⎭-sinϕ⎫⎛-sinθ⎫⎛-sin(θ+ϕ)⎫⎪⎪⎪=⎪⎪cosϕ⎭⎝cosθ⎭⎝cos(θ+ϕ)⎪⎭C⎪⎪=θ+ϕ)-sin(θ+ϕ)⎫⎛0⎫⎛-sin(θ+ϕ)⎫⎛cos(⎛0⎫⎛0⎫⎪⎪⎪⎪⎪=BA=C，于是⎪⎪⎪⎪⎪θ+ϕ)⎭⎝1⎭⎝cos(θ+ϕ)⎭⎝sin(θ+ϕ)cos(⎝1⎭⎝1⎭于是BA=C，得证。范例2：两次压缩0⎫⎛10⎫⎪1，B=1⎪（m,n≠1），它们分别表示把平面上的任意⎪0⎪m⎭n⎭⎝11向量α向x方向垂直压缩为和的变换。利用A和B，对平面上的任意向量连续做两次mn⎛1我们取矩阵A=0⎝11，再向x方向垂直压缩为。实际上连续完成这两个变mn1换，变换的效果可以用一个变换来表示，即使向量α向x方向垂直压缩为的变换，矩mn⎛10⎫1⎪恰好表示了这一变换。让我们来看下面的几何解释：阵C=（其中向量α就是向0⎪mn⎭⎝量OD）变换，即先向x方向垂直压缩为y矩阵变换A下面我们再从代数上验证，通过前面的学习我们知道所谓两个矩阵变换相等就是指它们作用在任意平面向量上的结果一样，并且为了验证的简便，根据矩阵的线性性质，我们只需要验证它们在单位向量再用等式y⎪⎪=x0⎪⎪+y1⎪⎪即可说明任意向0⎪⎪和1⎪⎪上的作用，⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭量的正确性。⎛1⎫⎛0⎫⎛x⎫⎛1⎫⎛0⎫⎛⎛1⎛1⎫⎛1⎫BA0⎪⎪=B(A0⎪⎪)=B0⎝⎭⎝⎭⎝⎝0⎫⎛1⎫⎫1⎛1⎫⎛⎪⎪1⎪⎪⎪=B0⎪⎪=0⎪0⎝⎭⎝m⎭⎝⎭⎭0⎫⎛1⎫⎛1⎫1⎪⎪⎪=0⎪⎪⎪0n⎭⎝⎭⎝⎭1⎛1⎫⎛C0⎪⎪=0⎝⎭⎝0⎫⎛1⎫⎛1⎫⎛1⎫⎛1⎫1⎪⎪⎪⎪⎪BA=C，于是，=⎪⎪⎪⎪⎪00⎝0⎭⎝0⎭mn⎭⎝⎭⎝⎭0⎫⎛0⎫⎫⎛0⎫⎛10⎫⎛0⎫⎛0⎫⎪⎪1⎪=B1⎪=1⎪1⎪=1⎪⎪⎪⎝1⎭⎪⎪1⎪⎪⎪m⎭⎝m⎭⎝n⎭⎝m⎭⎝mn⎭⎭⎛⎛1⎛0⎫⎛0⎫BA1⎪⎪=B(A1⎪⎪)=B0⎝⎭⎝⎭⎝⎝10⎫⎛0⎫⎛0⎫0⎫⎛0⎫⎛0⎫⎛⎪⎪=1⎪，于是BA⎛1⎪⎪=CC=1⎪1⎪⎪1⎪01⎪⎪⎪⎝⎭⎝⎭⎝⎭⎝mn⎭⎝⎭⎝m⎭于是BA=C，得证。范例3：先切变，再反射⎛1a⎫⎛10⎫⎪⎪我们取矩阵A=，B=⎪⎪，它们分别表示把平面上的任意向量α作切010-1⎝⎭⎝⎭变变换和关于x轴对称的变换。利用A和B，对平面上的任意向量连续做两次变换，即先作切变变换，再关于x轴对称。实际上连续完成这两个变换，变换的效果可以用一个变换来表示，矩阵C=0⎝⎛1a⎫⎪恰好表示了这一变换。让我们来看下面的几何解释：（其中向量α-1⎪⎭就是向量OD）下面我们再从代数上验证，通过前面的学习我们知道所谓两个矩阵变换相等就是指它们作用在任意平面向量上的结果一样，并且为了验证的简便，根据矩阵的线性性质，我们只需⎛1⎫⎛0⎫⎛x⎫要验证它们在单位向量再用等式y⎪⎪=0⎪⎪和1⎪⎪上的作用，⎝⎭⎝⎭⎝⎭量的正确性。⎛1⎫x0⎪⎪+⎝⎭⎛0⎫y1⎪⎪即可说明任意向⎝⎭⎛⎛1a⎫⎛1⎫⎫⎛1⎫⎛1⎫⎛1⎫⎛10⎫⎛1⎫⎛1⎫⎪⎪⎪⎪⎪BA=B(A)=B=B0⎪0⎪0⎪⎪=0-1⎪⎪0⎪⎪=0⎪⎪01⎪0⎪⎪⎝⎭⎝⎭⎭⎝⎭⎭⎝⎭⎝⎭⎝⎭⎝⎭⎝⎝⎛1⎫⎛1a⎫⎛1⎫⎛1⎫⎛1⎫⎛1⎫⎪⎪⎪⎪⎪⎪C⎪==⎪，于是BA⎪=C⎪⎪⎪，00-10000⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭⎛⎛1a⎫⎛0⎫⎫⎛0⎫⎛0⎫⎛a⎫⎛10⎫⎛a⎫⎛a⎫⎪⎪⎪⎪⎪BA=B(A)=B=B1⎪1⎪1⎪⎪=0-1⎪⎪1⎪⎪=-1⎪⎪01⎪1⎪⎪⎝⎭⎝⎭⎭⎝⎭⎭⎝⎭⎝⎭⎝⎭⎝⎭⎝⎝⎛0⎫⎛1a⎫⎛0⎫⎛a⎫⎛0⎫⎛0⎫C1⎪⎪=0-1⎪⎪1⎪⎪=-1⎪⎪，于是BA1⎪⎪=C1⎪⎪⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭⎝⎭于是BA=C，得证。请同学们按照上述证明的办法求下面矩阵变换的合成变换矩阵，并用矩阵的乘法公式验证：⎛0-1⎫⎛-10⎫⎪⎪1、A=，B=⎪⎪100-1⎝⎭⎝⎭⎛1⎫⎛10⎫0⎪1⎪2、A=2，B=⎪0⎪013⎝⎭⎝⎭3、A=⎛1a⎫⎛1b⎫⎪⎪B=，⎪⎪实验一测试矩阵乘法的性能（计算机性能评测）实验目的计算机系统的性能指标体现于时间和空间两个方面，在系统上程序实际运行的时间应该是衡量机器时间(速度)性能是最可靠的标准。通过本实验可以帮助学生加深对计算机性能评测问题的理解，掌握利用时间衡量某种系统结构某种算法的方法。实验环境开发工具使用windows平台下的vc++6.0。实验内容利用随机函数产生出两个随机矩阵，计算两个矩阵相乘。在此过程中，记录下两个矩阵相乘之前之后的两个系统时间，求出时间差。绘制出矩阵维数和时间的坐标图。访问网站和,了解专业的测试方式和测试报告。注意此网站提供的测试程序和数据均是针对服务器和机架式计算机的。运行PC机测试软件测试实验室计算机的硬件参数和性能(侧重中央处理器)。实验结果出现的问题及解决方案问题：计算机性能的主要评价指标都有哪些？解决方案：运算速度；字长；内存储器的容量；外存储器的容量。问题：计算机性能的基准测试程序都有哪些？解决方案：1.综合型（如Dhrystone,Whetstone等）；2.核心型（如LivemoreFortranKernals,NASA之NAS等）；3.数学库（如Linpack,FFT等）；4.应用型（如SPEC，Perfect,Splash等）；5.并行型（如NAS之NPB,PARKBENCK等）。思考题除了编写算法程序测量计算机性能外，还有什么其它的方法来评测计算机的性能？游戏性能测试买电脑的朋友很少有不玩游戏的，而且游戏可以说是对电脑性能的综合测试，包含了对CPU、内存、显卡、主板、显示器、光驱、键盘鼠标、声卡、音箱等的测试。所以，电脑首先应该进行的就是游戏测试。我们可以选择几款常见的游戏来测试爱机。例如：极品飞车、古墓丽影、QUAKE、CS、虚幻竞技场、魔兽争霸。不一定要把这些游戏都试用一下，可以选择其中的几款来测试电脑性能。电脑配置高一些的朋友可以选择高一些的游戏版本来测试，配置低一些的朋友可以选择版本低一些的游戏来测试。测试主要应该注意游戏安装速度、游戏运行速度、游戏画质、游戏流畅程度、游戏音质等几方面。可以更改显示器设置、显卡设置、BIOS设置、系统设置、游戏设置来感受不同设置下电脑的不同表现。例如改变显示器的亮度、对比度，改变游戏的分辨率，改变显卡的频率，改变内存的延时，改变CPU频率，改变系统硬件加速比例，改变系统缓存设置等等。大家要注意的是在测试以前最好把所有的补丁程序安装齐全，改变设置测试完成以后要把设置改回来(或者改到最佳状态)。有条件的朋友可以和配置相近的电脑对比一下，相信能感受出自己爱机的性能。游戏性能方面的测试现在主要以Fraps为主，这个软件主要用于游戏运行过程中的实时帧速测试，并可以记录测试过程中的平均、最高以及最低帧速，帮助用户考量本身配置的性能。另外，也有一些在线的测试系统（比如INMARK），可以帮助用户查看配置在某些游戏中的运行流畅度。播放电影测试接下来可以考虑播放一段电影来测试自己的电脑。建议选择常用的播放器和比较熟悉的电影这样可能不用和其他电脑对比就能看出自己爱机的“优势”。这时候应该注意的是播放有没有异常、画面的鲜艳程度、调整显示器亮度后的画面变化情况、电影画面的清晰程度等等。对于软解性能，可以参考H264Encoder测试软件。图片处理测试再下来可以考虑测试一下电脑的图片处理能力。笔者推荐用常用的图形处理软件来测试，例如PHOTOSHOP、FIREWORKS、AUTOCAD、3DMAX等等。可以试着打开多个图片文件、更改图片或者编辑图片来测试电脑图片处理速度、观察画质。或者通过渲染性能测试软件测试，比如POV-Ray（根据管线追踪绘制3D图像）、CineBench（重负载渲染性能测试软件），这些软件均可测试机器本身在3D/2D工具软件中的性能表现。拷贝文件测试此操作比较简单，尽量选择大一些的文件拷贝，大家可以选择拷贝VCD，DVD或者包含大量文件的文件夹。压缩测试可以选择我们常用的WINZIP或者WINRAR来压缩大一些的文件。也可以通过压缩CD、VCD来测试电脑，选择我们常用的超级解霸软件来测试。以上测试重点查看速度。网络性能测试网络测试主要检查网络是否能正常连接、连接速度是否正常。多线程性能测试多线程性能，主要指CPU核心的处理性能，通过Wprime2.0（主要针对多线程处理性能的测试）或者FritzChess（IBM“深蓝”进行过的测试项目，国际象棋预测计算）两款软件，可以测试用户现有配置中CPU多线程以及基础运算的性能。其他测试软件工具3DMARK系列软件：主要针对显卡，以及部分CPU性能的专业测试软件；PCMARK系列软件：主要针对整体平台运算性能的专业测试软件；MemTest软件：对内存进行监测以及性能测试的辅助软件；Furmark软件：针对显卡图形能力的重负载测试软件，同时可监测重负载时的显卡温度；HDTuch软件：针对硬盘的磁盘性能测试软件。实验总结通过这次的实验，一些系统设计上的缺陷。如在测试中经常会遇到系统I/O的问题,必须在程序编译时选取合适的"开关",程序才能计算出正确的结果。再如高级语言编译软件设计标准问题,有些系统对标准的Fortran77程序都不能顺利通过。又如有的系统上常需要降低编译优化的级别才能正常计算。在对一些系统进行测试的过程中,还由于显示出的计算结果不正确或计算结果的精度达不到要求,从而帮助厂家查出了计算机系统硬件或软件设计上的一些问题。好像目前还没有这方面题目的总结。这几天连续看到四个问这类题目的人，今天在这里简单写一下。这里我们不介绍其它有关矩阵的知识，只介绍矩阵乘法和相关性质。不要以为数学中的矩阵也是黑色屏幕上不断变化的绿色字符。在数学中，一个矩阵说穿了就是一个二维数组。一个n行m列的矩阵可以乘以一个m行p列的矩阵，得到的结果是一个n行p列的矩阵，其中的第i行第j列位置上的数等于前一个矩阵第i行上的m个数与后一个矩阵第j列上的m个数对应相乘后所有m个乘积的和。比如，下面的算式表示一个2行2列的矩阵乘以2行3列的矩阵，其结果是一个2行3列的矩阵。其中，结果的那个4等于2*2+0*1：下面的算式则是一个1x3的矩阵乘以3x2的矩阵，得到一个1x2的矩阵：矩阵乘法的两个重要性质：一，矩阵乘法不满足交换律；二，矩阵乘法满足结合律。为什么矩阵乘法不满足交换律呢？废话，交换过来后两个矩阵有可能根本不能相乘。为什么它又满足结合律呢？仔细想想你会发现这也是废话。假设你有三个矩阵A、B、C，那么(AB)C和A(BC)的结果的第i行第j列上的数都等于所有A(ik)*B(kl)*C(lj)的和（枚举所有的k和l）。经典题目1给定n个点，m个操作，构造O(m+n)的算法输出m个操作后各点的位置。操作有平移、缩放、翻转和旋转这里的操作是对所有点同时进行的。其中翻转是以坐标轴为对称轴进行翻转（两种情况），旋转则以原点为中心。如果对每个点分别进行模拟，那么m个操作总共耗时O(mn)。利用矩阵乘法可以在O(m)的时间里把所有操作合并为一个矩阵，然后每个点与该矩阵相乘即可直接得出最终该点的位置，总共耗时O(m+n)。假设初始时某个点的坐标为x和y，下面5个矩阵可以分别对其进行平移、旋转、翻转和旋转操作。预先把所有m个操作所对应的矩阵全部乘起来，再乘以(x,y,1)，即可一步得出最终点的位置。经典题目2给定矩阵A，请快速计算出A^n（n个A相乘）的结果，输出的每个数都modp。由于矩阵乘法具有结合律，因此A^4=A*A*A*A=(A*A)*(A*A)=A^2*A^2。我们可以得到这样的结论：当n为偶数时，A^n=A^(n/2)*A^(n/2)；当n为奇数时，A^n=A^(n/2)*A^(n/2)*A（其中n/2取整）。这就告诉我们，计算A^n也可以使用二分快速求幂的方法。例如，为了算出A^25的值，我们只需要递归地计算出A^12、A^6、A^3的值即可。根据这里的一些结果，我们可以在计算过程中不断取模，避免高精度运算。经典题目3POJ3233(感谢rmq)题目大意：给定矩阵A，求A+A^2+A^3+...+A^k的结果（两个矩阵相加就是对应位置分别相加）。输出的数据modm。k<=10^9。这道题两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：A+A^2+A^3+A^4+A^5+A^6=(A+A^2+A^3)+A^3*(A+A^2+A^3)应用这个式子后，规模k减小了一半。我们二分求出A^3后再递归地计算A+A^2+A^3，即可得到原问题的答案。经典题目4VOJ1049题目大意：顺次给出m个置换，反复使用这m个置换对初始序列进行操作，问k次置换后的序列。m<=10,k<2^31。首先将这m个置换“合并”起来（算出这m个置换的乘积），然后接下来我们需要执行这个置换k/m次（取整，若有余数则剩下几步模拟即可）。注意任意一个置换都可以表示成矩阵的形式。例如，将1234置换为3124，相当于下面的矩阵乘法：置换k/m次就相当于在前面乘以k/m个这样的矩阵。我们可以二分计算出该矩阵的k/m次方，再乘以初始序列即可。做出来了别忙着高兴，得意之时就是你灭亡之日，别忘了最后可能还有几个置换需要模拟。经典题目5《算法艺术与信息学竞赛》207页（2.1代数方法和模型，[例题5]细菌，版次不同可能页码有偏差）大家自己去看看吧，书上讲得很详细。解题方法和上一题类似，都是用矩阵来表示操作，然后二分求最终状态。经典题目6给定n和p，求第n个Fibonacci数modp的值，n不超过2^31根据前面的一些思路，现在我们需要构造一个2x2的矩阵，使得它乘以(a,b)得到的结果是(b,a+b)。每多乘一次这个矩阵，这两个数就会多迭代一次。那么，我们把这个2x2的矩阵自乘n次，再乘以(0,1)就可以得到第n个Fibonacci数了。不用多想，这个2x2的矩阵很容易构造出来：经典题目7VOJ1067我们可以用上面的方法二分求出任何一个线性递推式的第n项，其对应矩阵的构造方法为：在右上角的(n-1)*(n-1)的小矩阵中的主对角线上填1，矩阵第n行填对应的系数，其它地方都填0。例如，我们可以用下面的矩阵乘法来二分计算f(n)=4f(n-1)-3f(n-2)+2f(n-4)的第k项：利用矩阵乘法求解线性递推关系的题目我能编出一卡车来。这里给出的例题是系数全为1的情况。经典题目8给定一个有向图，问从A点恰好走k步（允许重复经过边）到达B点的方案数modp的值把给定的图转为邻接矩阵，即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即可。经典题目9用1x2的多米诺骨牌填满MxN的矩形有多少种方案，M<=5，N<2^31，输出答案modp的结果我们以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。后面我写了一份此题的源代码。你可以再次看到位运算的相关应用。经典题目10POJ2778题目大意是，检测所有可能的n位DNA串有多少个DNA串中不含有指定的病毒片段。合法的DNA只能由ACTG四个字符构成。题目将给出10个以内的病毒片段，每个片段长度不超过10。数据规模n<=2000000000。下面的讲解中我们以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原创，转贴请注明出处。#include<cstdio>#defineSIZE(1<<m)#defineMAX_SIZE32usingnamespacestd;classCMatrix{public:longelement[MAX_SIZE][MAX_SIZE];voidsetSize(int);voidsetModulo(int);CMatrixoperator*(CMatrix);CMatrixpower(int);private:intsize;longmodulo;};voidCMatrix::setSize(inta){for(inti=0;i<a;i++)for(intj=0;j<a;j++)element[i][j]=0;size=a;}voidCMatrix::setModulo(inta){modulo=a;}CMatrixCMatrix::operator*(CMatrixparam){CMatrixproduct;product.setSize(size);product.setModulo(modulo);for(inti=0;i<size;i++)for(intj=0;j<size;j++)for(intk=0;k<size;k++){product.element[i][j]+=element[i][k]*param.element[k][j];product.element[i][j]%=modulo;}returnproduct;}CMatrixCMatrix::power(intexp){CMatrixtmp=(*this)*(*this);if(exp==1)return*this;elseif(exp&1)returntmp.power(exp/2)*(*this);elsereturntmp.power(exp/2);}intmain(){constintvalidSet[]={0,3,6,12,15,24,27,30};longn,m,p;CMatrixunit;scanf("%d%d%d",&n,&m,&p);unit.setSize(SIZE);for(inti=0;i<SIZE;i++)for(intj=0;j<SIZE;j++)if(((~i)&j)==((~i)&(SIZE-1))){boolisValid=false;for(intk=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]);return0;}#defineR1#defineG2#defineVS3#defineCS4#defineVCCS5#defineVCVS6#defineCCCS7#defineCCVS8#defineOPAMP9typedefstructtagBranch{ intnumber;//支路号 chartype;//类型号 intnfrom;//始节点 intnto;//终节点 floatvalue;//参数值 intncfrom;//控制始节点/控制终节点；G/R支路的“设置为电压定义支路”标记 intncto;//控制终节点；G/R支路的“控制用电压定义支路”标记}BRANCH;BRANCHBranch0;*BRANCH;Branch=newBRANCH[BranchNumber+1];intBranchNumber;intNodeNumber;intvBranchNumber;intNN;doubleval[NN];introw[NN];intcol[NN];intup[NN];intdown[NN];intleft[NN];intright[NN];intrp[NN];intcp[NN];intdel[NN];intnza;intn;typedefstructtagLINKED_LIST{ intnza;//非零元素个数 intn;//矩阵阶数 int*row;//元素行号 int*col;//元素列号 int*up;//元素的上邻元素号 int*down;//元素的下邻元素号 int*left;//元素的左邻元素号 int*right;//元素的右邻元素号 int*rp;//行链指针 int*cp;//列链指针 char*del;//元素删除标志；1=删除}LINKED_LIST;typedefstructtagTRIANGULAR_TABLE{ intnza;//非零元素个数 intn;//矩阵阶数 double*val;//矩阵及右端向量/中间解向量元素值 int*roco;//对角元素行列号，U元素列号，L元素行号 int*urp;//U阵行指针 int*lcp;//L阵列指针 intlun;//由符号LU分解确定的需修正的元素个数 intfen;//由符号前消确定的需修正的向量元素个数 int*lup;//由符号LU分解确定的需修正的元素号 int*fep;//前消需修正的元素号}TRIANGULAR_TABLE;typedefstructtagTRIANGULAR_TABLE_C{ intnza; intn; intlun; intfen; double*vre; double*vim; int*roco; int*urp; int*lcp; int*lup; int*fep;}tagTRIANGULAR_TABLE_C;#include"XISHU.h"//创建双重链接表voidnew_list(intrank_max,intnonzero,LINKED_LIST*list){ row=newint[nonzero]; col=newint[nonzero]; up=newint[nonzero]; down=newint[nonzero]; left=newint[nonzero]; right=newint[nonzero]; rp=newint[rank_max]; cp=newint[rank_max]; del=newint[nonzero]; nza=0;//创建时尚无元素 inti; for(i=0;i<rank_max;i++){//初始化 rp[i]=0; cp[i]=0; } for(i=0;i<nonzero+1;i++){ row[i]=0; col[i]=0; up[i]=0; down[i]=0; left[i]=0; right[i]=0; del[i]=0; } }//-----------------------------------------------------------------charinsert_row(intnza,intx,inty,LINKED_LIST*list){intk,kl,z;k=(*list).rp[x];if(k==0){//第x行还未有非零元素，a[x][y]为第一个元素 (*list)->rp[x]=nza; (*list)->left[nza]=0; (*list)->right[nza]=0; return(0); }//成功插入for(;;){//搜索第x行非零元素 z=(*list).col[k];//a[x][y]为非零元素 if(y<z){//nza在k左 kl=(*list).left[k]; if(kl==0){//k为原行首，nza在其左，为新行首 (*list)->left[nza]=0; (*list)->right[nza]=k; (*list)->left[k]=nza; (*list)->rp[x]=nza; } else{//nza在kl和k之间 (*list)->left[nza]=kl; (*list)->right[nza]=k; (*list)->left[k]=nza; (*list)->right[k]=nza; } return(0); } elseif(y>z){//nza在k右 kl=(*list).right[k]; if(kl==0){//k为原行尾，nza在其中，为新行尾 (*list)->right[nza]=0; (*list)->right[k]=nza; (*list)->left[nza]=k; return(0); } k=kl;//右边还有非零元素，继续向右搜索 continue; } else{//y=z,a[x][y]已存在 return(1);//已有，未插入 }}}//----------------------------------------------------------------charinsert_col(intnza,intx,inty,LINKED_LIST*list){ intk,kl,z; k=(*list).cp[y]; if(k==0){//第y列还未有非零元素，a[x][y]为第一个元素 (*list)->cp[y]=nza; (*list)->up[nza]=0; (*list)->down[nza]=0; return(0); }//成功插入 for(;;){//搜索第y列非零元素 z=(*list).row[k];//a[x][y]为非零元素 if(x<z){//nza在k左 kl=(*list).up[k]; if(kl==0){//k为原列首，nza在其上，为新列首 (*list)->up[nza]=0; (*list)->down[nza]=k; (*list)->up[k]=nza; (*list)->cp[x]=nza; } else{//nza在kl和k之间 (*list)->up[nza]=kl; (*list)->down[nza]=k; (*list)->up[k]=nza; (*list)->down[k]=nza; } return(0); } elseif(x>z){//nza在k下 kl