yanprog @ ウィキ
yanprog @ ウィキ

線形漸化式一般項推定エンジン

行列

線形漸化式で定義される数列の一般項を確定する。数列の最初のほうの項をaとしたとき
est.build(a)でaを与え、est.get(n)でn項目を取得する。

TODO:Codeforces 230 Div1C 

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;
	ll d=extgcd(a,b, x, y);
	assert(d==1);
	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;
	}*/
};
 
struct EST{
	EST(){};
	vector<mint> q_init;	
	vector<mint> q;
	void build(vector<mint> given){
		q_init=given;
		int n=given.size()/2;
		for(int k=0; k<=n; k++){
			vector<mint> a(given.begin(), given.begin()+2*k);
			bool b=build_partial(a, q);
			//deb(k);deb(b);debl;
		}
		//debug(q);
		//debug(q.size());
	}
	bool build_partial(vector<mint> a, vector<mint> &ret){
		assert(a.size()%2==0);
		int n=a.size()/2;
		vector<vector<mint> > m(n, vector<mint>(n+1));		
		REP(i,n) REP(j,n+1) m[i][j]=a[i+j];
		REP(y,n){
			FOR(x,y,n) if(m[x][y]!=0){
				REP(y2,n+1) swap(m[y][y2], m[x][y2]);
				break;				
			}
			if(m[y][y]==0)return false;
			{
				mint z=m[y][y];
				REP(y2,n+1) m[y][y2]/=z; 
			}
			assert(m[y][y]==1);
			REP(x,n){
				if(x==y)continue;
				mint z=m[x][y];
				REP(y2,n+1)m[x][y2]-=z*m[y][y2];
			} 
		}
		//debug(m);
		ret.resize(n);
		REP(x,n) ret[x]=m[x][n];
		return true;
	}
	mint get(ll n){
		int k=q.size();
		matrix<mint> m(k,k);
		REP(x,k) REP(y,k) if(y==x+1) m[x][y]=1;
		REP(y,k) m[k-1][y]=q[y];
		m=m.pow(n);
		mint ans=0;
		REP(y,k) ans+=m[0][y]*q_init[y];
		return ans;
	}
};
 

使用例　SRM521 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;
	ll d=extgcd(a,b, x, y);
	assert(d==1);
	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;
	}*/
};
 
struct EST{
	EST(){};
	vector<mint> q_init;	
	vector<mint> q;
	void build(vector<mint> given){
		q_init=given;
		int n=given.size()/2;
		for(int k=0; k<=n; k++){
			vector<mint> a(given.begin(), given.begin()+2*k);
			bool b=build_partial(a, q);
			//deb(k);deb(b);debl;
		}
		//debug(q);
		//debug(q.size());
	}
	bool build_partial(vector<mint> a, vector<mint> &ret){
		assert(a.size()%2==0);
		int n=a.size()/2;
		vector<vector<mint> > m(n, vector<mint>(n+1));		
		REP(i,n) REP(j,n+1) m[i][j]=a[i+j];
		REP(y,n){
			FOR(x,y,n) if(m[x][y]!=0){
				REP(y2,n+1) swap(m[y][y2], m[x][y2]);
				break;				
			}
			if(m[y][y]==0)return false;
			{
				mint z=m[y][y];
				REP(y2,n+1) m[y][y2]/=z; 
			}
			assert(m[y][y]==1);
			REP(x,n){
				if(x==y)continue;
				mint z=m[x][y];
				REP(y2,n+1)m[x][y2]-=z*m[y][y2];
			} 
		}
		//debug(m);
		ret.resize(n);
		REP(x,n) ret[x]=m[x][n];
		return true;
	}
	mint get(ll n){
		int k=q.size();
		matrix<mint> m(k,k);
		REP(x,k) REP(y,k) if(y==x+1) m[x][y]=1;
		REP(y,k) m[k-1][y]=q[y];
		m=m.pow(n);
		mint ans=0;
		REP(y,k) ans+=m[0][y]*q_init[y];
		return ans;
	}
};
 
