1.矩陣是個啥
一個m×n的矩陣就是m×n個數(shù)排成m行n列的一個數(shù)陣。矩陣可以用于批量解決一些線性問題,例如遞推方程。
另,1*m 的矩陣又叫做向量(沒錯就是vector)
2.啥是矩陣乘法
蒟陣乘法,顧名思義,就是兩個矩陣(蒟蒻)乘起來,不過和普通的乘法不一樣。
A.兩個矩陣能相乘的條件:
當一個矩陣的行數(shù)等于另一個矩陣的列數(shù)時,兩個相乘才有意義。
B.矩陣咋相乘。
三個矩陣a,b,c
如果a等于bc的乘積,那么寫作a=bc
假設b,c都是n*n的矩陣那么有

寫成矩陣的形式就是

這也就是為什么必須要求一個矩陣的行數(shù)等于另一個矩陣的列數(shù),否則元素數(shù)不相等
C.注意事項
如果兩個大小為a*n,n*m
的矩陣相乘,最后乘出來的矩陣就是是a*m
的
而矩陣乘法的次數(shù)就是a*n*m
例如2*3
和3*2
的矩陣相乘

如果問原因的話,我們可以發(fā)現(xiàn),新矩陣的行數(shù)和左邊矩陣的行數(shù)相等,而新矩陣的列數(shù)和后邊矩陣的列數(shù)相等。
D.這東西有啥性質啊
a.矩陣乘法滿足結合律,即a(bc)=(ab)c,可以通過這個性質減少乘法次數(shù)
b.而且,矩陣乘法沒有交換律
因為普通的加法和乘法滿足交換律的原因是因為加法和乘法的結果都只是一個數(shù),而矩陣的本質是數(shù)的排列,所以運算結果還是矩陣,并不是一個數(shù),自然和數(shù)有區(qū)別
c.矩陣的單位元..…這東西也不知道放在哪,就寫這吧。

就是對角線上全為1
3.這玩意有啥用
md我咋知道啊
在某些題中可以加速遞推公式,可以把一個O(n)的式子優(yōu)化到O(logn)
A.應用前提:矩陣快速冪(md這又是啥啊)
因為矩陣的乘法也有結合律,所以也可以用快速冪加速,也是矩陣優(yōu)化遞推的原理,即把一個不能用快速冪優(yōu)化的式子轉換成快速冪的形式
4.來個栗子啊
出門左轉炒鍋
f(n) = f(n-1) + f(n-2) (n ≥ 2 且n為整數(shù))
請你求出$f(n) \mod 10^9+7$
的值(n在long long(INT64)范圍內)。
顯然線性爆炸,所以考慮加速。那么很容易想到矩陣乘法(我也不會別的)
那么考慮把這個遞推式子變成可以矩陣的形式,那么就考慮把答案弄進一個矩陣里,(那么就設一個1*1
的吧!(啪))顯然,設一個1*1
的矩陣是一個很不明智的行為,因為沒法轉移??紤]到當前的狀態(tài)與前兩項有關。
所以構造一個1*2
向量表示[F_n&F_{n-1}]
那么轉移也就很顯然了
我們要達到這樣一個目的:

變個形

那么就可以構造中間的矩陣了
顯然就是

帶回去驗證一下qwq

發(fā)現(xiàn)成立了誒。
5.更大的栗子
出門右轉大炒鍋
廣義的斐波那契數(shù)列是指形如$a_n=p\times a_{n-1}+q\times a_{n-2}$
的式子,給出$p,q,a_1,a_2,n,m$
求$a_n\mod m$
。
除了和上一題的系數(shù)不變沒區(qū)別啊喂~~
考慮到只有轉移時的系數(shù)變了,那么只要改變你構造的矩陣就可以了,那么我們再考慮轉移是如何進行的

那么注意到改變的實際上只有這一項$F_{n-1}*1+F_n*1$
,考慮這一項是怎么來的,根據(jù)矩陣乘法的定義我們知道這一項實際上是
向量的第一行乘矩陣的第一列,也就是

那么我們只需要改動這一列即可,即把這一列改為

那么最后的轉移矩陣就是

