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
きたまさ法  (0) 2016.11.24
Discrete Fourier Transform  (0) 2016.11.17
Number Theoretic Fast Fourier Transform  (1) 2016.11.02

Comment +0

Discrete Fourier Transform

수학2016. 11. 17. 14:27

이산 푸리에 변환이란 주기가 \(N\)인 복소수열 \(x_{0},x_{1},\dots,x_{N-1}\)에 대하여 다음과 같이 정의되는 주기가 \(N\)인 복소수열 \(X_{0},X_{1},\dots,X_{N-1}\)로 바꾸는 변환을 말한다.


\[X_{k}=\sum_{n=0}^{N-1}e^{-2\pi ikn/N}x_{n}\]


여기서 \(e^{iz}=\cos z+i\sin z\)이며, 이로부터 \(e^{2\pi i}=e^{-2\pi i}=1\)임을 쉽게 알 수 있다.


~


이산 푸리에 역변환


또한 이렇게 변환된 복소수열 \(X_{k}\)는 다음과 같은 과정을 통하여 복소수열 \(x_{n}\)으로 역변환될 수 있다.


\[x_{n}=\frac{1}{N}\sum_{k=0}^{N-1}e^{2\pi ikn/N}X_{k}\]


~


이산 푸리에 역변환의 증명


\(X_{k}\)의 정의대로 식을 풀어보면 다음과 같다.


\[\sum_{k=0}^{N-1}e^{2\pi ikn/N}X_{k}=\sum_{k=0}^{N-1}\sum_{j=0}^{N-1}e^{2\pi ik(n-j)/N}x_{j}=\sum_{j=0}^{N-1}x_{j}\sum_{k=0}^{N-1}e^{2\pi ik(n-j)/N}\]


\(j=n\)인 경우에는 다음처럼 된다.


\[\sum_{k=0}^{N-1}e^{2\pi ik(n-j)/N}=\sum_{k=0}^{N-1}1=N\]


\(j\neq n\)인 경우에는 다음처럼 된다.


\[\sum_{k=0}^{N-1}e^{2\pi ik(n-j)/N}=\frac{e^{2\pi i(n-j)}-1}{e^{2\pi i(n-j)/N}-1}=0\]


따라서 다음이 성립한다.


\[\frac{1}{N}\sum_{k=0}^{N-1}e^{2\pi ikn/N}X_{k}=\frac{1}{N}Nx_{n}=x_{n}\]


~


이산 푸리에 변환과 Convolution


주기가 \(N\)인 수열 \(a_{n}\)과 \(b_{n}\)을 각각 이산 푸리에 변환한 수열을 \(A_{n}\)과 \(B_{n}\)이라 하자.


또한 새로운 수열 \(C_{n}=A_{n}B_{n}\)을 정의하고, 이를 이산 푸리에 역변환한 수열을 \(c_{n}\)이라 하면 다음과 같은 식이 성립한다.


\[c_{n}=\sum_{k=0}^{N-1}a_{k}b_{n-k}\]


~


이산 푸리에 변환과 Convolution의 증명


\(C_{n}\)의 정의대로 식을 풀어보면 다음과 같다.


\[C_{n}=A_{n}B_{n}=\sum_{j=0}^{N-1}\sum_{k=0}^{N-1}e^{-2\pi i(j+k)n/N}a_{k}b_{j}=\sum_{j=0}^{N-1}e^{-2\pi ijn/N}\sum_{k=0}^{N-1}a_{k}b_{j-k}\]


따라서 이를 역변환한 수열 \(c_{n}\)에 대한 식이 성립한다.


~


고속 푸리에 변환


이산 푸리에 변환을 단순한 반복문으로 계산하면 \(O(N^{2})\)의 시간복잡도를 가지지만, 이를 \(O(N\log N)\)의 시간복잡도에 행할 수 있는 방법이 고속 푸리에 변환이다.


