(color{red}{ ext{多项式全集之二 任长任模的FFT:}})
三模NTT实现任模FFT：
(1.)为什么要用MTT：当(p)不是NTT模数或者多项式长度大于模数限制时，就要使用MTT。
(2.)MTT的使用原理：我们对初始多项式取模，那么如果在不取模卷积情况下，答案(x)不会超过(N imes p^2)。我们取三个NTT模数(p_1,p_2,p_3)，分别做多项式乘法，得到(x)分别(mod~p_1,p_2,p_3)的答案，通过CRT合并可以得到(x~mod~p_1p_2p_3)的答案，如果(x<p_1p_2p_3)那么就可以得到准确的答案，再对(p)取模即可。
(3.)CRT合并的小优化：
(step~0:)初始式子
[{egin{cases}xequiv c_1(mod~p_1)\xequiv c_2(mod~p_2)\xequiv c_3(mod~p_3)end{cases}}
]
(step~1:)把一式二式合并(LL范围内)。
[{egin{cases}xequiv a(mod~p_1p_2)\xequiv c_3(mod~p_3)end{cases}}
]
(step~2:)再次合并(不需要(long~double) 快速乘)。
(4.)常用NTT模数：
以下模数的共同(g=3189)
| (p=r imes 2^k+1) | (k) | (g) |
|---|---|---|
| (104857601) | (22) | (3) |
| (167772161) | (25) | (3) |
| (469762049) | (26) | (3) |
| (950009857) | (21) | (7) |
| (998244353) | (23) | (3) |
| (1004535809) | (21) | (3) |
| (2013265921) | (27) | (31) |
| (2281701377) | (27) | (3) |
| (3221225473) | (30) | (5) |
拆系数FFT(CFFT)实现任模FFT：
(1.)实现原理：运用实数FFT不取模做乘法，然后取模回归到整数。但是由于误差较大（值域是(10^{23})），我们令(t=sqrt{m})把系数(a_i=k_it+b_i)，对(k_i,t_i)交叉做四遍卷积，求出答案按系数贡献取模加入。
(2.)可按合并DFT的方法优化DFT次数。
(bluestein)算法实现任长FFT：
当(m)不是(2)的幂次的时候，我们从式子入手：
令(X_i=a_iw_m^{frac {i^2} 2},Y_i=w_m^{frac{-i^2}2})
(egin{align*} A_k & = sum_{j=0}^{m-1}a_jw_m^{jk}\ & = sum_{j=0}^{m-1}a_jw_m^{frac{j^2+k^2-{(k-j)}^2}{2}}\ &=w_m^{frac {k^2} 2}sum_{j=0}^{m-1}a_jw_m^{frac{j^2} 2}w_m^{frac{-{(k-j)}^2}{2}}\ &=w_m^{frac {k^2} 2}sum_{j=0}^{m-1}X_jY_{k-j} end{align*})
喜闻乐见的模板：
三模NTT模板（注意：不可以MTT回来，因为系数会取模）
namespace MTT{
typedef long long LL;
int n, m;
LL p, mod;
const LL p1 = 998244353;
const LL p2 = 1004535809;
const LL p3 = 104857601;
const int g = 3189;
LL a[300005], b[300005], c[300005], cpa[300005], cpb[300005];
LL c3[300005], c1[300005], c2[300005];
LL qpow(LL a, LL b, LL mod) {
LL ans = 1;
while(b) {
if(b & 1) ans = ans * a % mod;
a = a * a % mod;
b >>= 1;
}
return ans;
}
const LL inv12 = qpow(p1, p2 - 2, p2);
const LL inv123 = qpow(p1 * p2 % p3, p3 - 2, p3);
struct p_l_e{
int wz[300005];
void MTT(LL *a, int N, int op) {
for(int i = 0; i < N; i++)
if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int le = 2; le <= N; le <<= 1) {
int mid = le >> 1;
LL wn = qpow(g, (mod - 1) / le, mod);
if(op == -1) wn = qpow(wn, mod - 2, mod);
for(int i = 0; i < N ;i += le) {
LL w = 1, x, y;
for(int j = 0; j < mid; j++) {
x = a[i + j];
y = a[i + j + mid] * w % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = w * wn % mod;
}
}
}
}
void mult(LL *a, LL *b, LL *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
MTT(a, N, 1); MTT(b, N, 1);
for(int i = 0; i < N; i++) c[i] = a[i] * b[i] % mod;
MTT(c, N, -1);
LL t = qpow(N, mod - 2, mod);
for(int i = 0; i < N; i++) c[i] = c[i] * t % mod;
}
}PLE;
LL CRT(LL c1, LL c2, LL c3) {
LL x = (c1 + p1 * ((c2 - c1 + p2) % p2 * inv12 % p2));
LL y = (x % p + p1 * p2 % p * ((c3 - x % p3 + p3) % p3 * inv123 % p3) % p) % p;
return y;
}
void merge(LL *c1, LL *c2, LL *c3, LL *c, int N) {
for(int i = 0; i < N; i++)
c[i] = CRT(c1[i], c2[i], c3[i]);
return;
}
void main() {
scanf("%d%d%lld", &n, &m, &p); n++; m++;
for(int i = 0; i < n; i++) scanf("%lld", &a[i]);
for(int i = 0; i < m; i++) scanf("%lld", &b[i]);
mod = p1; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c1, n + m - 1);
mod = p2; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c2, n + m - 1);
mod = p3; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c3, n + m - 1);
merge(c1, c2, c3, c, n + m - 1);
for(int i = 0; i < n + m - 1; i++) printf("%lld ", (c[i] % p + p) % p);
return;
}
}
拆系数FFT模板（注意：相同系数的两项可以合并一起IDFT。采用共轭优化法，只进行四次DFT）
namespace CFFT{
typedef long long LL;
int n, m, p ,sqrp;
int a[300005], b[300005];
const long double pi = acos(-1);
struct cp{
long double x, y;
cp() {x = y = 0;}
cp(long double X,long double Y) {x = X; y = Y; }
cp conj() {return (cp) {x, -y};}
}ka[300005], kb[300005], ta[300005], tb[300005], kk[300005], kt[300005], tt[300005], c[300005], I(0, 1), d[300005];
cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; }
cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; }
cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};}
cp operator* (const cp &a, long double b) {return (cp){a.x * b, a.y * b};}
cp operator/ (const cp &a, long double b) {return (cp){a.x / b, a.y / b};}
struct p_l_e{
int wz[300005];
void FFT(cp *a, int N, int op){
for(int i = 0; i < N; i++)
if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int le = 2; le <= N; le <<= 1){
int mid = le >> 1;
cp x, y, w, wn = (cp){cos(op * 2 * pi / le), sin(op * 2 * pi / le)};
for(int i = 0; i < N; i += le){
w = (cp){1, 0};
for(int j = 0; j < mid; j++){
x = a[i + j];
y = a[i + j + mid] * w;
a[i + j] = x + y;
a[i + j + mid] = x - y;
w = w * wn;
}
}
}
}
void D_FFT(cp *a, cp *b, int N, int op){
for(int i = 0; i < N; i++) d[i] = a[i] + I * b[i];
FFT(d, N, op);
d[N] = d[0];
if(op == 1){
for(int i = 0; i < N; i++){
a[i] = (d[i] + d[N - i].conj()) / 2;
b[i] = I * (-1) * (d[i] - d[N - i].conj()) / 2;
}
} else {
for(int i = 0; i < N; i++){
a[i] = cp(d[i].x, 0);
b[i] = cp(d[i].y, 0);
}
}
d[N] = cp(0, 0);
}
void mult(int *a, int *b, int M){
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
for(int i = 0; i < N; i++){
ka[i].x = a[i] >> 15;
kb[i].x = b[i] >> 15;
ta[i].x = a[i] & 32767;
tb[i].x = b[i] & 32767;
}
D_FFT(ta, ka, N, 1); D_FFT(tb, kb, N, 1);
for(int i = 0; i < N; i++){
kk[i] = ka[i] * kb[i];
kt[i] = ka[i] * tb[i] + ta[i] * kb[i];
tt[i] = ta[i] * tb[i];
}
D_FFT(tt, kk, N, -1); FFT(kt, N, -1);
for(int i = 0; i < N; i++){
tt[i] = tt[i] / N;
kt[i] = kt[i] / N;
kk[i] = kk[i] / N;
}
}
}PLE;
void main() {
scanf("%d%d%d", &n, &m, &p); n++; m++;
for(int i = 0; i < n; i++) scanf("%d", &a[i]),a[i] = a[i] % p;
for(int i = 0; i < m; i++) scanf("%d", &b[i]),b[i] = b[i] % p;
PLE.mult(a, b, n + m - 1);
for(int i = 0; i < n + m - 1; i++)
printf("%lld ",(((((LL)round(kk[i].x)) % p) << 30) + ((((LL)round(kt[i].x)) % p) << 15) + ((LL)round(tt[i].x)) % p) % p);
}
}
(blue\_stein)模板：
struct polynie {
CP getw(int m, int k, int op) {
return CP(cos(2 * pi * k / m), op * sin(2 * pi * k / m));
}
int wz[MAXN];
CP A[MAXN], B[MAXN], C[MAXN];
void FFT(CP *a, int N, int op) {
rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int l = 2; l <= N; l <<= 1) {
int mid = l >> 1;
CP x, y, w, wn = CP(cos(pi / mid), sin(op * pi / mid));
for(int i = 0; i < N; i += l) {
w = CP(1, 0);
rop(j, 0, mid) {
x = a[i + j];
y = w * a[i + j + mid];
a[i + j] = x + y;
a[i + j + mid] = x - y;
w = w * wn;
}
}
}
}
void mult(CP *a, CP *b, CP *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
FFT(a, N, 1); FFT(b, N, 1);
rop(i, 0, N) c[i] = a[i] * b[i];
FFT(c, N, -1);
rop(i, 0, N) c[i].x = c[i].x / N, c[i].y = c[i].y / N;
}
void blue_stein(CP *a, int M, int op) {
int M2 = M << 1;
memset(A, 0, sizeof(A));
memset(B, 0, sizeof(B));
memset(C, 0, sizeof(C));
rop(i, 0, M) A[i] = a[i] * getw(M2, 1ll * i * i % M2, op);
rop(i, 0, M2) B[i] = getw(M2, 1ll * (i - M) * (i - M) % M2, -op);
mult(A, B, C, M2 + M - 1);
rop(i, 0, M) a[i] = C[i + M] * getw(M2, 1ll * i * i % M2, op);
if(op == -1) rop(i, 0, M) a[i].x = a[i].x / M, a[i].y = a[i].y / M;
}
}PLE;
(color{red}{ ext{多项式全集之三 多项式求逆与除法:}})
多项式求逆：
(1.)问题描述：
已知(F(x))，且(F(x)G(x)equiv 1 (mod~x^n))，求(G(x))
(2.)推导过程：
egin{align}
B(x)&equiv F(x){-1}&(mod~x{lceil frac n 2
ceil})
由于
G(x)&equiv F(x)^{-1}& (mod~x^n)
所以
G(x)&equiv F(x)^{-1}& (mod~x^{lceil frac n 2
ceil})
G(x)& equiv B(x)&(mod~x^{lceil frac n 2
ceil})
G(x)-B(x)& equiv 0&(mod~x^{lceil frac n 2
ceil})
end{align}
两边平方，得：
由于([G(x)-B(x)]^2)的第(k<n)项为
[sum_{i=0}^k[g_ix^i-b_ix^i][g_{k-i}x^{k-i}-b_{k-i}x^{k-i}]
]
(i,k-i)一定有一项(<frac n 2)，所以
egin{align}
[G(x)-B(x)]^2&equiv 0 &(mod~x^n)
G2(x)+B2(x)-2G(x)B(x)&equiv 0&(mod~x^n)
end{align}
两边同乘(A(x))，得：
egin{align}
A(x)G2(x)+A(x)B2(x)-2A(x)G(x)B(x)&equiv 0& (mod~x^n)
G(x)+A(x)B^2(x)-2B(x)&equiv 0& (mod~x^n)
G(x)&equiv 2B(x)-A(x)B2(x)&(mod~xn)
G(x)&equiv B(x)[2-A(x)B(x)]&(mod~x^n)
end{align}
多项式除法：
(1.)问题描述：
已知一个(n)次多项式(A(x))，一个(m)次多项式(B(x))，且(A(x)=B(x)C(x)+D(x))，求(n-m)次多项式(C(x))，(<m)次多项式(D(x))。
(2.)推导过程：
由(A(x) = B(x)C(x)+D(x))得：
egin{align}
A(frac 1 x)&=B(frac 1 x)C(frac 1 x)+D(frac 1 x)
x^nA(frac 1 x)&=x^nB(frac 1 x)C(frac 1 x)+x^nD(frac 1 x)
x^nA(frac 1 x)&=x^mB(frac 1 x)x^{n-m}C(frac 1 x)+x{m-1}x{n-m+1}D(frac 1 x)
A_r(x)&=B_r(x)C_r(x)+x^{n-m+1}D(x)
A_r(x)&=B_r(x)C_r(x)&(mod ;x^{n-m+1})
B_r(x)&=A_r^{-1}(x)C_r(x)&(mod ;x^{n-m+1})
end{align}
求逆可得(B_r(x))，再反转得(B(x))，然后乘(C(x))去减(A(x))得(D(x)).
喜闻乐见的模板：
namespace INV{
typedef long long LL;
int n, a[300005], b[300005];
const int mod = 998244353;
const int g = 3189;
int qpow(int a, int b){
int ans = 1;
while(b){
if(b & 1) ans = 1ll * ans * a % mod;
a = 1ll * a * a % mod;
b >>= 1;
}
return ans;
}
struct p_l_e{
int wz[300005], i_c[300005];
void NTT(int *a, int N, int op){
for(int i = 0; i < N; i++)
if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int le = 2; le <= N; le <<= 1){
int mid = le >> 1, wn = qpow(g, (mod - 1) / le);
if(op == -1) wn = qpow(wn, mod - 2);
for(int i = 0; i < N; i += le){
LL w = 1; int x, y;
for(int j = 0; j < mid; j++){
x = a[i + j];
y = w * a[i + j + mid] % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = w * wn % mod;
}
}
}
}
int init(int M){
int N = 1, len = 0;
while(N < M) N <<= 1, len++;
for(int i = 0; i < N; i++)
wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
return N;
}
void INV(int *a, int *b, int deg){
if(deg == 1){b[0] = qpow(a[0], mod - 2); return;}
INV(a, b, (deg + 1) >> 1);
int N = init(deg + deg - 1);
for(int i = 0; i < deg; i++) i_c[i] = a[i];
for(int i = deg; i < N; i++) i_c[i] = 0;
NTT(b, N, 1);NTT(i_c, N, 1);
for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * (2 - 1ll * b[i] * i_c[i] % mod + mod) % mod;
NTT(b, N, -1);
int t = qpow(N, mod - 2);
for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod;
for(int i = deg; i < N; i++) b[i] = 0;
}
}PLE;
void main(){
scanf("%d", &n);
for(int i = 0; i < n; i++) scanf("%d", &a[i]);
PLE.INV(a, b, n);
for(int i = 0; i < n; i++) printf("%d ",b[i]);
}
}
(color{red}{ ext{多项式全集之四 多项式ln与exp:}})
多项式ln：
(1.)做法：
设
[G(x) =ln F(x)
]
两边求导得
[G'(x)=frac {F'(x)} {F(x)}
]
积分回去即可。
(2.)应用：
[e^{F(x)}=sum_{k ge 0}frac{F^k(x)}{k!}=G(x)
]
[F(x) = ln G(x)
]
这个的组合意义是：无序组合。
设(F(x))，(f_i)表示一些东西，那么这些东西有序组合的方案数为
[F^0(x) + F^1(x)+F^2(x)+cdots=frac 1 {1 - F(x)}
]
而无序组成的方案数为：
[frac{F^0(x)}{0!} + frac {F^1(x)}{1!}+frac{F^2(x)}{2!}+cdots=e^{F(x)}
]
如果无序组合方案数好求，那么求(ln)就能得到(F(x))。
多项式(exp)：
喜闻乐见的代码：
多项式(ln):
namespace PLE_ln{
struct polyme {
int li[SZ], wz[SZ];
void NTT(int *a, int N, int op) {
rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]);
for(int l = 2; l <= N; l <<= 1) {
int mid = l >> 1;
int x, y, w, wn = qpow(g, (mod - 1) / l);
if(op) wn = qpow(wn, mod - 2);
for(int i = 0; i < N; i += l) {
w = 1;
for(int j = 0; j < mid; ++j) {
x = a[i + j]; y = 1ll * w * a[i + j + mid] % mod;
a[i + j] = (x + y) % mod;
a[i + j + mid] = (x - y + mod) % mod;
w = 1ll * w * wn % mod;
}
}
}
}
void qd(int *a, int *b, int n) {
rop(i, 0, n) b[i] = 1ll * a[i + 1] * (i + 1) % mod;
}
void jf(int *a, int *b, int n) {
rop(i, 1, n) b[i] = 1ll * a[i - 1] * qpow(i, mod - 2) % mod;
}
void mult(int *a, int *b, int *c, int M) {
int N = 1, len = 0;
while(N < M) N <<= 1, len ++;
rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
NTT(a, N, 0); NTT(b, N, 0);
rop(i, 0, N) c[i] = 1ll * a[i] * b[i] % mod;
NTT(c, N, 1);
int t = qpow(N, mod - 2);
rop(i, 0, N) c[i] = 1ll * c[i] * t % mod;
}
void inv(int *a, int *b, int deg) {
if(deg == 1) {b[0] = qpow(a[0], mod - 2) % mod; return;}
inv(a, b, (deg + 1) >> 1);
rop(i, 0, deg) li[i] = a[i];
int N = 1, len = 0;
while(N < deg + deg - 1) N <<= 1, len ++;
rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1));
rop(i, deg, N) li[i] = 0;
NTT(li, N, 0); NTT(b, N, 0);
rop(i, 0, N) b[i] = 1ll * b[i] * (2 - 1ll * li[i] * b[i] % mod + mod) % mod;
NTT(b, N, 1);
int t = qpow(N, mod - 2);
for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod;
rop(i, deg, N) b[i] = 0;
}
}PLE;
int a[SZ], da[SZ], ia[SZ], dla[SZ], la[SZ], n;
void main() {
scanf("%d", &n);
rop(i, 0, n) scanf("%d", &a[i]);
PLE.qd(a, da, n);
PLE.inv(a, ia, n);
PLE.mult(ia, da, dla, n + n - 1);
PLE.jf(dla, la, n);
rop(i, 0, n) printf("%d ", la[i]);
}
}