$FFT$

快速傅里叶变换 ($\text{fast Fourier transform}$), 即利用计算机计算离散傅里叶变换($DFT$)的高效、快速计算方法的统称,简称FFT。

(上面是从百科抄的

FFT能在$\Theta(nlog_n)$的时间内求出一个多项式的点值表达的算法

1
注:假设我们有一个n-1次的多项式,它的点值表达即用n个不同的x代入多项式所得到的n个y,这n对(x,y)唯一确定了该多项式。

然后FFT就是用一些特别的x的取值来快速的得出这n个y。

前置芝士:复数

一个复数可以写成$a+bi$的形式

其中$i$叫虚数单位,$i=\sqrt{-1}$,可以把一个复数想象成一个向量$(a,bi)$

复数相乘的规则:模长相乘,辐角相加。

但是$FFT$只用到了后半句(

如果写成向量那么复数的运算就可以这样写:

1
2
3
4
5
假设我们现在有两个复数 x(a,bi) 和 y(c,di)
x + y = (a + c,(b + d)i)
x - y = (a - c,(b - d)i)
x * y = (a * c - b * d,(a * c + b * d)i)
其中乘法就是写成(a + bi) * (c + di)的形式再展开

$STL$提供了复数$complex$的模板,可以直接用,但是常数不小,建议手写。

复数类的模板:

1
2
3
4
5
6
7
struct complex {
double x, y;
complex(double x = 0, double y = 0) :x(x), y(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); }
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); }
}

然后就是重头戏了:

单位根

看上去非常高大上对不对

由于我们在上面说过,复数可以看做一个二维平面的向量,因此横轴表示实部,纵轴表示虚部。

那么我们在这上面画一个单位圆

从点$(1,0)$开始,逆时针将这$n$个点从$0$开始编号,第$k$个点记为$\omega^k_n$,由复数辐角相加的运算法则可知,$\omega^k_n$是$\omega^1_n$的$k$次方,因此$\omega_n^1$被称为$n$次单位根。

显然$\omega^k_n$所对应的向量为$(cos\dfrac {k2\pi} n,sin\dfrac {k2\pi} n)$

单位根的性质:

首先由定义可知:

  • $\omega^{2k}_{2n}=\omega^k_n$

然后由三角函数的性质可知

  • $\omega^{k+\frac n2}=-\omega^k_n$

由单位根的性质可得

把多项式A(x)的离散傅里叶变换结果作为另一个多项式B(x)的系数,取单位根的倒数即$\omega^0_n$,$ω^{−1}_n$,$ω^{−2}_n$,…,$ω^{−(n−1)}_n$作为x代入B(x),得到的每个数再除以n,得到的是A(x)的各项系数。这实现了傅里叶变换的逆变换——把点值表示转换成多项式系数表示,这就是离散傅里叶变换神奇的特殊性质。

证明略,想看的可以看这个讲的非常的详细。

FFT是基于分治的原理,利用单位根的奇偶性来使复杂度达到优秀的$\Theta(nlog_n)$,详细证明可以看上面给出的链接,这里不再赘述。

其实FFT并不需要背过原理也不需要会证明,因为过一段时间就忘了,所以只有记住代码就可以了,代码其实理解了原理再稍微看看就刻在$DNA$里了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
const int N = 4e6 + 1;
struct complex {
double x, y;
complex(double x = 0, double y = 0) :x(x), y(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); }
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); }
}a[N], b[N];
const double Pi = acos(-1.0);
int rev[N], lim = 1, lg;
inl void FFT(complex *a, double op) {
for (re i = 0; i < lim; i++)if (i < rev[i])swap(a[i], a[rev[i]]);
for (re mid = 1; mid < lim; mid <<= 1) {
complex W(cos(Pi / mid), op * sin(Pi / mid));
for (re r = mid << 1, j = 0; j < lim; j += r) {
complex w(1, 0);
for (re k = 0; k < mid; k++, w = w * W) {
complex x = a[j + k], y = w * a[j + k + mid];
a[j + k] = x + y;
a[j + k + mid] = x - y;
}
}
}
}

NTT

