行列演算 - yanprog @ ウィキ

行列

高速行列冪演算系の問題を解く場合、次を貼り付ける。

ll mod=1000000007;
struct mint{
	ll value;
	mint():value(0){}
	mint(ll v):value((v%mod+mod)%mod){}
};
mint& operator+=(mint&a, mint b){return a=a.value+b.value;}
mint& operator-=(mint&a, mint b){return a=a.value-b.value;}
mint& operator*=(mint&a, mint b){return a=a.value*b.value;}
mint operator+(mint a, mint b){return a+=b;}
mint operator-(mint a, mint b){return a-=b;}
mint operator*(mint a, mint b){return a*=b;}
mint operator-(mint a){return 0-a;}
/*bool operator==(mint a, mint b){return a.value==b.value;}
bool operator!=(mint a, mint b){return a.value!=b.value;}
 
 
std::ostream& operator<<(std::ostream& os, const mint& m){
return ( os << m.value );}
ll extgcd(ll a, ll b, ll &x, ll &y){
	ll d=a;
	if(b!=0){
		d=extgcd(b, a%b, y, x);
		y-=(a/b)*x;
	}
	else{
		x=1,y=0;
	}
	return d;
}
ll modinverse(ll a, ll b){
	ll x,y;
	extgcd(a,b, x, y);
	return (x%b+b)%b;
}
mint& operator/=(mint&a, mint b){return a=a.value*modinverse(b.value,mod);}
mint operator/(mint a, mint b){return a/=b;}
*/
 
 
template <typename T> struct matrix : vector<vector<T > >{
	matrix(int n){
		vector<vector<T> >::clear();
		resize(n, vector<T>(n) );
		REP(i,n) (*this)[i][i]=1;
	}
	matrix(int n, int m){
		vector<vector<T> >::clear();
		resize(n, vector<T>(m) ); 
	}
	matrix<T> operator*(const matrix<T> &b){
		int n1=this->size(), n2=b.size(), n3=b[0].size();
		assert(this->size()==n2);
		matrix<T> ret(n1,n3);
		REP(i,n1) REP(j,n3) REP(k,n2) ret[i][j]=ret[i][j]+(*this)[i][k]*b[k][j];
		return ret;
	}
	matrix<T> pow(ll n){
		int m=this->size();
		assert((*this)[0].size()==m);
		return pow2(matrix<T>(m) ,n);
	}
	matrix<T> pow2(matrix t, ll n){
		if(n==0)return t;
		if(n%2==0)return ((*this)*(*this)).pow2( t, n/2);
		else return ((*this)*(*this)).pow2 ((*this)*t, n/2);
	}
	/*T det(){
		int n=this->size();
		vector<vector<T> >x;
		REP(i,n) x.pb((*this)[i]);
		T h=1;
		REP(i,n){
			FOR(j,i,n){
				if(x[j][i]!=0 && j!=i){
					REP(k,n)swap(x[j][k], x[i][k]);
					h*=-1;
					break;
				}
			}
			if(x[i][i]==0)return 0;
			T v=1/x[i][i];
			h=h*x[i][i];
			FOR(j,i+1,n){
				T v2=v*x[j][i];
				REP(k,n) x[j][k]=x[j][k]-x[i][k]*v2;
			}
		}
		return h;
	}*/
};
 

使用例　SRM544 Div1HARD

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <map>
#include <utility>
#include <set>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <sstream>
#include <complex>
#include <stack>
#include <queue>
#include <cstring>
#include <assert.h>
 
#define FOR(i,a,b) for(int i=(a);i<(b);++i)
#define REP(i,n)  FOR(i,0,n)
#define EACH(i,c) for(typeof((c).begin()) i=(c).begin(); i!=(c).end(); ++i)
#define EXIST(s,e) ((s).find(e)!=(s).end())
#define ALL(s) (s).begin(),(s).end
 
#define clr(a) memset((a),0,sizeof(a))
#define nclr(a) memset((a),-1,sizeof(a))
#define pb push_back
#define INRANGE(x,s,e) ((s)<=(x) && (x)<(e))
#define MP(a,b) make_pair((a),(b))
 
using namespace std;
 
// BEGIN CUT HERE
#define debug(x) cerr << #x << " = " << (x) << " (L" << __LINE__ << ")" << " " << __FILE__ << endl;
#define deb(x) cerr << #x << " = " << (x) << " , ";
#define debl cerr << " (L" << __LINE__ << ")"<< endl;
template<typename T> std::ostream& operator<<(std::ostream& os, const vector<T>& z){
os << "[ ";REP(i,z.size())os << z[i] << ", " ;return ( os << "]" << endl);}
template<typename T> std::ostream& operator<<(std::ostream& os, const set<T>& z){
os << "set( "; EACH(p,z)os << (*p) << ", " ;return ( os << ")" << endl);}
template<typename T,typename U> std::ostream& operator<<(std::ostream& os, const map<T,U>& z){
os << "{ "; EACH(p,z)os << (p->first) << ": " << (p->second) << ", " ;return ( os << "}" << endl);}
template<typename T,typename U> std::ostream& operator<<(std::ostream& os, const pair<T,U>& z){
return ( os << "(" << z.first << ", " << z.second << ",)" );}
// END CUT HERE
 