고속 푸리에 변환에는 다양한 알고리즘이 있지만, 여기서는 \(N=2^{M}\) 꼴의 푸리에 변환에 최적화된 Cooley-Tukey 알고리즘에 대해 설명한다.


\(X_{k}\)를 다음처럼 고쳐 적을 수 있다.


\[X_{k}=\sum_{n=0}^{N-1}e^{-2\pi ikn/N}x_{n}=\sum_{n=0}^{N/2-1}e^{-2\pi ik(2n)/N}x_{2n}+\sum_{n=0}^{N/2-1}e^{-2\pi ik(2n+1)/N}x_{2n+1}=\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n}+e^{-2\pi ik/N}\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n+1}\]


이제 수열 \(Y_{k}\)와 \(Z_{k}\)를 다음처럼 정의하자.


\[Y_{k}=\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n}\]

\[Z_{k}=\sum_{n=0}^{N/2-1}e^{-2\pi ikn/(N/2)}x_{2n+1}\]


수열 \(Y_{k}\)와 \(Z_{k}\)를 계산할 수 있다면, 수열 \(Z_{k}\)는 다음과 같은 방법으로 \(O(N)\)의 시간복잡도로 계산할 수 있다.


\[X_{k}=Y_{k}+e^{-2\pi ik/N}Z_{k}\]

\[X_{k+N/2}=Y_{k}-e^{-2\pi ik/N}Z_{k}\]


그런데, 수열 \(Y_{k}\)와 \(Z_{k}\)는 수열 \(x_{2n}\)과 \(x_{2n+1}\)에 대한 이산 푸리에 변환이므로, 분할정복을 통해 총 \(O(N\log N)\)의 시간복잡도로 이산 푸리에 변환을 구할 수 있다.


#include <cstdio>
#include <cmath>
#include <complex>

const int SZ = 20, N = 1 << SZ;
const double PI = acos(-1);

using namespace std;

int Rev(int x) {
	int i, r = 0;
	for (i = 0; i < SZ; i++) {
		r = r << 1 | x & 1;
		x >>= 1;
	}
	return r;
}

void FFT(complex<double> *a, bool f) {
	complex<double> x, y, z;
	double w;
	int i, j, k;
	for (i = 0; i < N; i++) {
		j = Rev(i);
		if (i < j) {
			z = a[i];
			a[i] = a[j];
			a[j] = z;
		}
	}
	for (i = 1; i < N; i <<= 1) {
		w = PI / i;
		if (f) w = -w;
		x = complex<double>(cos(w), sin(w));
		for (j = 0; j < N; j += i << 1) {
			y = complex<double>(1);
			for (k = 0; k < i; k++) {
				z = a[i + j + k] * y;
				a[i + j + k] = a[j + k] - z;
				a[j + k] += z;
				y *= x;
			}
		}
	}
	if (f) for (i = 0; i < N; i++) a[i] /= N;
}

complex<double> X[N];

int main() {
	double x;
	int i, n;
	scanf("%d", &n);
	for (i = 0; i <= n; i++) {
		scanf("%lf", &x);
		X[i] = complex<double>(x, 0);
	}
	FFT(X, false);
	for (i = 0; i < N; i++) X[i] *= X[i];
	FFT(X, true);
	for (i = 0; i <= n + n; i++) printf("%.0f ", X[i].real());
}


~


일반화된 이산 푸리에 변환


지금까지는 함수 \(\mathbb{Z}_{N}\rightarrow\mathbb{C}\)에 대한 이산 푸리에 변환에 대해서만 살펴보았다.


한편 이산 푸리에 변환에 필요한 연산은 덧셈과 복소수곱셈뿐이므로 이 둘이 정의된 모든 집합에 대해서 이산 푸리에 변환이 가능하다.


예를 들어, 함수 \(\mathbb{Z}_{M}\rightarrow\mathbb{C}\)는 덧셈과 복소수곱셈이 가능하므로 \(\mathbb{Z}_{N}\rightarrow\mathbb{Z}_{M}\rightarrow\mathbb{C}\)에 대한 이산 푸리에 변환이 가능하다.


