"Triangular Toeplitz"(三角テプリッツ行列)と"Circulant"(巡回行列)とそのBlock Matrix

参照

Matrix Reference Manual "Special Matrices"
Autumn Fest 2012 解説スライド "Ninja of Train"
あとwikipedia。(wikipediaをreferenceにするのって…)

このエントリについて

行列累乗の高速化としての問題が作れる、行列の特殊形を3つメモする。
実際の問題例と実装したコードも貼っておく。

Toeplitz行列

正方行列Aが"Toeplitz行列"であるとは、
ある関数fが存在して(A[i][j] = f(i-j))と表せることを言う。
直感的に言うと、「斜め(左上から右下へ)の要素が同じ」。
Toeplitzであるだけでは乗算の元で閉じていない(¬(∀x y, Toeplitz x ∧ Toeplitz y → Toeplitz (x * y)))ので、あまり使えない。

上三角Toeplitz行列

"上三角Toeplitz行列"とは、Toeplitz行列であって上三角行列(i > j ならば A[i][j] = 0)であるもの。
下三角行列も同様の性質を持っているが、簡単に変換ができるのでここでは上三角のほうだけを書く。
重要なことに、これは加算・乗算の元で閉じている。
さらに重要なことに、この行列は乗算がO(n^2)でできる。
この行列はサイズnのベクトルによって表すことが出来る。乗算はそのベクトルの畳込みになっている(ならFFTでO(n log n)でできたりするのかな？FFTはわかってないけど！)
式で実際に書くと(C[i] = Σ j=0..i (A[j] * B[i - j]))。
典型的には、ある状態からそれ(以上/以下)の状態にのみ遷移するような場合が三角行列になる。

巡回行列

あるサイズnの正方行列Aが巡回行列であるとは、
ある関数fが存在して(A[i][j] = f( (i-j) mod n))と表せることを言う。
(任意のkに対してA[i][j] = A[(i+k) mod n][(j+k) mod n])であることと同値(たぶん)。
巡回行列ならばToeplitz行列である。
この巡回行列も上三角Toeplitz行列と同じように、
加算・乗算の元で閉じていて、乗算がO(n^2)でできる。
またサイズnのベクトルによって表すことができ、
乗算の式は(C[i] = Σ j=0..n-1 (A[j] * B[(i-j) mod n]))と書ける。
典型的には、円状のボード上で相対的な動きをする場合に巡回行列になる。

Block Matrix

任意の行列は縦横いくつかの"Block"に分割して表現できる。
Blockで表現した時の乗算は(C[i][j] = Σk (A[i][k] * B[k][j]))と、普通の行列の乗算と同じ式で計算できる。
重要なことに、ある性質が乗算と加算のもとで閉じているなら、乗算・加算のもとで各ブロックの性質が閉じている。
つまり巡回行列や上三角Toeplitz行列を並べた"Block Matrix"も「使える」。
上三角Toeplitz行列などの"Block Matrix"の中である項の総和を取ったりしたい時がある。これは対角行列の形にしてやればいい(下の"Ninja of Train"のコード参照)。

SRM 566 DIV1 Medium PenguinEmperor

http://community.topcoder.com/stat?c=problem_statement&pm=12338
この問題の動きを(iからjの都市へ行く場合の数 = A[i][j])行列にするとそれは巡回行列になっている。
あとは行列累乗する。
時間制限がきついので剰余を減らす高速化をしないと通らない。

コード(下にあるライブラリを引いたもの)

#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
struct PenguinEmperor {
	int countJourneys(int n, long long daysPassed) {
		vector<CirculantMatrix> C(n, CirculantMatrix(n, n));
		rep(k, n) {
			rep(i, n)
				C[k].at((- k + n) % n) = C[k].at(k) = 1;
		}
		++ daysPassed;
		CirculantMatrix B = CirculantMatrix::identity(n);
		CirculantMatrix D = CirculantMatrix::identity(n);
		rep(k, n) {
			B *= C[k];
			if(k+1 == daysPassed % n) D = B;
		}
		B = B ^ (daysPassed / n);
		B *= D;
		return B.at(0);
	}
};

"Autumn Fest 2012" J "Ninja of Train"

