先介绍一篇矩阵好的博文 Matrix 67：http://www.matrix67.com/blog/archives/276
一．高斯消元
我觉得不错的模板
// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解，但无整数解，
//-1表示无解，0表示唯一解，大于0表示无穷解，并返回自由变元的个数)
//有equ个方程，var个变元。增广矩阵行数为equ,分别为0到equ-1,列数为var+1,分别为0到var.
int Gauss(int equ,int var){
int i,j,k;
int max_r;// 当前这列绝对值最大的行.
int col;//当前处理的列
int ta,tb;
int LCM;
int temp;
int free_x_num;
int free_index;
for(int i=0;i<=var;i++){
x[i]=0;
free_x[i]=true;
}
//转换为阶梯阵.
col=0; // 当前处理的列
for(k = 0;k < equ && col < var;k++,col++){// 枚举当前处理的行.
// 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
max_r=k;
for(i=k+1;i<equ;i++)
if(abs(a[i][col])>abs(a[max_r][col]))
max_r=i;
if(max_r!=k)// 与第k行交换
for(j=k;j<var+1;j++)
swap(a[k][j],a[max_r][j]);
if(a[k][col]==0){// 说明该col列第k行以下全是0了，则处理当前行的下一列.
k--;
continue;
}
for(i=k+1;i<equ;i++){// 枚举要删去的行.
if(a[i][col]!=0){
LCM = lcm(abs(a[i][col]),abs(a[k][col]));
ta = LCM/abs(a[i][col]);
tb = LCM/abs(a[k][col]);
if(a[i][col]*a[k][col]<0)tb=-tb;//异号的情况是相加
for(j=col;j<var+1;j++)
a[i][j] = a[i][j]*ta-a[k][j]*tb;
}
}
}
// 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
// 对于无穷解来说，如果要判断哪些是自由变元，那么初等行变换中的交换就会影响，则要记录交换.
for (i = k; i < equ; i++)
if (a[i][col] != 0) return -1;
// 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行，即说明没有形成严格的上三角阵.
// 且出现的行数即为自由变元的个数.
if (k < var){
// 首先，自由变元有var - k个，即不确定的变元至少有var - k个.
for (i = k - 1; i >= 0; i--){
// 第i行一定不会是(0, 0, ..., 0)的情况，因为这样的行是在第k行到第equ行.
// 同样，第i行一定不会是(0, 0, ..., a), a != 0的情况，这样的无解的.
free_x_num = 0; // 用于判断该行中的不确定的变元的个数，如果超过1个，则无法求解，它们仍然为不确定的变元.
for (j = 0; j < var; j++)
if (a[i][j] != 0 && free_x[j])
free_x_num++, free_index = j;
if (free_x_num > 1) continue; // 无法求解出确定的变元.
// 说明就只有一个不确定的变元free_index，那么可以求解出该变元，且该变元是确定的.
temp = a[i][var];
for (j = 0; j < var; j++)
if (a[i][j] != 0 && j != free_index)
temp -= a[i][j] * x[j];
x[free_index] = temp / a[i][free_index]; // 求出该变元.
free_x[free_index] = 0; // 该变元是确定的.
}
return var - k; // 自由变元有var - k个.
}
// 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
// 计算出Xn-1, Xn-2 ... X0.
for (i = var - 1; i >= 0; i--){
temp = a[i][var];
for (j = i + 1; j < var; j++)
if (a[i][j] != 0)
temp -= a[i][j] * x[j];
if (temp % a[i][i] != 0) return -2; // 说明有浮点数解，但无整数解.
x[i] = temp / a[i][i];
}
return 0;
}
1.POJ 1222 EXTENDED LIGHTS OUT
题目：关灯问题，按下开关，该灯和相邻的灯都会变化灯的状态。现在给出初始状态，问如何设置开关，使得灯的最终状态全为关闭的。。
分析：我们发现对于i灯，必有ai ^ xi ^ (A) = 0，相当于xi ^ (A) = ai
A为相邻的几个变量xj ^ xk ^ xl & xm(假设有四个)
所以我们可以列出30个方程组出来。然后解方程组就行了
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define rep(i,n) for(int i=0;i<n;i++)
#define foreach(i,vec) for(unsigned i=0;i<vec.size();i++)
#define pb push_back
const int X = 35;
const int n = 30;
int a[X][X];
int dirx[4] = {-1,0,0,1};
int diry[4] = {0,-1,1,0};
void build(){
rep(i,5){
rep(j,6){
int x = i*6+j;
a[x][x] = 1;
rep(k,4){
int dx = dirx[k]+i;
int dy = diry[k]+j;
if(dx>=0&&dx<5 && dy>=0&&dy<6){
int y = dx*6+dy;
a[x][y] = 1;
}
}
}
}
}
void debug(){
rep(i,30){
rep(j,31)
cout<<a[i][j]<<" ";
cout<<endl;
}
}
void gauss(){
int i = 0,j = 0;
while(i<n&&j<n){
int r = i;
for(int k=i;k<n;k++)
if(a[k][j]){
r = k;
break;
}
if(a[r][j]){
if(r!=i)
rep(k,n+1)
swap(a[r][k],a[i][k]);
for(int u=i+1;u<n;u++)
if(a[u][j])
for(int k=j;k<n+1;k++)
a[u][k] ^= a[i][k];
i ++;
}
j ++;
}
for(int i=n-2;i>=0;i--)
for(int j=n-1;j>i;j--)
a[i][n] ^= (a[i][j]&&a[j][n]);
}
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
int ncase;
cin>>ncase;
rep(cnt,ncase){
printf("PUZZLE #%d\n",cnt+1);
memset(a,0,sizeof(a));
rep(i,n)
scanf("%d",&a[i][n]);
build();
gauss();
rep(i,5){
printf("%d",a[i*6][30]);
for(int j=1;j<6;j++)
printf(" %d",a[i*6+j][30]);
puts("");
}
}
return 0;
}
相似的几题
POJ 1681 Painter's Problem
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define rep(i,n) for(int i=0;i<n;i++)
#define foreach(i,vec) for(unsigned i=0;i<vec.size();i++)
#define pb push_back
const int X = 230;
int n,m;
int dirx[] = {0,0,-1,1};
int diry[] = {-1,1,0,0};
int a[X][X];
void debug(){
rep(i,n){
rep(j,n+1)
cout<<a[i][j]<<" ";
cout<<endl;
}
}
void build(){
rep(i,m){
rep(j,m){
int x = i*m+j;
a[x][x] = 1;
rep(k,4){
int dx = dirx[k]+i;
int dy = diry[k]+j;
if(dx>=0&&dx<m && dy>=0&&dy<m){
int y = dx*m+dy;
a[x][y] = 1;
}
}
}
}
}
int gauss(){
int i = 0,j = 0;
while(i<n&&j<n){
int r = i;
for(int k=i;k<n;k++)
if(a[k][j]){
r = k;
break;
}
if(a[r][j]){
if(r!=i)
rep(k,n+1)
swap(a[r][k],a[i][k]);
for(int u=i+1;u<n;u++)
if(a[u][j])
for(int k=j;k<n+1;k++)
a[u][k] ^= a[i][k];
i ++;
}
j ++;
}
//cout<<"dsaaaaaa "<<i<<endl;
for(;i<n;i++)
if(a[i][n])
return -1;
for(int i=n-2;i>=0;i--)
for(int j=n-1;j>i;j--)
a[i][n] ^= (a[i][j]&&a[j][n]);
int cnt = 0;
for(int i=0;i<n;i++)
cnt += a[i][n];
return cnt;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
int ncase;
cin>>ncase;
while(ncase--){
memset(a,0,sizeof(a));
cin>>m;
char s[16];
n = m*m;
rep(i,m){
scanf("%s",s);
rep(j,m)
a[i*m+j][n] = (s[j]=='w');
}
build();
int ans = gauss();
if(ans==-1)
puts("inf");
else
printf("%d\n",ans);
}
return 0;
}
POJ 1830 开关问题
/*
题目：
同关灯问题，高斯消元
*/
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define rep(i,n) for(int i=0;i<n;i++)
#define foreach(i,vec) for(unsigned i=0;i<vec.size();i++)
#define pb push_back
const int X = 32;
int n;
int dirx[] = {0,0,-1,1};
int diry[] = {-1,1,0,0};
int a[X][X];
void debug(){
rep(i,n){
rep(j,n+1)
cout<<a[i][j]<<" ";
cout<<endl;
}
}
int gauss(){
int i = 0,j = 0;
while(i<n&&j<n){
int r = i;
for(int k=i;k<n;k++)
if(a[k][j]){
r = k;
break;
}
if(a[r][j]){
if(r!=i)
rep(k,n+1)
swap(a[r][k],a[i][k]);
for(int u=i+1;u<n;u++)
if(a[u][j])
for(int k=j;k<n+1;k++)
a[u][k] ^= a[i][k];
i ++;
}
j ++;
}
int cnt = n-i;
for(;i<n;i++)
if(a[i][n])
return -1;
return 1<<cnt;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
int ncase;
cin>>ncase;
while(ncase--){
cin>>n;
memset(a,0,sizeof(a));
int x[X],y;
rep(i,n)
scanf("%d",&x[i]);
rep(i,n){
scanf("%d",&y);
a[i][n] = x[i] ^ y;
a[i][i] = 1;
}
int p,q;
while(scanf("%d%d",&p,&q),p||q)
a[--q][--p] = 1;
int ans = gauss();
if(ans==-1)
puts("Oh,it's impossible~!!");
else
printf("%d\n",ans);
}
return 0;
}
URAL 1042 Central Heating
/*
题目：
同关灯问题，高斯消元
*/
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define rep(i,n) for(int i=0;i<n;i++)
#define foreach(i,vec) for(unsigned i=0;i<vec.size();i++)
#define pb push_back
const int X = 300;
int n;
int a[X][X];
void debug(){
rep(i,n){
rep(j,n+1)
cout<<a[i][j]<<" ";
cout<<endl;
}
}
int gauss(){
int i = 0,j = 0;
while(i<n&&j<n){
int r = i;
for(int k=i;k<n;k++)
if(a[k][j]){
r = k;
break;
}
if(a[r][j]){
if(r!=i)
rep(k,n+1)
swap(a[r][k],a[i][k]);
for(int u=i+1;u<n;u++)
if(a[u][j])
for(int k=j;k<n+1;k++)
a[u][k] ^= a[i][k];
i ++;
}
j ++;
}
//cout<<"dsaaaaaa "<<i<<endl;
for(;i<n;i++)
if(a[i][n])
return -1;
for(int i=n-2;i>=0;i--)
for(int j=n-1;j>i;j--)
a[i][n] ^= (a[i][j]&&a[j][n]);
int cnt = 0;
for(int i=0;i<n;i++)
cnt += a[i][n];
return cnt;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
while(cin>>n){
memset(a,0,sizeof(a));
rep(i,n){
int y;
a[i][n] = 1;
while(scanf("%d",&y),y!=-1)
a[--y][i] = 1;
}
int ans = gauss();
if(ans==-1)
puts("No solution");
else{
bool ok = false;
rep(i,n)
if(a[i][n]){
ok?putchar(' '):ok = true;
printf("%d",i+1);
}
puts("");
}
}
return 0;
}
2.BZOJ 1013 [JSOI2008]球形空间产生器sphere
题目：中文题。。
分析：裸的高斯消元题
二维平面上的圆上的点与圆心的距离有(x-a)^2+(y-b)^2 = r^2
三维空间上的球上的点与球心的距离有(x-a)^2+(y-b)^2+(z-c)^2 = r^2
同理：在n维空间上的球上的点与球心的距离有sigma((xi-ai)^2) = r^2，圆心为(a1,a2,...,an)
另外，在二维平面上，可由三点（不共线）确定一个园，在三维上四点（不共线）确定一个球，同理，在n维平面上，可由n+1个点（不共线）确定一个n维的球。
这样，题目就可以转化为n+1个方程组，但是是平方级别的，如何转化为一维的？
我们不妨对于相邻的两个方程组左右分别相减，可以发现：
2*(x21 - x11)*x1 + 2*(x22 - x12)*x2 +...+2*(x2n - x1n) = (x21^2 - x11^1)+...+(x2n^2 - x1n^2)
这样，由n+1个方程组就可以转化为了n个一维的方程组了。下面，直接用高斯消元法即 可解决该问题
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
using namespace std;
const int X = 20;
#define esp 1e-8
double dp[X][X];
double a[X][X],b[X][X];
double f[X];
int n;
double gauss(){
for(int i=1;i<=n;i++){
int x = i;
for(int j=i+1;j<=n;j++)
if(fabs(a[j][i]-a[x][i])>esp)
x = j;
if(x!=i)
for(int j=1;j<=n+1;j++)
swap(a[i][j],a[x][j]);
for(int j=i+1;j<=n;j++)
if(fabs(a[j][i])>esp){
double temp = a[j][i] / a[i][i];
for(int k=i;k<=n+1;k++)
a[j][k] -= temp*a[i][k];
}
}
for(int i=n;i;i--){
double temp = a[i][n+1];
for(int j=i+1;j<=n;j++)
temp -= f[j]*a[i][j];
f[i] = temp / a[i][i];
}
}
int main(){
cin>>n;
for(int i=1;i<n+2;i++)
for(int j=1;j<=n;j++)
scanf("%lf",&b[i][j]);
for(int i=1;i<=n;i++){
double temp = 0;
for(int j=1;j<=n;j++){
temp += b[i+1][j]*b[i+1][j]-b[i][j]*b[i][j];
a[i][j] = 2*(b[i+1][j]-b[i][j]);
}
a[i][n+1] = temp;
}
gauss();
for(int i=1;i<n;i++)
printf("%.3lf ",f[i]);
printf("%.3lf\n",f[n]);
return 0;
}
二．矩阵快速幂
1.poj 3070 Fibonacci
题目：求斐波那契数列第n项的后四位数
分析：由于n很大，我们可以构造矩阵A，B
A = 1 1 B = f1
1 0 f0
则 A^(n-1) * B = fn
fn-1
所以我们可以直接利用矩阵快速幂来求fn%10000
/*
题目：
斐波那契数列的矩阵算法
分析：
裸的矩阵快速幂
*/
#include <iostream>
#include <cstring>
#include <cstdio>
using namespace std;
const int X = 5;
const int mod = 10000;
int n;
class matrix{
public:
int a[X][X];
int size;
int mod;
matrix(){
memset(a,0,sizeof(a));
}
matrix(int _size,int _mod):size(_size),mod(_mod){
memset(a,0,sizeof(a));
}
void setE(){
for(int i=0;i<size;i++)
a[i][i] = 1;
}
matrix operator * (matrix p){
matrix c(size,mod);
for(int i=0;i<size;i++)
for(int j=0;j<size;j++)
for(int k=0;k<size;k++){
c.a[i][j] += a[i][k]*p.a[k][j];
c.a[i][j] %= mod;
}
return c;
}
matrix pow(int exp){
matrix cur = *this;
matrix c(size,mod);
c.setE();
while(exp){
if(exp&1)
c = c*cur;
cur = cur*cur;
exp = exp>>1;
}
return c;
}
void display(){
for(int i=0;i<size;i++){
for(int j=0;j<size;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
}
};
int main(){
freopen("sum.in","r",stdin);
while(scanf("%d",&n),n!=-1){
if(!n){
puts("0");
continue;
}
matrix ans(2,mod);
ans.a[0][0] = ans.a[0][1] = ans.a[1][0] = 1;
ans.a[1][1] = 0;
ans = ans.pow(n-1);
printf("%d\n",ans.a[0][0]);
}
return 0;
}
相似的题目:
hdu 1005 number sequence
/*
题目：
f(1) = 1, f(2) = 1, f(n) = (A * f(n - 1) + B * f(n - 2)) mod 7.
给出A,B求f(n)模7
分析：
我们可以构造一个矩阵
[ a b ] * [ fn-1 ] = [ fn ]
[ 1 0 ] [ fn-2 ] [ fn-1 ]
最后发现最要求左边的矩阵的(n-2)次幂后所得的上面两项的和值就是
fn，所以用到了矩阵的快速幂可以做~~
*/
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
const int X = 3;
class matrix{
public:
int a[X][X];
matrix(){
memset(a,0,sizeof(a));
}
matrix(int _size,int _mod):size(_size),mod(_mod){
memset(a,0,sizeof(a));
}
matrix operator * (matrix p){
matrix c(size,mod);
for(int i=0;i<size;i++)
for(int j=0;j<size;j++)
for(int k=0;k<size;k++){
c.a[i][j] += a[i][k]*p.a[k][j];
c.a[i][j] %= mod;
}
return c;
}
void setE(){
for(int i=0;i<size;i++)
a[i][i] = 1;
}
matrix pow(int exp){
matrix temp(size,mod);
temp.setE();
matrix cur = *this;
while(exp){
if(exp&1)
temp = temp*cur;
cur = cur*cur;
exp = exp>>1;
}
return temp;
}
private:
int size;
int mod;
};
int main(){
freopen("sum.in","r",stdin);
int a,b,n;
while(cin>>a>>b>>n,a||b||n){
if(n==1){
cout<<1<<endl;
continue;
}
else if(n==2){
cout<<1<<endl;
continue;
}
matrix ans(2,7);
ans.a[0][0] = a;
ans.a[0][1] = b;
ans.a[1][0] = 1;
ans.a[1][1] = 0;
ans = ans.pow(n-2);
cout<<(ans.a[0][0]+ans.a[0][1])%7<<endl;
}
return 0;
}
2.HOJ 1991 Happy 2005
题目：
给出n，问2005^n的各个因子数之和对29取模
分析：
2005 = 5*401，我们可以对于401进行分类：
401^0 : 1 5 5^2 ... 5^n
401^1 : 401 401*5 401*5^2 ... 401*5^n
.
.
.
401^n : 401^n 401^n*5 401^n*5^2 ... 401^n*5^n
由此我们可以发现，问题可以转换为
(1+401+401^2+...401^n)*(1+5+5^2+...+5^n)%29
方法一：
二分再二分。首先，a^n用一次二分，求和的时候再用一次二分。
a^n二分的时候就是快速幂。
求和二分：
A+A^2+A...+A^(2k+1)= A+A^2+...+A^k+A^(k+1)+A^(k+1)*(A+A^2+...+A^k).
A+A^2+...+A^2k = A+A^2+...+A^k+A^k*(A+A^2+...+A^k).
方法二：
构造矩阵matrix如下：
A 1
0 1
我们发现matrix^(n+1)项的时候，第一行第二列就是问题所求
所以在求A+A^2+A^3+...+A^k % 29的时候，我们可以直接转化为对矩阵进行
快速幂取模。
我下面的代码为构造矩阵求解。。。
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
typedef long long ll;
#define debug puts("here");
const int X = 3;
const int MOD = 29;
struct Matrix{
ll a[X][X];
int n;
int mod;
Matrix(){}
Matrix(int _n,int _mod):n(_n),mod(_mod){
memset(a,0,sizeof(a));
}
void setE(){
memset(a,0,sizeof(a));
for(int i=1;i<=n;i++)
a[i][i] = 1;
}
Matrix operator * (Matrix p){
Matrix c(n,mod);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++){
c.a[i][j] += a[i][k]*p.a[k][j];
c.a[i][j] %= mod;
}
return c;
}
Matrix pw(int exp){
Matrix cur = *this;
Matrix c(n,MOD);
c.setE();
while(exp>0){
if(exp&1)
c = c*cur;
cur = cur*cur;
exp = exp>>1;
}
return c;
}
void di(){
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
}
void init(int x){
a[1][2] = a[2][2] = 1;
a[2][1] = 0;
a[1][1] = x;
}
}matrix;
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
ll n;
while(cin>>n,n){
Matrix a = Matrix(2,MOD);
a.init(5);
a = a.pw(n+1);
ll ans = a.a[1][2]%MOD;
a.init(401);
a = a.pw(n+1);
ans = ans*a.a[1][2]%MOD;
cout<<ans<<endl;
}
return 0;
}
相似的题目：
hoj 2635 Weights
/*
题目：
直接求1+3^1+...+3^n的和
分析：
构造矩阵matrix如下：
A 1
0 1
我们发现matrix^(n+1)项的时候，第一行第二列就是问题所求
所以在求A+A^2+A^3+...+A^k % p的时候，我们可以直接转化为对矩阵进行
快速幂取模。
我下面的代码为构造矩阵求解。。。
*/
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
typedef long long ll;
#define debug puts("here");
const int X = 3;
const int MOD = 9999997;
struct Matrix{
ll a[X][X];
int n;
int mod;
Matrix(){}
Matrix(int _n,int _mod):n(_n),mod(_mod){
memset(a,0,sizeof(a));
}
void setE(){
memset(a,0,sizeof(a));
for(int i=1;i<=n;i++)
a[i][i] = 1;
}
Matrix operator * (Matrix p){
Matrix c(n,mod);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++){
c.a[i][j] += a[i][k]*p.a[k][j];
c.a[i][j] %= mod;
}
return c;
}
Matrix pw(int exp){
Matrix cur = *this;
Matrix c(n,MOD);
c.setE();
while(exp>0){
if(exp&1)
c = c*cur;
cur = cur*cur;
exp = exp>>1;
}
return c;
}
void di(){
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
}
void init(int x){
a[1][2] = a[2][2] = 1;
a[2][1] = 0;
a[1][1] = x;
}
}matrix;
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
ll n;
while(cin>>n){
Matrix a = Matrix(2,MOD);
a.init(3);
a = a.pw(n);
//a.di();
ll ans = a.a[1][2]%MOD;
cout<<ans<<endl;
}
return 0;
}
poj 3233 Matrix Power Series
/*
题目:
求出 S = A + A^2 + A^3 + … + A^k.
分析：
解法一
Let B= A I
0 I
B^(k+1) = A^k I+A+...+A^k
0 I
解法二
设f[n]=A^1+A^2+....A^n;
当n是偶数,f[n]=f[n/2]+f[n/2]*A^(n/2);
但n是奇数,f[n]=f[n-1]+A^(n);
*/
#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
const int X = 31<<2;
int n,m,k;
class matrix{
public:
int size;
int mod;
int a[X][X];
matrix(){
memset(a,0,sizeof(a));
}
matrix(int _size,int _mod):size(_size),mod(_mod){
memset(a,0,sizeof(a));
}
void setE(){
for(int i=0;i<2*size;i++)
a[i][i] = 1;
}
void print(){
for(int i=0;i<size;i++){
printf("%d",a[i][size]);
for(int j=1;j<size;j++)
printf(" %d",a[i][j+size]);
puts("");
}
}
matrix operator * (matrix p){
matrix c(size,mod);
for(int i=0;i<2*size;i++)
for(int j=0;j<2*size;j++)
for(int k=0;k<2*size;k++){
c.a[i][j] += a[i][k]*p.a[k][j];
c.a[i][j] %= mod;
}
return c;
}
void operator -- (){
for(int i=0;i<size;i++)
a[i][i+size] = (--a[i][i+size]+mod)%mod;
}
matrix pow(int exp){
matrix temp(size,mod);
temp.setE();
matrix cur = *this;
while(exp){
if(exp&1)
temp = temp*cur;
cur = cur*cur;
exp = exp>>1;
}
return temp;
}
void display(){
for(int i=0;i<2*size;i++){
for(int j=0;j<2*size;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
}
};
int main(){
freopen("sum.in","r",stdin);
while(cin>>n>>k>>m){
matrix ans(n,m);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
scanf("%d",&ans.a[i][j]);
for(int i=n;i<2*n;i++)
ans.a[i-n][i] = ans.a[i][i] = 1;
ans = ans.pow(k+1);
--ans;
ans.print();
}
return 0;
}
3.线性递推
对于线性递推f[k] = a*f[k-1]+b*f[k-2]+c*f[k-3]...z*f[0]
我们可以构造矩阵A,B:
A = 0 1 0... B = f[k-1]
0 0 1... f[k-2]
.......... ........
a b c... f[0]
我们可以发现
A*B = f[k]
f[k-1]
........
f[1]
所以A^(n-k) * B = f[n]
f[n-1]
....
f[n-k+1]
HOJ 1790 Firepersons
题目:
线性递推关系，an=Σ1<=i<=k*an-i*bi，问ai
分析：
可以构造矩阵A如下
0 1 0 ...0
0 0 1 ...0
...
0 0 0 ...1
bk bk-1 bk-2...b0
矩阵B为
a0
a1
a2
...
ak-1
则ai = ai(i<k)
= A^(i-k+1)*B最后一行的元素
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define debug puts("here")
const int MOD = 10000;
const int X = 102;
class Matrix{
public:
int n,m;
int a[X][X];
Matrix(){
memset(a,0,sizeof(a));
}
Matrix(int _n,int _m):n(_n),m(_m){
memset(a,0,sizeof(a));
}
Matrix operator * (Matrix p){
int q = p.m;
Matrix c(n,q);
for(int i=0;i<n;i++)
for(int j=0;j<q;j++)
for(int k=0;k<m;k++)
c.a[i][j] = (c.a[i][j]+a[i][k]*p.a[k][j])%MOD;
return c;
}
void getE(){
memset(a,0,sizeof(a));
for(int i=0;i<n;i++)
a[i][i] = 1;
}
Matrix bin(int exp){
Matrix temp(n,n);
temp.getE();
Matrix cur = *this;
while(exp>0){
if(exp&1)
temp = temp*cur;
cur = cur*cur;
exp = exp >> 1;
}
return temp;
}
void di(){
for(int i=0;i<n;i++){
for(int j=0;j<m;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
};
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
int n;
while(cin>>n,n){
Matrix a = Matrix(n,n);
Matrix b = Matrix(n,1);
for(int i=0;i<n-1;i++)
a.a[i][i+1] = 1;
for(int i=0;i<n;i++)
scanf("%d",&b.a[i][0]);
for(int i=0;i<n;i++)
scanf("%d",&a.a[n-1][n-i-1]);
int exp;
cin>>exp;
if(exp<n){
printf("%d\n",b.a[exp][0]);
continue;
}
exp ++;
exp -= n;
a = a.bin(exp);
a = a*b;
printf("%d\n",a.a[n-1][0]);
}
return 0;
}
BZOJ 2875: [ NOI2012 ] 随机数生成器
http://www.lydsy.com/JudgeOnline/problem.php?id=2875
题目：f[n+1] = (f[n]+c)%m , 给出f[0],n,m,c,g，求f[n]%g
分析：很明显可以构造矩阵
A = a 1 B = f[0]
0 1 c
但是由于数据太大了，所以在矩阵乘法中间计算时会溢出，所以我们需要在乘的时候做一下处理。改为跟快速幂乘法相似的计算方式来防止溢出。具体看代码
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef unsigned long long ll;
#define debug puts("here")
#define rep(i,n) for(int i=0;i<n;i++)
#define foreach(i,vec) for(unsigned i=0;i<vec.size();i++)
#define pb push_back
const int X = 5;
ll MOD;
class Matrix{
public:
int n,m;
ll a[X][X];
ll cal(ll x,ll y){
ll sum = 0;
while(y>0){
if(y&1)
sum = (sum+x)%MOD;
x = (x<<1)%MOD;
y >>= 1;
}
return sum;
}
Matrix(){
memset(a,0,sizeof(a));
}
Matrix(int _n,int _m):n(_n),m(_m){
memset(a,0,sizeof(a));
}
Matrix operator * (Matrix p){
int q = p.m;
Matrix c(n,q);
for(int i=0;i<n;i++)
for(int j=0;j<q;j++)
for(int k=0;k<m;k++)
c.a[i][j] = (c.a[i][j]+cal(a[i][k],p.a[k][j]))%MOD;
return c;
}
void getE(){
memset(a,0,sizeof(a));
for(int i=0;i<n;i++)
a[i][i] = 1;
}
Matrix bin(ll exp){
Matrix temp(n,n);
temp.getE();
Matrix cur = *this;
while(exp>0){
if(exp&1)
temp = temp*cur;
cur = cur*cur;
exp = exp >> 1;
}
return temp;
}
void di(){
for(int i=0;i<n;i++){
for(int j=0;j<m;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
};
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
ll m,a,c,x0,n,g;
while(cin>>m>>a>>c>>x0>>n>>g){
MOD = m;
Matrix ans(2,2);
ans.a[0][0] = a;
ans.a[0][1] = ans.a[1][1] = 1;
ans = ans.bin(n);
//ans.di();
Matrix temp(2,1);
temp.a[0][0] = x0;
temp.a[1][0] = c;
ans = ans*temp;
cout<<ans.a[0][0]%g<<endl;
}
return 0;
}
HOJ 2060 Fibonacci Problem Again
题目：
计算斐波那契数列[a,b]的和值对于1e9取模
分析：
对于斐波那契求第n项，我们可以构造矩阵
A = 1 1 B = f[n]
1 0 f[n-1]
则f[n]为矩阵 A^n * B的第一项
对于这题，我们可以额外构造多一维的矩阵出来为
1 1 0 f[n]
A = 1 0 0 B = f[n-1]
1 0 1 sum[n-1]
我们同样可以算出 sum[n] = A^n * B 的第三项
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define debug puts("here")
#define rep(i,n) for(int i=0;i<n;i++)
#define foreach(i,vec) for(unsigned i=0;i<vec.size();i++)
#define pb push_back
const int X = 5;
const ll MOD = 1e9;
class Matrix{
public:
int n,m;
ll a[X][X];
Matrix(){
memset(a,0,sizeof(a));
}
Matrix(int _n,int _m):n(_n),m(_m){
memset(a,0,sizeof(a));
}
Matrix operator * (Matrix p){
int q = p.m;
Matrix c(n,q);
for(int i=0;i<n;i++)
for(int j=0;j<q;j++)
for(int k=0;k<m;k++)
c.a[i][j] = (c.a[i][j]+a[i][k]*p.a[k][j])%MOD;
return c;
}
void getE(){
memset(a,0,sizeof(a));
for(int i=0;i<n;i++)
a[i][i] = 1;
}
Matrix bin(int exp){
Matrix temp(n,n);
temp.getE();
Matrix cur = *this;
while(exp>0){
if(exp&1)
temp = temp*cur;
cur = cur*cur;
exp = exp >> 1;
}
return temp;
}
void di(){
for(int i=0;i<n;i++){
for(int j=0;j<m;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
};
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
int n,m;
while(scanf("%d%d",&n,&m),n||m){
Matrix a = Matrix(3,3);
a.a[0][0] = a.a[0][1] = a.a[1][0] = a.a[2][0] = a.a[2][2] = 1;
Matrix f = Matrix(3,1);
f.a[0][0] = 1;
f.a[1][0] = 1;
f.a[2][0] = 1;
ll pre = 0;
ll now = 0;
if(n){
Matrix x = a.bin(n-1);
x = x*f;
pre = x.a[2][0];
}
Matrix y = a.bin(m);
y = y*f;
now = y.a[2][0];
ll mod = ll(1000000000);
cout<<(now-pre+mod)%mod<<endl;
}
return 0;
}
（HOJ 2255 Not Fibonacci这题跟上面的题目基本一模一样，代码略）
HOJ 2930 Perfect Fill IIl
详情看上一篇博文 http://www.cnblogs.com/yejinru/archive/2013/02/01/2888368.html
4.HDU 2157 How many ways
题目：
给定一个有向图，问从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 <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define debug puts("here")
#define rep(i,n) for(int i=0;i<n;i++)
#define foreach(i,vec) for(unsigned i=0;i<vec.size();i++)
#define pb push_back
const int MOD = 1000;
const int X = 25;
class Matrix{
public:
int n,m;
int a[X][X];
Matrix(){
memset(a,0,sizeof(a));
}
Matrix(int _n,int _m):n(_n),m(_m){
memset(a,0,sizeof(a));
}
Matrix operator * (Matrix p){
int q = p.m;
Matrix c(n,q);
for(int i=0;i<n;i++)
for(int j=0;j<q;j++)
for(int k=0;k<m;k++)
c.a[i][j] = (c.a[i][j]+a[i][k]*p.a[k][j])%MOD;
return c;
}
void getE(){
memset(a,0,sizeof(a));
for(int i=0;i<n;i++)
a[i][i] = 1;
}
Matrix bin(int exp){
Matrix temp(n,n);
temp.getE();
Matrix cur = *this;
while(exp>0){
if(exp&1)
temp = temp*cur;
cur = cur*cur;
exp = exp >> 1;
}
return temp;
}
void di(){
for(int i=0;i<n;i++){
for(int j=0;j<m;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
};
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
int n,m,x,y,k;
while(scanf("%d%d",&n,&m),n||m){
Matrix a = Matrix(n,n);
while(m--){
scanf("%d%d",&x,&y);
a.a[x][y] = 1;
}
int cnt;
cin>>cnt;
while(cnt--){
scanf("%d%d%d",&x,&y,&k);
Matrix temp = a.bin(k);
printf("%d\n",temp.a[x][y]);
}
}
return 0;
}
不会分类的几题：
vijos 1049 送给圣诞夜的礼品
顺次给出m个置换，反复使用这m个置换对初始序列进行操作，问k次置换后的序列。m<=10, k<2^31。
首先将这m个置换“合并”起来（算出这m个置换的乘积），然后接下来我们需要执行这个置换k/m次（取整，若有余数则剩下几步模拟即可）。
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define debug puts("here")
const int MAXN = 102;
const int MAXM = 11;
int n,m;
class Matrix{
public:
int a[MAXN][MAXN];
Matrix(){
memset(a,0,sizeof(a));
}
void setE(){
memset(a,0,sizeof(a));
for(int i=0;i<n;i++)
a[i][i] = 1;
}
Matrix operator * (Matrix b){
Matrix c = Matrix();
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
for(int k=0;k<n;k++)
c.a[i][j] += a[i][k]*b.a[k][j];
return c;
}
Matrix bin(int exp){
Matrix ans;
Matrix cur = *this;
ans.setE();
while(exp>0){
if(exp&1)
ans = ans*cur;
cur = cur*cur;
exp >>= 1;
}
return ans;
}
void di(){
for(int i=0;i<n;i++){
for(int j=0;j<n;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
};
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
int k;
while(cin>>n>>m>>k){
int x;
Matrix ans;
ans.setE();
Matrix ret[12];
ret[0].setE();
for(int i=0;i<m;i++){
Matrix cur;
for(int j=0;j<n;j++){
scanf("%d",&x);
cur.a[j][x-1] = 1;
}
ret[i+1] = cur*ret[i];
//ret[i+1].di();
ans = cur*ans;
}
ans = ans.bin(k/m);
ans = ret[k%m]*ans;
//ans.di();
for(int i=0;i<n;i++)
for(int j=0;j<n;j++){
if(ans.a[i][j]){
if(i)
cout<<" ";
cout<<j+1;
break;
}
}
cout<<endl;
}
return 0;
}
HOJ 2446 Cellular Automaton
题目：
在一个环中有n个格子，每个格子的值为ai，距离该格子不足d的所有格子的和对于
m取余为新的值，问第k次变换后的所有n个格子的值
分析：
很容易可以构造出一个循环的矩阵出来，但是如果是O(n^3*logn)会TLE。
我们可以注意到循环矩阵a * b只需要计算a的第一行*b，然后下面的移位均可以得
到。时间为O(n^2*logn)
#include <set>
#include <map>
#include <cmath>
#include <queue>
#include <stack>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long ll;
#define debug puts("here")
#define rep(i,n) for(int i=0;i<n;i++)
#define foreach(i,vec) for(unsigned i=0;i<vec.size();i++)
#define pb push_back
const int X = 505;
ll a[X][X],c[X][X];
ll ans[X];
int main(){
#ifndef ONLINE_JUDGE
freopen("sum.in","r",stdin);
#endif
int n,m,k,d;
while(~scanf("%d%d%d%d",&n,&m,&d,&k)){
rep(i,n)
scanf("%lld",&ans[i]);
rep(i,n)
rep(j,n)
if( (i-j+n)%n<=d || (j-i+n)%n<=d )
a[i][j] = 1;
else
a[i][j] = 0;
while(k>0){
if(k&1){
rep(i,n){
c[0][i] = 0;
rep(j,n)
c[0][i] += ans[j]*a[j][i];
}
rep(i,n)
ans[i] = c[0][i]%m;
}
rep(i,n){
c[0][i] = 0;
rep(j,n)
c[0][i] += a[0][j]*a[j][i];
}
rep(i,n)
rep(j,n)
if(i==0)
a[i][j] = c[i][j]%m;
else
a[i][j] = a[i-1][(j-1+n)%n];
k >>= 1;
}
rep(i,n-1)
printf("%lld ",ans[i]);
printf("%lld\n",ans[n-1]);
}
return 0;
}
压缩矩阵
poj 3318 Matrix Multiplication
题目：
判断矩阵a * b == c
分析：
方法一：
O(n^3)算法提示会TLE,但是原矩阵是一个稀疏矩阵，所以可
以在相乘的时候判断是否为0，这样同样不会TLE~~
方法二：
压缩矩阵，左乘1*n的矩阵，使得左边以及右边都变成1*n的
矩阵，然后两边判断是否相等就行了~~但是如果这样压缩的话，
不保证每个元素的特性，所以这个1*n的矩阵得要体现特性，构
造的时候可以取随机数，或者令(1,2...n)
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
const int X = 505;
#define debug puts("here");
typedef __int64 ll;
int n;
void mul(ll a[X][X],ll *b){
ll temp[X] = {0};
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
temp[i] += b[j]*a[j][i];
for(int i=0;i<n;i++)
b[i] = temp[i];
}
bool eq(ll a[],ll b[]){
for(int i=0;i<n;i++)
if(a[i]!=b[i])
return false;
return true;
}
ll a[X][X],b[X][X],c[X][X];
ll q[X],p[X];
int main(){
freopen("sum.in","r",stdin);
while(cin>>n){
for(int i=0;i<n;i++)
q[i] = p[i] = i+1;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
scanf("%I64d",&a[i][j]);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
scanf("%I64d",&b[i][j]);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
scanf("%I64d",&c[i][j]);
mul(a,p);
mul(b,p);
mul(c,q);
if(eq(p,q))
puts("YES");
else
puts("NO");
}
return 0;
}
以后若还有的话，会更新一下的。。。