码迷，mamicode.com
首页 > 其他好文 > 详细

矩阵小结

时间：2015-05-28 02:01:42      阅读：324      评论：0      收藏：0      [点我收藏+]

标签：矩阵   数学   

1.矩阵快速幂，用倍增来加速(O(n^3*logk))

2.矩阵求解递推关系第n项（n很大）可以构造矩阵，用矩阵快速幂迅速求出。

3.给定起点和终点求从起点到终点恰好进过k步的方案数可以直接对可达矩阵相乘k次得到结果

4.矩阵乘法的顺序对时间影响比较大（提高Cache命中率），kij最快而且还可以进行稀疏矩阵加速（当a[i][k]为0时没必要进行运算）。

因为最近在搞矩阵，所以准备写一个矩阵模板类。结果遇到不少坑，毕竟平时没怎么使用动态分配内存，只好先用静态数组水过，搞完之后调试很久才写出矩阵的动态版本，主要是开始写模板矩阵的时候没有创建复制构造和重载=，导致类的浅复制，出现很多bug（值很随机）。最后查了好多遍才想到是那里的问题，不过搞完之后对这方面了解更多点了，嘿嘿！

题目：

1.稀疏矩阵乘法

?题意：两个n*n的矩阵相乘(n<=800)，结果对3取模（hdu4920）
?思路：先对3取模，所以两个矩阵里面会出现很多0，所以可以先枚举一个矩阵，只有当该位置不是0的时候才和另一个矩阵做乘法,我还用了一下输入挂。
#include <iostream>
#include <cstdio>
#include <cstdlib>
typedef long long ll;
#define foreach(it,v) for(__typeof((v).begin()) it = (v).begin(); it != (v).end(); ++it)
using namespace std;
int mod = 3;
const int buf_len = 4096;
char buf[buf_len], *bufb(buf), *bufe(buf + 1);
int E = 1;
#define readBuf() {     if (++ bufb == bufe)         bufe = (bufb = buf) + (E=fread(buf, 1, sizeof(buf), stdin));     if(!E)return 0; }
#define readInt(_y_) {     register int _s_(0);     do {         readBuf();     } while (!isdigit(*bufb));     do {         _s_ = (_s_<<1) + (_s_<<3) + *bufb - '0';         readBuf();     } while (isdigit(*bufb));     _y_ = _s_; }
template <typename T> struct Matrix
{
	T ** base;
	int row,colnum;
	Matrix(int n = 0,int m = 0):row(n),colnum(m) {
		base = new T * [n];
		for(int i = 0; i < n; i++)
			base[i] = new T [m];
	}
	Matrix (const Matrix<T> & A) {
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	Matrix operator + (const Matrix <T> &A)const{
		Matrix<T> res(row,colnum);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++) {
				res.base[i][j] = base[i][j] + A.base[i][j];
				if(mod)res.base[i][j] %= mod;
			}
		}
		return res;
	}
	void operator = (const Matrix<T> & A) {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	T * operator [] (const int & i) {
		return base[i];
	}
	void setv(const T & val) {
		for(int i = 0; i < row; i++)
			for(int j = 0; j < colnum; j++)
				base[i][j] = val;
	}
	Matrix ones(int n) const{
		Matrix<T> res(n,n);
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++)
				res.base[i][j] = (T)(i==j);
		return res;
	}
	Matrix operator * (const Matrix<T> & rhs) const {
		if(colnum != rhs.row) {
			cerr<<"worng size of two Matrix"<<endl;
			exit(-1);
		}
		Matrix<T> c(row,rhs.colnum);
		c.setv(T(0));
		for(int k = 0; k < colnum; k++) {
			for(int i = 0; i < row; i++) {
				T r = base[i][k];
				if(!r)continue;
				for(int j = 0; j < c.colnum; j++) {
					c[i][j] += r*rhs.base[k][j];
					if(mod)c[i][j] %= mod;
				}
			}
		}
		return c;
	}
	Matrix exp(int n) const{
		if(row!=colnum) {
			cerr<<"can't exp on different row and colnum"<<endl;
			exit(-1);
		}
		Matrix<T>res = ones(row),b(*this);
		while(n > 0) {
			if(n & 1)res = res * b;
			b = (b * b);
			n >>= 1;
		}
		return res;
	}
	void debug()const{
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++)
				cout<<base[i][j]<<" \n"[j+1==colnum];
		}
	}
	~Matrix() {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
	}
};
int main(int argc, char const *argv[])
{
	while(1) {
		int n;readInt(n);
		Matrix<int> A(n,n),B(n,n);
		for(int i = 0; i < n; i++) {
			for(int j = 0; j < n; j++) {
				readInt(A[i][j]);
				A[i][j] %= 3;
			}
		}
		for(int i = 0; i < n; i++) {
			for(int j = 0; j < n; j++) {
				readInt(B[i][j]);
				B[i][j] %= 3;
			}
		}
		(A*B).debug();
	}
	return 0;
}

