yukicoder : No.426 往復漸化式

本番うまく式変形できなかった。
http://yukicoder.me/problems/no/426

問題

3次正方行列の列A[N]、2次正方行列の列B[N]があり、これらは初期状態では単位行列である。
a[i],b[i]は以下のように定義される。

$a_{i+1} = A_i a_i$
$b_{i-1} = B_i b_i + \begin{pmatrix} 6i & 6i+1 & 6i+2 \\ 6i+3 & 6i+4 & 6i+5\end{pmatrix} a_i$

3要素のベクトルa[0]、2要素のベクトルb[N]が与えられる。
この状態で以下のクエリを順次処理せよ。

A[i]を与えられた値に更新する。
B[i]を与えられた値に更新する。
a[i]の値を答える。
b[i]の値を答える。

解法

＊以下、TeXの式を書くのが面倒だったので、公式解説のTeX表記を一部コピーペーストしています。

行列の積の形に持ち込み、SegTreeで答える。
まず、A,aのクエリを処理しよう。
こちらは簡単で、 $a_{j} = A_{j-1} A_{j-2} ... A_0 a_0$ に答えられるよう連続したA[i]の積を求められるようにしておくとよい。

問題はB,bの方である。
$S_i = \begin{pmatrix}6i& 6i + 1& 6i + 2 \\ 6i + 3& 6i + 4& 6i + 5\end{pmatrix}$ とおくと、 $b_{i-1} = B_i b_i + S_i a_i$ となる。
この漸化式を展開すると以下のようになる。
$b_i = B_{i+1} B_{i+2} \cdots B_{n} b_n + \\ (B_{i+1} B_{i+2} \cdots B_{n-1} S_n A_{n-1} A_{n-2} \cdots A_{i+1} + \\ B_{i+1} B_{i+2} \cdots B_{n-2} S_{n-1} A_{n-1} A_{n-2} \cdots A_{i+1} + \\ B_{i+1} S_{i+2} A_{i+1} +\\ S_{i+1})\\ a_{i+1}$

最初の項は、A同様に連続したB[i]の積を求められるようにすれば良い。また、a[i+1]はA[i]の連続した積で求められる。
問題は括弧内の部分である。A,B,Sそれぞれ、区間に対する以下の値を考える。

$A_{[i, j]}=A_{j}A_{j-1}...A_{i}$
$B_{[i, j]}=B_{i}B_{i+1}...B_{j}$
$S_{[i, j]}=S_i + B_i S_{i+1} A_i + B_i B_{i+1} S_{i+2} A_{i+1} A_i + \cdots + B_i B_{i+1} \cdots B_{j-1} S_{j} A_{j-1} \cdots A_{i+1} A_i$

するとそれぞれ以下の通り区間の結合ができるので、SegTreeで高速に区間に対する値を求めることができる。

$A_{[i, h]} = A_{[j+1, h]} A_{[i, j]}$
$B_{[i, h]} = B_{[i, j]} B_{[j+1, h]}$
$S_{[i, h]} = S_{[i, j]} + B_{[i, j]} S_{[j+1, h]} A_{[i, j]}$

あとはA,Bを更新するクエリに対してはSegTreeを更新しつつ、bを求めるクエリに対しては $b_i = B_{[i+1, n$ } b_n + S_{[i+1, n]} A_{[0, i]} a_0]を答えればよい。
以下のコードはBに関する処理も３次正方行列で行っている。

ll mo=1000000007;
const int MAT=3;
struct Mat { ll v[MAT][MAT]; };
Mat mulmat(Mat& a,Mat& b,int n=MAT) {
	int x,y,z; Mat r;
	FOR(x,n) FOR(y,n) r.v[x][y]=0;
	FOR(x,n) FOR(z,n) FOR(y,n) r.v[x][y] += (a.v[x][z]*b.v[z][y]) % mo;
	FOR(x,n) FOR(y,n) r.v[x][y]%=mo;
	return r;
}

Mat def={.v={{1,0,0},{0,1,0},{0,0,1}}};
Mat zero={.v={{0,0,0},{0,0,0},{0,0,0}}};

template<class V,int NV> class SegTree_A {
public:
	vector<V> val;
	virtual V comp(V l,V r){ return mulmat(r,l);};
	SegTree_A(){val=vector<V>(NV*2,def);};
	V getval(int x,int y,int l=0,int r=NV,int k=1) { // x<=i<y
		if(r<=x || y<=l) return def;
		if(x<=l && r<=y) return val[k];
		return comp(getval(x,y,l,(l+r)/2,k*2),getval(x,y,(l+r)/2,r,k*2+1));
	}
	void update(int entry, V v) {
		entry += NV;
		val[entry]=v;
		while(entry>1) entry>>=1, val[entry]=comp(val[entry*2],val[entry*2+1]);
	}
};
template<class V,int NV> class SegTree_B : public SegTree_A<V,NV>{
public:
	V comp(V l,V r){ return mulmat(l,r);};
};

