Harigami
FFTによる多項式の乗算
Public Domain
C
#include <complex.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
typedef double complex cmplx;
int exp_table[] = {
1,
2,
4,
8,
16,
32,
64,
128,
256,
512,
1024,
2048,
4096,
8192,
16384,
32768,
65536,
131072,
262144,
524288,
1048576,
2097152,
4194304,
8388608,
16777216,
33554432,
67108864,
134217728,
268435456,
536870912,
1073741824,
};
cmplx omega(int k, int v) {
return cexp(-v * 2 * M_PI * I / exp_table[k]);
}
int find_next_power(int n) {
// ループカウンタ
int i;
for(i = 0; i < 32; i++) {
if(n <= exp_table[i]) return i;
}
return -1;
}
// 2^k個の領域を確保
cmplx *alloc_fourier(int k) {
cmplx *result;
result = malloc(sizeof(cmplx) * exp_table[k]);
return result;
}
void output_cmplx_array(const cmplx *eA, int k) {
// ループカウンタ
int i;
for(i = 0; i < exp_table[k]; i++) {
printf("%f+%f\n", creal(eA[i]), cimag(eA[i]));
}
}
// A: 変換対象の配列
// k: 次数
// 次数を1ずつ下げながらDFTを行なう
void fftr(cmplx *A, int k) {
// feとfo
cmplx *fe, *fo;
// ループカウンタ
int i;
if(k <= 0) return;
// feとfoを確保する
fe = alloc_fourier(k-1);
fo = alloc_fourier(k-1);
// feとfoに値を詰め込む
for(i = 0; i < exp_table[k-1]; i++) {
fe[i] = A[i*2];
fo[i] = A[i*2+1];
}
// 再帰フーリエ変換
fftr(fe, k-1);
fftr(fo, k-1);
// 上書き
for(i = 0; i < exp_table[k-1]; i++) {
A[i] = fe[i] + omega(k, i) * fo[i];
A[i + exp_table[k-1]] = fe[i] - omega(k, i) * fo[i];
}
// 解放する
free(fe);
free(fo);
}
void conj_array(cmplx *A, int k) {
// ループカウンタ
int i;
// 共役をとる
for(i = 0; i < exp_table[k]; i++) {
A[i] = conj(A[i]);
}
}
void mult_array(const cmplx *A, const cmplx *B, cmplx *C, int k) {
// ループカウンタ
int i;
// 乗算する
for(i = 0; i < exp_table[k]; i++) {
C[i] = A[i] * B[i];
}
}
void ifft(cmplx *A, int k) {
// ループカウンタ
int i;
// 共役をとる
conj_array(A, k);
// fft
fftr(A, k);
// 共役をとる
conj_array(A, k);
// わる
for(i = 0; i < exp_table[k]; i++) A[i] /= exp_table[k];
}
void copy_int_to_cmplx(const int *A, cmplx *eA, int n) {
// ループカウンタ
int i;
for(i = 0; i < n; i++) eA[i] = A[i];
}
void convolve(const int *A, const int *B, int *C, int n) {
// 拡張されたA, B
cmplx *eA, *eB;
// 畳み込みの結果
cmplx *eC;
// n * 2のパワー
int power = find_next_power(n * 2);
// 拡張された長さ
int eLen = exp_table[power];
// ループカウンタ
int i;
// eA, eB, eCを確保,初期化
eA = malloc(sizeof(cmplx) * eLen);
eB = malloc(sizeof(cmplx) * eLen);
eC = malloc(sizeof(cmplx) * eLen);
memset(eA, 0, sizeof(cmplx) * eLen);
memset(eB, 0, sizeof(cmplx) * eLen);
copy_int_to_cmplx(A, eA, n);
copy_int_to_cmplx(B, eB, n);
// 巡回畳込みを行なう
fftr(eA, power);
fftr(eB, power);
mult_array(eA, eB, eC, power);
ifft(eC, power);
// 結果をCに書き戻す
for(i = 0; i < n * 2 - 1; i++) {
C[i] = round(creal(eC[i]));
}
// eA, eB, eCを解放
free(eA);
free(eB);
free(eC);
}
void input_array(int *A, int n) {
// ループカウンタ
int i;
for(i = 0; i < n; i++) scanf("%d", A + i);
}
void output_array(const int *C, int n) {
// ループカウンタ
int i;
for(i = 0; i < n; i++) printf("%d\n", C[i]);
}
int main(void) {
// 多項式の次数
int n;
// 入力、出力の配列サイズ
int isize, osize;
// 入力配列
int *A, *B;
// 出力配列
int *C;
// ループカウンタ
int i;
scanf("%d", &n);
isize = n + 1;
osize = (n + 1) * 2 - 1;
// nに応じてA, Bを確保
A = malloc(sizeof(int) * isize);
B = malloc(sizeof(int) * isize);
// 値を入力
for(i = 0; i < isize; i++) {
scanf("%d", A + i);
scanf("%d", B + i);
}
// nに応じてCを確保
C = malloc(sizeof(int) * osize);
// 畳み込みを行なう
convolve(A, B, C, isize);
// 出力
output_array(C, osize);
// 全て解放
free(A);
free(B);
free(C);
return 0;
}