多项式多点求值和快速插值学习笔记

前言

常数好大


多点求值

描述

给一个 \(n-1\) 次的多项式 \(f(x)\)\(m\) 个位置 \(x_0,x_1,..,x_{m-1}\),求\(f(x_0),f(x_1),..,f(x_{m-1})\)

暴力

做法

不难发现

\[f(x_0)\equiv f(x)\pmod{(x-x_0)}\]

因为设 \(f(x)=d(x)(x-x_0)+r\),显然有\(f(x_0)=r\)

于是我们可以分治,设

\[ \begin{align} L(x)&=\prod_{i=0}^{\lfloor\frac{m}{2}\rfloor}(x-x_i) \\ R(x)&=\prod_{i=\lfloor\frac{m}{2}\rfloor+1}^{m-1}(x-x_i) \end{align} \]

对于 \(0\le i\le \lfloor\frac{m}{2}\rfloor\),有 \(f(x_i)=(f\ mod\ L)(x_i)\)

对于 \(\lfloor\frac{m}{2}\rfloor<i<m\),有 \(f(x_i)=(f\ mod\ R)(x_i)\)

可以把次数减少一半

求出每个 \(L(x)\)\(R(x)\) 可以用分治FFT,复杂度和后面的取模相同

\(n,m\) 同阶,总复杂度 \(\mathcal O(m\log^2m)\)


插值

描述

\(n\) 个点 \((x_0,y_0),(x_1,y_1),..,(x_{n-1},y_{n-1})\),求一个 \(n-1\) 次多项式 \(f(x)\) 满足 \(f(x_i)=y_i\)

暴力

\(\mathcal O(n^2)\) 拉格朗日插值,这里不展开

快速做法

还是考虑拉格朗日插值

\[ f(x)= \sum_{i=0}^{n-1} \frac {\prod_{j\ne i}(x-x_j)} {\prod_{j\ne i}(x_i-x_j)} y_i \]

先计算每个 \(\prod_{j\ne i}(x_i-x_j)\)

\[ \begin{align} M(x)&=\prod_{i=0}^{n-1} (x-x_i) \\ g_i(x)&=\frac{M(x)}{x-x_i} \end{align} \]

那么

\[\prod_{j\ne i}(x_i-x_j)=g_i(x_i)\]

分子分母都是 \(0\),根据洛必达法则

\[\prod_{j\ne i}(x_i-x_j)=g_i(x_i)=M'(x_i)\]

于是可以使用多点求值在 \(\mathcal O(n\log^2n)\) 的时间内求得

\(z_i=\frac{y_i}{\prod_{j\ne i}(x_i-x_j)}\),问题就是求 \(f(x)=\sum_{i=0}^{n-1}z_i\prod_{j\ne i}(x-x_j)\)

可以分治FFT,同样地设

\[ \begin{align} L(x)&=\prod_{i=0}^{\lfloor\frac{n}{2}\rfloor}(x-x_i) \\ R(x)&=\prod_{i=\lfloor\frac{n}{2}\rfloor+1}^{n-1}(x-x_i) \end{align} \]

于是

\[ f(x)= R(x) \sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}z_i \prod_{0\le j\le \lfloor\frac{n}{2}\rfloor,j\ne i} (x-x_j) + L(x) \sum_{i=\lfloor\frac{n}{2}\rfloor+1}^{n-1} \prod_{\lfloor\frac{n}{2}\rfloor<j<n,j\ne i} (x-x_j) \]

总复杂度还是 \(\mathcal O(n\log^2n)\)


例题

NFLSOJ #53. 多项式多点求值和多点插值

模板

这是一个不完整的多项式板子.

之后会放一个比较全的

在这里挑战多项式

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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
#include<cstdio>
#include<algorithm>
#include<ctype.h>
#include<string.h>
#include<math.h>
#include<vector>

using namespace std;
#define ll long long