?uva10655 给三个非负整数p,q,n，求a^n+b^n，其中a+b=p,ab=q(保证结果在10^18以内)
?解决思路：从a^n+b^n入手构造递推关系而且只用p,q表示，手算几项之后会发现f[n] = f[n-1]*p - f[n-2]*q,构造下矩阵然后快速幂。
#include <bits/stdc++.h>
typedef long long ll;
#define foreach(it,v) for(__typeof((v).begin()) it = (v).begin(); it != (v).end(); ++it)
using namespace std;
const int maxn = 2e5 + 10;
ll mod = 0;
template <typename T> struct Matrix
{
	T ** base;
	int row,colnum;
	Matrix(int n = 0,int m = 0):row(n),colnum(m) {
		base = new T * [n];
		for(int i = 0; i < n; i++)
			base[i] = new T [m];
	}
	Matrix (const Matrix<T> & A) {
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	void operator = (const Matrix<T> & A) {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	T * operator [] (const int & i) {
		return base[i];
	}
	void setv(const T & val) {
		for(int i = 0; i < row; i++)
			for(int j = 0; j < colnum; j++)
				base[i][j] = val;
	}
	Matrix operator * (const Matrix<T> & rhs) const {
		if(colnum != rhs.row) {
			cerr<<"worng size of two Matrix"<<endl;
			exit(-1);
		}
		Matrix<T> c(row,rhs.colnum);
		c.setv(T(0));
		for(int k = 0; k < colnum; k++) {
			for(int i = 0; i < row; i++) {
				T r = base[i][k];
				for(int j = 0; j < c.colnum; j++) {
					c[i][j] += r*rhs.base[k][j];
					if(mod)c[i][j] %= mod;
				}
			}
		}
		return c;
	}
	Matrix exp(int n) {
		if(row!=colnum) {
			cerr<<"can't exp on different row and colnum"<<endl;
			exit(-1);
		}
		Matrix<T>res(row,colnum),b(*this);
		res.setv(T(0));
		for(int i = 0; i < row; i++) res.base[i][i] = T(1);
		while(n > 0) {
			if(n & 1)res = res * b;
			b = (b * b);
			n >>= 1;
		}
		return res;
	}
	void debug() {
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++)
				cout<<base[i][j]<<" \n"[j+1==colnum];
		}
	}
	~Matrix() {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
	}
};
int main(int argc, char const *argv[])
{
	int p,q,n;
	while(scanf("%d%d%d",&p,&q,&n)==3) {
		if(n==0)puts("2");
		else {
			Matrix<ll> A(1,2);
			A[0][0] = p;A[0][1] = 2;
			Matrix<ll> B(2,2);
			B[0][0] = p;B[0][1] = 1;
			B[1][0] = -q; B[1][1] = 0;
			printf("%lld\n", (A*B.exp(n-1))[0][0]);
		}
	}
	return 0;
}