map<vi, mint> memo;
mint rec(vi z){
 
	if(EXIST(memo,z)) return memo[z];
	mint&ret=memo[z];
	int tot=0;
	REP(i,z.size()) tot+=z[i];
	if(tot==0) return ret=1;
	mint ans=0;
	REP(i1,8){
		int i2=(i1+1)%8;
		if(z[i1]==z[i2] && z[i1]>=1 && z[i1]%2==i1%2){
			vi nz=z;
			nz[i1]--;
			nz[i2]--;
			ans+=rec(nz);
		}
	}
	return ret=ans;
}
 
class Chimney {
public:
	int countWays(long long n) {
		vector<mint> a;
		REP(i,10) a.pb(rec(vi(8,i)));
		//debug(a);
		EST est;
		est.build(a);
		int ans=est.get(n).value;
 
		return ans;
	}
 
 
};
 


使用例その２(SRM544Div1HAARD) 
#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;
	ll d=extgcd(a,b, x, y);
	assert(d==1);
	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;
	}*/
};
 
struct EST{
	EST(){};
	vector<mint> q_init;	
	vector<mint> q;
	void build(vector<mint> given){
		q_init=given;
		int n=given.size()/2;
		for(int k=0; k<=n; k++){
			vector<mint> a(given.begin(), given.begin()+2*k);
			bool b=build_partial(a, q);
			//deb(k);deb(b);debl;
		}
		//debug(q);
		//debug(q.size());
		if(q.size()==0){
			q.pb(0);
		}
	}
	bool build_partial(vector<mint> a, vector<mint> &ret){
		assert(a.size()%2==0);
		int n=a.size()/2;
		vector<vector<mint> > m(n, vector<mint>(n+1));		
		REP(i,n) REP(j,n+1) m[i][j]=a[i+j];
		REP(y,n){
			FOR(x,y,n) if(m[x][y]!=0){
				REP(y2,n+1) swap(m[y][y2], m[x][y2]);
				break;				
			}
			if(m[y][y]==0)return false;
			{
				mint z=m[y][y];
				REP(y2,n+1) m[y][y2]/=z; 
			}
			assert(m[y][y]==1);
			REP(x,n){
				if(x==y)continue;
				mint z=m[x][y];
				REP(y2,n+1)m[x][y2]-=z*m[y][y2];
			} 
		}
		//debug(m);
		ret.resize(n);
		REP(x,n) ret[x]=m[x][n];
		return true;
	}
	mint get(ll n){
		int k=q.size();
		matrix<mint> m(k,k);
		REP(x,k) REP(y,k) if(y==x+1) m[x][y]=1;
		REP(y,k) m[k-1][y]=q[y];
		m=m.pow(n);
		mint ans=0;
		REP(y,k) ans+=m[0][y]*q_init[y];
		return ans;
	}
};
 
int dx[4]={1,0,-1,0};
int dy[4]={0,1,0,-1};
class SplittingFoxes {
public:
	int sum(long long n, int S, int L, int R) {
		map<vi, mint> z;
		vi init(3,0);
		z[init]=1;
		vector<mint> a;
		REP(i,10){
			//debug(z);
			mint tot=0;
			EACH(pp,z){
				int x=pp->first[0];
				int y=pp->first[1];
				int r=pp->first[2];
				tot+=x*y*pp->second;		
			}
			a.pb(tot);
			map<vi, mint> nz;
			EACH(pp,z){
				int x=pp->first[0];
				int y=pp->first[1];
				int r=pp->first[2];
				mint v=pp->second;
				REP(k,3){
					int nx=x,ny=y,nr=r;
					if(k==0){
						nx+=dx[r];
						ny+=dy[r];
					}
					if(k==1)nr++;
					if(k==2)nr+=3;
					nr%=4;
					vi p(3);
					p[0]=nx;
					p[1]=ny;
					p[2]=nr;
					int weight[3]={S,L,R};
					nz[p]+=weight[k]*v;
				}
			}
			z=nz;
		}
		//debug(a);
		//debug(n);
		EST est;
		est.build(a);
 
		int ans=est.get(n).value;
 
		return ans;
	}	
};
 

タグ：

+ タグ編集
  • タグ：
タグの更新に失敗しました
エラーが発生しました。ページを更新してください。
ページを更新
「線形漸化式一般項推定エンジン」をウィキ内検索
最終更新：2014年02月19日 00:47