http://autumn_fest.contest.atcoder.jp/tasks/autumn_fest_10
この問題では、(残り滞空時間 = ブロックのインデックス), (場所 = ブロック中のインデックス)とすると"Block (Upper/Lower) Toeplitz Matrix"になる。
それを累乗する。数え上げ部分はコード参照

コード(これもライブラリ部分を引いたもの, ただしそのMOD部分は1000003)

#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))
typedef vector<int> vi;
int main() {
	int h, t;
	scanf("%d%d", &h, &t);
	int s = (int)sqrt((double)h)+1;
	BlockMatrix V(vi(1, h), vi(s+1, h));
	V.at(0, 0) = BlockMatrix::Block::identity(h);
	vector<int> hs; rep(i, s+1) hs.push_back(h);
	BlockMatrix A(hs, hs);
	rep(i, h) rep(j, s) {
		if(j == 0) {
			A.at(0, s).at(i, i) = 1;
			rer(k, 1, h*2) {
				int t = (int)sqrt((double)k);
				int u = i+k-t;
				if(!(0 <= u && u < h)) continue;
				A.at(0, t-1).at(i, u) = 1;
			}
		}else {
			A.at(j, j-1).at(i, i) = 1;
		}
	}

	A.at(s, s) = BlockMatrix::Block::identity(h);
	BlockMatrix B = A ^ t+1;
	BlockMatrix C = V; C *= B;
	ll r = 0;
	rep(i, h) (r += C.at(0, s).at(i, h-1)) %= MOD;
	cout << r << endl;
	return 0;
}

行列とその特殊形ライブラリ

一時的なオブジェクトを確保するのが嫌なのでダブルバッファ的にやっている。
assert多めだったり冗長な部分が多いので使う時に消すべき。特にassertは提出時に"#define assert(_) ((void)0)"するべき
剰余を減らす高速化をしているが、MODが√(2^64-0xf*2^60)くらいより大きい場合はこれだと駄目。書き換えないといけない。
高速化は素人なのでもっとできるとは思う。ループはアンロールしたらいいかもしれない。

//このコードはいかなることも自由に使っていいです
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))
#define reu(i,l,u) for(int (i)=(int)(l);(i)<(int)(u);++(i))

typedef unsigned long long Num;
#define MOD 1000000007ULL	//！！！！！！！MODが大きい場合、"*="を書き換えること！

struct Matrix {
	vector<vector<Num> > v, w;
	Matrix() {}
	Matrix(int n, int m): v(n, vector<Num>(m)) { assert(n > 0 && m > 0); }
	inline int height() const { return (int)v.size(); }
	inline int width() const { return (int)v[0].size(); }
	inline Num& at(int i, int j) { assert(0 <= i && i < height() && 0 <= j && j < width()); return v[i][j]; }
	inline const Num& at(int i, int j) const { assert(0 <= i && i < height() && 0 <= j && j < width()); return v[i][j]; }
	static Matrix identity(int n) {
		Matrix A(n, n);
		rep(i, n) A.at(i, i) = 1;
		return A;
	}
	inline static Matrix identity(const Matrix& A) {
		assert(A.height() == A.width()); return identity(A.height()); }
	Matrix& operator*=(const Matrix& B) {
		int n = height(), m = B.width(), p = B.height();
		assert(p == width());
		w.resize(n, vector<Num>(m, 0));
		rep(i, n) rep(j, m) {
			//MOD = 1000000007の場合、1000000012000000036となり、
			//2^64にオーバーフローするのは2^64-1-1000000006^2くらいで、それは最上位4bitが立っているから
			//！MODが違う場合、書き換えること！
			Num x = 0;
			rep(k, p) {
				x += at(i, k) * B.at(k, j);
				if((x >> 60) == 0xf) x %= MOD;
			}
			x %= MOD;
			w[i][j] = x;
		}
		v.swap(w);
		return *this;
	}
	Matrix& operator+=(const Matrix& B) {
		int n = height(), m = width();
		assert(n == B.height() && m == B.width());
		rep(i, n) rep(j, m) {
			at(i, j) += B.at(i, j);
			if(at(i, j) >= MOD) at(i, j) -= MOD;
		}
		return *this;
	}
	void undomult() { v.swap(w); }
};
ostream& operator<<(ostream& o, const Matrix& A) {
	int n = A.height(), m = A.width();
	rep(i, n) {
		o << "["; rep(j, m) o << A.at(i, j) << (j+1 == m ? "]\n" : ",");
	}
	return o;
}