이를 통해 두 \(N\times M\) 행렬 \(A_{i,j}\)와 \(B_{i,j}\)에 대하여 다음을 만족하는 행렬 \(C_{i,j}\)를 구할 수 있다.


\[C_{i,j}=\sum_{k=0}^{N-1}\sum_{l=0}^{M-1}A_{k,l}B_{i-k,j-l}\]


template을 사용하여 구현하면 다음과 같다.


#include <cstdio>
#include <cmath>
#include <complex>

const int SZ = 10, N = 1 << SZ;
const double PI = acos(-1);

using namespace std;

int Rev(int x) {
	int i, r = 0;
	for (i = 0; i < SZ; i++) {
		r = r << 1 | x & 1;
		x >>= 1;
	}
	return r;
}

template <typename T>
struct Func {
	T a[N];

	T &operator [](int x) {
		return a[x];
	}

	Func<T> operator +(Func<T> &x) {
		Func<T> r;
		int i;
		for (i = 0; i < N; i++) r.a[i] = a[i] + x[i];
		return r;
	}

	Func<T> operator -(Func<T> &x) {
		Func<T> r;
		int i;
		for (i = 0; i < N; i++) r.a[i] = a[i] - x[i];
		return r;
	}

	Func<T> operator *(complex<double> &x) {
		Func<T> r;
		int i;
		for (i = 0; i < N; i++) r.a[i] = a[i] * x;
		return r;
	}

	Func<T> operator *=(complex<double> &x) {
		int i;
		for (i = 0; i < N; i++) a[i] *= x;
		return *this;
	}

	Func<T> operator *=(Func<T> &x) {
		int i;
		for (i = 0; i < N; i++) a[i] *= x[i];
		return *this;
	}

	void FFT(bool f) {
		complex<double> x, y;
		T z;
		double w;
		int i, j, k;
		for (i = 0; i < N; i++) {
			j = Rev(i);
			if (i < j) {
				z = a[i];
				a[i] = a[j];
				a[j] = z;
			}
		}
		for (i = 1; i < N; i <<= 1) {
			w = PI / i;
			if (f) w = -w;
			x = complex<double>(cos(w), sin(w));
			for (j = 0; j < N; j += i << 1) {
				y = complex<double>(1);
				for (k = 0; k < i; k++) {
					z = a[i + j + k] * y;
					a[i + j + k] = a[j + k] - z;
					a[j + k] = a[j + k] + z;
					y *= x;
				}
			}
		}
		if (f) {
			x = complex<double>((double)1 / N);
			for (i = 0; i < N; i++) a[i] *= x;
		}
	}
};

Func<Func<complex<double> > > X;

int main() {
	double x;
	int i, j, n, m;
	scanf("%d%d", &n, &m);
	for (i = 0; i <= n; i++) for (j = 0; j <= m; j++) {
		scanf("%lf", &x);
		X[i][j] = complex<double>(x, 0);
	}
	for (j = 0; j < N; j++) X[j].FFT(false);
	X.FFT(false);
	X *= X;
	X.FFT(true);
	for (j = 0; j < N; j++) X[j].FFT(true);
	for (i = 0; i <= n + n; i++) {
		for (j = 0; j <= m + m; j++) printf("%.0f ", X[i][j].real());
		printf("\n");
	}
}


~


XOR Convolution


\(2^{M}\) 크기의 배열 \(A_{i}\)와 \(B_{i}\)에 대하여 다음처럼 정의되는 배열 \(C_{i}\)를 구하는 문제를 생각해 보자. (\(\oplus\)는 Bitwise XOR 연산을 뜻함)


\[C_{i}=\sum_{k=0}^{2^{M}-1}A_{k}B_{i\oplus k}\]