?poj3233 A是一个n*n矩阵，求S = （A + A^2 + A^3 +… + A^k）%m,n<30,k<10^9
思路：矩阵快速幂加二分。如果我们求出了x = A+A^2+...A^(k/2),那么y = A^(k/2+1)+A^(k/2+2)+...A^k，可以由x的乘A^
(k/2) + {A^k|(当k为奇数时才需要}；
#include <iostream>
#include <cstdio>
#include <cstdlib>
typedef long long ll;
#define foreach(it,v) for(__typeof((v).begin()) it = (v).begin(); it != (v).end(); ++it)
using namespace std;
const int maxn = 2e5 + 10;
int mod = 0;
template <typename T> struct Matrix
{
	T ** base;
	int row,colnum;
	Matrix(int n = 0,int m = 0):row(n),colnum(m) {
		base = new T * [n];
		for(int i = 0; i < n; i++)
			base[i] = new T [m];
	}
	Matrix (const Matrix<T> & A) {
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	Matrix operator + (const Matrix <T> &A)const{
		Matrix<T> res(row,colnum);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++) {
				res.base[i][j] = base[i][j] + A.base[i][j];
				if(mod)res.base[i][j] %= mod;
			}
		}
		return res;
	}
	void operator = (const Matrix<T> & A) {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	T * operator [] (const int & i) {
		return base[i];
	}
	void setv(const T & val) {
		for(int i = 0; i < row; i++)
			for(int j = 0; j < colnum; j++)
				base[i][j] = val;
	}
	Matrix ones(int n) const{
		Matrix<T> res(n,n);
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++)
				res.base[i][j] = (T)(i==j);
		return res;
	}
	Matrix operator * (const Matrix<T> & rhs) const {
		if(colnum != rhs.row) {
			cerr<<"worng size of two Matrix"<<endl;
			exit(-1);
		}
		Matrix<T> c(row,rhs.colnum);
		c.setv(T(0));
		for(int k = 0; k < colnum; k++) {
			for(int i = 0; i < row; i++) {
				T r = base[i][k];
				if(!r)continue;
				for(int j = 0; j < c.colnum; j++) {
					c[i][j] += r*rhs.base[k][j];
					if(mod)c[i][j] %= mod;
				}
			}
		}
		return c;
	}
	Matrix exp(int n) const{
		if(row!=colnum) {
			cerr<<"can't exp on different row and colnum"<<endl;
			exit(-1);
		}
		Matrix<T>res = ones(row),b(*this);
		while(n > 0) {
			if(n & 1)res = res * b;
			b = (b * b);
			n >>= 1;
		}
		return res;
	}
	void debug()const{
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++)
				cout<<base[i][j]<<" \n"[j+1==colnum];
		}
	}
	~Matrix() {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
	}
};
template<typename T>Matrix<T> gao(const Matrix<T> &A,int k) {
	if(k==0)return A.ones(A.row);
	if(k==1)return A;
	Matrix<T> c = gao(A,k>>1);
	Matrix<T> t = A.exp(k>>1);
	c = c + c*t;
	if(k&1)c = c + t*t*A;
	return c;
}
int main(int argc, char const *argv[])
{
	int n,k;
	while(scanf("%d%d%d",&n,&k,&mod)==3) {
		Matrix<int>A(n,n);
		for(int i = 0; i < A.row; i++)
			for(int j = 0; j < A.colnum; j++)
				scanf("%d",&A[i][j]);
		(gao(A,k)).debug();
	}
	return 0;
}
?hdu1588 f(0)=0  f(1)=1  f(n)=f(n-1)+f(n-2) (n>=2),g(i)=k*i+b ，给定k,b,n,求sigma(f(g(i))(0<=i<=n)
?(k,b,n<=1,000,000,000,有取模）
因为g是一个等差数列，所以f是个f(g)是个等比矩阵。和上题思路一样，对公比的幂求和之后乘上初始值
#include <iostream>
#include <cstdio>
#include <cstdlib>
typedef long long ll;
#define foreach(it,v) for(__typeof((v).begin()) it = (v).begin(); it != (v).end(); ++it)
using namespace std;
const int maxn = 2e5 + 10;
ll mod = 0;
template <typename T> struct Matrix
{
	T ** base;
	int row,colnum;
	Matrix(int n = 0,int m = 0):row(n),colnum(m) {
		base = new T * [n];
		for(int i = 0; i < n; i++)
			base[i] = new T [m];
	}
	Matrix (const Matrix<T> & A) {
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	Matrix operator + (const Matrix <T> &A)const{
		Matrix<T> res(row,colnum);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++) {
				res.base[i][j] = base[i][j] + A.base[i][j];
				if(mod)res.base[i][j] %= mod;
			}
		}
		return res;
	}
	void operator = (const Matrix<T> & A) {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	T * operator [] (const int & i) {
		return base[i];
	}
	void setv(const T & val) {
		for(int i = 0; i < row; i++)
			for(int j = 0; j < colnum; j++)
				base[i][j] = val;
	}
	Matrix ones(int n) const{
		Matrix<T> res(n,n);
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++)
				res.base[i][j] = (T)(i==j);
		return res;
	}
	Matrix operator * (const Matrix<T> & rhs) const {
		if(colnum != rhs.row) {
			cerr<<"worng size of two Matrix"<<endl;
			exit(-1);
		}
		Matrix<T> c(row,rhs.colnum);
		c.setv((T)0);
		for(int k = 0; k < colnum; k++) {
			for(int i = 0; i < row; i++) {
				T r = base[i][k];
				if(!r)continue;
				for(int j = 0; j < c.colnum; j++) {
					c[i][j] += r*rhs.base[k][j];
					if(mod)c[i][j] %= mod;
				}
			}
		}
		return c;
	}
	Matrix exp(int n) const{
		if(row!=colnum) {
			cerr<<"can't exp on different row and colnum"<<endl;
			exit(-1);
		}
		Matrix<T>res = ones(row),b(*this);
		while(n > 0) {
			if(n & 1)res = res * b;
			b = (b * b);
			n >>= 1;
		}
		return res;
	}
	void debug()const{
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++)
				cout<<base[i][j]<<" \n"[j+1==colnum];
		}
	}
	~Matrix() {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
	}
};
template<typename T>Matrix<T> gao(const Matrix<T> &A,int k) {
	if(k==0)return A.ones(A.row);
	if(k==1)return A;
	Matrix<T> c = gao(A,k>>1);
	Matrix<T> t = A.exp(k>>1);
	c = c + c*t;
	if(k&1)c = c + t*t*A;
	return c;
}
int main(int argc, char const *argv[])
{
	int n,k,b;
	Matrix<ll> a0(1,2);
	a0[0][0] = 1; a0[0][1] = 0;
	Matrix<ll> q0(2,2); q0.setv(1);
	q0[1][1] = 0;
	while(scanf("%d%d%d%I64d",&k,&b,&n,&mod)==4) {
		Matrix<ll> A = a0*(q0.exp(b));
		Matrix<ll> Q = q0.exp(k);
		if(n==1) printf("%I64d\n", A[0][1]);
		else {
			printf("%I64d\n",(A*(Q.ones(Q.row) + gao(Q,n-1)))[0][1]);
		}
	}
	return 0;
}