/* ジェネリックなべき乗算 */
template<typename Mat>
Mat operator^(const Mat& t, ll k) {
	Mat A = t, B = Mat::identity(t);
	while(k) {
		if(k & 1) {
			B *= A;
		}
		A *= A;
		k >>= 1;
	}
	return B;
}

struct CirculantMatrix {
	vector<Num> v, w;
	CirculantMatrix() {}
	CirculantMatrix(int n, int m): v(n, 0) { (void)m; assert(n == m); assert(n > 0); }
	CirculantMatrix(const Matrix& A) {
		int n = A.height();
		assert(n > 0);
		assert(n == A.width());
		rep(i, n) rep(j, n) assert(A.at(i, j) == A.at((i+1)%n, (j+1)%n));
		v.assign(n, 0);
		rep(i, n) v[i] = A.at(i, 0);
	}
	inline int size() const { return (int)v.size(); }
	inline int height() const { return size(); }
	inline int width() const { return size(); }
	inline Num& at(int i) { assert(0 <= i && i < size()); return v[i]; }
	inline const Num& at(int i) const { assert(0 <= i && i < size()); return v[i]; }
	inline Num& at(int i, int j) { assert(0 <= i && i < size() && 0 <= j && j < size()); return v[(i - j + size()) % size()]; }
	inline const Num& at(int i, int j) const { assert(0 <= i && i < size() && 0 <= j && j < size()); return v[(i - j + size()) % size()]; }

	inline static CirculantMatrix identity(int n) {
		CirculantMatrix A(n, n);
		A.at(0) = 1;
		return A;
	}
	inline static CirculantMatrix identity(const CirculantMatrix& A) { return identity(A.size()); }
	CirculantMatrix& operator*=(const CirculantMatrix& B) {
		int n = size(); assert(n == B.size()); w.resize(n);
		Num *vv = &v[0], *ww = &w[0]; const Num *zz = &B.v[0];
		rep(j, n) {
			Num x = 0;
			rer(k, 0, j) {
				x += vv[k] * zz[j-k];
				if((x >> 60) == 0xf) x %= MOD;
				//MOD = 1000000007の場合、1000000012000000036となり、
				//2^64にオーバーフローするのは2^64-1-1000000006^2くらいで、それは最上位4bitが立っているから
				//！MODが違う場合、ここと下の部分も書き換えること！
			}
			rer(k, j+1, n-1) {
				x += vv[k] * zz[j-k+n];
				if((x >> 60) == 0xf) x %= MOD;
			}
			x %= MOD;
			ww[j] = x;
		}
		v.swap(w);
		return *this;
	}
	CirculantMatrix& operator+=(const CirculantMatrix& B) {
		int n = size();
		assert(n == B.size());
		rep(i, n) {
			v[i] += B.v[i];
			if(v[i] >= MOD) v[i] -= MOD;
		}
		return *this;
	}
	void undomult() { v.swap(w); }
};

ostream& operator<<(ostream& o, const CirculantMatrix& A) {
	int n = A.size();
	//普通の行列のように表示する
//	rep(i, n) {
//		o << "["; rep(j, n) o << A[(i - j + n) % n] << (j+1 == n ? "]\n" : ",");
//	}
	o << "["; rep(j, n) o << A.v[j] << (j+1 == n ? "]" : ",");
	return o;
}

