# FFT&NTT小结

## Process

##### 1.DFT

$F(x) = (a_0 + a_2 \cdot x^2 + \cdots + a_{n - 2} \cdot x^{n - 2}) + (a_1 \cdot x + a_3 \cdot x^3 + \cdots + a_{n - 1} \cdot x^{n - 1})$

$G(x) = a_0 + a_2 \cdot x + \cdots + a_{n - 2} \cdot x^{\frac{n}{2} - 1} \\ G'(x) = a_1 + a_3 \cdot x + \cdots + a_{n - 1} \cdot x^{\frac{n}{2} - 1}$

$F(x) = G(x^2) + x \cdot G'(x^2)$

$\omega$的本质是一个复数，且满足$\omega^n = 1$，所以显然$\omega$只能在单位圆上

• $\omega_n^k = \omega_n^{k + a \cdot n}$
• $\omega_n^{k_1} \cdot \omega_n^{k_2}= \omega_n^{k_1 + k_2}$
• $\omega_{d \cdot n}^{d \cdot k} = \omega_n^k$
• $\omega_n^{k + \frac{n}{2}} = - \omega_n^{k}$

$F(\omega_n^k) = G((\omega_n^k)^2) + \omega_n^k \cdot G'((\omega_n^k)^2) \\ = G(\omega_n^{2k}) + \omega_n^k \cdot G'(\omega_n^{2k}) \\ = G(\omega_{n / 2}^k) + \omega_n^k \cdot G'(\omega_{n / 2}^k)$

$F(\omega_n^{k + n / 2}) = G(\omega_{n / 2}^k) - \omega_n^k \cdot G'(\omega_{n / 2}^k)$

## Code

##### FFT
#include <bits/stdc++.h>

using namespace std;

#define fst first
#define snd second
#define mp make_pair
#define squ(x) ((LL)(x) * (x))
#define debug(...) fprintf(stderr, __VA_ARGS__)

typedef long long LL;
typedef pair<int, int> pii;

template<typename T> inline bool chkmax(T &a, const T &b) { return a < b ? a = b, 1 : 0; }
template<typename T> inline bool chkmin(T &a, const T &b) { return a > b ? a = b, 1 : 0; }

int sum = 0, fg = 1; char c = getchar();
for (; !isdigit(c); c = getchar()) if (c == '-') fg = -1;
for (; isdigit(c); c = getchar()) sum = (sum << 3) + (sum << 1) + (c ^ 0x30);
return fg * sum;
}

namespace FFT {

const int MAX_LEN = (1 << 21) + 5;
const double PI = acos(-1.0);

struct com {
double a, b;
com (double _a = 0.0, double _b = 0.0): a(_a), b(_b) { }
com operator + (const com &t) const { return com(a + t.a, b + t.b); }
com operator - (const com &t) const { return com(a - t.a, b - t.b); }
com operator * (const com &t) const { return com(a * t.a - b * t.b, a * t.b + b * t.a); }
};

int len, cnt, rev[MAX_LEN];
com g[MAX_LEN];

inline void init(int N) {
for (cnt = -1, len = 1; len <= N; len <<= 1) ++cnt;
for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << cnt);
g[0] = com(1.0, 0.0);
for (int i = 1; i <= len; i++) g[i] = com(cos(PI * 2 / len * i), sin(PI * 2 / len * i));
}

inline void DFT(com *x, int op) {
for (int i = 0; i < len; i++) if (i < rev[i]) swap(x[i], x[rev[i]]);
for (int k = 2; k <= len; k <<= 1)
for (int j = 0; j < len; j += k)
for (int i = 0; i < k / 2; i++) {
com X = x[j + i], Y = x[j + i + k / 2] * g[len / k * (~op ? i : k - i)];
x[j + i] = X + Y, x[j + i + k / 2] = X - Y;
}
if (op == -1) for (int i = 0; i < len; i++) x[i].a /= len, x[i].b /= len;
}

inline void mul(int *a, int n, int *b, int m, int *c) {
if (n + m == 0) { c[0] = a[0] * b[0]; return; }
init(n + m);
static com F[MAX_LEN], G[MAX_LEN], S[MAX_LEN];
for (int i = 0; i < len; i++) F[i] = com(i <= n ? a[i] : 0.0, 0.0);
for (int i = 0; i < len; i++) G[i] = com(i <= m ? b[i] : 0.0, 0.0);
DFT(F, 1), DFT(G, 1);
for (int i = 0; i < len; i++) S[i] = F[i] * G[i];
DFT(S, -1);
for (int i = 0; i <= n + m; i++) c[i] = round(S[i].a);
}

}

