首页 > 其他 > 详细

矩阵小结

时间：2015-05-28 01:59:23      阅读：310      评论：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年09月23日 (328)
2021年09月24日 (313)
2021年09月17日 (191)
2021年09月15日 (369)
2021年09月16日 (411)
2021年09月13日 (439)
2021年09月11日 (398)
2021年09月12日 (393)
2021年09月10日 (160)
2021年09月08日 (222)
最新文章
更多>
教程昨日排行
更多>
友情链接
汇智网   PHP教程   插件网  
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣，分享程序人生！