就是把FFT从复数域拉回了整数域,避免了精度的损失,同时能更方便的处理取模,常数也更小,就是把单位根通过数论的知识换成了原根,具体证明我不会,同时NTT对模数有很苛刻的要求,基本常用的就那么几个

  • 998244353

  • 1004535809

  • 469762049

这些模数必须是满足可以写成$r*2^k+1$的形式的质数

原根g的定义:就是一个在$\text{mod p}$意义下,g的所有次幂正好对应了0到$p-1$内所有的数。

这个性质非常的棒,可以在特定情况下把乘变成g的指数相加

求原根的方法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
inl void ordm(int m) {
re n = m - 1;
for (re i = 2; i*i < n; i++) {
if (n%i == 0) {
fac[++cnt] = i;
while (n%i == 0)n /= i;
}
}
if (n - 1)fac[++cnt] = n;
for (re i = 2; i < m; i++) {
bool flag = true;
for (re j = 1; j <= cnt; j++) {
if (qpow(i, (m - 1) / fac[j], m) == 1) {
flag = false; break;
}
}
if (flag) { G = i; break; }
}
}

NTT就是把FFT单位根换成了原根,这样理解就够用了。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
inl void NTT(int *a, int lim, int op = 1) {
for (re i = 1; i < lim; i++)if (i < rev[i]) swap(a[i], a[rev[i]]);
for (re mid = 1, r = 2; mid < lim; mid <<= 1, r <<= 1) {
re g = qpow(op == 1 ? G : iG, (p - 1) / r);//iG为G在模P意义下的逆元
for (re j = 0; j < lim; j += r) {
for (re k = 0, gn = 1, x, y; k < mid; k++, gn = 1ll * g * gn % p) {
x = a[j + k], y = 1ll * a[j + k + mid] * gn % p;
a[j + k] = (x + y) % p, a[j + k + mid] = (x - y + p) % p;
}
}
}
if (op != 1) {
re inv = qpow(lim, p - 2);
for (re i = 0; i < lim; i++) a[i] = 1ll * a[i] * inv % p;
}
}

例题:

​ 已经学会了多项式乘法,我们就可以开始做题了,多项式乘法的题目最大的特点和难点就是要把题目给的信息抽象成一个卷积的形式,主要难点在化式子。

luogu P3338

luogu P3321

luogu P5488

luogu P3723

$\large upd: 2020.1.9:$鸽子回来补笔记了QAQ

求逆:给定多项式$F(x)$,求出$G(x)$满足$F(x)*G(x)\equiv 1 \ (mod\ x^n)$

注:如无特殊说明,下文证明过程中的$F$均指$F(x)$

求逆运用了经典的倍增思想,假设我们已经求出了$F_0*G\equiv 1\ (mod\ x^{\frac{n}{2}})$

显然$F*G\equiv 1\ (mod\ x^{\frac{n}{2}})$

那么可以得到$F-F_0 \equiv 0\ (mod\ x^{\frac{n}{2}})$

两边平方得 $F^2 -2F_0F+F_0^2\equiv 0\ (mod\ x^n)$

两边同乘以$G(x)$得 $F-2F_0\times G+F_0^2\equiv 0(mod\ x^n)$

移项得 $F\equiv F_0(2G-F_0)(mod\ x^n)$

然后就可以倍增得出F,边界条件显然为$F(0)=inv(G(0))$

求逆复杂度为$\Theta (nlogn)$,证明可通过主定理:$T(n)=T(n/2)+\Theta(n\log n)=\Theta (n \log n)$

$luogu P4238$:模板题

$\large Code:$

1
2
3
4
5
6
7
8
9
10
11
12
inl void getinv(int *f, int *g, int n) {
if (n == 1) return (void)(g[0] = qpow(f[0], p - 2));
getinv(f, g, n >> 1);
re lim = 1, lg = -1, *a = tmp[++cnt];
while (lim <= n) lim <<= 1, lg++;
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
for (re i = 0; i < n; i++) a[i] = f[i];
NTT(a, lim, 1), NTT(g, lim, 1);
for (re i = 0; i < lim; i++) g[i] = (2 - 1ll * a[i] * g[i] % p + p) % p * g[i] % p, a[i] = 0;
NTT(g, lim, -1), cnt--;
for (re i = n; i < lim; i++) g[i] = 0;
}