inline char read() {
static const int IN_LEN = 1000000;
static char buf[IN_LEN], *s, *t;
return (s==t?t=(s=buf)+fread(buf,1,IN_LEN,stdin),(s==t?-1:*s++):*s++);
}
template<class T>
inline void read(T &x) {
static bool iosig;
static char c;
for (iosig=false, c=read(); !isdigit(c); c=read()) {
if (c == '-') iosig=true;
if (c == -1) return;
}
for (x=0; isdigit(c); c=read()) x=((x+(x<<2))<<1)+(c^'0');
if (iosig) x=-x;
}
const int OUT_LEN = 10000000;
char obuf[OUT_LEN], *ooh=obuf;
inline void print(char c) {
if (ooh==obuf+OUT_LEN) fwrite(obuf, 1, OUT_LEN, stdout), ooh=obuf;
*ooh++=c;
}
template<class T>
inline void print(T x) {
static int buf[30], cnt;
if (x==0) print('0');
else {
if (x<0) print('-'), x=-x;
for (cnt=0; x; x/=10) buf[++cnt]=x%10+48;
while(cnt) print((char)buf[cnt--]);
}
}
inline void flush() { fwrite(obuf, 1, ooh - obuf, stdout); }

int n, m;
vector<int> x, y, z;

namespace Poly{
const int P = 998244353;

vector<int> ans;// for Evaluate()
vector<vector<int>> p;// for Evaluate() & Interpolate()

inline int Pow(ll x, int y=P-2){ // x^y
int ans=1;
for(; y; y>>=1, x=x*x%P) if(y&1) ans=ans*x%P;
return ans;
}
inline int Ge(int x){ int n=1; while(n<=x) n<<=1; return n;}
inline int Mod(int x){ return x<P?x:x-P;}
inline void NTT(vector<int> &f, int g, int n){
f.resize(n);
for(int i=0, j=0; i<n; ++i){
if(i>j) swap(f[i], f[j]);
for(int k=n>>1; (j^=k)<k; k>>=1);
}
vector<int> w(n>>1);
for(int i=1; i<n; i<<=1){
for(int j=w[0]=1, w0=(g==1?Pow(3, (P-1)/i/2):Pow(Pow(3, (P-1)/i/2))); j<i; ++j) w[j]=(ll)w[j-1]*w0%P;
for(int j=0; j<n; j+=i<<1){
for(int k=j; k<j+i; ++k){
int t=(ll)f[k+i]*w[k-j]%P;
f[k+i]=Mod(f[k]-t+P);
f[k]=Mod(f[k]+t);
}
}
}
if(g==-1) for(int i=0, I=Pow(n); i<n; ++i) f[i]=(ll)f[i]*I%P;
}
inline vector<int> Add(const vector<int> &f, const vector<int> &g){
vector<int> ans=f;
for(unsigned i=0; i<f.size(); ++i) (ans[i]+=g[i])%=P;
return ans;
}
inline vector<int> Mul(const vector<int> &f, const vector<int> &g){// f*g
vector<int> F=f, G=g;
int p=Ge(f.size()+g.size()-2);
NTT(F, 1, p), NTT(G, 1, p);
for(int i=0; i<p; ++i) F[i]=(ll)F[i]*G[i]%P;
NTT(F, -1, p);
return F.resize(f.size()+g.size()-1), F;
}
inline vector<int> PolyInv(const vector<int> &f, int n=-1){// 1/f
if(n==-1) n=f.size();
vector<int> ans;
if(n==1) return ans.push_back(Pow(f[0])), ans;
ans=PolyInv(f, (n+1)/2);
vector<int> tmp(&f[0], &f[0]+n);
int p=Ge(n*2-2);
NTT(tmp, 1, p), NTT(ans, 1, p);
for(int i=0; i<p; ++i) ans[i]=(2-(ll)ans[i]*tmp[i]%P+P)*ans[i]%P;
NTT(ans, -1, p);
return ans.resize(n), ans;
}
inline void PolyDiv(const vector<int> &a, const vector<int> &b, vector<int> &d, vector<int> &r){//a=d*b+r
if(b.size()>a.size()) return (void)(d.clear(), r=a);

vector<int> A=a, B=b, iB;
int n=a.size(), m=b.size();
reverse(A.begin(), A.end()), reverse(B.begin(), B.end());
B.resize(n-m+1), iB=PolyInv(B, n-m+1);
d=Mul(A, iB);
d.resize(n-m+1), reverse(d.begin(), d.end());

r=Mul(b, d);
for(int i=0; i<m-1; ++i) r[i]=(P+a[i]-r[i])%P;
r.resize(m-1);
}
inline vector<int> Derivative(const vector<int> &a){ // a'
vector<int> ans;
ans.resize(a.size()-1);
for(unsigned i=1; i<a.size(); ++i) ans[i-1]=(ll)a[i]*i%P;
return ans;
}
void Evaluate_Interpolate_Init(int l, int r, int t, const vector<int> &a){
if(l==r) return p[t].clear(), p[t].push_back(P-a[l]), p[t].push_back(1);
int mid=(l+r)/2, k=t<<1;
Evaluate_Interpolate_Init(l, mid, k, a), Evaluate_Interpolate_Init(mid+1, r, k|1, a);
p[t]=Mul(p[k], p[k|1]);
}
inline void Evaluate(int l, int r, int t, const vector<int> &f, const vector<int> &a){
if(r-l+1<=512){
for(int i=l; i<=r; ++i){
int x=0, j=f.size(), a1=a[i], a2=(ll)a[i]*a[i]%P, a3=(ll)a[i]*a2%P, a4=(ll)a[i]*a3%P, a5=(ll)a[i]*a4%P, a6=(ll)a[i]*a5%P, a7=(ll)a[i]*a6%P, a8=(ll)a[i]*a7%P;
while(j>=8)
x=((ll)x*a8+(ll)f[j-1]*a7+(ll)f[j-2]*a6+(ll)f[j-3]*a5+(ll)f[j-4]*a4+(ll)f[j-5]*a3+(ll)f[j-6]*a2+(ll)f[j-7]*a1+f[j-8])%P, j-=8;
while(j--) x=((ll)x*a[i]+f[j])%P;
ans.push_back(x);
}
return;
}
vector<int> tmp;
PolyDiv(f, p[t], tmp, tmp);
Evaluate(l, (l+r)/2, t<<1, tmp, a), Evaluate((l+r)/2+1, r, t<<1|1, tmp, a);
}
inline vector<int> Evaluate(const vector<int> &f, const vector<int> &a, int flag=-1){// f(a_i)
if(flag==-1) p.resize(a.size()<<2), Evaluate_Interpolate_Init(0, a.size()-1, 1, a);
ans.clear(), Evaluate(0, a.size()-1, 1, f, a);
return ans;
}
vector<int> Interpolate(int l, int r, int t, const vector<int> &x, const vector<int> &f){
if(l==r){
vector<int> ans;
return ans.push_back(f[l]), ans;
}
int mid=(l+r)/2, k=t<<1;
return Add(Mul(Interpolate(l, mid, k, x, f), p[k|1]), Mul(Interpolate(mid+1, r, k|1, x, f), p[k]));
}
inline vector<int> Interpolate(const vector<int> &x, const vector<int> &y){// (x_i,y_i)
int n=x.size();
p.resize(n<<2), Evaluate_Interpolate_Init(0, n-1, 1, x);
vector<int> f=Evaluate(Derivative(p[1]), x, 0);
for(int i=0; i<n; ++i) f[i]=(ll)y[i]*Pow(f[i])%P;
return Interpolate(0, n-1, 1, x, f);
}
}
using namespace Poly;

int main() {
read(n), x.resize(n), y.resize(n);
for(int i=0; i<n; ++i) read(x[i]), read(y[i]);
read(m), z.resize(m);
for(int i=0; i<m; ++i) read(z[i]);

x=Evaluate(Interpolate(x, y), z);
for(int i:x) print(i), print(' ');
return flush(), 0;
}
0%