cubelover의 블로그

きたまさ法

수학2016. 11. 24. 01:51

키타마사법은 다음과 같은 점화식을 가지는 수열의 \(n\)번째 항을 빠르게 계산하는 데 사용된다.


\[a_{n}=\sum_{k=1}^{m}a_{n-k}c_{k}\]


이 점화식을 단순 행렬곱으로 풀면 \(O(m^{3}\log n)\)의 시간이 필요하지만, 키타마사법을 사용하면 \(O(m^{2}\log n)\)의 시간이면 충분하다.


~


임의의 \(a_{n}\)을 다음처럼 \(a_{1}, a_{2}, \dots, a_{m}\)에 대한 식으로 표현할 수 있다.


\[a_{n}=\sum_{k=1}^{m}a_{k}d_{k}\]


키타마사법의 아이디어는, 위 식을 만족하는 \(n\), \(d_{k}\)와 임의의 정수 \(l\)에 대하여 다음이 성립한다는 것이다.


\[a_{n+l}=\sum_{k=1}^{m}a_{k+l}d_{k}\]


이 식이 성립하는 이유는 수학적 귀납법으로 쉽게 보일 수 있으므로 생략한다.


이를 통해 \(a_{2n}\)에 대한 식을 정리하면 다음과 같다.


\[a_{2n}=\sum_{k=1}^{m}a_{n+k}d_{k}=\sum_{k=1}^{m}\sum_{j=1}^{m}a_{j+k}d_{j}d_{k}\]


따라서 \(a_{2n}\)을 \(a_{2}, a_{3}, \dots, a_{2m}\)에 관한 선형식으로 표현할 수 있고, 이들의 계수는 \(O(m^{2})\) 시간 안에 계산할 수 있다.


또한 \(a_{m+1}, a_{m+2}, \dots, a_{2m}\)을 \(a_{1}, a_{2}, \dots, a_{m}\)에 대한 선형식으로 나타내는 데 \(O(m^{2})\)의 시간이면 충분하다.


이를 종합하면, \(a_{n}\)에 대한 식을 \(a_{2n}\)에 대한 식으로 바꾸는 데 \(O(m^{2})\)의 시간밖에 소요되지 않으므로 \(O(m^{2}\log n)\)의 시간복잡도로 \(n\)번째 항을 계산할 수 있다.


#include <cstdio>

const int P = 1000000007;

int a[1001];
int c[1001];
int d[1001];
int t[2002];

void Kitamasa(int n, int m) {
	int i, j;
	if (n == 1) {
		d[1] = 1;
		for (i = 2; i <= m; i++) d[i] = 0;
		return;
	}
	if (n & 1) {
		Kitamasa(n ^ 1, m);
		j = d[m];
		for (i = m; i >= 1; i--) d[i] = (d[i - 1] + (long long)c[m - i + 1] * j) % P;
		return;
	}
	Kitamasa(n >> 1, m);
	for (i = 1; i <= m + m; i++) t[i] = 0;
	for (i = 1; i <= m; i++) for (j = 1; j <= m; j++) t[i + j] = (t[i + j] + (long long)d[i] * d[j]) % P;
	for (i = m + m; i > m; i--) for (j = 1; j <= m; j++) t[i - j] = (t[i - j] + (long long)c[j] * t[i]) % P;
	for (i = 1; i <= m; i++) d[i] = t[i];
}

int main() {
	int i, n, m, r;
	scanf("%d%d", &n, &m);
	for (i = 1; i <= m; i++) scanf("%d", &a[i]);
	for (i = 1; i <= m; i++) scanf("%d", &c[i]);
	Kitamasa(n, m);
	r = 0;
	for (i = 1; i <= m; i++) r = (r + (long long)a[i] * d[i]) % P;
	printf("%d", r);
}

'수학' 카테고리의 다른 글

Dividing Polynomials  (0) 2017.01.12
Discrete Fourier Transform  (0) 2016.11.17
Number Theoretic Fast Fourier Transform  (1) 2016.11.02