好像又很久没写博客了
总结成一个小札还是香啊 awa
把 @tly 的博客学了一遍然后翻译成自己看得懂的备忘
# 基础多项式版子
你可能会在下面看到 BASIC_POLY 这么一个命名空间，代码如下。
点击展开/折叠代码
1 |
|
# 多项式多点求值
- 问题1
给定一个最高次项为 $n-1$ 的多项式 $f(x)$，求 $f(a_0),f(a_1)\cdots f(a_{m-1})$。
- 新的解法
之前的做法要用到多项式取模这种大常数计算，而且代码还很复杂……
先定义一种所谓的“差卷积”，记作
我们可以构造出 $g(x)=\frac{1}{1-ax}$ 使得 $h=f\otimes g$ 满足 $[x^0]h(x)=f(a)$。
点击展开/折叠证明
把 $g(x)$ 的闭形式展开，则 $(f\otimes g)$ 相当于
$$(f_0+f_1x+f_2x^2+\cdots)\otimes(1+ax+a^2x^2+\cdots)$$观察其常数项，一定是 $f_ix^i\otimes a^ix^i$，也就是 $\sum a^if_i$。
这样我们可以把一个点的求值问题转化成卷积问题。那怎么扩展到多点求值呢？
关于差卷积，我们还需要一个性质——
根据定义很容易理解，即 $f_ig_jh_kx^{i-j-k}=f_ix^i(g_jh_k)x^{-(j+k)}$。
有了这个性质，就可以考虑用线段树的方式分治。我们希望能在线段树的第 $i$ 个叶子处得到 $x^0$。那么线段树上区间 $[l,r]$ 就维护
要从 $[l,r]$ 递推到子区间 $[l,m]$，我们需要计算：
代入即可证明。于是我们还需要对线段树的每个节点 $[l,r]$ 先预处理出
$T_{l,r}$ 的次数和节点的大小是一致的，因为线段树上每一层节点大小减半，可以直接 $O(n\log n)$ 预处理。
似乎现在就可以解决多点求值了？但还剩了一个非常棘手的问题——$S_{l,r}$ 是一个形式幂级数，我们应该保留多少项才能计算出叶子处的常数项？
实际上我们只需要保留 $S_{l,r}$ 的前 $r-l+1$ 项，归纳证明如下：
- 归纳边界：叶子处只需要保留常数项；
- 考虑 $[l,r]$ 的子节点 $[l,m]$，$S_{l,m}$ 需要保留前 $m-l+1$ 项（最高项为 $m-l$ 次），而 $S_{l,m}$ 的计算方法是：注意到 $T_{m+1,r}$ 最高项是 $r-m-1$ 次的，根据差卷积的计算方式，$S_{l,r}$ 的最高项应为 $r-l$ 次，即只需要保留前 $r-l+1$ 项。
这样的话复杂度就有保证了，每层项数减半，仍然可以 $O(n\log n)$ 解决。
- 源代码1
点击展开/折叠代码
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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59//BASIC_POLY 中主要是NTT这些东西
namespace UPPER_POLY{
typedef vector<int> vint;
//普通卷积
vint polyMul(con(vint)a,con(vint)b){
static int tmp1[N+10],tmp2[N+10];
int na=(int)a.size(),nb=(int)b.size(),nc=na+nb-1;
vint c(nc);
int l=BASIC_POLY::modelize(nc);
for(int i=0;i<na;i++) tmp1[i]=a[i];fill(tmp1+na,tmp1+l,0);
for(int i=0;i<nb;i++) tmp2[i]=b[i];fill(tmp2+nb,tmp2+l,0);
BASIC_POLY::ntt(tmp1,l,1),BASIC_POLY::ntt(tmp2,l,1);
for(int i=0;i<l;i++) tmp1[i]=mul(tmp1[i],tmp2[i]);
BASIC_POLY::ntt(tmp1,l,-1);
for(int i=0;i<nc;i++) c[i]=tmp1[i];
return c;
}
//差卷积，卷积结果保留前 nc项
vint polySubMul(con(vint)a,con(vint)b,con(int)nc){
vint c=a;reverse(c.begin(),c.end());
c=polyMul(c,b);
c.resize(a.size()),reverse(c.begin(),c.end());
c.resize(nc);
return c;
}
vint seg[N];
//预处理出 T[l,r]
void build(int *a,con(int)le,con(int)ri){
if(le==ri){
seg[idx(le,ri)].clear();
seg[idx(le,ri)].push_back(1);
seg[idx(le,ri)].push_back(sub(0,a[le]));
return;
}
int mi=(le+ri)>>1;
build(a,le,mi),build(a,mi+1,ri);
seg[idx(le,ri)]=polyMul(seg[idx(le,mi)],seg[idx(mi+1,ri)]);
}
//p即 S[l,r]
void solve(int *res,con(int)le,con(int)ri,con(vint)p){
if(le==ri){res[le]=p[0];return;}
int mi=(le+ri)>>1;
solve(res,le,mi,polySubMul(p,seg[idx(mi+1,ri)],mi-le+1));
solve(res,mi+1,ri,polySubMul(p,seg[idx(le,mi)],ri-mi));
}
void multiVal(int *f,con(int)n,int *pos,con(int)m,int *r){
static int tmp1[N+10],tmp2[N+10];
build(pos,0,m-1);
int rt=idx(0,m-1);
for(int i=0,ii=min(n,(int)seg[rt].size());i<ii;i++)
tmp1[i]=seg[rt][i];
BASIC_POLY::polyInv(tmp1,tmp2,n);
//T[0,n-1] 求个逆再和 f差卷积得到 S[0,n-1]
vint v1(n),v2(n);
for(int i=0;i<n;i++) v1[i]=f[i],v2[i]=tmp2[i];
solve(r,0,m-1,polySubMul(v1,v2,m));
}
}
# 多项式快速插值
- 问题2
给定 $n$ 组 $(x_i,y_i)$，$n-1$ 次多项式 $f(x)$ 满足 $\forall i,f(x_i)=y_i$，求 $f(x)$。
- 拉格朗日插值
由若干个 $n-1$ 次多项式叠加得到 $f(x)$。具体构造如下：
这样构造的目的是 $g_i(x_j)$ 当且仅当 $i=j$ 时 $g_i(x_j)\neq0$，于是可以将 $g_0(x),g_1(x)\cdots g_{n-1}(x)$ 直接相加得到 $f(x)$。
- 快速插值
考虑对拉格朗日插值进行优化（优化计算方式以加速）。观察 $f(x)$ 的表达式：
先看看怎么计算分式的分母部分，该部分对于固定的 $i$ 来说是一个常数，记为 $k_i$。
根据极限的相关知识，不难证明下面这个式子是成立的：
便于书写，记 $g(x)=\prod\limits_{i=0}^{n-1}(x-x_i)$。分式上下都趋近于 $0$，符合洛必达法则的适用条件——由此可得 $k_i=g’(x_i)$
。
可以先用类似线段树的分治（本质是启发式合并）在 $O(n\log n)$ 的复杂度内算出 $g(x)$，然后 $O(n)$ 多项式求导得到 $g’(x)$，再多点求值就可以得到 $k_i$ 了。
再看 $f(x)$，我们现在可以把它写成
则需要解决后面这个式子。这其实可以用到多点求值的旧方法中的一个技巧——仍然是分治：
分治
记 $h_{l,r}(x)$：
$$h_{l,r}(x)=\sum_{i=l}^r\frac{y_i}{k_i}\prod_{j\in[l,r]}^{j\neq i}(x-x_j)$$
仍然用类似线段树的方法分治，考虑如何从 $h_{l,m}(x)$ 和 $h_{m+1,r}(x)$ 合并到 $h_{l,r}(x)$。
$$ \begin{aligned} S_{l,r}&=\prod_{i=l}^r(x-x_i)\\ h_{l,r}(x)&=h_{l,m}S_{m+1,r}+S_{l,m}h_{m+1,r} \end{aligned} $$
这样就可以 $O(n\log n)$ 求出答案。
总的复杂度就是 $O(n^2+n\log n)$
- 源代码2
点击展开/折叠代码
1 | namespace UPPER_POLY{ |
THE END
Thank for reading!
若陪伴是我所有
那痕迹 藏匿了整个宇宙
——《听风捕梦》By (×28)双笙/封茗囧菌/司南
> Link 听风捕梦-网易云