?题意：给出n，求 技术分享 (n<=10^9)
?
?F(k)是标准fibonacci数列，C(n,k)是组合数
?结果对1e9+7取模
打表之后发现结果是F(2n)。证明如下
技术分享
#include <iostream>
#include <cstdio>
#include <cstdlib>
typedef long long ll;
#define foreach(it,v) for(__typeof((v).begin()) it = (v).begin(); it != (v).end(); ++it)
using namespace std;
const int maxn = 2e5 + 10;
ll mod = 0;
template <typename T> struct Matrix
{
	T ** base;
	int row,colnum;
	Matrix(int n = 0,int m = 0):row(n),colnum(m) {
		base = new T * [n];
		for(int i = 0; i < n; i++)
			base[i] = new T [m];
	}
	Matrix (const Matrix<T> & A) {
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	Matrix operator + (const Matrix <T> &A)const{
		Matrix<T> res(row,colnum);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++) {
				res.base[i][j] = base[i][j] + A.base[i][j];
				if(mod)res.base[i][j] %= mod;
			}
		}
		return res;
	}
	void operator = (const Matrix<T> & A) {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	T * operator [] (const int & i) {
		return base[i];
	}
	void setv(const T & val) {
		for(int i = 0; i < row; i++)
			for(int j = 0; j < colnum; j++)
				base[i][j] = val;
	}
	Matrix ones(int n) const{
		Matrix<T> res(n,n);
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++)
				res.base[i][j] = (T)(i==j);
		return res;
	}
	Matrix operator * (const Matrix<T> & rhs) const {
		if(colnum != rhs.row) {
			cerr<<"worng size of two Matrix"<<endl;
			exit(-1);
		}
		Matrix<T> c(row,rhs.colnum);
		c.setv((T)0);
		for(int k = 0; k < colnum; k++) {
			for(int i = 0; i < row; i++) {
				T r = base[i][k];
				if(!r)continue;
				for(int j = 0; j < c.colnum; j++) {
					c[i][j] += r*rhs.base[k][j];
					if(mod)c[i][j] %= mod;
				}
			}
		}
		return c;
	}
	Matrix exp(int n) const{
		if(row!=colnum) {
			cerr<<"can't exp on different row and colnum"<<endl;
			exit(-1);
		}
		Matrix<T>res = ones(row),b(*this);
		while(n > 0) {
			if(n & 1)res = res * b;
			b = (b * b);
			n >>= 1;
		}
		return res;
	}
	void debug()const{
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++)
				cout<<base[i][j]<<" \n"[j+1==colnum];
		}
	}
	~Matrix() {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
	}
};
template<typename T>Matrix<T> gao(const Matrix<T> &A,int k) {
	if(k==0)return A.ones(A.row);
	if(k==1)return A;
	Matrix<T> c = gao(A,k>>1);
	Matrix<T> t = A.exp(k>>1);
	c = c + c*t;
	if(k&1)c = c + t*t*A;
	return c;
}
int main(int argc, char const *argv[])
{
	int n,k,b;
	Matrix<ll> a0(1,2);
	a0[0][0] = 1; a0[0][1] = 0;
	Matrix<ll> q0(2,2); q0.setv(1);
	q0[1][1] = 0;
	while(scanf("%d%d%d%I64d",&k,&b,&n,&mod)==4) {
		Matrix<ll> A = a0*(q0.exp(b));
		Matrix<ll> Q = q0.exp(k);
		if(n==1) printf("%I64d\n", A[0][1]);
		else {
			printf("%I64d\n",(A*(Q.ones(Q.row) + gao(Q,n-1)))[0][1]);
		}
	}
	return 0;
}