$ln$:即自然对数,给出函数$G(x)$,我们要求出满足$ln(G(x))\equiv F(x) \ (mod\ x^n)$的$F(x)$

前置芝士:基础求导和积分,百度百科和数学选修$2-2$课本都可以用来入门

符号:$F’(x)$表示函数$F(x)$的一阶导,$\int F(x)$为函数$F(x)$的积分

对于ln来说我们需要知道$ln’(x)=\dfrac 1 x$,多项式函数$F(x)=a_0+a_1x+a_2x^2+…+a_nx^n$的导数为$F’(x)=a_1+2a_2x+3a_3x^2+…+na_nx^{n-1}$,积分即为求导的逆运算$\int F’(x)+C=F(x)$其中C为常数,于是我们发现求导再积分会丢掉常数项,在某些情况下需要特殊处理。

现在我们开始化式子

$lnG\equiv F\ (mod\ x^n)$

复合函数求导法则:$F(G(x))’=F’(G(x))*G’(x)$

证明:设$u=G(x)$,则$F(u)’=\dfrac{dy}{du}\dfrac{du}{dx}=F’(G(x))G’(x)$

两边求导得 $\dfrac{G’}{G}\equiv F’\ (mod\ x^n)$

即$F\equiv \int (G’*G^{-1})\ (mod\ x^n)$

之后一遍求导,一遍求逆,一遍积分即可。

复杂度还是$\Theta(n\log n)$

$\large Code:$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
inl void der(int *a, int *b, int n) {
for (re i = 1; i < n; i++) b[i - 1] = 1ll * a[i] * i % p;
b[n - 1] = 0;
}//求导
inl void inter(int *a, int *b, int n) {
for (re i = 1; i < n; i++) b[i] = 1ll * a[i - 1] * qpow(i, p - 2) % p;
b[0] = 0;
}//积分
inl void getln(int *f, int *g, int n) {
re *a = tmp[++cnt], *b = tmp[++cnt], lim = 1, lg = -1;
while (lim <= n) lim <<= 1, lg++;
der(f, a, n), getinv(f, b, n);
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(a, lim, 1), NTT(b, lim, 1);
for (re i = 0; i < lim; i++) a[i] = 1ll * a[i] * b[i] % p, b[i] = 0;
NTT(a, lim, -1), inter(a, g, n);
for (re i = n; i < lim; i++) g[i] = a[i] = 0;
for (re i = 0; i < n; i++) a[i] = 0;
cnt -= 2;
}

exp:多项式指数函数,即给出$F(x)$,求出$exp(G(x))\equiv F(x)\ (mod\ x^n)$的$F(x)$

前置芝士:泰勒展开,牛顿迭代

泰勒展开是用多项式函数拟合一个复杂的函数

泰勒公式:$f(x)≈\sum\limits_{n=0}^{∞}\dfrac{f^{(n)}(x_0)}{n!}(x-x_0)^n$

其中$f^{(n)}$为$f$的$n$阶导,通常$x_0$取$0$,此时也叫麦克劳林公式

$f(x)≈\sum\limits_{n=0}^{∞}\dfrac{f^{(n)}(0)}{n!}x^n$

牛顿迭代:已知多项式函数$G(x)$,用于求解复合函数$G(F(x))\equiv 0\ (mod\ x^n)$的$F(x)$

同样使用与求逆类似的倍增思想,假设我们现在已经求出了 $G(F_0(x))\equiv 0\ (mod\ x^{\frac{n}{2}})$

将函数 $G(F(x))$在$F_0(x)$处泰勒展开得$G(F(x))=G(F_0(x))+\dfrac{G(F_0(x))’}{1!}(F(x)-F_0(x))+\dfrac{G(F_0(x))’’}{2!}(F(x)-F_0(x))^2+…+\dfrac{G(F_0(x))^{(n)}}{n!}(F(x)-F_0(x))^n$

因为 $F_0$与$F$的前$\dfrac{n}{2}$项相同,所以$(F(x)-F_0(x))^2$的最小非零项的次数大于$n$

