文章目录
基于NTT的多项式操作
因为实在是太多东西啦，所以就全部整理了以下，持续更新ing~
前置知识：NTT
顺便说一句，代码采用重载vector封装的形式（因为懒得自己开结构体）
下面是几个已经封装的基础代码，为了方便浏览，先贴出来：
typedef std::vector<int> VI;
int fix(int x) {return (x >> 31 & P) + x;}
int Pow(int x, int k) {
int r = 1;
for(;k; k >>= 1, x = 1LL * x * x % P)
if(k & 1)
r = 1LL * r * x % P;
return r;
}
int Inv(int x) {return Pow(x, P - 2);}
void Pre(int m) {
int x = 0; L = 1;
for(;(L <<= 1) < m; ++x) ;
for(int i = 1;i < L; ++i)
R[i] = R[i >> 1] >> 1 | (i & 1) << x;
int wn = Pow(3, (P - 1) / L); w[0] = 1;
for(int i = 1;i < L; ++i)
w[i] = 1LL * w[i - 1] * wn % P;
InvL = Inv(L);
}
void NTT(int *F) {
for(int i = 0;i < L; ++i)
if(R[i] > i)
std::swap(F[i], F[R[i]]);
for(int i = 1, d = L >> 1; i < L; i <<= 1, d >>= 1)
for(int j = 0;j < L; j += i << 1) {
int *l = F + j, *r = F + i + j, *p = w, tp;
for(int k = i; k--; ++l, ++r, p += d)
tp = 1LL * *p * *r % P, *r = (*l - tp) % P, *l = (*l + tp) % P;
}
}
void Fill(const VI &a, int *A, int m) {
m = std::min(m, (int)a.size());
for(int i = 0;i < m; ++i)
A[i] = a[i];
for(int i = m; i < L; ++i)
A[i] = 0;
}
void Fill(int *A, int *B, int m) {
for(int i = 0;i < m; ++i)
B[i] = A[i];
for(int i = m; i < L; ++i)
B[i] = 0;
}
VI operator * (const VI &a, const VI &b) {
const int Lim = 3000;
int m = a.size() + b.size() - 1;
static VI c; static int A[N], B[N];
c.resize(m, 0);
if(a.size() * b.size() <= Lim) { //小范围暴力
for(int i = 0;i < m; ++i)
c[i] = 0;
for(int i = 0;i < a.size(); ++i)
for(int j = 0;j < b.size(); ++j)
c[i + j] = (c[i + j] + 1LL * a[i] * b[j]) % P;
return c;
}
Pre(m);
Fill(a, A, a.size());
Fill(b, B, b.size());
NTT(A); NTT(B);
for(int i = 0;i < L; ++i)
A[i] = 1LL * A[i] * B[i] % P;
NTT(A);
for(int i = 0;i < m; ++i)
c[i] = fix(1LL * A[L - i & L - 1] * InvL % P);
return c;
}
VI operator - (const VI &a, const VI &b) {
int n = std::max(a.size(), b.size());
static VI c; c.resize(n, 0);
for(int i = 0;i < a.size(); ++i)
c[i] = a[i];
for(int i = 0;i < b.size(); ++i)
c[i] = fix(c[i] - b[i]);
return c;
}
VI operator + (const VI &a, const VI &b) {
int n = std::max(a.size(), b.size());
static VI c; c.resize(n, 0);
for(int i = 0;i < a.size(); ++i)
c[i] = a[i];
for(int i = 0;i < b.size(); ++i)
c[i] = (c[i] + b[i]) % P;
return c;
}
迭代相关：多项式求逆，开根，求exp
这些我们熟悉的多项式操作都是采用迭代法解决，既然都是迭代法，必有其共性，因此下面便脱离传统的证明方法，采用牛顿迭代法来进行一波推导
知识点：泰勒展开
下面首先介绍一下泰勒展开：
对于一个函数
f
(
x
)
f(x)
f(x)，其在
x
=
x
0
x=x_0
x=x0处的泰勒公式如下：
f
(
x
)
=
f
(
x
0
)
0
!
+
f
′
(
x
0
)
1
!
(
x
−
x
0
)
+
f
′
′
(
x
0
)
2
!
(
x
−
x
0
)
2
⋯
f
(
n
)
(
x
0
)
n
!
(
x
−
x
0
)
n
f(x)=\frac{f(x_0)}{0!}+\frac{f'(x_0)}{1!}(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2\cdots\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n
f(x)=0!f(x0)+1!f′(x0)(x−x0)+2!f′′(x0)(x−x0)2⋯n!f(n)(x0)(x−x0)n
知识点：牛顿迭代
这是一种逼近零点的科技。（高中数学课本上有，文化课选手应该了解）
对于一个函数
f
(
x
)
f(x)
f(x)，我们考虑以
x
i
x_i
xi作为起点逼近其零点。
方法就是，作
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi)处的切线，设其交
x
x
x轴于
x
i
+
1
x_{i+1}
xi+1，以
x
i
+
1
x_{i+1}
xi+1作为下一次的起点重复上述过程不断逼近零点。
推一推式子：
y
−
f
(
x
i
)
=
f
′
(
x
i
)
(
x
−
x
i
)
y-f(x_i)=f'(x_i)(x-x_i)
y−f(xi)=f′(xi)(x−xi)
令
y
=
0
y=0
y=0，得到
x
i
+
1
=
x
i
−
f
(
x
i
)
f
′
(
x
i
)
x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)}
xi+1=xi−f′(xi)f(xi)
有没有发现式子其实挺熟悉的？
考虑直线方程：
y
=
f
(
x
i
)
+
f
′
(
x
i
)
(
x
−
x
i
)
y=f(x_i)+f'(x_i)(x-x_i)
y=f(xi)+f′(xi)(x−xi)
发现直线方程就是
f
(
x
)
f(x)
f(x)在
x
i
x_i
xi处的一阶泰勒展开。
形式化地，我们可以这样描述牛顿迭代法：
用函数
f
(
x
)
f(x)
f(x)在
f
(
x
i
)
f(x_i)
f(xi)处的一阶展开方程的零点作为下一次迭代的起点
x
i
x_i
xi，这样的迭代法就是牛顿迭代法。
这个东西对于多项式有什么好处呢？
我们不难发现定义一个以多项式为自变量的函数
G
(
x
)
G(x)
G(x)，泰勒公式仍然成立。
我们首先将题目的形式转化为求
F
(
x
)
∣
G
(
F
(
x
)
)
≡
0
(
m
o
d
  
x
n
)
F(x)|G(F(x))\equiv 0(\mod x^n)
F(x)∣G(F(x))≡0(modxn)
迭代法的一般考虑步骤是：
已知
G
(
F
0
(
x
)
)
≡
0
(
m
o
d
  
x
n
)
G(F_0(x))\equiv 0(\mod x^n)
G(F0(x))≡0(modxn)
求
G
(
F
(
x
)
)
≡
0
(
m
o
d
  
x
2
n
)
G(F(x))\equiv 0 (\mod x^{2n})
G(F(x))≡0(modx2n)
考虑
G
(
x
)
G(x)
G(x)在
F
0
(
x
)
F_0(x)
F0(x)处的泰勒展开逼近
F
(
x
)
F(x)
F(x)
G
(
F
(
x
)
)
≡
G
(
F
0
(
x
)
)
+
G
′
(
F
0
(
x
)
)
(
F
(
x
)
−
F
0
(
x
)
)
+
G
′
′
(
F
0
(
x
)
)
(
F
(
x
)
−
F
0
(
x
)
)
2
⋯
(
m
o
d
  
x
n
)
G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))+G''(F_0(x))(F(x)-F_0(x))^2\cdots(\mod x^n)
G(F(x))≡G(F0(x))+G′(F0(x))(F(x)−F0(x))+G′′(F0(x))(F(x)−F0(x))2⋯(modxn)
事实上，我们有
F
(
x
)
≡
F
0
(
x
)
(
m
o
d
  
x
n
)
F(x)\equiv F_0(x) (\mod x^n)
F(x)≡F0(x)(modxn)
那么
F
(
x
)
−
F
0
(
x
)
≡
0
(
m
o
d
  
x
n
)
⇒
(
F
(
x
)
−
F
0
(
x
)
)
2
≡
0
(
m
o
d
  
x
2
n
)
F(x)-F_0(x)\equiv 0 (\mod x^n)\Rightarrow (F(x)-F_0(x))^2\equiv 0 (\mod x^{2n})
F(x)−F0(x)≡0(modxn)⇒(F(x)−F0(x))2≡0(modx2n)
这一步推导是因为
F
(
x
)
−
F
0
(
x
)
F(x)-F_0(x)
F(x)−F0(x)的最低次幂是
x
n
x^n
xn，平方之后最低次幂就变成了
x
2
n
x^{2n}
x2n
于是我们有
G
(
F
(
x
)
)
≡
G
(
F
0
(
x
)
)
+
G
′
(
F
0
(
x
)
)
(
F
(
x
)
−
F
0
(
x
)
)
(
m
o
d
  
x
2
n
)
G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x)) (\mod x^{2n})
G(F(x))≡G(F0(x))+G′(F0(x))(F(x)−F0(x))(modx2n)
也就是说，
F
(
x
)
F(x)
F(x)为
F
0
(
x
)
F_0(x)
F0(x)牛顿迭代的结果。
根据牛顿迭代的结果，我们可以直接得到：
F
(
x
)
=
F
0
(
x
)
−
G
(
F
0
(
x
)
)
G
′
(
F
0
(
x
)
)
F(x)=F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}
F(x)=F0(x)−G′(F0(x))G(F0(x))
基于这个结论，我们可以采用类似倍增的方法来求得最终需要的
F
(
x
)
)
F(x))
F(x))
多项式求逆
已知
F
(
x
)
F(x)
F(x)，求
R
(
x
)
∣
F
(
x
)
R
(
x
)
≡
1
(
m
o
d
  
x
n
)
R(x)|F(x)R(x)\equiv 1(\mod x^n)
R(x)∣F(x)R(x)≡1(modxn)
原式等价于
F
(
x
)
−
R
−
1
(
x
)
≡
0
(
m
o
d
  
x
n
)
F(x)-R^{-1}(x)\equiv 0(\mod x^n)
F(x)−R−1(x)≡0(modxn)
注意构造的时候一般让
F
(
x
)
F(x)
F(x)在常数位置上，这样求导的时候直接消掉
构造函数
G
(
R
(
x
)
)
=
F
(
x
)
−
R
−
1
(
x
)
G(R(x))=F(x)-R^{-1}(x)
G(R(x))=F(x)−R−1(x)
求导
G
′
(
R
(
x
)
)
=
R
−
2
(
x
)
G'(R(x))=R^{-2}(x)
G′(R(x))=R−2(x)
假设求得
G
(
R
0
(
x
)
)
≡
0
(
m
o
d
  
x
n
)
G(R_0(x))\equiv 0(\mod x^n)
G(R0(x))≡0(modxn)
由刚才的结论
R
(
x
)
≡
R
0
(
x
)
−
G
(
R
0
(
x
)
)
G
′
(
R
0
(
x
)
)
≡
R
0
(
x
)
−
F
(
x
)
−
R
0
−
1
(
x
)
R
0
−
2
(
x
0
)
≡
2
R
0
(
x
)
−
F
(
x
)
R
0
2
(
x
)
(
m
o
d
  
x
2
n
)
R(x)\equiv R_0(x)-\frac{G(R_0(x))}{G'(R_0(x))}\equiv R_0(x)-\frac{F(x)-R_0^{-1}(x)}{R_0^{-2}(x_0)}\equiv 2R_0(x)-F(x)R_0^2(x) (\mod x^{2n})
R(x)≡R0(x)−G′(R0(x))G(R0(x))≡R0(x)−R0−2(x0)F(x)−R0−1(x)≡2R0(x)−F(x)R02(x)(modx2n)
初值直接费马小定理。
结合多项式乘法，结束了!
复杂度由主定理得到
T
(
n
)
=
T
(
n
2
)
+
O
(
n
l
o
g
n
)
=
O
(
n
l
o
g
n
)
T(n)=T(\frac{n}{2})+O(nlogn)=O(nlogn)
T(n)=T(2n)+O(nlogn)=O(nlogn)
VI Inv(const VI &a, const int m) { //多项式求逆
static int A[N], B[N], C[N];
for(int i = 0;i < m; ++i)
A[i] = 0;
A[0] = Inv(a[0]);
int n = 1;
for(;n < m;) {
Pre(n << 2);
Fill(a, B, n << 1);
Fill(A, C, n);
NTT(B); NTT(C);
for(int i = 0;i < L; ++i)
B[i] = 1LL * B[i] * C[i] % P * C[i] % P;
NTT(B);
n <<= 1;
for(int i = 0;i < n; ++i)
A[i] = ((A[i] << 1) - 1LL * B[L - i & L - 1] * InvL) % P;
}
static VI c; c.resize(m);
for(int i = 0;i < m; ++i)
c[i] = fix(A[i]);
return c;
}
多项式开根
已知
F
(
x
)
F(x)
F(x)，求
B
(
x
)
∣
B
(
x
)
≡
F
(
x
)
(
m
o
d
  
x
n
)
B(x)|B(x)\equiv \sqrt {F(x)}(\mod x^n)
B(x)∣B(x)≡F(x)(modxn)
同样地让
F
(
x
)
F(x)
F(x)在常数位上
B
2
(
x
)
−
F
(
x
)
≡
0
(
m
o
d
  
x
n
)
B^2(x)-F(x)\equiv 0(\mod x^n)
B2(x)−F(x)≡0(modxn)
构造函数
G
(
B
(
x
)
)
=
B
2
(
x
)
−
F
(
x
)
G(B(x))=B^2(x)-F(x)
G(B(x))=B2(x)−F(x)
求导
G
′
(
B
(
x
)
)
=
2
B
(
x
)
G'(B(x))=2B(x)
G′(B(x))=2B(x)
迭代
B
(
x
)
=
B
0
(
x
)
−
G
(
B
0
(
x
)
)
G
′
(
B
0
(
x
)
)
=
B
0
(
x
)
−
B
0
2
(
x
)
−
F
(
x
)
2
B
0
(
x
)
=
F
(
x
)
+
B
0
2
(
x
)
2
B
0
(
x
)
B(x)=B_0(x)-\frac{G(B_0(x))}{G'(B_0(x))}=B_0(x)-\frac{B_0^2(x)-F(x)}{2B_0(x)}=\frac{F(x)+B_0^2(x)}{2B_0(x)}
B(x)=B0(x)−G′(B0(x))G(B0(x))=B0(x)−2B0(x)B02(x)−F(x)=2B0(x)F(x)+B02(x)
初值一般为1，否则的话要采用二次剩余或者原根BSGS
结合多项式求逆，多项式乘法可以解决，复杂度同上。
VI Sqrt(const VI &a) { //多项式开根
static VI b, c;
b.resize(1, 1); //仅仅在a0=1的时候成立，否则的话要用BSGS
int n = 1, m = a.size();
for(;n < m;) {
n <<= 1;
c = Inv(b, n);
b = b * b; b.resize(n);
for(int i = 0;i < std::min(n, m); ++i)
b[i] = 1LL * (b[i] + a[i]) * inv2 % P;
b = b * c;
b.resize(n);
}
b.resize(m);
return b;
}
多项式求exp
前置知识：多项式求ln
这个东西简单很多
A
(
x
)
≡
L
n
(
F
(
x
)
)
(
m
o
d
  
x
n
)
A(x)\equiv Ln(F(x))(\mod x^n)
A(x)≡Ln(F(x))(modxn)
两边同时对
x
x
x求导(注意和之前对多项式求导不一样！)
A
′
(
x
)
≡
F
′
(
x
)
F
(
x
)
A'(x)\equiv \frac{F'(x)}{F(x)}
A′(x)≡F(x)F′(x)
结合求逆再积分回去即可，复杂度同。
求导积分就不多说了吧！
VI deri(const VI &a) { //多项式求导
static VI c;
if(a.size() == 1) return VI(1, 0);
c.resize(a.size() - 1);
for(int i = 0;i < a.size() - 1; ++i)
c[i] = 1LL * a[i + 1] * (i + 1) % P;
return c;
}
VI inter(const VI &a) { //多项式积分
static VI c; c.resize(a.size() + 1);
c[0] = 0;
for(int i = 0;i < a.size(); ++i)
c[i + 1] = 1LL * a[i] * Inv(i + 1) % P;
return c;
}
VI Ln(const VI &a, const int m) { //多项式求Ln
static VI c;
c = deri(a) * Inv(a, m - 1);
c.resize(m - 1);
return inter(c);
}
接下来是求exp
已知
F
(
x
)
F(x)
F(x)，求
A
(
x
)
∣
A
(
x
)
≡
e
F
(
x
)
(
m
o
d
  
x
n
)
A(x)|A(x)\equiv e^ {F(x)}(\mod x^n)
A(x)∣A(x)≡eF(x)(modxn)
F
(
x
)
−
l
n
A
(
x
)
≡
(
m
o
d
  
x
n
)
F(x) - lnA(x) \equiv(\mod x^n)
F(x)−lnA(x)≡(modxn)
构造函数
G
(
A
(
x
)
)
=
F
(
x
)
−
l
n
(
A
(
x
)
)
G(A(x))=F(x)-ln(A(x))
G(A(x))=F(x)−ln(A(x))
求导
G
′
(
A
(
x
)
)
=
−
A
−
1
(
x
)
G'(A(x))=-A^{-1}(x)
G′(A(x))=−A−1(x)
迭代
A
(
x
)
=
A
0
(
x
)
−
G
(
A
0
(
x
)
)
G
′
(
A
0
(
x
)
)
=
A
0
(
x
)
+
F
(
x
)
−
l
n
(
A
(
x
)
)
A
0
−
1
(
x
)
=
A
0
(
x
)
(
1
+
F
(
x
)
−
l
n
(
A
(
x
)
)
)
A(x)=A_0(x)-\frac{G(A_0(x))}{G'(A_0(x))}=A_0(x)+\frac{F(x)-ln(A(x))}{A_0^{-1}(x)}=A_0(x)(1+F(x)-ln(A(x)))
A(x)=A0(x)−G′(A0(x))G(A0(x))=A0(x)+A0−1(x)F(x)−ln(A(x))=A0(x)(1+F(x)−ln(A(x)))
初值一般为1，否则无意义。
结合多项式求逆，多项式求Ln，多项式乘法可以解决，复杂度同上。
VI Exp(const VI &a, const int m) { //多项式求exp
static VI f, g; f.resize(1, 1);
int n = 1, sz = a.size();
for(;n < m;) {
n <<= 1; g = Ln(f, n);
for(int i = 0;i < n; ++i)
g[i] = i < sz ? fix(a[i] - g[i]) : fix(-g[i]);
(++g[0]) %= P;
f = f * g; f.resize(n);
}
return f;
}
各种黑科技
多项式快速幂
给定
k
,
A
(
x
)
k,A(x)
k,A(x)，求
B
(
x
)
≡
A
(
x
)
k
(
m
o
d
  
x
n
)
B(x)\equiv A(x)^k(\mod x^n)
B(x)≡A(x)k(modxn)
k
≤
1
0
1
0
5
k \le 10^{10^5}
k≤10105
B
(
x
)
≡
e
k
l
n
(
A
(
x
)
)
(
m
o
d
  
x
n
)
B(x)\equiv e^{kln(A(x))}(\mod x^n)
B(x)≡ekln(A(x))(modxn)
VI Pow(const VI &a, int k, int n) { //多项式快速幂
static VI f; f = Ln(a, n);
for(int i = 0;i < n; ++i)
f[i] = 1LL * f[i] * k % P;
return Exp(f, n);
}
多项式除法
给定
F
(
x
)
,
G
(
x
)
F(x),G(x)
F(x),G(x)，阶数分别为
n
,
m
n,m
n,m，求
F
(
x
)
=
Q
(
x
)
G
(
x
)
+
R
(
x
)
F(x)=Q(x)G(x)+R(x)
F(x)=Q(x)G(x)+R(x)，其中
R
(
x
)
R(x)
R(x)的阶数小于
m
m
m。
F
(
x
)
=
Q
(
x
)
G
(
x
)
+
R
(
x
)
⇒
F
(
1
x
)
=
Q
(
1
x
)
G
(
1
x
)
+
R
(
1
x
)
F(x)=Q(x)G(x)+R(x)\Rightarrow F(\frac{1}{x})=Q(\frac{1}{x})G(\frac{1}{x})+R(\frac{1}{x})
F(x)=Q(x)G(x)+R(x)⇒F(x1)=Q(x1)G(x1)+R(x1)
⇒
x
n
F
(
1
x
)
=
x
n
−
m
Q
(
1
x
)
x
m
G
(
1
x
)
+
x
n
R
(
1
x
)
⇒
F
R
(
x
)
=
Q
R
(
x
)
G
R
(
x
)
+
x
n
−
m
+
1
R
R
(
x
)
\Rightarrow x^nF(\frac{1}{x})=x^{n-m}Q(\frac{1}{x})x^mG(\frac{1}{x})+x^nR(\frac{1}{x})\Rightarrow F_R(x)=Q_R(x)G_R(x)+x^{n-m+1}R_R(x)
⇒xnF(x1)=xn−mQ(x1)xmG(x1)+xnR(x1)⇒FR(x)=QR(x)GR(x)+xn−m+1RR(x)
其中
F
R
(
x
)
F_R(x)
FR(x)表示系数反转后的多项式。
显然
Q
(
x
)
Q(x)
Q(x)的阶数为
n
−
m
n-m
n−m，于是
F
R
(
x
)
≡
Q
R
(
x
)
G
R
(
x
)
(
m
o
d
  
x
n
−
m
+
1
)
F_R(x)\equiv Q_R(x)G_R(x) (\mod x^{n-m+1})
FR(x)≡QR(x)GR(x)(modxn−m+1)
⇒
Q
R
(
x
)
≡
F
R
(
x
)
G
R
−
1
(
m
o
d
  
x
n
−
m
+
1
)
\Rightarrow Q_R(x)\equiv F_R(x)G_R^{-1} (\mod x^{n-m+1})
⇒QR(x)≡FR(x)GR−1(modxn−m+1)
结合多项式求逆，即可求得
Q
(
x
)
Q(x)
Q(x)，然后有
R
(
x
)
=
F
(
x
)
−
Q
(
x
)
G
(
x
)
R(x)=F(x)-Q(x)G(x)
R(x)=F(x)−Q(x)G(x)
VI operator / (const VI &a, const VI &b) { //多项式除法
static VI f, g, q; f = a; g = b;
std::reverse(f.begin(), f.end());
std::reverse(g.begin(), g.end());
int n = f.size(), m = g.size();
q = Inv(g, n - m + 1); f.resize(n - m + 1);
q = q * f; q.resize(n - m + 1);
std::reverse(q.begin(), q.end());
return q;
}
VI operator % (const VI &a, const VI &b) { //多项式取余
if(a.size() < b.size()) return a;
static VI c; c = a - (a / b) * b;
c.resize(b.size() - 1);
return c;
}
线性递推
给定
a
0
,
a
1
⋯
a
k
−
1
,
f
0
,
f
1
⋯
f
k
−
1
a_0,a_1\cdots a_{k-1},f_0,f_1\cdots f_{k-1}
a0,a1⋯ak−1,f0,f1⋯fk−1，求前
k
k
k项为
a
a
a的递推式
F
n
=
∑
i
=
0
k
−
1
f
i
F
n
−
k
+
i
F_n=\sum_{i=0}^{k-1}f_iF_{n-k+i}
Fn=∑i=0k−1fiFn−k+i的第
n
n
n项。
n
≤
1
0
9
n\le 10^9
n≤109
推导比较麻烦，这里直接写结论。
求
G
(
x
)
≡
x
n
(
m
o
d
  
x
k
−
f
k
−
1
x
k
−
1
⋯
−
f
0
)
G(x)\equiv x^n(\mod x^k-f_{k-1}x^{k-1}\cdots -f_0)
G(x)≡xn(modxk−fk−1xk−1⋯−f0)
答案就是
A
n
s
=
∑
i
=
0
k
−
1
G
i
a
i
Ans=\sum_{i=0}^{k-1} G_ia_i
Ans=∑i=0k−1Giai
取模的时候快速幂即可。
int Push(const VI &f, const VI &a, int n) { //线性递推
static VI mod, r, b; mod = f;
std::reverse(mod.begin(), mod.end());
for(int i = 0;i < mod.size(); ++i)
mod[i] *= -1;
mod.push_back(1);
r.resize(1, 1); b.resize(2, 0); b[1] = 1;
for(;n; b = b * b % mod, n >>= 1)
if(n & 1)
r = r * b % mod;
long long ans = 0;
for(int i = 0;i < a.size(); ++i)
ans += 1LL * r[i] * a[i] % P;
return fix(ans % P);
}
多项式多点求值
给定
F
(
x
)
F(x)
F(x)和
a
0
,
a
1
⋯
a
m
a_0,a_1\cdots a_m
a0,a1⋯am，对每个
a
i
a_i
ai求
F
(
a
i
)
F(a_i)
F(ai)。
有一个结论是若
R
(
x
)
≡
F
(
x
)
m
o
d
  
(
x
−
a
i
)
R(x)\equiv F(x)\mod (x-a_i)
R(x)≡F(x)mod(x−ai)，那么
F
(
x
)
=
R
(
x
)
F(x)=R(x)
F(x)=R(x)
这个结论是显然的，因为
F
(
x
)
=
(
x
−
a
i
)
Q
(
x
)
+
R
(
x
)
F(x)=(x-a_i)Q(x)+R(x)
F(x)=(x−ai)Q(x)+R(x)，带进去前面消掉了。
考虑分治，构造多项式
P
0
(
x
)
=
∏
i
=
0
m
2
x
−
a
i
P_0(x)=\prod_{i=0}^{\frac{m}{2}}x-a_i
P0(x)=∏i=02mx−ai,
P
1
(
x
)
=
∏
i
=
0
m
2
+
1
x
−
a
i
P_1(x)=\prod_{i=0}^{\frac{m}{2}+1}x-a_i
P1(x)=∏i=02m+1x−ai
然后用
F
(
x
)
F(x)
F(x)对
P
0
,
P
1
P_0,P_1
P0,P1取模，得到
F
0
,
F
1
F_0,F_1
F0,F1，这样我们得到了两个规模减半的相同子问题。递归即可。
P
P
P要事先预处理。
namespace Value { //多项式多点求值
VI res, p[N];
void Pre(int i, int L, int R, const VI &a) {
if(L == R) {
p[i].resize(2);
p[i][0] = fix(-a[L]);
p[i][1] = 1;
return ;
}
int m = L + R >> 1;
Pre(i << 1, L, m, a); Pre(i << 1 | 1, m + 1, R, a);
p[i] = p[i << 1] * p[i << 1 | 1];
}
void Solve(int i, int L, int R, VI f, const VI &a) {
if(L == R)
return res[L] = f[0], void();
int m = L + R >> 1;
Solve(i << 1, L, m, f % p[i << 1], a);
Solve(i << 1 | 1, m + 1, R, f % p[i << 1 | 1], a);
}
VI Query(const VI &f, const VI &a) {
int m = a.size(); res.resize(m);
Pre(1, 0, m - 1, a);
Solve(1, 0, m - 1, f % p[1], a);
return res;
}
}
多项式快速插值
给定
N
N
N个点
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi)，求
N
−
1
N-1
N−1阶拟合多项式。
先复习一下拉格朗日插值的式子吧！
L
(
x
)
=
∑
i
=
0
n
−
1
y
i
ℓ
(
i
)
L(x)=\sum_{i=0}^{n-1} y_i \ell(i)
L(x)=∑i=0n−1yiℓ(i)
其中插值基函数
ℓ
(
i
)
=
∏
j
=
0
,
j
≠
i
n
−
1
x
−
x
j
x
i
−
x
j
\ell(i)=\prod\limits_{j=0,j\neq i}^{n-1}\frac{x-x_j}{x_i-x_j}
ℓ(i)=j=0,j̸=i∏n−1xi−xjx−xj
我们先来考虑求
g
i
=
∏
j
=
0
,
j
≠
i
n
−
1
x
i
−
x
j
g_i=\prod\limits_{j=0,j\neq i}^{n-1}x_i-x_j
gi=j=0,j̸=i∏n−1xi−xj
令
π
(
x
)
=
∏
i
=
0
n
−
1
(
x
−
x
i
)
\pi(x)=\prod_{i=0}^{n-1}(x-x_i)
π(x)=∏i=0n−1(x−xi)
g
i
=
lim
x
→
x
i
π
(
x
i
)
x
−
x
i
=
π
′
(
x
i
)
g_i=\lim_{x\to x_i}\frac{\pi(x_i)}{x-x_i}=\pi'(x_i)
gi=limx→xix−xiπ(xi)=π′(xi)
结合多点求值可以得到。
接下来同样考虑分治，设
L
l
,
r
(
x
)
L_{l,r}(x)
Ll,r(x)表示第
l
l
l到
r
r
r个点插值得到的结果，
m
=
l
+
r
2
m=\frac{l+r}{2}
m=2l+r
有
L
l
,
r
(
x
)
=
∑
i
=
l
r
y
i
g
i
∏
j
=
l
,
j
≠
i
r
x
−
x
j
L_{l,r}(x)=\sum\limits_{i=l}^r\frac{y_i}{g_i}\prod\limits_{j=l,j\neq i}^rx-x_j
Ll,r(x)=i=l∑rgiyij=l,j̸=i∏rx−xj
=
∑
i
=
l
m
y
i
g
i
∏
j
=
l
,
j
≠
i
m
x
−
x
j
∏
j
=
m
+
1
r
x
−
x
j
+
∑
i
=
m
+
1
r
y
i
g
i
∏
j
=
m
+
1
,
j
≠
i
r
x
−
x
j
∏
j
=
l
m
x
−
x
j
=\sum\limits_{i=l}^m\frac{y_i}{g_i}\prod\limits_{j=l,j\neq i}^mx-x_j\prod\limits_{j=m+1}^rx-x_j+\sum\limits_{i=m+1}^r\frac{y_i}{g_i}\prod\limits_{j=m+1,j\neq i}^rx-x_j\prod\limits_{j=l}^mx-x_j
=i=l∑mgiyij=l,j̸=i∏mx−xjj=m+1∏rx−xj+i=m+1∑rgiyij=m+1,j̸=i∏rx−xjj=l∏mx−xj
=
L
l
,
m
(
x
)
∏
j
=
m
+
1
r
x
−
x
j
+
L
m
+
1
,
r
(
x
)
∏
j
=
l
m
x
−
x
j
=L_{l,m}(x)\prod\limits_{j=m+1}^rx-x_j+L_{m+1,r}(x)\prod\limits_{j=l}^mx-x_j
=Ll,m(x)j=m+1∏rx−xj+Lm+1,r(x)j=l∏mx−xj
同前面的套路，预处理出
∏
j
=
l
r
x
−
x
j
\prod\limits_{j=l}^rx-x_j
j=l∏rx−xj，分治即可。
namespace Lag { //多项式快速插值
VI g;
VI Solve(int i, int L, int R, const VI &y) {
if(L == R) return VI(1, 1LL * y[L] * g[L] % P);
int m = L + R >> 1;
return Solve(i << 1, L, m, y) * Value::p[i << 1 | 1] +
Value::p[i << 1] * Solve(i << 1 | 1, m + 1, R, y);
}
VI Calc(const VI &x, const VI &y) {
int n = x.size();
Value::Pre(1, 0, n - 1, x);
g = deri(Value::p[1]);
Value::res.resize(n);
Value::Solve(1, 0, n - 1, g, x);
for(int i = 0;i < n; ++i)
g[i] = Inv(Value::res[i]);
return Solve(1, 0, n - 1, y);
}
}
代码汇总
luogu模板区改改主程序开O2都能过，但是因为是封装的vector所以常数巨大。
#include<bits/stdc++.h>
const int N = 524288, P = 998244353, inv2 = P + 1 >> 1;
typedef std::vector<int> VI;
int fix(int x) {return (x >> 31 & P) + x;}
int ri() {
char c = getchar(); int x = 0, f = 1; for(;c < '0' || c > '9'; c = getchar()) if(c == '-') f = -1;
for(;c >= '0' && c <= '9'; c = getchar()) x = (x << 1) + (x << 3) - '0' + c; return x * f;
}
int mri() {
char c = getchar(); int x = 0, f = 1; for(;c < '0' || c > '9'; c = getchar()) if(c == '-') f = -1;
for(;c >= '0' && c <= '9'; c = getchar()) x = (1LL * x * 10 - '0' + c) % P; return fix(x * f);
}
int InvL, L, R[N], w[N];
int Pow(int x, int k) {
int r = 1;
for(;k; k >>= 1, x = 1LL * x * x % P)
if(k & 1)
r = 1LL * r * x % P;
return r;
}
int Inv(int x) {return Pow(x, P - 2);}
void Pre(int m) {
int x = 0; L = 1;
for(;(L <<= 1) < m; ++x) ;
for(int i = 1;i < L; ++i)
R[i] = R[i >> 1] >> 1 | (i & 1) << x;
int wn = Pow(3, (P - 1) / L); w[0] = 1;
for(int i = 1;i < L; ++i)
w[i] = 1LL * w[i - 1] * wn % P;
InvL = Inv(L);
}
void NTT(int *F) {
for(int i = 0;i < L; ++i)
if(R[i] > i)
std::swap(F[i], F[R[i]]);
for(int i = 1, d = L >> 1; i < L; i <<= 1, d >>= 1)
for(int j = 0;j < L; j += i << 1) {
int *l = F + j, *r = F + i + j, *p = w, tp;
for(int k = i; k--; ++l, ++r, p += d)
tp = 1LL * *p * *r % P, *r = (*l - tp) % P, *l = (*l + tp) % P;
}
}
void Fill(const VI &a, int *A, int m) {
m = std::min(m, (int)a.size());
for(int i = 0;i < m; ++i)
A[i] = a[i];
for(int i = m; i < L; ++i)
A[i] = 0;
}
void Fill(int *A, int *B, int m) {
for(int i = 0;i < m; ++i)
B[i] = A[i];
for(int i = m; i < L; ++i)
B[i] = 0;
}
VI operator * (const VI &a, const VI &b) {
const int Lim = 3000;
int asz = a.size(), bsz = b.size(), m = asz + bsz - 1;
static VI c; c.resize(m);
if(1LL * asz * bsz <= Lim) {
for(int i = 0;i < m; ++i)
c[i] = 0;
for(int i = 0;i < asz; ++i)
for(int j = 0;j < bsz; ++j)
c[i + j] = (c[i + j] + 1LL * a[i] * b[j]) % P;
return c;
}
Pre(m); static int A[N], B[N];
Fill(a, A, asz); Fill(b, B, bsz);
NTT(A); NTT(B);
for(int i = 0;i < L; ++i)
A[i] = (1LL * A[i] * B[i]) % P;
NTT(A);
for(int i = 0;i < m; ++i)
c[i] = fix(1LL * A[L - i & L - 1] * InvL % P);
return c;
}
VI operator - (const VI &a, const VI &b) {
int n = std::max(a.size(), b.size());
static VI c; c.resize(n, 0);
for(int i = 0;i < a.size(); ++i)
c[i] = a[i];
for(int i = 0;i < b.size(); ++i)
c[i] = fix(c[i] - b[i]);
return c;
}
VI operator + (const VI &a, const VI &b) {
int n = std::max(a.size(), b.size());
static VI c; c.resize(n, 0);
for(int i = 0;i < a.size(); ++i)
c[i] = a[i];
for(int i = 0;i < b.size(); ++i)
c[i] = (c[i] + b[i]) % P;
return c;
}
VI Inv(const VI &a, const int m) { //多项式求逆
static int A[N], B[N], C[N];
for(int i = 0;i < m; ++i)
A[i] = 0;
A[0] = Inv(a[0]);
int n = 1;
for(;n < m;) {
Pre(n << 2);
Fill(a, B, n << 1);
Fill(A, C, n);
NTT(B); NTT(C);
for(int i = 0;i < L; ++i)
B[i] = 1LL * B[i] * C[i] % P * C[i] % P;
NTT(B);
n <<= 1;
for(int i = 0;i < n; ++i)
A[i] = ((A[i] << 1) - 1LL * B[L - i & L - 1] * InvL) % P;
}
static VI c; c.resize(m);
for(int i = 0;i < m; ++i)
c[i] = fix(A[i]);
return c;
}
VI deri(const VI &a) { //多项式求导
static VI c;
if(a.size() == 1) return VI(1, 0);
c.resize(a.size() - 1);
for(int i = 0;i < a.size() - 1; ++i)
c[i] = 1LL * a[i + 1] * (i + 1) % P;
return c;
}
VI inter(const VI &a) { //多项式积分
static VI c; c.resize(a.size() + 1);
c[0] = 0;
for(int i = 0;i < a.size(); ++i)
c[i + 1] = 1LL * a[i] * Inv(i + 1) % P;
return c;
}
VI Ln(const VI &a, const int m) { //多项式求Ln
static VI c;
c = deri(a) * Inv(a, m - 1);
c.resize(m - 1);
return inter(c);
}
VI Sqrt(const VI &a) { //多项式开根
static VI b, c;
b.resize(1);
b[0] = 1; //仅仅在a0=1的时候成立，否则的话要用BSGS
int n = 1, m = a.size();
for(;n < m;) {
n <<= 1;
c = Inv(b, n);
b = b * b; b.resize(n);
for(int i = 0;i < std::min(n, m); ++i)
b[i] = 1LL * (b[i] + a[i]) * inv2 % P;
b = b * c;
b.resize(n);
}
b.resize(m);
return b;
}
VI Exp(const VI &a, const int m) { //多项式求exp
static VI f, g; f.resize(1); f[0] = 1;
int n = 1, sz = a.size();
for(;n < m;) {
n <<= 1; g = Ln(f, n);
for(int i = 0;i < n; ++i)
g[i] = i < sz ? fix(a[i] - g[i]) : fix(-g[i]);
(++g[0]) %= P;
f = f * g; f.resize(n);
}
f.resize(m);
return f;
}
VI Pow(const VI &a, int k, int n) { //多项式快速幂
static VI f; f = Ln(a, n);
for(int i = 0;i < n; ++i)
f[i] = 1LL * f[i] * k % P;
return Exp(f, n);
}
VI operator / (const VI &a, const VI &b) { //多项式除法
static VI f, g, q; f = a; g = b;
std::reverse(f.begin(), f.end());
std::reverse(g.begin(), g.end());
int n = f.size(), m = g.size();
q = Inv(g, n - m + 1); f.resize(n - m + 1);
q = q * f; q.resize(n - m + 1);
std::reverse(q.begin(), q.end());
return q;
}
VI operator % (const VI &a, const VI &b) { //多项式取余
if(a.size() < b.size()) return a;
static VI c; c = a - (a / b) * b;
c.resize(b.size() - 1);
return c;
}
int Push(const VI &f, const VI &a, int n) { //线性递推
static VI mod, r, b; mod = f;
std::reverse(mod.begin(), mod.end());
for(int i = 0;i < mod.size(); ++i)
mod[i] *= -1;
mod.push_back(1);
r.resize(1, 1); b.resize(2, 0); b[1] = 1;
for(;n; b = b * b % mod, n >>= 1)
if(n & 1)
r = r * b % mod;
long long ans = 0;
for(int i = 0;i < a.size(); ++i)
ans += 1LL * r[i] * a[i] % P;
return fix(ans % P);
}
namespace Value { //多项式多点求值
VI res, p[N];
void Pre(int i, int L, int R, const VI &a) {
if(L == R) {
p[i].resize(2);
p[i][0] = fix(-a[L]);
p[i][1] = 1;
return ;
}
int m = L + R >> 1;
Pre(i << 1, L, m, a); Pre(i << 1 | 1, m + 1, R, a);
p[i] = p[i << 1] * p[i << 1 | 1];
}
void Solve(int i, int L, int R, VI f, const VI &a) {
if(L == R)
return res[L] = f[0], void();
int m = L + R >> 1;
Solve(i << 1, L, m, f % p[i << 1], a);
Solve(i << 1 | 1, m + 1, R, f % p[i << 1 | 1], a);
}
VI Query(const VI &f, const VI &a) {
int m = a.size(); res.resize(m);
Pre(1, 0, m - 1, a);
Solve(1, 0, m - 1, f % p[1], a);
return res;
}
}
namespace Lag { //多项式快速插值
VI g;
VI Solve(int i, int L, int R, const VI &y) {
if(L == R) return VI(1, 1LL * y[L] * g[L] % P);
int m = L + R >> 1;
return Solve(i << 1, L, m, y) * Value::p[i << 1 | 1] +
Value::p[i << 1] * Solve(i << 1 | 1, m + 1, R, y);
}
VI Calc(const VI &x, const VI &y) {
int n = x.size();
Value::Pre(1, 0, n - 1, x);
g = deri(Value::p[1]);
Value::res.resize(n);
Value::Solve(1, 0, n - 1, g, x);
for(int i = 0;i < n; ++i)
g[i] = Inv(Value::res[i]);
return Solve(1, 0, n - 1, y);
}
}
int main() {
return 0;
}
扫一扫
291
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