hdu 2157给出一个n个点，m条边的有向图，T次询问，每次输入a，b，k问从a到b经过k个点的方案是多少，直接快速幂。输出

#include <iostream>
#include <cstdio>
#include <cstdlib>
typedef long long ll;
#define foreach(it,v) for(__typeof((v).begin()) it = (v).begin(); it != (v).end(); ++it)
using namespace std;
const int maxn = 2e5 + 10;
ll mod = 0;
template <typename T> struct Matrix
{
	T ** base;
	int row,colnum;
	Matrix(int n = 0,int m = 0):row(n),colnum(m) {
		base = new T * [n];
		for(int i = 0; i < n; i++)
			base[i] = new T [m];
	}
	Matrix (const Matrix<T> & A) {
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	Matrix operator + (const Matrix <T> &A)const{
		Matrix<T> res(row,colnum);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++) {
				res.base[i][j] = base[i][j] + A.base[i][j];
				if(mod)res.base[i][j] %= mod;
			}
		}
		return res;
	}
	void operator = (const Matrix<T> & A) {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
		row = A.row;
		colnum = A.colnum;
		base = new T *[row];
		for(int i = 0; i < row; i++) {
			base[i] = new T [colnum];
			for(int j = 0; j < colnum; j++)
				base[i][j] = A.base[i][j];
		}
	}
	T * operator [] (const int & i) {
		return base[i];
	}
	void setv(const T & val) {
		for(int i = 0; i < row; i++)
			for(int j = 0; j < colnum; j++)
				base[i][j] = val;
	}
	Matrix ones(int n) const{
		Matrix<T> res(n,n);
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++)
				res.base[i][j] = (T)(i==j);
		return res;
	}
	Matrix operator * (const Matrix<T> & rhs) const {
		if(colnum != rhs.row) {
			cerr<<"worng size of two Matrix"<<endl;
			exit(-1);
		}
		Matrix<T> c(row,rhs.colnum);
		c.setv((T)0);
		for(int k = 0; k < colnum; k++) {
			for(int i = 0; i < row; i++) {
				T r = base[i][k];
				if(!r)continue;
				for(int j = 0; j < c.colnum; j++) {
					c[i][j] += r*rhs.base[k][j];
					if(mod)c[i][j] %= mod;
				}
			}
		}
		return c;
	}
	Matrix exp(int n) const{
		if(row!=colnum) {
			cerr<<"can't exp on different row and colnum"<<endl;
			exit(-1);
		}
		Matrix<T>res = ones(row),b(*this);
		while(n > 0) {
			if(n & 1)res = res * b;
			b = (b * b);
			n >>= 1;
		}
		return res;
	}
	void debug()const{
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < colnum; j++)
				cout<<base[i][j]<<" \n"[j+1==colnum];
		}
	}
	~Matrix() {
		for(int i = 0; i < row; i++) delete [] base[i];
		delete [] base;
	}
};
int main(int argc, char const *argv[])
{
	int n,m;
	mod = 1000;
	while(scanf("%d%d",&n,&m)==2) {
		if(n+m==0)return 0;
		Matrix<int> A(n,n);
		A.setv(0);
		while(m--) {
			int u,v;scanf("%d%d",&u,&v);
			A[u][v] = 1;
		}
		int T;scanf("%d",&T);
		while(T--) {
			int s,t,k; scanf("%d%d%d",&s,&t,&k);
			printf("%d\n",A.exp(k)[s][t]);
		}
	}
	return 0;
}

