拆系数FFT

拆系数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;
}