最后放上這個題的核心代碼
ll mod;
struct mar
{
ll a[4][4];
}ans,base;
mar mul(mar x,mar y)
{
mar t;
memset(t.a,0,sizeof(t.a));
for(int i=1;i<>2;++i)
for(int j=1;j<>2;++j)
for(int k=1;k<>2;++k)
t.a[i][j]+=(x.a[i][k])%mod*(y.a[k][j])%mod,t.a[i][j]%=mod;
return t;
}
ans.a[1][1]=a2,ans.a[1][2]=a1;
base.a[1][2]=1;
base.a[2][1]=q;
base.a[1][1]=p;
ksm(n-2);
6.關于矩陣乘法的優(yōu)化
A.緩存優(yōu)化
先普及一點知識: CPU cache
通俗點說就是平時說的緩存,是位于CPU 和內存之間的臨時存儲器,具有體積小,交換速率快的特點,存儲的是 CPU 即將訪問的內存。
所以緩存存在的意義就是為了加快內存訪問的速度。
那么對于OIer來說,這東西唯一的用處就是卡常,可能會有人聽過“數(shù)組不要開偶偶數(shù)大小 之類的玄學”。但這可能真的是有用的。
程序在一段時間內訪問的數(shù)據(jù)通常具有局部性,比如對一維數(shù)組來說,訪問了地址 x 上的元素,那么以后訪問地址 x+1 、 x+2 上元素的可能性就比較高;這也稱為 cache 的空間局限性。
如果說人話的話。就是數(shù)組的訪問順序是水平連續(xù)的。
讓我們再來一個栗子
1
for(int i=1;i<>
for(int j=1;j<>
for(int k=1;k<>
c[i][j]=a[i][k]*b[k][j];
2
for(int i=1;i<>
for(int k=1;k<>
for(int j=1;j<>
c[i][j]=a[i][k]*b[k][j];
眼尖的同學可能會發(fā)現(xiàn),這兩段程序唯一的不同只有二三層的循環(huán)循順序啊
那么實際上的效率可能會在數(shù)據(jù)很大的時候產生極大的差異。
因為根據(jù)我們上邊所了解到的信息,對于代碼2中的 a[i] 和 b[k] 這兩行的訪問都是連續(xù)的,所以緩存會被充分利用。這樣就達到了卡常的目的。
這種做法的術語是 loop permutation(循環(huán)置換)
測試結果如下
為了驗證相等,輸出了第5行,可以看到,效率還是有很大差異的
另附測試代碼
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#define file(a) freopen(a'.in','r',stdin);freopen(a'.out','w',stdout);
using namespace std;
char IN[117],*SS=IN,*TT=IN;
inline char gc(){return (SS==TT&&(TT=(SS=IN)+fread(IN,1,117,stdin),SS==TT)?EOF:*SS++);}
inline int read()
{
int now=0,f=1;register char c=gc();
for(;!isdigit(c);c=='-'&&(f=-1),c=gc());
for(;isdigit(c);now=now*10+c-'0',c=gc());
return now*f;
}
const int mod=1e9+7;
long long a[1100][1100];
long long temp[1010][1010];
long long b[1010][1010];
long long c[1010][1010];
int main()
{
freopen('a.out','r',stdin);
freopen('x1.out','w',stdout);
int n=read();
for(int i=1;i<>
for(int j=1;j<>
for(int i=1;i<>
for(int j=1;j<>
for(int i=1;i<>
for(int k=1;k<>
for(int j=1;j<>
(a[i][j]+=b[i][k]*c[k][j])%=mod;
// for(int i=1;i<>
// for(int j=1;j<>
// for(int k=1;k<>
// (a[i][j]+=(b[i][k]*c[k][j]))%=mod;
for(int i=1;i<>printf('%lld ',a[5][i]);
printf('\n');
printf('Elapsed time:%.2lf secs.\n',1.0*clock()/1000.0);
return 0;
}
有興趣的同學可以去嘗試做一下HUD4920。
經典的矩陣乘法(du)常(liu)題.
B.矩陣分塊
當矩陣特別大時,可能會發(fā)生緩存放不下的情況,那么之前的ikj
優(yōu)化就會收效甚微,所以我們考慮把一個矩陣分為幾塊子矩陣,不過需要注意的是A的行的劃分與B的列的劃分必須一致。
int size=n/25;
for(int it=1;it<>
for(int jt=1;jt<>
for(int kt=1;kt<>
for(int i=it;i<>
for(int j=jt;j<>
for(int k=kt;k<>
(a[i][j]+=b[i][k]*c[k][j])%=mod;
C.次數(shù)優(yōu)化
大家可能做過最優(yōu)矩陣鏈乘問題問題。
本質就是利用結合律來減少乘法次數(shù)。
例如
A是一個50×10的矩陣,B是10×20的矩陣,C是20×5的矩陣
那么要把ABC乘起來的話,有兩種順序,一種是(A(BC)),另一種是((AB)C)
而這兩種順序的乘法次數(shù)一個是3500次,而另一個卻是15000次,所以當數(shù)據(jù)非常大的時候,我們可以考慮構造一個乘法的順序,使得乘法次數(shù)最小,但只適用于矩陣少而大的情況。
D.取模優(yōu)化
眾所周知,計算機中的取模運算非常之慢,而long long的取模是最慢的,所以很多人程序的常數(shù)大的原因也是這個。所以我們可能會在某些人的程序中看到這樣的函數(shù)
ll mo(ll x,ll y){return x+y>=mod?x+y-mod:x+y;}
這就是利用減法代替取模從而達到卡常的效果。
對于矩陣乘法也是這樣。我們可以設定一個ull上界$Q=16ull*mod$
,當超過這個上界時就可以直接用減法代替取模,因為當mod在10^9
級別時ull最多只能存下16個$mod*mod$
相加。
當然中間過程需要long long
來存儲,而最后的結果最好還是int
類型。
所以取模次數(shù)由O(n^3)
降到O(n^2)
次
但實際效果略遜于ikj
,可能是本機的問題
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#define file(a) freopen(a'.in','r',stdin);freopen(a'.out','w',stdout);
using namespace std;
char IN[117],*SS=IN,*TT=IN;
inline char gc(){return (SS==TT&&(TT=(SS=IN)+fread(IN,1,117,stdin),SS==TT)?EOF:*SS++);}
inline int read()
{
int now=0,f=1;register char c=gc();
for(;!isdigit(c);c=='-'&&(f=-1),c=gc());
for(;isdigit(c);now=now*10+c-'0',c=gc());
return now*f;
}
const int maxn=100005;
typedef unsigned long long ull;
const int mod=1e9+7;
const ull Q=16ull*mod*mod;
int a[1100][1100];
long long temp;
ull b[1010][1010];
ull c[1010][1010];
ull mo(ull x,ull y){return x+y>=Q?x+y-Q:x+y;}
int main()
{
freopen('a.out','r',stdin);
freopen('x2.out','w',stdout);
int n=read();
for(int i=1;i<>
for(int j=1;j<>
for(int i=1;i<>
for(int j=1;j<>
for(int i=1;i<>
for(int j=1;j<>
{
temp=0;
for(int k=1,s=1;k<>
{
if(s==16)temp=mo(temp,b[i][k]*c[k][j]),s=1;
else temp+=b[i][k]*c[k][j];
}
a[i][j]=temp%mod;
}
for(int i=1;i<>printf('%d ',a[5][i]);
printf('\n');
printf('Elapsed time:%.2lf secs.\n',1.0*clock()/1000.0);
return 0;
}
和ikj
優(yōu)化一起用的效果

E.特判優(yōu)化
特判矩陣的元素是否為0,如果是0
可以直接continue,在稀疏矩陣里有奇效,由于數(shù)據(jù)比較難生成,所以不實際運行.還請見諒
感謝審核的@ComeIntoPower提出的寶貴意見。
如果還有不足請指出,感謝。
另:測試數(shù)據(jù)均為以下代碼生成
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#define file(a) freopen(a'.in','r',stdin);freopen(a'.out','w',stdout);
using namespace std;
char IN[117],*SS=IN,*TT=IN;
inline char gc(){return (SS==TT&&(TT=(SS=IN)+fread(IN,1,117,stdin),SS==TT)?EOF:*SS++);}
inline int read()
{
int now=0,f=1;register char c=gc();
for(;!isdigit(c);c=='-'&&(f=-1),c=gc());
for(;isdigit(c);now=now*10+c-'0',c=gc());
return now*f;
}
const int maxn=100005;
int main()
{
srand(time(NULL)*20020307%19260817);
freopen('a.out','w',stdout);
int n=1000;
printf('1000\n');
for(int i=1;i<>
{
for(int j=1;j<>
printf('%d ',rand()*10+50000);
printf('\n');
}
for(int i=1;i<>
{
for(int j=1;j<>
printf('%d ',rand()*10+50000);
printf('\n');
}
// printf('Elapsed time:%.2lf secs.\n',1.0*clock()/1000.0);
return 0;
}