typedef long long ll;
typedef pair<int,int> pii;
typedef pair<long long, long long> pll;
typedef vector<int> vi;
typedef vector<vector<int> > vvi;
typedef vector<long long> vl;
typedef vector<vector<long long> > vvl;
typedef vector<double> vd;
typedef vector<vector<double> > vvd;
typedef vector<string> vs;
 
 
ll mod=1000000007;
struct mint{
	ll value;
	mint():value(0){}
	mint(ll v):value((v%mod+mod)%mod){}
};
mint& operator+=(mint&a, mint b){return a=a.value+b.value;}
mint& operator-=(mint&a, mint b){return a=a.value-b.value;}
mint& operator*=(mint&a, mint b){return a=a.value*b.value;}
mint operator+(mint a, mint b){return a+=b;}
mint operator-(mint a, mint b){return a-=b;}
mint operator*(mint a, mint b){return a*=b;}
mint operator-(mint a){return 0-a;}
//bool operator==(mint a, mint b){return a.value==b.value;}
//bool operator!=(mint a, mint b){return a.value!=b.value;}
 
/*
std::ostream& operator<<(std::ostream& os, const mint& m){
return ( os << m.value );}
ll extgcd(ll a, ll b, ll &x, ll &y){
	ll d=a;
	if(b!=0){
		d=extgcd(b, a%b, y, x);
		y-=(a/b)*x;
	}
	else{
		x=1,y=0;
	}
	return d;
}
ll modinverse(ll a, ll b){
	ll x,y;
	extgcd(a,b, x, y);
	return (x%b+b)%b;
}
mint& operator/=(mint&a, mint b){return a=a.value*modinverse(b.value,mod);}
mint operator/(mint a, mint b){return a/=b;}
*/
 
 
template <typename T> struct matrix : vector<vector<T > >{
	matrix(int n){
		vector<vector<T> >::clear();
		resize(n, vector<T>(n) );
		REP(i,n) (*this)[i][i]=1;
	}
	matrix(int n, int m){
		vector<vector<T> >::clear();
		resize(n, vector<T>(m) ); 
	}
	matrix<T> operator*(const matrix<T> &b){
		int n1=this->size(), n2=b.size(), n3=b[0].size();
		assert(this->size()==n2);
		matrix<T> ret(n1,n3);
		REP(i,n1) REP(j,n3) REP(k,n2) ret[i][j]=ret[i][j]+(*this)[i][k]*b[k][j];
		return ret;
	}
	matrix<T> pow(ll n){
		int m=this->size();
		assert((*this)[0].size()==m);
		return pow2(matrix<T>(m) ,n);
	}
	matrix<T> pow2(matrix t, ll n){
		if(n==0)return t;
		if(n%2==0)return ((*this)*(*this)).pow2( t, n/2);
		else return ((*this)*(*this)).pow2 ((*this)*t, n/2);
	}
	/*T det(){
		int n=this->size();
		vector<vector<T> >x;
		REP(i,n) x.pb((*this)[i]);
		T h=1;
		REP(i,n){
			FOR(j,i,n){
				if(x[j][i]!=0 && j!=i){
					REP(k,n)swap(x[j][k], x[i][k]);
					h*=-1;
					break;
				}
			}
			if(x[i][i]==0)return 0;
			T v=1/x[i][i];
			h=h*x[i][i];
			FOR(j,i+1,n){
				T v2=v*x[j][i];
				REP(k,n) x[j][k]=x[j][k]-x[i][k]*v2;
			}
		}
		return h;
	}*/
};
 
class SplittingFoxes {
public:
	int sum(long long n, ll _S, ll _L, ll _R) {
		mint S=_S, L=_L, R=_R;
		matrix<mint> a(4,4);
		enum{N,X,Y,XY};
		a[N][N]=S-L-R;
 
		a[XY][XY]=S+L+R;
		a[XY][Y]=S;
 
		a[X][N]=S;
		a[X][X]=S;
		a[X][Y]=-L+R;
 
		a[Y][Y]=S;
		a[Y][X]=L-R;
		//debug(a);
		//debug(a.pow(n));
		int ans=a.pow(n)[XY][N].value;
 
		return ans;
	}
};