即 $G(F(x))\equiv G(F_0(x))+G(F_0(x))’*(F(x)-F_0(x))\ (mod \ x^n)$

因为 $G(F(x))\equiv 0\ (mod\ x^{n})$

所以 $F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G(F_0(x))’}$

然后考虑求$exp$

$e^{G(x)}\equiv F(x)\ (mod\ x^n)$

两边取$ln$得

$ln(F(x))-G(x)\equiv 0\ (mod\ x^n)$

我们不妨把$G(x)$看做常数,则可以设函数$H(F(x))=ln(F(x))-G(x)$

代入牛顿迭代得$F(x)\equiv F_0(x)-\dfrac{ln(F_0(x))-G(x)}{\frac{1}{F_0(x)}}\ (mod\ x^n)$

则$F(x)\equiv F_0(x)*(1-ln(F_0(x))+G(x))\equiv(mod\ x^n)$

递归求解即可,边界条件因为$G(0)=0$,则$F(0)=1$

同样由主定理分析可知,复杂度为$\Theta(n \log n)$

$\large Code:$

1
2
3
4
5
6
7
8
9
10
11
12
13
inl void getexp(int *f, int *g, int n) {
if (n == 1) return (void)(g[0] = 1);
getexp(f, g, n >> 1);
re lim = 1, lg = -1, *a = tmp[++cnt], *b = tmp[++cnt];
while (lim <= n) lim <<= 1, lg++;
for (re i = 0; i < n; i++) b[i] = f[i];
getln(g, a, lim >> 1);
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(g, lim, 1), NTT(a, lim, 1), NTT(b, lim, 1);
for (re i = 0; i < lim; i++) g[i] = 1ll * g[i] * ((1ll - a[i] + b[i] + p) % p) % p, a[i] = b[i] = 0;
NTT(g, lim, -1), cnt -= 2;
for (re i = n; i < lim; i++) g[i] = 0;
}

带余除法:给定一个 $n$ 次多项式 $F(x)$ 和一个 $m$ 次多项式 $G(x)$,求出多项式 $Q(x)$, $R(x)$,满足以下条件:

  • $Q(x)$ 次数为 $n−m$,$R(x)$ 次数小于 $m$
  • $F(x)=Q(x)∗G(x)+R(x)$

显然有 $F(x^{-1})=G(x^{-1})*Q(x^{-1})+R(x^{-1})$

定义$F_R=\sum\limits_{i=0}^nf_{i}x^{n-i}$

两边同乘$x^n$得 $F_R(x)=G_R(x)*Q_R(x)+x^{n-m+1}R_R(x)$

  • 因为 $F(x^{-1})=\sum\limits_{i=0}^nf_ix^{-i}$

  • 所以$x^n*F(x^{-1})=F_R(x)$,其中$n$为$F$的项数

两边模$x^{n-m+1}$得$F_R(x)\equiv G_R(x)*Q_R(x)\ (mod\ x^{n-m+1})$

即 $Q(x)\equiv reverse(F_R(x)*G_R^{-1}(x)) (mod\ x^{n-m+1})$

得到$Q(x)$以后就可以得到 $R(x)\equiv F(x)-Q(x)*G(x)(mod\ x^{m-1})$

复杂度$\Theta(n\log n)$

$\large Code:$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
//F / G = Q + R
inl void division(int *F, int *G, int *Q, int *R, int n, int m) {
static int tmp1[N], tmp2[N];
reverse(F, F + 1 + n), reverse(G, G + 1 + m);
re lim = 1, lg = -1;
while (lim <= (n << 1)) lim <<= 1, lg++;
for (re i = 0; i <= n; i++) tmp1[i] = F[i];
getinv(G, tmp2, lim >> 1);
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(tmp1, lim, 1), NTT(tmp2, lim, 1);
for (re i = 0; i < lim; i++) Q[i] = 1ll * tmp1[i] * tmp2[i] % p, tmp1[i] = tmp2[i] = 0;
NTT(Q, lim, -1);
for (re i = n - m + 1; i < lim; i++) Q[i] = 0;
reverse(Q, Q + n - m + 1), reverse(F, F + 1 + n), reverse(G, G + 1 + m);
for (re i = 0; i <= n - m; i++) tmp1[i] = Q[i];
for (re i = 0; i <= m; i++) tmp2[i] = G[i];
NTT(tmp1, lim, 1), NTT(tmp2, lim, 1);
for (re i = 0; i < lim; i++) R[i] = 1ll * tmp1[i] * tmp2[i] % p, tmp1[i] = tmp2[i] = 0;
NTT(R, lim, -1);
for (re i = m; i < lim; i++) R[i] = 0;
for (re i = 0; i < m; i++) R[i] = (F[i] - R[i] + p) % p;
}

