矩阵快速幂
我们先介绍矩阵乘法，a为一个m*p矩阵，b为一个p*n矩阵，矩阵a*矩阵b=矩阵c(m*n矩阵)，c[i][j]=sum{a[i][k]*b[k][j]},(1<=k<=p)。矩阵快速幂即是将矩阵乘法与普通快速幂结合。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
using namespace std;
inline ll read()
{
ll l=0,w=1;
char ch=0;
while(ch<'0'||ch>'9')
{
if(ch=='-')
w=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9')
{
l=(l<<3)+(l<<1)+(ch^'0');
ch=getchar();
}
return w*l;
}
struct pp
{
ll matrix[120][120];
void clear(){memset(matrix,0,sizeof(matrix));}
}a,b,ans;
ll n;
inline void copy(pp &x,pp &y)
{
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
x.matrix[i][j]=y.matrix[i][j];
}
inline void mul(pp &x,pp &y,pp &z)
{
z.clear();
for(int i=1;i<=n;i++)
for(int k=1;k<=n;k++)
if(x.matrix[i][k])
for(int j=1;j<=n;j++)
{
z.matrix[i][j]+=x.matrix[i][k]*y.matrix[k][j];
z.matrix[i][j]%=1000000007;
}
}
void kasumi(ll k)
{
if(k==1)
{
copy(ans,a);
return;
}
kasumi(k>>1);
mul(ans,ans,b);
if(k&1)
{
mul(a,b,ans);
return;
}
copy(ans,b);
}
int main()
{
ll k;
n=read(),k=read();
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
a.matrix[i][j]=read();
kasumi(k);
for(int i=1;i<=n;i++)
{
for(int j=1;j<=n;j++)
printf("%d ",ans.matrix[i][j]);
printf("\n");
}
return 0;
}
矩阵加速
矩阵可以加速形如F[n]=F[n-1]+F[n-2]+…+F[1]等数列的计算。以F[n]=F[n-1]+F[n-3](n>3),F[1]=F[2]=F[3]=1为例，F[n]=F[n-1]*1+F[n-2]*0+F[n-3]*1,F[n-1]=F[n-1]*1+F[n-2]*0+F[n-3]*0,F[n-2]=F[n-1]*0+F[n-2]*1+F[n-3]*0。于是我们构造出一个3x3的矩阵a[{1,0,1},{1,0,0},{0,1,0}]与一个3x1矩阵b{1,1,1},矩阵{F[n],F[n-1],F[n-2]}=a^(n-3)*b，其中a^(n-3)的计算可以用矩阵快速幂加速。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
using namespace std;
inline ll read()
{
ll ans=0;
char ch=0;
while(ch<'0'||ch>'9')
ch=getchar();
while(ch>='0'&&ch<='9')
{
ans=(ans<<3)+(ans<<1)+(ch^'0');
ch=getchar();
}
return ans;
}
struct pp
{
ll matrix[4][4];
void init(){matrix[1][1]=1,matrix[1][2]=1,matrix[2][3]=1,matrix[3][1]=1;}
void clear(){memset(matrix,0,sizeof(matrix));}
}a,b,ans;
inline void copy(pp &x,pp &y)
{
for(int i=1;i<=3;i++)
for(int j=1;j<=3;j++)
x.matrix[i][j]=y.matrix[i][j];
}
inline void mul(pp &x,pp &y,pp &z)
{
z.clear();
for(int i=1;i<=3;i++)
for(int k=1;k<=3;k++)
if(x.matrix[i][k])
for(int j=1;j<=3;j++)
{
z.matrix[i][j]+=x.matrix[i][k]*y.matrix[k][j];
z.matrix[i][j]%=mod;
}
}
void kasumi(ll k)
{
if(k==1)
{
copy(ans,a);
return;
}
kasumi(k>>1);
mul(ans,ans,b);
if(k&1)
{
mul(a,b,ans);
return;
}
copy(ans,b);
}
int main()
{
ll T,n,k,f[4]={0,1,1,1};
a.init();
T=read();
while(T--)
{
k=0;
n=read();
if(n<=3)
{
printf("%lld\n",f[n]);
continue;
}
kasumi(n-3);
for(int i=1;i<=3;i++)
k+=f[i]*ans.matrix[i][1];
printf("%lld\n",k%mod);
}
return 0;
}
高斯消元法
先将多元方程组和解一起列成矩阵形式，从第一行开始将a[i][i]化为1，该方程乘以相同倍数化为等价式子，再用该式子把所有a[j][i]消去(j>i),由此得到一个倒三角，倒三角顶角即为该未知数的解，再将该未知数代入上层方程，即可求出所有未知数。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
using namespace std;
double a[120][120],ans[120];
int n;
inline void mul(int x,double y)
{
for(int i=x;i<=n+1;i++)
a[x][i]*=y;
}
inline void Minus(int x,int y)
{
for(int i=1;i<=n+1;i++)
a[x][i]-=a[y][i];
}
int main()
{
cin>>n;
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)
cin>>a[i][j];
for(int i=1;i<=n;i++)
{
if(!a[i][i])
{
printf("No Solution\n");
return 0;
}
for(int j=i+1;j<=n;j++)
{
double temp=a[j][i]/a[i][i];
mul(i,temp);
Minus(j,i);
}
}
ans[n]=a[n][n+1]/a[n][n];
for(int i=n-1;i>0;i--)
{
for(int j=i+1;j<=n;j++)
a[i][n+1]-=ans[j]*a[i][j];
ans[i]=a[i][n+1]/a[i][i];
}
for(int i=1;i<=n;i++)
printf("%.2lf\n",ans[i]);
return 0;
}