TOP > 数学
任意mod畳み込み(Arbitrary-Mod-Convolution)
説明
高速フーリエ変換による任意 mod 畳み込みを行う。
計算量
- $O((n + m) \log (n + m))$
実装例
依存ライブラリ Mod-Int Fast-Fourier-Transform
テンプレート引数としてMod-Intが渡されることを想定している。
- multiply($a$, $b$):= 配列 $a$ と配列 $b$ を畳み込みした結果を返す。
以下は mod $10^9+7$ で verified。ふつうはこっち。
template< typename T >
struct ArbitraryModConvolution {
using real = FastFourierTransform::real;
using C = FastFourierTransform::C;
ArbitraryModConvolution() = default;
vector< T > multiply(const vector< T > &a, const vector< T > &b, int need = -1) {
if(need == -1) need = a.size() + b.size() - 1;
int nbase = 0;
while((1 << nbase) < need) nbase++;
FastFourierTransform::ensure_base(nbase);
int sz = 1 << nbase;
vector< C > fa(sz);
for(int i = 0; i < a.size(); i++) {
fa[i] = C(a[i].x & ((1 << 15) - 1), a[i].x >> 15);
}
fft(fa, sz);
vector< C > fb(sz);
if(a == b) {
fb = fa;
} else {
for(int i = 0; i < b.size(); i++) {
fb[i] = C(b[i].x & ((1 << 15) - 1), b[i].x >> 15);
}
fft(fb, sz);
}
real ratio = 0.25 / sz;
C r2(0, -1), r3(ratio, 0), r4(0, -ratio), r5(0, 1);
for(int i = 0; i <= (sz >> 1); i++) {
int j = (sz - i) & (sz - 1);
C a1 = (fa[i] + fa[j].conj());
C a2 = (fa[i] - fa[j].conj()) * r2;
C b1 = (fb[i] + fb[j].conj()) * r3;
C b2 = (fb[i] - fb[j].conj()) * r4;
if(i != j) {
C c1 = (fa[j] + fa[i].conj());
C c2 = (fa[j] - fa[i].conj()) * r2;
C d1 = (fb[j] + fb[i].conj()) * r3;
C d2 = (fb[j] - fb[i].conj()) * r4;
fa[i] = c1 * d1 + c2 * d2 * r5;
fb[i] = c1 * d2 + c2 * d1;
}
fa[j] = a1 * b1 + a2 * b2 * r5;
fb[j] = a1 * b2 + a2 * b1;
}
fft(fa, sz);
fft(fb, sz);
vector< T > ret(need);
for(int i = 0; i < need; i++) {
int64_t aa = llround(fa[i].x);
int64_t bb = llround(fb[i].x);
int64_t cc = llround(fa[i].y);
aa = T(aa).x, bb = T(bb).x, cc = T(cc).x;
ret[i] = aa + (bb << 15) + (cc << 30);
}
return ret;
}
};
以下は $10^{11}$ 以下の素数 mod で verified。FFT の精度を long double にすること。
template< typename T >
struct ArbitraryModConvolutionLong {
using real = FastFourierTransform::real;
using C = FastFourierTransform::C;
ArbitraryModConvolutionLong() = default;
vector< T > multiply(const vector< T > &a, const vector< T > &b, int need = -1) {
if(need == -1) need = a.size() + b.size() - 1;
int nbase = 0;
while((1 << nbase) < need) nbase++;
FastFourierTransform::ensure_base(nbase);
int sz = 1 << nbase;
vector< C > fa(sz);
for(int i = 0; i < a.size(); i++) {
fa[i] = C(a[i].x & ((1 << 19) - 1), a[i].x >> 19);
}
fft(fa, sz);
vector< C > fb(sz);
if(a == b) {
fb = fa;
} else {
for(int i = 0; i < b.size(); i++) {
fb[i] = C(b[i].x & ((1 << 19) - 1), b[i].x >> 19);
}
fft(fb, sz);
}
real ratio = 0.25 / sz;
C r2(0, -1), r3(ratio, 0), r4(0, -ratio), r5(0, 1);
for(int i = 0; i <= (sz >> 1); i++) {
int j = (sz - i) & (sz - 1);
C a1 = (fa[i] + fa[j].conj());
C a2 = (fa[i] - fa[j].conj()) * r2;
C b1 = (fb[i] + fb[j].conj()) * r3;
C b2 = (fb[i] - fb[j].conj()) * r4;
if(i != j) {
C c1 = (fa[j] + fa[i].conj());
C c2 = (fa[j] - fa[i].conj()) * r2;
C d1 = (fb[j] + fb[i].conj()) * r3;
C d2 = (fb[j] - fb[i].conj()) * r4;
fa[i] = c1 * d1 + c2 * d2 * r5;
fb[i] = c1 * d2 + c2 * d1;
}
fa[j] = a1 * b1 + a2 * b2 * r5;
fb[j] = a1 * b2 + a2 * b1;
}
fft(fa, sz);
fft(fb, sz);
vector< T > ret(need);
auto mul1 = T(2).pow(19);
auto mul2 = T(2).pow(38);
for(int i = 0; i < need; i++) {
int64_t aa = llround(fa[i].x);
int64_t bb = llround(fb[i].x);
int64_t cc = llround(fa[i].y);
aa = T(aa).x, bb = T(bb).x, cc = T(cc).x;
ret[i] = (mul1 * bb) + (mul2 * cc) + aa;
}
return ret;
}
};
検証
int main() {
int N;
cin >> N;
vector< modint > A(N + 1), B(N + 1);
for(int i = 1; i <= N; i++) cin >> A[i] >> B[i];
ArbitraryModConvolution< modint > fft;
auto C = fft.multiply(A, B);
for(int i = 1; i <= 2 * N; i++) cout << C[i] << "\n";
}