快速幂:给定一个$n−1$次多项式$A(x)$,求一个在 $\bmod\ x^n$意义下的多项式$B(x)$,使得$B(x) \equiv A^k(x) \ (\bmod\ x^n)$,保证$a_0=1$

  • 做法一:显然有$B(x)=exp(k*ln(A(x)))$,复杂度$\Theta(n\log n)$,常数极大
  • 做法二:像正常快速幂一样倍增,不过乘法换成NTT,复杂度$\Theta(n \log^2n )$

这里只给出做法一的代码

$\large Code:$

1
2
3
getln(f, h, n);
for (re i = 0; i < n; i++)h[i] = (1ll * h[i] * k) % p;
getexp(h, g, n);

当$a_0!=1$时该怎么做呢?当然是把它变成$a_0=1$啦

我们考虑把式子写成$A(x)=(\dfrac{A(x)}{a_0})^k*a_0^k$

这样就完事了…吗?

$a_0$是可以等于$0$的,所以当$a_0$等于$0$时便没法再这样做。

假设$a_t$是第一个非零项,那么我们可以把$A$左移$t$位然后再按上面的做法去处理,只不过做完以后要再右移$min(t*k,n)$位
$\large Code:$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
inl void getpow(int *f, int *g, int n, int k) {
re pos = 0;
if (f[0] == 0) {
for (re i = 1; i < n; i++)
if (f[i]) {
pos = i, n -= pos;
break;
}
}
re *a = tmp[++cnt], *b = tmp[++cnt], innv = qpow(f[pos], p - 2), poi = qpow(f[pos], k), m = 1;
while (m <= n) m <<= 1;
for (re i = 0; i < n; i++) a[i] = 1ll * f[i + pos] * innv % p;
getln(a, b, m);
for (re i = 0; i < n; i++) b[i] = 1ll * b[i] * k % p;
for (re i = n; i < m; i++) b[i] = 0;
getexp(b, g, m);
for (re i = 0; i < n; i++) a[i] = b[i] = 0;
cnt -= 2, n += pos, pos = min(1ll * n, 1ll * pos * k);
for (re i = n - 1; ~i; i--) {
if (i + pos < n) g[i + pos] = g[i];
if (pos) g[i] = 0;
}
for (re i = pos; i < n; i++) g[i] = 1ll * g[i] * poi % p;
}

开根:给定一个$n−1$次多项式$A(x)$,求一个在 $\bmod\ x^n$意义下的多项式$B(x)$,使得$B^2(x) \equiv A(x) \ (\bmod\ x^n)$。若有多解,取零次项系数较小的作为答案,保证$a_0=1$

  • 做法1:直接当成快速幂做,做法与上面相同,取$k$为$2$的逆元即可,复杂度$\Theta(n\log n)$
  • 做法2:化式子

