拆系数FFT
表示才发现自己没有掌握这个似乎烂大街了的科技了……
概念：
应对那种模数比较恶心人的多项式乘法，大概就是吧一个多项式拆成两个，然后让乘法不会爆掉，最后再进行取模。既然拆成了两个多项式，\(DFT\)和\(IDFT\)次数自然就会变多，一共有\(7\)次的和\(4\)次的两种写法，自然是后面的快一些啦，但是后边的精度要求比较高，并且一般也不会卡的这么严重的，这里就只介绍第一种吧。
具体实现：
我们回忆\(FFT\)的过程，是将多项式转化为点值表示，相乘之后再插值回来。如果模数较好的话，模意义下\(NTT\)比较好，如果模数不好，但是数值范围较小的时候，\(FFT\)也是可以的。然而如果没有可以做\(NTT\)的模数，并且直接\(FFT\)强上会爆掉的时候，就需要拆系数\(FFT\)了。 我们想将多项式\(A(x)\)进行拆分，得到两个新的多项式\(B(x), C(x)\)。其中\(A_i = B_i * x^{\frac{n}{2}} + C_i\)，如此处理，让我们直接进行乘法的时候不会爆精。具体过程就是这样，假设我们做\(A * B\)的卷积，那么拆分后\[A(x) = C(x) * x^{\frac{n}{2}} + D(x)，B(x) = E(x) * x^{\frac{n}{2}} + F(x)\] 原来的卷积式变成了\[(C(x) * x^{\frac{n}{2}} + D(x)) * (E(x) * x^{\frac{n}{2}} + F(x))\] 暴力展开可得\[C(x) * E(x) * x^n + (C(x) * F(x) + D(x) * E(x)) * x^{\frac{n}{2}} + D(x) * F(x)\] 这样总计\(4\)次\(DFT\)，\(3\)次\(IDFT\)，常数变得超级大……
\(4\)次的坑以后填吧……
Code：
#include <bits/stdc++.h>
using namespace std;
const int N = 1e5 + 500;
typedef long long ll;
const double PI = acos(-1);
typedef vector<int> Vec;
int Md;
inline int Add(const int &x, const int &y) { return (x + y >= Md) ? (x + y - Md) : (x + y);}
inline int Sub(const int &x, const int &y) { return (x - y < 0) ? (x - y + Md) : (x - y);}
inline int Mul(const int &x, const int &y) { return (ll)x * y % Md;}
int Powe(int x, int y) {
int ans = 1;
while(y) {
if(y & 1) ans = Mul(ans, x);
x = Mul(x, x);
y >>= 1;
}
return ans;
}
namespace Poly {
struct Complex {
double x, y;
void operator = (int a) {
x = a;
y = 0;
}
friend Complex operator + (Complex A, Complex B) {
return (Complex) { A.x + B.x, A.y + B.y};
}
friend Complex operator - (Complex A, Complex B) {
return (Complex) { A.x - B.x, A.y - B.y};
}
friend Complex operator * (Complex A, Complex B) {
return (Complex) { A.x * B.x - A.y * B.y, A.x * B.y + A.y * B.x};
}
friend Complex operator / (Complex A, int B) {
return (Complex) { A.x / B, A.y / B};
}
} w[N << 2], A[N << 2], B[N << 2], C[N << 2], D[N << 2];
int rev[N << 2 | 1], l = 1;
void Pre(int len) {
for(; l < len; l <<= 1) {
for(int i = l; i < (l << 1); i++) {
w[i] = (Complex) { cos(PI / l * (i - l)), sin(PI / l * (i - l))};
}
}
for(int i = 0; i < len; i++) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? len >> 1 : 0);
}
}
void DFT(Complex *A, int len) {
for(int i = 0; i < len; i++) if(i < rev[i]) swap(A[i], A[rev[i]]);
for(int i = 1; i < len; i <<= 1) {
for(int j = 0; j < len; j += i << 1) {
Complex x, y;
for(int k = 0; k < i; k++) {
x = A[j + k], y = A[i + j + k] * w[i + k];
A[j + k] = x + y, A[i + j + k] = x - y;
}
}
}
}
void IDFT(Complex *A, int len) {
reverse(A + 1, A + len);
DFT(A, len);
for(int i = 0; i < len; i++) A[i] = A[i] / len;
}
Vec MUL(Vec a, Vec b) {
int n = a.size(), m = b.size(), len;
for(len = 1; len < n + m - 1; len <<= 1);
a.resize(len), b.resize(len);
Pre(len);
for(int i = 0; i < len; i++) {
A[i] = a[i] >> 15;
B[i] = a[i] & 32767;
C[i] = b[i] >> 15;
D[i] = b[i] & 32767;
}
DFT(A, len), DFT(B, len), DFT(C, len), DFT(D, len);
for(int i = 0; i < len; i++) {
Complex AA = A[i] * C[i], BB = A[i] * D[i] + B[i] * C[i], CC = B[i] * D[i];
A[i] = AA; B[i] = BB; C[i] = CC;
}
IDFT(A, len); IDFT(B, len); IDFT(C, len);
for(int i = 0; i < len; i++) {
ll x = llround(A[i].x) % Md, y =llround(B[i].x) % Md, z = llround(C[i].x) % Md;
a[i] = ((x << 30) % Md + (y << 15) % Md + z ) % Md;
}
a.resize(n + m - 1);
return a;
}
}
int n, m;
int main() {
scanf("%d%d%d", &n, &m, &Md);
Vec a(n + 1), b(m + 1);
for(int i = 0; i <= n; i++) scanf("%d", &a[i]);
for(int i = 0; i <= m; i++) scanf("%d", &b[i]);
a = Poly::MUL(a, b);
for(int i = 0; i < (int)a.size(); i++) printf("%d ", a[i]);
puts("");
return 0;
}
转载于:https://www.cnblogs.com/Apocrypha/p/10507645.html
相关资源：FFT_XILINX实现
