矩阵小结

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;
}
时间： 2024-10-07 22:17:44

矩阵小结的相关文章

稀疏矩阵的三元组顺序表存储及矩阵相乘算法小结

稀疏矩阵的三元组顺序表存储及矩阵相乘算法小结 巧若拙(欢迎转载,但请注明出处:http://blog.csdn.net/qiaoruozhuo) 一:稀疏矩阵的三元组顺序表数据结构 typedef int ElemType; typedef struct { intx, y;  //该非零元素的行下标和列下标 ElemTypee; //该非零元素的值 } Triple; typedef struct { Tripledata[MAXSIZE]; //非零元素三元组顺序表 intmu, nu, t

小结：矩阵乘法

概要: 在一些递推式中,我们发现好像不能在优化了(例如斐波那契数列普通递推是O(n)的),但是这个特殊的递推式我们可以用矩阵来实现O(logn)(忽略了矩阵自身计算的O(n^3)).而矩阵乘法运算是a[i, k]*b[k, j]=c[i, j],从这个式子可看出朴素是n^3的(当然那些神算法我不会),而fib数列递推式可以用矩阵乘法来表示,即 ┏ 1,1 ┓[f(n-1), f(n-2)]*┃        ┃=[f(n), f(n-1)] ┗ 1, 0 ┛ 而矩阵乘法是可以做快速幂的!也就是说

矩阵快速幂的一份小结

矩阵真是个好东西!虽然矩乘的复杂度有点难看... ... 这几天也做了不少矩阵题目,还是有几道好题目的.不过我打算从入门开始. 矩阵乘法:A[i][k]*B[k][j]=C[i][j];(A的第i行的每项依次乘以B的第j列的每项的和) 很显然这是一个n^3的算法,还是比较难看的. 代码就差不多是这样了. struct Matrix{int T[51][51];}; Matrix Mul(Matrix a,Matrix b,int I,int K,int J) { Matrix S=S0; for

矩阵快速幂小结

矩阵 并不想扯什么高端线代的内容因为我也不会 定义 由$n \times m$个数$a_{ij}$排成的$n$行$m$列的数表称为$n$行$m$列的矩阵,简称$n \times m$矩阵. $$A =\begin{bmatrix}a_{11} & a_{12} & \dots a_{1m} \\ a_{21}, & \dots & \dots \\ a_{31}, & \dots & \dots \\ a_{41} & \dots & a_{

矩阵乘法的Strassen算法详解

题目描述 请编程实现矩阵乘法,并考虑当矩阵规模较大时的优化方法. 思路分析 根据wikipedia上的介绍:两个矩阵的乘法仅当第一个矩阵B的列数和另一个矩阵A的行数相等时才能定义.如A是m×n矩阵和B是n×p矩阵,它们的乘积AB是一个m×p矩阵,它的一个元素其中 1 ≤ i ≤ m, 1 ≤ j ≤ p. 值得一提的是,矩阵乘法满足结合律和分配率,但并不满足交换律,如下图所示的这个例子,两个矩阵交换相乘后,结果变了: 下面咱们来具体解决这个矩阵相乘的问题. 解法一.暴力解法 其实,通过前面的分析

基于OpenMP的矩阵乘法实现及效率提升分析

一．  矩阵乘法串行实现 例子选择两个1024*1024的矩阵相乘,根据矩阵乘法运算得到运算结果.其中,两个矩阵中的数为double类型,初值由随机数函数产生.代码如下: #include <iostream> #include <omp.h> // OpenMP编程需要包含的头文件 #include <time.h> #include <stdlib.h> using namespace std; #define MatrixOrder 1024 #def

逻辑回归原理小结

逻辑回归是一个分类算法,它可以处理二元分类以及多元分类.虽然它名字里面有"回归"两个字,却不是一个回归算法.那为什么有"回归"这个误导性的词呢?个人认为,虽然逻辑回归是分类模型,但是它的原理里面却残留着回归模型的影子,本文对逻辑回归原理做一个总结. 1. 从线性回归到逻辑回归 我们知道,线性回归的模型是求出输出特征向量Y和输入样本矩阵X之间的线性关系系数\(\theta\),满足\(\mathbf{Y = X\theta}\).此时我们的Y是连续的,所以是回归模型.

SVD小结

1.矩阵分解 假设一个矩阵Data是m行n列,SVD(奇异值分解)将Data分解为U,E,VT 三个矩阵: Datam*n=Um*kEk*kVTk*n E是一个对角矩阵,对角元素为奇异值,对应Data的奇异值,即Data*DataT特征值的平方 2.选取特征 下面确定选取哪几维特征实现降维,去除噪声和冗余信息,用低维数据集表示原数据集. 典型做法是保留矩阵90%能量信息,公式如下,先选一个值h: 奇异阵的平方 sig=ETE 如果奇异阵的平方中前i项的和大于奇异阵的平方总和,即sum(sig[:

线性回归原理小结

线性回归可以说是机器学习中最基本的问题类型了,这里就对线性回归的原理和算法做一个小结. 1. 线性回归的模型函数和损失函数 线性回归遇到的问题一般是这样的.我们有m个样本,每个样本对应于n维特征和一个结果输出,如下: (x(0)1,x(0)2,...x(0)n,y0),(x(1)1,x(1)2,...x(1)n,y1),...(x(m)1,x(m)2,...x(m)n,yn)(x1(0),x2(0),...xn(0),y0),(x1(1),x2(1),...xn(1),y1),...(x1(m)