着重讲解一下做法2,因为做法1没有什么好讲的(

$B^2-A\equiv 0(\bmod\ x^n)$

这个式子很像上面牛顿迭代的式子,同样把A看做常数,直接代入。

$B\equiv B_0-\dfrac{B_0^2-A}{2B_0}(\bmod\ x^n)$

即$B\equiv (B_0^2+A)*2B_0^{-1} (\bmod\ x^n)$

倍增即可,复杂度$\Theta(n\log n)$,实测常数大概比做法1小一倍。

$\large Code:$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
inl void getsqrt(int *f, int *g, int n) {
if (n == 1) return (void)(g[0] = 1);
getsqrt(f, g, n >> 1);
re *a = tmp[++cnt], *b = tmp[++cnt], lim = 1, lg = -1;
while (lim <= n) lim <<= 1, lg++;
for (re i = 0; i < n; i++) a[i] = f[i], b[i] = mod(g[i] << 1);
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(a, lim, 1), NTT(g, lim, 1);
for (re i = 0; i < lim; i++) g[i] = (1ll * g[i] * g[i] % p + a[i]) % p, a[i] = 0;
NTT(g, lim, -1);
getinv(b, a, n);
for (re i = n; i < lim; i++) g[i] = 0;
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(a, lim, 1), NTT(g, lim, 1);
for (re i = 0; i < lim; i++) g[i] = 1ll * a[i] * g[i] % p, a[i] = b[i] = 0;
NTT(g, lim, -1), cnt -= 2;
for (re i = n; i < lim; i++) g[i] = 0;
}

然后同理,我们考虑$a_0!=1$但保证$a_0$存在二次剩余时该怎么做。

$\sqrt a_0$在模$P$意义下是什么?我们可以通过求二次剩余找出$\sqrt a_0$,不过显然不能$\Theta(p)$的枚举,我们需要一些高效的方法。

设$p$的原根为$g$,由原根的性质可得,在模$p$意义下任何一个数$a$都可以被唯一的表示为$g^x\equiv a(\bmod\ p)$

那么我们设$\sqrt{a_0}=3^x(\bmod\ 998244353)$

则有$a_0=3^{2*x}(\bmod\ 998244353)$

那么我们可以通过$BSGS$求解出$2x$,之后再一遍快速幂就能得出$\sqrt{a_0}$,复杂度为$\Theta(\sqrt p)$

$\large Code:$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
inl int BSGS(int a, int b) {
unordered_map<int, int> mp;
re t = (int)sqrt(p) + 1;
for (re j = 0; j < t; j++) {
re w = 1ll * b * qpow(a, j) % p;
mp[w] = j;
}
a = qpow(a, t);
if (!a) return b == 0 ? 1 : -1;
for (re i = 0; i <= t; i++) {
re w = qpow(a, i), j = mp.find(w) == mp.end() ? -1 : mp[w];
if (j >= 0 && i * t - j >= 0) return i * t - j;
}
return -1;
}
inl int getkrt(int x, int k) {
re ret = qpow(3, BSGS(3, x) / k);
return min(ret, p - ret);
}
inl void getsqrt(int *f, int *g, int n) {
if (n == 1) return (void)(g[0] = getkrt(f[0], 2));
getsqrt(f, g, n >> 1);
re *a = tmp[++cnt], *b = tmp[++cnt], lim = 1, lg = -1;
while (lim <= n) lim <<= 1, lg++;
for (re i = 0; i < n; i++) a[i] = f[i], b[i] = mod(g[i] << 1);
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(a, lim, 1), NTT(g, lim, 1);
for (re i = 0; i < lim; i++) g[i] = (1ll * g[i] * g[i] % p + a[i]) % p, a[i] = 0;
NTT(g, lim, -1);
getinv(b, a, n);
for (re i = n; i < lim; i++) g[i] = 0;
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(a, lim, 1), NTT(g, lim, 1);
for (re i = 0; i < lim; i++) g[i] = 1ll * a[i] * g[i] % p, a[i] = b[i] = 0;
NTT(g, lim, -1), cnt -= 2;
for (re i = n; i < lim; i++) g[i] = 0;
}

多项式三角函数:给定一个$n-1$次多项式$A(x)$,求一个$\bmod{:x^n}$下的多项式$F(x)$,满足$F(x)\equiv\sin{A(x)}$或$F(x)\equiv\cos{A(x)})$。

对于这种函数图像复杂的函数,我们考虑用泰勒展开去拟合它。

直接套麦克劳林公式能发现一些奇妙的性质:

  • $sin(x)≈\sum\limits_{n=0}^∞(-1)^n\dfrac{x^{2n+1}}{(2n+1)!}$
  • $cos(x)≈\sum\limits_{n=0}^∞(-1)^n\dfrac{x^{2n}}{(2n)!}$
  • $e(x)≈\sum\limits_{n=0}^∞\dfrac{x^n}{n!}$