const int maxn = 2e6 + 10;

int main() {
#ifdef xunzhen
freopen("FFT.in", "r", stdin);
freopen("FFT.out", "w", stdout);
#endif

static int a[maxn], b[maxn], c[maxn];
for (int i = 0; i <= n; i++) a[i] = read();
for (int i = 0; i <= m; i++) b[i] = read();

FFT::mul(a, n, b, m, c);

for (int i = 0; i <= n + m; i++) printf("%d%c", c[i], i < n + m ? ' ' : '\n');

return 0;
}
##### NTT
#include <bits/stdc++.h>

using namespace std;

#define fst first
#define snd second
#define mp make_pair
#define squ(x) ((LL)(x) * (x))
#define debug(...) fprintf(stderr, __VA_ARGS__)

typedef long long LL;
typedef pair<int, int> pii;

template<typename T> inline bool chkmax(T &a, const T &b) { return a < b ? a = b, 1 : 0; }
template<typename T> inline bool chkmin(T &a, const T &b) { return a > b ? a = b, 1 : 0; }

int sum = 0, fg = 1; char c = getchar();
for (; !isdigit(c); c = getchar()) if (c == '-') fg = -1;
for (; isdigit(c); c = getchar()) sum = (sum << 3) + (sum << 1) + (c ^ 0x30);
return fg * sum;
}

namespace NTT {

const int MAX_LEN = 1 << 21;
const int mod = 998244353, g0 = 3;

int len, cnt, rev[MAX_LEN], g[MAX_LEN];

inline int add(int x, int y) { return (x += y) < mod ? (x >= 0 ? x : x + mod) : x - mod; }
inline int mul(int x, int y) { return (LL)x * y % mod; }
inline int Pow(int x, int y) {
if (y < 0) y = -1LL * y * (mod - 2) % (mod - 1);
int res = 1;
for (; y; y >>= 1, x = mul(x, x)) if (y & 1) res = mul(res, x);
return res;
}

void init(int N) {
for (cnt = -1, len = 1; len <= N; len <<= 1) ++cnt;
for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << cnt);
g[0] = 1;
for (int G = Pow(g0, (mod - 1) / len), i = 1; i < len; i++) g[i] = mul(g[i - 1], G);
}

void DFT(int *x, int op) {
for (int i = 0; i < len; i++) if (i < rev[i]) swap(x[i], x[rev[i]]);
for (int k = 2; k <= len; k <<= 1)
for (int j = 0; j < len; j += k)
for (int i = 0; i < k / 2; i++) {
int X = x[j + i], Y = mul(x[j + i + k / 2], g[~op ? len / k * i : len / k * (i ? k - i : i)]);
x[j + i] = add(X, Y), x[j + i + k / 2] = add(X, -Y);
}
if (op == -1) for (int inv = Pow(len, -1), i = 0; i < len; i++) x[i] = mul(x[i], inv);
}

void mul(int *a, int n, int *b, int m, int *c) {
init(n + m);
static int F[MAX_LEN], G[MAX_LEN], S[MAX_LEN];
for (int i = 0; i < len; i++) F[i] = i <= n ? a[i] : 0;
for (int i = 0; i < len; i++) G[i] = i <= m ? b[i] : 0;
DFT(F, 1), DFT(G, 1);
for (int i = 0; i < len; i++) S[i] = mul(F[i], G[i]);
DFT(S, -1);
for (int i = 0; i <= n + m; i++) c[i] = S[i];
}

}

const int maxn = 2e6 + 10;

int main() {
#ifdef xunzhen
freopen("NTT.in", "r", stdin);
freopen("NTT.out", "w", stdout);
#endif

static int a[maxn], b[maxn], c[maxn];
for (int i = 0; i <= n; i++) a[i] = read();
for (int i = 0; i <= m; i++) b[i] = read();

NTT::mul(a, n, b, m, c);

for (int i = 0; i <= n + m; i++) printf("%d%c", c[i], i < n + m ? ' ' : '\n');

return 0;
}

## Summary

NTT可以用来避免浮点数的缓慢运算 但好像取模运算更满(雾

IDFT就先留个坑，等以后再来填算了

##### 19.2.14upd

posted @ 2019-02-03 18:20 xunzhen 阅读(...) 评论(...) 编辑 收藏