FZU 2173 给出一个n个点，H条边的有向图，还有输入的k，经过一条边需要1s，有一个花费w，问从1到n恰好k秒的最小花费是多少，其实就是跑k遍floyd，可以用矩阵快速幂加速一下，乘法换加法就ok，这题没用模板写
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
typedef long long ll;
const ll inf = 1e17 + 10;
#define Type ll
#define SIZER 50
#define SIZEC 50
#define foreach(it,v) for(__typeof((v).begin()) it = (v).begin(); it != (v).end(); ++it)
using namespace std;
int mod = 0;
const int maxn = 5000000 + 4;
struct Matrix
{
	Type a[SIZER][SIZEC];
	int row,colunm;
	Matrix(int n = 1,int m = 1):row(n),colunm(m){}
	bool samerc(const Matrix & rhs) const{
		return rhs.row == row && rhs.colunm == colunm;
	}
	bool canmul(const Matrix & rhs) const{
		return colunm == rhs.row;
	}
	Matrix operator + (const Matrix & rhs) const{
		if(!samerc(rhs)) {
			puts("size Dont match");
			exit(-1);
		}
		Matrix c(row,colunm);
		for(int i = 0; i < row; i++)
			for(int j = 0; j < colunm; j++) {
				c.a[i][j] = a[i][j] + rhs.a[i][j];
				if(mod) c.a[i][j] %= mod;
			}
		return c;
	}
	void one() {
		for(int i = 0; i < row; i++)
			for(int j = 0; j < colunm; j++)
				a[i][j] = (i==j);
	}
	void debug() {
		for(int i = 0; i < row; i++)
			for(int j = 0; j < colunm; j++)
				cout<<a[i][j]<<" \n"[j+1==colunm];
	}
	void Zero() { memset(a,0,sizeof a); }
	Matrix operator * (const Matrix & rhs) const{
		if(!canmul(rhs)) {
			puts("can‘t mul Matrix with wrong size");
			exit(-1);
		}
		Matrix c(row,rhs.colunm);
		for(int i = 0; i < c.row; i++) {
			for(int j = 0; j < c.colunm; j++) {
				c.a[i][j] = inf;
			}
		}
		for(int i = 0; i < row; i++)
			for(int j = 0; j < rhs.colunm; j++)
				for(int k = 0; k < colunm; k++)
					c.a[i][j] = min(c.a[i][j],a[i][k] + rhs.a[k][j]);
		return c;
	}
};
Matrix quick_exp(const Matrix & A,int k)
{
	Matrix res = A,b = A;
	k--;
	while(k > 0) {
		if(k&1)res = res * b;
		k >>= 1;
		b = b*b;
	}
	return res;
}
int main()
{
	int T;scanf("%d",&T);
	while(T--) {
		ll n,h,k;scanf("%I64d%I64d%I64d",&n,&h,&k);
		Matrix A(n,n);
		for(int i = 0; i < n; i++) {
			for(int j = 0; j < n; j++) {
				A.a[i][j] = inf;
			}
		}
		while(h--) {
			ll u,v,w;scanf("%I64d%I64d%I64d",&u,&v,&w);
			u--,v--;
			A.a[u][v] = min(A.a[u][v],w);
		}
		ll ans = quick_exp(A,k).a[0][n-1];
		if(ans>= inf) ans = -1;
		printf("%I64d\n",ans);
	}
	return 0;
}

矩阵小结

标签：矩阵   数学   

原文地址：http://blog.csdn.net/acvcla/article/details/46064845
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！