struct UTToeplitzMatrix {
	vector<Num> v, w;	//A[0][i]を保持するイメージ
	static const Num Zero = 0;
	UTToeplitzMatrix() {}
	UTToeplitzMatrix(int n, int m): v(n, 0) { (void)m; assert(n == m); assert(n > 0); }
	inline int size() const { return (int)v.size(); }
	inline int height() const { return size(); }
	inline int width() const { return size(); }
	inline Num& at(int i) { assert(0 <= i && i < size()); return v[i]; }
	inline const Num& at(int i) const { assert(0 <= i && i < size()); return v[i]; }
	inline Num& at(int i, int j) {
		assert(0 <= i && i < size() && 0 <= j && j < size());
		assert(i <= j);	//0の部分を変更することはできない
		return v[j - i];
	}
	inline const Num& at(int i, int j) const {
		assert(0 <= i && i < size() && 0 <= j && j < size());
		if(i > j) return Zero;	//0の部分
		return v[j - i];
	}
	static UTToeplitzMatrix identity(int n) {
		UTToeplitzMatrix A(n, n);
		A.at(0) = 1;
		return A;
	}
	inline static UTToeplitzMatrix identity(const UTToeplitzMatrix& A) { return identity(A.size()); }
	UTToeplitzMatrix& operator*=(const UTToeplitzMatrix& B) {
		int n = size(); assert(n == B.size()); w.resize(n);
		Num *vv = &v[0], *ww = &w[0]; const Num *zz = &B.v[0];
		rep(i, n) {
			Num x = 0;
			rer(j, 0, i) {
				x += v[j] * B.v[i - j];
				if((x >> 60) == 0xf) x %= MOD;
			}
			x %= MOD;
			w[i] = x;
		}
		v.swap(w);
		return *this;
	}
	UTToeplitzMatrix& operator+=(const UTToeplitzMatrix& B) {
		int n = size();
		assert(n == B.size());
		rep(i, n) {
			v[i] += B.v[i];
			if(v[i] >= MOD) v[i] -= MOD;
		}
		return *this;
	}
	void undomult() { v.swap(w); }
};

struct BlockMatrix {
	typedef UTToeplitzMatrix Block;
	vector<vector<Block> > blocks, blocks2;
	vector<int> heights, widths;
	BlockMatrix(const vector<int>& hs, const vector<int>& ws)
		: heights(hs), widths(ws) {
		int n = heights.size(), m = widths.size();
		assert(n > 0 && m > 0);
		blocks.resize(n);
		rep(i, n) rep(j, m)
			blocks[i].push_back(Block(heights[i], widths[j]));
	}
	inline int height() const { return heights.size(); }
	inline int width() const { return widths.size(); }
	inline int blockHeight(int i) const { assert(0 <= i && i < height()); return blocks[i][0].height(); }
	inline int blockWidth(int j) const { assert(0 <= j && j < width()); return blocks[0][j].width(); }
	inline Block& at(int i, int j) { assert(0 <= i && i < height() && 0 <= j && j < width()); return blocks[i][j]; }
	inline const Block& at(int i, int j) const { assert(0 <= i && i < height() && 0 <= j && j < width()); return blocks[i][j]; }
	static BlockMatrix identity(int n, int m) {
		vector<int> v;
		rep(i, n) v.push_back(m);
		BlockMatrix A(v, v);
		rep(i, n) A.at(i, i) = Block::identity(m);
		return A;
	}
	inline static BlockMatrix identity(const BlockMatrix& A) {
		int n = A.height(); assert(n == A.width());
		int m = A.heights[0];
		assert(count(A.heights.begin(), A.heights.end(), m) == n);
		assert(count(A.widths.begin(), A.widths.end(), m) == n);
		return identity(n, m);
	}
	BlockMatrix& operator*=(const BlockMatrix& B) {
		int n = height(), m = B.width(), p = B.height();
		assert(p == width());
		blocks2.resize(n, vector<Block>(m, Block()));
		rep(i, n) rep(j, m) {
			Block x = Block(blockHeight(i), B.blockWidth(j));
			rep(k, p) {
				at(i, k) *= B.at(k, j);
				x += at(i, k);
				at(i, k).undomult();
			}
			blocks2[i][j] = x;
		}
		blocks.swap(blocks2);
		return *this;
	}
};
ostream& operator<<(ostream& o, const BlockMatrix& A) {
	int n = A.height(), m = A.width();
	rep(i, n) rep(i2, A.blockHeight(i)){
		o << "[";
		rep(j, m) rep(j2, A.blockWidth(j))
			o << A.at(i, j).at(i2, j2) << (j+1 == m && j2+1 == A.blockWidth(j) ? "]\n" : ",");
	}
	return o;
}