FFT
FFT能幹嘛？
多項式乘法
多項式
表示多項式
\(f(x) = a + bx + cx^2 + dx^3 + ex^4\)
\(f(x) = (x_1, y_1)\\(x_2, y_2)\\(x_3, y_3)\\(x_4, y_4)\\(x_5, y_5)\)
多項式乘法
\(f(x) = ax^2 + bx + c\)
\(g(x) = dx^2 + ex + f\)
\(f(x)g(x) = adx^4 + (ae + bd)x^3 + (af + be + cd)x^2 + (bf + ce)x + cf\)
\(f(x) = (x_1, y_1), (x_2, y_2), (x_3, y_3)\)
\(g(x) = (x_1, y_4), (x_2, y_5), (x_3, y_6)\)
\(f(x)g(x) = (x_1, y_1y_4), (x_2, y_2y_5), (x_3, y_3y_6)\)
複數平面
座標表示
複數運算
\((a + bi)+ (c+ di) = (a + c) + (b + d)i \\ (a + bi) - (c + di) = (a - c) + (b - d)i \\ (a + bi)(c + di) = (ac - bd)+ (ad + bc)i \\ \frac{a + bi}{c + di} = \frac{(a + bi)(c - di)}{c^2 + d^2} = \frac{ac + bd}{c^2 + d^2} + \frac{bc - ad}{c^2 + d^2}i\)
很有名的東西
\(e^{i\pi} + 1 = 0\)
\(e = \lim\limits_{n \to \infty}(1 + \frac{1}{n})^n\)
\(e^x = \lim\limits_{n \to \infty}((1 + \frac{1}{n})^x) ^ n\\ = \lim\limits_{n \to \infty}(\binom{x}{0}1^x(\frac{1}{n})^0 + \binom{x}{1}1^{x - 1}(\frac{1}{n})^1 + \binom{x}{2}1^{x - 2}(\frac{1}{n})^2 + \dots + \binom{x}{x}1^{0}(\frac{1}{n})^x)^ n\\ = \lim\limits_{n \to \infty}(1 + \frac{x}{n}) ^ n\)
\(e^{ix} = \lim\limits_{n \to \infty}(1 + \frac{ix}{n})^n \\ = \lim\limits_{n \to \infty}\binom{n}{0}1^n(\frac{ix}{n})^0 + \binom{n}{1}1^{n - 1}(\frac{ix}{n})^1 + \binom{n}{2}1^{n - 2}(\frac{ix}{n})^2 + \binom{n}{3}1^{n - 3}(\frac{ix}{n})^3 + \dots + \binom{n}{n}1^0(\frac{ix}{n})^n \\ = \lim\limits_{n \to \infty} \frac{n!}{n!0!} \cdot 1 + \frac{n!}{(n - 1)!1!} \cdot \frac{ix}{n} + \frac{n!}{(n - 2)!2!} \cdot \frac{-x}{n^2} + \frac{n!}{(n - 3)!3!} \cdot \frac{-ix}{n^3} + \dots + \frac{n!}{0!n!} \times \frac{(ix)^n}{n^n} \\ = \lim\limits_{n \to \infty} \frac{1}{0!} + \frac{ix}{1!} + \frac{-x}{2!} + \frac{-ix}{3!} + \dots + \frac{(ix)^n}{n!} \\ = \lim\limits_{n \to \infty} (\frac{1}{0!} - \frac{x^2}{2!} + \frac{x^4}{4!} \dots) + i(\frac{x}{1!} - \frac{x^3}{3!} + \frac{x^5}{5!} \dots) \\ = cosx + i \cdot sinx\)
泰勒展開式
\(f(x) = \int_0^x f'(x) \cdot dx + f(0) \\ f'(x) = \int_0^x f''(x) \cdot dx + f'(0) \\ f(x) = \int_0^x (\int_0^x f''(x) \cdot dx + f'(0)) \cdot dx + f(0) \\ = \int_0^x (\int_0^x f''(x) \cdot dx) \cdot dx + \int_0^x f'(0) \cdot dx + f(0) \\ = \int_0^x(\int_0^x(\int_0^x f'''(x) \cdot dx) \cdot dx) \cdot dx + f''(0) \cdot \int_0^x(\int_0^x dx) \cdot dx + f'(0) \cdot \int_0^x dx + f(0) \\ = \int_0^x(\int_0^x(\int_0^x f'''(x) \cdot dx) \cdot dx) \cdot dx + f''(0) \cdot \frac{x^2}{2!} + f'(0) \cdot \frac{x}{1!} + f(0) \cdot \frac{x^0}{0!} \\ = \sum\limits_{n = 0}^{\infty}{f^{(n)}(0) \cdot \frac{x^n}{n!}} = \sum\limits_{n = 0}^{\infty}{f^{(n)}(a) \cdot \frac{(x - a)^n}{n!}}\)
\(cos(x) = \sum\limits_{n = 0}^{\infty}{cos^{(n)}(0) \cdot \frac{x^n}{n!}} \\ = cos(0) \cdot \frac{1}{0!} - sin(0) \cdot \frac{x}{1!} - cos(0) \cdot \frac{x^2}{2!} + sin(0) \cdot \frac{x^3}{3!} + cos(0) \cdot \frac{x^4}{4!} \dots \\ = \frac{1}{0!} - \frac{x^2}{2!} + \frac{x^4}{4!} - \frac{x^6}{6!} + \frac{x^8}{8!} \dots\)
\(sin(x) = \sum\limits_{n = 0}^{\infty}{sin^{(n)}(0) \cdot \frac{x^n}{n!}} \\ = sin(0) \cdot \frac{1}{0!} + cos(0) \cdot \frac{x}{1!} - sin(0) \cdot \frac{x^2}{2!} - cos(0) \cdot \frac{x^3}{3!} + sin(0) \cdot \frac{x^4}{4!} \dots \\ = \frac{x}{1!} - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} \dots\)
DFT
\(\begin{bmatrix} x_0^0 & x_0^1 & x_0^2 & \dots & x_0^{n - 1} \\ x_1^0 & x_1^1 & x_1^2 & \dots & x_1^{n - 1} \\ x_2^0 & x_2^1 & x_2^2 & \dots & x_2^{n - 1}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_{n - 1}^0 & x_{n - 1}^1 & x_{n - 1}^2 & \dots & x_{n - 1}^{n - 1}\end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n -1}\end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n - 1}\end{bmatrix}\)
\(\begin{bmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & w_n^1 & w_n^2 & \dots & w_n^{n - 1} \\ 1 & w_n^2 & w_n^4 & \dots & w_n^{2(n - 1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & w_n^{n - 1} & w_n^{2(n - 1)} & \dots & w_n^{(n - 1)^2}\end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n -1}\end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n - 1}\end{bmatrix}\)
\(w_n = e^{\frac{2i\pi}{n}}\)
IDFT
\(XA = Y \\ X^{-1}Y = A\)
\(X^{-1} = \begin{bmatrix} \frac{1}{n} & \frac{1}{n} & \frac{1}{n} & \dots & \frac{1}{n} \\ \frac{1}{n} & \frac{w_n^{-1}}{n} & \frac{w_n^{-2}}{n} & \dots & \frac{w_n^{-(n - 1)}}{n} \\ \frac{1}{n} & \frac{w_n^{-2}}{n} & \frac{w_n^{-4}}{n} & \dots & \frac{w_n^{-2(n - 1)}}{n}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{1}{n} & \frac{w_n^{-(n - 1)}}{n} & \frac{w_n^{-2(n - 1)}}{n} & \dots & \frac{w_n^{-(n - 1)^2}}{n} \end{bmatrix}\)
\(w_n^k = e^{\frac{2ki\pi}{n}} = cos(\frac{2k\pi}{n}) + i \times sin(\frac{2k\pi}{n}) \\ w_n^{-k} = e^{\frac{-2ki\pi}{n}} = cos(\frac{-2k\pi}{n}) + i \times sin(\frac{-2k\pi}{n}) = cos(\frac{2k\pi}{n}) - i \times sin(\frac{2k\pi}{n})\)
證明
\(n = 8, 第2列 \cdot 第3行為例\)
FFT
\(f(x) = f_{even}(x^2) + x \cdot f_{odd}(x^2)\)
砸分治
時間複雜度
\(分治 \ logN \ 層 \\ 每層共 \ N \ 個係數 / x \\ O(logN) \cdot O(N) = O(NlogN)\)
Recursive FFT
#include <bits/stdc++.h>
#define cpx complex<double>
#define pb push_back
using namespace std;
int N = 1 << 19;
const double pi = acos(-1);
vector<cpx> A, B, C;
cpx ei(double x){
return cpx(cos(x), sin(x));
}
void FFT(vector<cpx> &F, int inv){
int n = F.size(), h = n >> 1;
cpx w = ei(2 * inv * pi / n), x = 1;
if(n == 1) return;
vector<cpx> L, R;
for(int i = 0; i < n; i++){
if(i & 1) R.pb(F[i]);
else L.pb(F[i]);
}
FFT(L, inv), FFT(R, inv);
for(int i = 0; i < n; i++, x *= w){
F[i] = L[i % h] + x * R[i % h];
}
}
signed main(){
//input
FFT(A, 1), FFT(B, 1);
for(int i = 0; i < N; i++){
C.pb(A[i] * B[i]);
}
FFT(C, -1);
//output
return 0;
}Iterative FFT
#include <bits/stdc++.h>
#define cpx complex<double>
using namespace std;
int N = 1 << 19;
const double pi = acos(-1);
array<cpx, 1 << 19> A, B, C;
cpx ei(double x){
return cpx(cos(x), sin(x));
}
void FFT(array<cpx, 1 << 19> &F, int inv){
cpx x, w, f;
for(int i = 0, k = 0; i < N; i++, k = 0){
for(int j = 1; j < N; j <<= 1) k <<= 1, k |= !!(i & j);
if(i < k) swap(F[i], F[k]);
}
for(int k = 1; k < N; k <<= 1){
w = ei(inv * pi / k);
for(int j = 0; j < N; j += k << 1){
x = 1;
for(int i = j; i < j + k; i++, x *= w){
f = x * F[i + k];
F[i + k] = F[i] - f;
F[i] += f;
}
}
}
}
signed main(){
//input
FFT(A, 1), FFT(B, 1);
for(int i = 0; i < N; i++){
C[i] = A[i] * B[i];
}
FFT(C, -1);
//output
return 0;
}Recursive : Iterative
NTT
有模的FFT
所以不會有精度問題
原根
\(p \ 是一個質數 \\ x^{p - 1} \equiv 1 \ (mod \ p) \\ g \ 是 \ p \ 的原根 \\ 則 \ 0 \le i \not= j < p - 1,\ g^i \not= g^j\ (mod \ p)\)
\(拿 \ g \ 當 \ x \\ 大小為 \ n ,\ p - 1 = t \cdot n\)
\(\begin{bmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & w_n^1 & w_n^2 & \dots & w_n^{n - 1} \\ 1 & w_n^2 & w_n^4 & \dots & w_n^{2(n - 1)}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & w_n^{n - 1} & w_n^{2(n - 1)} & \dots & w_n^{(n - 1)^2}\end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n -1}\end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n - 1}\end{bmatrix}\)
\(w_n = g^{\frac{p - 1}{n}} = g^t\)
INTT
\(X^{-1} = \begin{bmatrix} \frac{1}{n} & \frac{1}{n} & \frac{1}{n} & \dots & \frac{1}{n} \\ \frac{1}{n} & \frac{w_n^{-1}}{n} & \frac{w_n^{-2}}{n} & \dots & \frac{w_n^{-(n - 1)}}{n} \\ \frac{1}{n} & \frac{w_n^{-2}}{n} & \frac{w_n^{-4}}{n} & \dots & \frac{w_n^{-2(n - 1)}}{n}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{1}{n} & \frac{w_n^{-(n - 1)}}{n} & \frac{w_n^{-2(n - 1)}}{n} & \dots & \frac{w_n^{-(n - 1)^2}}{n} \end{bmatrix}\)
證明
\(X \begin{bmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & w_n^{-1} & w_n^{-2} & \dots & w_n^{-(n - 1)} \\ 1 & w_n^{-2} & w_n^{-4} & \dots & w_n^{-2(n - 1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & w_n^{-(n - 1)} & w_n^{-2(n - 1)} & \dots & w_n^{-(n - 1)^2} \end{bmatrix} = P \\ P_{i, j} = \sum\limits_{k = 0}^{n - 1}{w_n^{ki} \cdot w_n^{-kj}} = \sum\limits_{k = 0}^{n - 1}{g^{tki} \cdot g^{-tkj}} = \sum\limits_{k = 0}^{n - 1}{g^{tk(i - j)}} \\ if(i = j) \ P_{i, j} = \sum\limits_{k = 0}^{n - 1}{g^{tk(0)}} = \sum\limits_{k = 0}^{n - 1}{1} = n \\ else \ (g^{t(i - j)} - 1)P_{i, j} = \sum\limits_{k = 0}^{n - 1}{g^{t(k + 1)(i- j)} - g^{tk(i - j)}} = g^{tn(i - j)} - 1 = g^{(p - 1)(i - j)} - 1 = 1 - 1 = 0 \\ i \not= j \implies g^{t(i - j)} \not= 1 \implies P_{i, j} = 0\)
扣
#include <bits/stdc++.h>
#define int long long
using namespace std;
int N = 1 << 20, iN;
const int g = 3, p = 998244353;
array<int, 1 << 21> A, B, C;
int exp(int x, int k){
int s = 1;
for(int i = 1; i <= k; i <<= 1, x = x * x % p){
if(i & k) s = s * x % p;
}
return s;
}
void NTT(array<int, 1 << 20> &F, int inv){
int w, f;
for(int i = 0, k = 0; i < N; i++, k = 0){
for(int j = 1; j < N; j <<= 1) k <<= 1, k |= !!(i & j);
if(i < k) swap(F[i], F[k]);
}
for(int k = 1; k < N; k <<= 1){
w = exp(g, (p - 1) / (k << 1) * (p - 1 + inv));
for(int j = 0; j < N; j += k << 1){
for(int i = j, x = 1; i < j + k; i++, x = x * w % p){
f = x * F[i + k] % p;
F[i + k] = (F[i] - f + p) % p;
F[i] = (F[i] + f) % p;
}
}
}
}
signed main(){
//input
iN = exp(N, p - 2);
NTT(A, 1), NTT(B, 1);
for(int i = 0; i < N; i++) C[i] = A[i] * B[i] % p;
NTT(C, -1);
for(int i = 0; i < N; i++) C[i] = C[i] * iN % p;
//output
return 0;
}題目
扣要戳中間才會出來
扣
#include <bits/stdc++.h>
#define int long long
#define cpx complex<double>
using namespace std;
int N = 1 << 19;
const double pi = acos(-1);
array<cpx, 1 << 19> A, B, C;
cpx ei(double x){
return cpx(cos(x), sin(x));
}
void FFT(array<cpx, 1 << 19> &F, int inv){
cpx x, w, f;
for(int i = 0, k = 0; i < N; i++, k = 0){
for(int j = 1; j < N; j <<= 1) k <<= 1, k |= !!(i & j);
if(i < k) swap(F[i], F[k]);
}
for(int k = 1; k < N; k <<= 1){
w = ei(inv * pi / k);
for(int j = 0; j < N; j += k << 1){
x = 1;
for(int i = j; i < j + k; i++, x *= w){
f = x * F[i + k];
F[i + k] = F[i] - f;
F[i] += f;
}
}
}
}
int rnd(double x){
double z = (int)x;
if(x - z >= 0.5) return (int)z + 1;
else return (int)z;
}
signed main(){
int k, n, m, a, b;
cin >> k >> n >> m;
for(int i = 0; i < n; i++){
cin >> a;
A[a] += 1;
}
for(int i = 0; i < m; i++){
cin >> b;
B[b] += 1;
}
FFT(A, 1), FFT(B, 1);
for(int i = 0; i < N; i++){
C[i] = A[i] * B[i];
}
FFT(C, -1);
for(int i = 2; i <= k << 1; i++){
cout << rnd(C[i].real() / N) << " ";
}
cout << "\n";
return 0;
}扣
#include <bits/stdc++.h>
#define cmp complex<double>
#define int long long
using namespace std;
const int N = 1 << 19;
const double pi = acos(-1);
array<cmp, 1 << 19> A, B, C;
cmp ei(double d){
return {cos(d), sin(d)};
}
void FFT(array<cmp, 1 << 19> &F, int inv){
cmp w, x, f;
for(int i = 0, k = 0; i < N; i++, k = 0){
for(int j = 1; j < N; j <<= 1) k <<= 1, k |= !!(i & j);
if(i < k) swap(F[i], F[k]);
}
for(int k = 1; k < N; k <<= 1){
w = ei(inv * pi / k);
for(int j = 0; j < N; j += k << 1){
x = 1;
for(int i = j; i < j + k; i++, x *= w){
f = x * F[i + k];
F[i + k] = F[i] - f;
F[i] += f;
}
}
}
}
int rnd(double x){
double z = (int)x;
if(x - z >= 0.5) return (int)z + 1;
else return (int)z;
}
signed main(){
int n;
string S;
cin >> S;
n = S.size();
for(int i = 0; i < n; i++){
A[i] = S[i] ^ '0';
B[n - i - 1] = A[i];
}
FFT(A, 1), FFT(B, 1);
for(int i = 0; i < N; i++) C[i] = A[i] * B[i];
FFT(C, -1);
for(int i = n; i < (n << 1) - 1; i++){
cout << rnd(C[i].real() / (double)N) << " ";
}
cout << "\n";
return 0;
}扣
#include <bits/stdc++.h>
#define cmp complex<double>
#define int long long
using namespace std;
const int N = 1 << 19;
const double pi = acos(-1);
array<cmp, 1 << 19> A, B, C;
cmp ei(double d){
return {cos(d), sin(d)};
}
void FFT(array<cmp, 1 << 19> &F, int inv){
cmp w, x, f;
for(int i = 0, k = 0; i < N; i++, k = 0){
for(int j = 1; j < N; j <<= 1) k <<= 1, k |= !!(i & j);
if(i < k) swap(F[i], F[k]);
}
for(int k = 1; k < N; k <<= 1){
w = ei(inv * pi / k);
for(int j = 0; j < N; j += k << 1){
x = 1;
for(int i = j; i < j + k; i++, x *= w){
f = x * F[i + k];
F[i + k] = F[i] - f;
F[i] += f;
}
}
}
}
int rnd(double x){
double z = (int)x;
if(x - z >= 0.5) return (int)z + 1;
else return (int)z;
}
signed main(){
int n, m;
cin >> n >> m;
for(int i = 0; i < n; i++) cin >> A[i];
for(int i = m - 1; i >= 0; i--) cin >> B[i];
FFT(A, 1), FFT(B, 1);
for(int i = 0; i < N; i++) C[i] = A[i] * B[i];
FFT(C, -1);
for(int i = 0; i < n + m - 1; i++){
cout << rnd(C[i].real() / (double)N) << " ";
}
cout << "\n";
return 0;
}扣
#include <bits/stdc++.h>
#define int long long
#define cpx complex<double>
using namespace std;
const int N = 1 << 19;
const double pi = acos(-1);
array<cpx, 1 << 19> A, B, C, X;
cpx ei(double x){
return {cos(x), sin(x)};
}
void FFT(array<cpx, 1 << 19> &F, int inv){
cpx w, x, f;
for(int i = 0, k = 0; i < N; i++, k = 0){
for(int j = 1; j < N; j <<= 1) k <<= 1, k |= !!(i & j);
if(i < k) swap(F[i], F[k]);
}
for(int k = 1; k < N; k <<= 1){
w = ei(inv * pi / k);
for(int j = 0; j < N; j += k << 1){
x = 1;
for(int i = j; i < j + k; i++, x *= w){
f = x * F[i + k];
F[i + k] = F[i] - f;
F[i] += f;
}
}
}
}
int rnd(double x){
double z = (int)x;
if(x - z >= 0.5) return (int)z + 1;
else return (int)z;
}
signed main(){
int n, pre = 0, zero = 0, cnt = 0;
string S;
cin >> S;
n = S.size();
B[n] = {1, 0};
for(char s : S){
if(s == '1'){
pre++;
zero += cnt * (cnt + 1) / 2;
cnt = 0;
}else cnt++;
A[pre] += 1;
B[n - pre] += 1;
}
zero += cnt * (cnt + 1) / 2;
FFT(A, 1), FFT(B, 1);
for(int i = 0; i < N; i++) C[i] = A[i] * B[i];
FFT(C, -1);
cout << zero << " ";
for(int i = n + 1; i <= 2 * n; i++){
cout << rnd(C[i].real() / (double)N) << " ";
}
cout << "\n";
return 0;
}扣
#include <bits/stdc++.h>
#define int long long
using namespace std;
int N, iN;
const int g1 = 3, g2 = 31, p1 = 1004535809, p2 = 2013265921, pp = p1 * p2;
array<int, 1 << 20> A, B, P, Q, C1, C2, C, rev;
void trans(string &s, array<int, 1 << 20> &S){
int n = s.size();
reverse(s.begin(), s.end());
for(int i = 0; i < n; i += 6){
for(int j = min(i + 6, n) - 1; j >= i; j--){
S[i / 6] *= 10;
S[i / 6] += s[j] ^ '0';
}
}
}
int find(int n){
for(int i = 1; i <= 1 << 20; i <<= 1){
if(i > 2 * n) return i;
}
}
int exp(int x, int k, bool t){
int s = 1;
for(int i = 1; i <= k; i <<= 1){
if(i & k){
if(t) s = s * x % p1;
else s = s * x % p2;
}
if(t) x = x * x % p1;
else x = x * x % p2;
}
return s;
}
int mul(int x, int k){
int s = 0;
for(int i = 1; i <= k; i <<= 1, x = (x << 1) % pp){
if(i & k) s = (s + x) % pp;
}
return s;
}
void NTT(array<int, 1 << 20> &F, int inv, bool t){
int w, f;
for(int i = 0; i < N; i++) if(i < rev[i]) swap(F[i], F[rev[i]]);
for(int k = 1; k < N; k <<= 1){
if(t) w = exp(g1, (p1 - 1) / (k << 1) * (p1 - 1 + inv), 1);
else w = exp(g2, (p2 - 1) / (k << 1) * (p2 - 1 + inv), 0);
for(int j = 0; j < N; j += k << 1){
for(int i = j, x = 1; i < j + k; i++){
if(t){
f = x * F[i + k] % p1;
F[i + k] = (F[i] - f + p1) % p1;
F[i] = (F[i] + f) % p1;
x = x * w % p1;
}else{
f = x * F[i + k] % p2;
F[i + k] = (F[i] - f + p2) % p2;
F[i] = (F[i] + f) % p2;
x = x * w % p2;
}
}
}
}
}
void out(array<int, 1 << 20> &S){
string s;
int car = 0, v;
for(int i = 0; i < N; i++){
S[i] += car;
car = S[i] / 1000000;
S[i] %= 1000000;
}
for(; car; car /= 10) s += ((car % 10) + '0');
reverse(s.begin(), s.end());
for(int i = N - 1; i >= 0; i--){
for(int k = 100000; k; k /= 10){
s += (S[i] / k + '0');
S[i] %= k;
}
}
for(v = 0; s[v] == '0' && v < s.size(); v++);
for(; v < s.size(); v++) cout << s[v];
cout << "\n";
}
void RUN(array<int, 1 << 20> &R, bool t){
if(t) iN = exp(N, p1 - 2, 1);
else iN = exp(N, p2 - 2, 0);
P = A, Q = B;
NTT(P, 1, t), NTT(Q, 1, t);
for(int i = 0; i < N; i++){
if(t) R[i] = P[i] * Q[i] % p1;
else R[i] = P[i] * Q[i] % p2;
}
NTT(R, -1, t);
for(int i = 0; i < N; i++){
if(t) R[i] = R[i] * iN % p1;
else R[i] = R[i] * iN % p2;
}
}
signed main(){
int ip1, ip2, p1p, p2p;
string a, b;
cin >> a >> b;
N = find(max(a.size(), b.size()) / 6 + 1);
trans(a, A);
trans(b, B);
for(int i = 0, k = 0; i < N; i++, k = 0){
for(int j = 1; j < N; j <<= 1) k <<= 1, k |= !!(i & j);
rev[i] = k;
}
ip2 = exp(p2, p1 - 2, 1), ip1 = exp(p1, p2 - 2, 0);
p1p = p1 * ip1, p2p = p2 * ip2;
RUN(C1, 1);
RUN(C2, 0);
for(int i = 0; i < N; i++) C[i] = (mul(C1[i], p2p) + mul(C2[i], p1p)) % pp;
out(C);
return 0;
}FFT
By thanksone