我们发现这三者有一些奇妙的相似,似乎我们将一个周期性引入$-1$的东西代入$e(x)$就可以找到三者的联系。

众所周知,虚数单位$i$的周期为$4$,$i^{4n}=1,i^{4n+1}=i,i^{4n+2}=-1,i^{4n+3}=-i$

那么我们不妨把$ix$代入$e(x)$,就得到了一个非常有趣的结果:

$e^{ix}≈ \sum\limits^∞_{n=0}\dfrac{(ix)^n}{n!}=(1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+…)+i(\frac{x}{1!}-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+…)$

两部分似乎和$cos$和$sin$的形式非常相似,事实上,他们是完 全 一 致的

$cos(x)≈\sum\limits_{n=0}^∞(-1)^n\dfrac{x^{2n}}{(2n)!}=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+…$

$sin(x)≈\sum\limits_{n=0}^∞(-1)^n\dfrac{x^{2n+1}}{(2n+1)!}=\frac{x}{1!}-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+…$

即$e^{ix}=cos(x)+isin(x)$,这也就是大名鼎鼎的欧拉公式。

同理我们代入$-ix$得到$e^{-ix}=cos(x)-isin(x)$,因为$cos$全为偶次幂,所以符号不变,$sin$全为奇数次幂,所以符号改变。

那么易知 $cos(x)=\dfrac{e^{ix}+e^{-ix}}{2}$,$sin(x)=\dfrac{e^{ix}-e^{ix}}{2i}$

然后通过二次剩余找出$i$在模$p$意义下的值,求$exp$即可,$i$在模$998244353$意义下为$86583718$

$\large Code:$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
inl void getcos(int *f, int *g, int n) {
re *a = tmp[++cnt], *b = tmp[++cnt], *c = tmp[++cnt], innv = 499122177;
for (re i = 0; i < n; i++) c[i] = 1ll * img * f[i] % p;
getexp(c, a, n), getinv(a, b, n);
for (re i = 0; i < n; i++) g[i] = 1ll * innv * (a[i] + b[i]) % p, a[i] = b[i] = c[i] = 0;
cnt -= 3;
}
inl void getsin(int *f, int *g, int n) {
re *a = tmp[++cnt], *b = tmp[++cnt], *c = tmp[++cnt], innv = 954952494;
for (re i = 0; i < n; i++) c[i] = 1ll * img * f[i] % p;
getexp(c, a, n), getinv(a, b, n);
for (re i = 0; i < n; i++) g[i] = 1ll * innv * ((a[i] - b[i]) % p + p) % p, a[i] = b[i] = c[i] = 0;
cnt -= 3;
}

反三角函数:给定一个$n−1$次多项式$A(x)$,求一个 $\bmod{:x^n}$下的多项式$F(x)$,满足$F(x)\equiv\text{asin}:A(x)$或$F(x)\equiv\text{atan}:A(x)$。

首先我们需要了解反函数求导法则:反函数的导数是原函数的导数的倒数(自变量为$y$)

例如 $y=\arcsin x$,显然可以转化为$\sin y=x$

两边求导得 $y’\cos y=1$,即$y’=\frac{1}{\cos y}$

由$\sin^2 x+\cos^2 x=1$得$y’=\dfrac{1}{\sqrt{1-\sin^2y}}=\dfrac{1}{\sqrt{1-x^2}}$

同理可得 $y=\arccos x$的导数为$\dfrac{1}{-\sqrt{1-\cos^2y}}=\dfrac{1}{-\sqrt{1-x^2}}$

即$F(x)\equiv\dfrac{1}{\sqrt{1-A^2(x)}}(\bmod:x^n)$

然后是$\arctan(x)$的导数,要求这个我们要先知道$\tan(x)$的导数。

导数的除法法则:$(\frac{u}{v})’=\dfrac{uv’-u’v}{u^2}$

证明的话套导数定义即可,可以先求出$\frac{1}{v}$的导数,然后再用复合函数求导得出$(\frac{u}{v})’$

知道了这个就可以很方便的求出$tan(x)$的导数。