int N,Q;
SegTree_A<Mat,1<<18> AS;
SegTree_B<Mat,1<<18> BS;

template<class V,int NV> class SegTree_S{
public:
	vector<V> val;
	SegTree_S(){val=vector<V>(NV*2,zero);};
	V getval(int x,int y,int l=0,int r=NV,int k=1) { // x<=i<y
		if(r<=x || y<=l) return zero;
		if(x<=l && r<=y) return val[k];
		
		if(y<=(l+r)/2) return getval(x,y,l,(l+r)/2,k*2);
		if(x>=(l+r)/2) return getval(x,y,(l+r)/2,r,k*2+1);
		
		Mat m1=getval(x,y,l,(l+r)/2,k*2);
		Mat m2=getval(x,y,(l+r)/2,r,k*2+1);
		Mat b=BS.getval(x,y,l,(l+r)/2,k*2);
		Mat a=AS.getval(x,y,l,(l+r)/2,k*2);
		
		m2 = mulmat(b,m2);
		m2 = mulmat(m2,a);
		int i,j;
		FOR(i,3) FOR(j,3) {
			m1.v[i][j]+=m2.v[i][j];
			if(m1.v[i][j]>=mo) m1.v[i][j]-=mo;
		}
		
		return m1;
	}
	void update(int x,int y,V v,int l=0,int r=NV,int k=1) { // x<=i<y
		if(r<=x || y<=l) return;
		if(x<=l && r<=y) {
			val[k]=v;
		}
		else {
			update(x,y,v,l,(l+r)/2,k*2);
			update(x,y,v,(l+r)/2,r,k*2+1);
			Mat m = mulmat(BS.val[k*2],val[k*2+1]);
			m = mulmat(m,AS.val[k*2]);
			int i,j;
			FOR(i,3) FOR(j,3) {
				val[k].v[i][j]=val[k*2].v[i][j]+m.v[i][j];
				if(val[k].v[i][j]>=mo) val[k].v[i][j]-=mo;
			}
		}
	}
};


SegTree_S<Mat,1<<18> SS;
Mat C[101010];
ll A[3],B[3];

void solve() {
	int i,j,k,l,r,x,y; string s;
	
	cin>>N;
	cin>>A[0]>>A[1]>>A[2];
	cin>>B[0]>>B[1];
	
	FOR(i,N+1) {
		C[i].v[0][0]=6*i+0;
		C[i].v[0][1]=6*i+1;
		C[i].v[0][2]=6*i+2;
		C[i].v[1][0]=6*i+3;
		C[i].v[1][1]=6*i+4;
		C[i].v[1][2]=6*i+5;
		SS.update(i,i+1,C[i]);
	}
	
	
	cin>>Q;
	while(Q--) {
		cin>>s>>i;
		if(s=="a") {
			Mat m;
			cin>>m.v[0][0]>>m.v[0][1]>>m.v[0][2];
			cin>>m.v[1][0]>>m.v[1][1]>>m.v[1][2];
			cin>>m.v[2][0]>>m.v[2][1]>>m.v[2][2];
			AS.update(i,m);
			SS.update(i,i+1,C[i]);
		}
		else if(s=="b") {
			Mat m;
			cin>>m.v[0][0]>>m.v[0][1];
			cin>>m.v[1][0]>>m.v[1][1];
			m.v[0][2]=m.v[1][2]=m.v[2][0]=m.v[2][1]=m.v[2][2]=0;
			BS.update(i,m);
			SS.update(i,i+1,C[i]);
		}
		else if(s=="ga") {
			Mat m=AS.getval(0,i);
			cout<<(m.v[0][0]*A[0]+m.v[0][1]*A[1]+m.v[0][2]*A[2])%mo<<" ";
			cout<<(m.v[1][0]*A[0]+m.v[1][1]*A[1]+m.v[1][2]*A[2])%mo<<" ";
			cout<<(m.v[2][0]*A[0]+m.v[2][1]*A[1]+m.v[2][2]*A[2])%mo<<endl;
		}
		else if(s=="gb") {
			Mat b=BS.getval(i+1,N+1);
			Mat s=SS.getval(i+1,N+1);
			Mat a=AS.getval(0,i+1);
			Mat sa=mulmat(s,a);
			cout<<((b.v[0][0]*B[0]+b.v[0][1]*B[1])%mo + (sa.v[0][0]*A[0]+sa.v[0][1]*A[1]+sa.v[0][2]*A[2])%mo)%mo<<" ";
			cout<<((b.v[1][0]*B[0]+b.v[1][1]*B[1])%mo + (sa.v[1][0]*A[0]+sa.v[1][1]*A[1]+sa.v[1][2]*A[2])%mo)%mo<<endl;
		}
		
	}
}

まとめ

なんでb,Bは2次にしたんだろ。

kmjp's blog

競技プログラミング参加記です

yukicoder : No.426 往復漸化式

問題

解法

まとめ