XOR 연산은 각 자리별로 더한 뒤 2로 나눈 나머지를 구하는 것과 같으므로, \(\mathbb{Z}_{2^{M}}\rightarrow\mathbb{C}\)를 다음처럼 펼쳐서 생각해 보자.


\[\mathbb{Z}_{2}\rightarrow\mathbb{Z}_{2}\rightarrow\cdots\rightarrow\mathbb{Z}_{2}\rightarrow\mathbb{C}\]


각각이 이산 푸리에 변환 가능하므로, \(\mathbb{Z}_{2}\)에 대한 이산 푸리에 변환을 \(M\)번 행하는 것임을 알 수 있다.


#include <cstdio>
#include <complex>

const int SZ = 20, N = 1 << SZ;

using namespace std;

int Rev(int x) {
	int i, r = 0;
	for (i = 0; i < SZ; i++) {
		r = r << 1 | x & 1;
		x >>= 1;
	}
	return r;
}

void FFT(int *a, bool f) {
	int i, j, k, z;
	for (i = 0; i < N; i++) {
		j = Rev(i);
		if (i < j) {
			z = a[i];
			a[i] = a[j];
			a[j] = z;
		}
	}
	for (i = 1; i < N; i <<= 1) for (j = 0; j < N; j += i << 1) for (k = 0; k < i; k++) {
		z = a[i + j + k];
		a[i + j + k] = a[j + k] - z;
		a[j + k] += z;
	}
	if (f) for (i = 0; i < N; i++) a[i] /= N;
}

int X[N];

int main() {
	int i, n;
	scanf("%d", &n);
	for (i = 0; i < 1 << n; i++) scanf("%d", &X[i]);
	FFT(X, false);
	for (i = 0; i < N; i++) X[i] *= X[i];
	FFT(X, true);
	for (i = 0; i < 1 << n; i++) printf("%d ", X[i]);
}



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

Dividing Polynomials  (0) 2017.01.12
きたまさ法  (0) 2016.11.24
Discrete Fourier Transform  (0) 2016.11.17
Number Theoretic Fast Fourier Transform  (1) 2016.11.02

Comment +0

Generate Connected Graph

DM2016. 11. 7. 02:03
#include "testlib.h"
#include <vector>
#include <unordered_set>

std::vector<int> P;

int Root(int X) {
	return X == P[X] ? X : P[X] = Root(P[X]);
}

void GenConnectedGraph(int N, int M, std::vector<std::pair<int, int> > &R) {

	if (N < 1 || M + 1 < N || M > (long long)N * (N - 1) / 2) return;

	P.resize(N);
	for (int i = 0; i < N; i++) P[i] = i;

	std::unordered_set<long long> S;

	R.resize(M);
	for (int i = 0; i < N - 1; i++) {
		int u, v;
		do {
			u = rnd.next(N);
			v = rnd.next(N);
		} while (Root(u) == Root(v));
		P[Root(u)] = Root(v);
		R[i].first = u;
		R[i].second = v;
		S.insert((long long)u << 32 | v);
	}

	for (int i = N - 1; i < M; i++) {
		int u, v;
		do {
			u = rnd.next(N);
			v = rnd.next(N);
		} while (u == v || S.find((long long)u << 32 | v) != S.end() || S.find((long long)v << 32 | u) != S.end());
		R[i].first = u;
		R[i].second = v;
		S.insert((long long)u << 32 | v);
	}

	shuffle(R.begin(), R.end());
}

int main(int argc, char *argv[]) {
	registerGen(argc, argv, 1);

	int N = 7;
	int M = 10;

	std::vector<std::pair<int, int> > E;
	GenConnectedGraph(N, M, E);

	printf("%d %d\n", N, M);
	for (int i = 0; i < M; i++) printf("%d %d\n", E[i].first, E[i].second);

	return 0;
}

'DM' 카테고리의 다른 글

Generate Connected Graph  (0) 2016.11.07
Generate Tree  (0) 2016.11.07
Testlib  (0) 2016.11.07

Comment +0