$tan’(x)=(\dfrac{\sin(x)}{\cos(x)})’=\dfrac{\sin^2(x)+\cos^2(x)}{\cos^2(x)}=\dfrac{1}{\cos^2(x)}=\dfrac{\sin^2(x)}{\cos^2(x)}+1=\tan^2(x)+1$

然后就可以求$\arctan$了,与上面相同,我们设$y=\arctan x$,即$\tan y=x$

两边求导得$y’=\dfrac{1}{1+\tan^2 y}=\dfrac{1}{1+x^2}$

即$F(x)\equiv\dfrac{1}{1+A^2(x)}(\bmod:x^n)$

直接上开根求逆即可,复杂度$\Theta(n\log n)$

$\large Code:$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
inl void arcsin(int *f, int *g, int n) {
re *a = tmp[++cnt], *b = tmp[++cnt], *c = tmp[++cnt], lim = 1, lg = -1;
for (re i = 0; i < n; i++) a[i] = f[i];
while (lim <= n) lim <<= 1, lg++;
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(a, lim, 1);
for (re i = 0; i < lim; i++) a[i] = (1ll - 1ll * a[i] * a[i] % p + p) % p;
NTT(a, lim, -1);
for (re i = n; i < lim; i++) a[i] = 0;
getsqrt(a, b, n), der(f, c, n);
for (re i = 0; i < n; i++) a[i] = 0;
getinv(b, a, n);
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(a, lim, 1), NTT(c, lim, 1);
for (re i = 0; i < lim; i++) b[i] = 1ll * a[i] * c[i] % p, a[i] = c[i] = 0;
NTT(b, lim, -1), inter(b, g, n);
for (re i = 0; i < lim; i++) b[i] = 0;
cnt -= 3;
}
inl void arctan(int *f, int *g, int n) {
re *a = tmp[++cnt], *b = tmp[++cnt], *c = tmp[++cnt], lim = 1, lg = -1;
for (re i = 0; i < n; i++) a[i] = f[i];
while (lim <= n) lim <<= 1, lg++;
for (re i = 1; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg);
NTT(a, lim, 1);
for (re i = 0; i < lim; i++) a[i] = (1ll + 1ll * a[i] * a[i] % p + p) % p;
NTT(a, lim, -1), der(f, c, n);
for (re i = n; i < lim; i++) a[i] = 0;
getinv(a, b, n);
for (re i = 0; i < n; i++) a[i] = 0;
NTT(b, lim, 1), NTT(c, lim, 1);
for (re i = 0; i < lim; i++) b[i] = 1ll * b[i] * c[i] % p, c[i] = 0;
NTT(b, lim, -1), inter(b, g, n);
for (re i = 0; i < lim; i++) b[i] = 0;
cnt -= 3;
}

分治NTT:给定长度为 $n−1$ 的数组 $g[1],g[2],..,g[n−1]$,求 $f[0],f[1],..,f[n−1]$,其中

$f[i]=\sum_{j=1}^if[i-j]g[j]$,边界为 $f[0]=1$,答案模 $998244353$。

除了分治$NTT$之外还可以用多项式求逆$\Theta(n\log n)$解决,先说明一下多项式求逆的做法:

设序列$f$的生成函数$F(x)=[n=0]+\sum\limits_{i=1}^∞ f_ix^i$,设$g[0]=0$,$g$的生成函数$G(x)=\sum\limits_{i=0}^∞ g_ix^i$

那么显然有$F(x)=f_0+F(x)*G(x)$

然后移项可得$F(x)=\dfrac{f_0}{1-G(x)}$,直接计算即可。

既然多项式求逆这么$dio$,为什么我们还要去学复杂度为$\Theta(n\log^2 n)$的分治$NTT$呢?

因为分治$NTT$可以用来优化一堆多项式的暴力卷积,还有一些其他的玄妙用途。

回到这题,对于这题来说,分治$NTT$就是运用了$cdq$分治优化$dp$的思路,每次把区间分为两段,先做左半段,然后用$NTT$计算左半段对右半段的贡献,再算右半段。

左边点对右边的贡献$w_x=\sum\limits_{i=l}^{mid}f_i*g_{x-i}$,这显然是个卷积的形式,直接$NTT$计算即可,然后加到右边的答案数组里。

至此,多项式的模板学习告一段落