多项式多点求值与快速插值

考虑给定多项式fff，如何求f(x1),f(x2),f(x3),...,f(xn)f(x_1),f(x_2),f(x_3),...,f_(x_n)f(x1),f(x2),f(x3),...,f(xn)。

构造g=∏i=1n2(x−xi),h=∏i=n2+1n(x−xi)g=\prod\limits_{i=1}^{\frac{n}{2}}(x-x_i),h=\prod\limits_{i=\frac{n}{2}+1}^n(x-x_i)g=i=1∏2n(x−xi),h=i=2n+1∏n(x−xi)，然后显然f(xi)=(f%g)(xi)(i≤n2)f(x_i)=(f\% g)(x_i)(i\le \frac{n}{2})f(xi)=(f%g)(xi)(i≤2n)，另一边同理。
然后就变成两个子问题了，时间复杂度T(n)=O(nlog⁡n)+2∗T(n/2)=O(nlog⁡2n)T(n) = O(n\log n) +2*T(n/2)=O(n \log^2 n)T(n)=O(nlogn)+2∗T(n/2)=O(nlog2n)。

考虑给定nnn个点值(x1,y1),(x2,y2),...,(xn,yn)(x_1,y_1),(x_2,y_2),...,(x_n,y_n)(x1,y1),(x2,y2),...,(xn,yn)，如何快速构造出其对应的多项式，考虑拉格朗日插值公式：
f=∑∏j̸=i(x−xj)∏j̸=i(xi−xj)yif=\sum\frac{\prod_{j \not = i}(x-x_j)}{\prod_{j \not = i}(x_i-x_j)} y_if=∑∏j̸=i(xi−xj)∏j̸=i(x−xj)yi

先考虑如何对于一个xix_ixi，求∏j̸=i(xi−xj)\prod_{j\not=i} (x_i-x_j)∏j̸=i(xi−xj)。

记M(x)=∏j(x−xj)M(x)=\prod_j(x-x_j)M(x)=∏j(x−xj)，等价于求c(x)=M(xi)(xi−xj)c(x)=\frac{M(x_i)}{(x_i-x_j)}c(x)=(xi−xj)M(xi)在x→xix \rightarrow x_ix→xi时的值，由洛必达法则，上下求导得其值等于M′(xi)M'(x_i)M′(xi)。于是分治FFT+多点求值即可，时间复杂度O(nlog⁡2n)O(n \log^2 n)O(nlog2n)。

求出这个值之后相当于还要求∑i∏j̸=i(x−xj)M′(xi)yi\sum_i\frac{\prod_{j \not=i}(x-x_j)}{M'(x_i)}y_i∑iM′(xi)∏j̸=i(x−xj)yi，分治FFT即可，时间复杂度O(nlog⁡2n)O(n \log^2 n)O(nlog2n)。

当然，这个常数巨大的算法在n=105n = 10^5n=105的时候需要跑20s+，于是就废了。

#include <bits/stdc++.h>
using namespace std;const int RLEN=1<<18|1;
inline char nc() {static char ibuf[RLEN],*ib,*ob;(ib==ob) && (ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));return (ib==ob) ? -1 : *ib++;
}
inline int rd() {char ch=nc(); int i=0,f=1;while(!isdigit(ch)) {if(ch=='-')f=-1; ch=nc();}while(isdigit(ch)) {i=(i<<1)+(i<<3)+ch-'0'; ch=nc();}return i*f;
}const int N=2e5+50, mod=998244353;
inline int add(int x,int y) {return (x+y>=mod) ? (x+y-mod) : (x+y);}
inline int dec(int x,int y) {return (x-y<0) ? (x-y+mod) : (x-y);}
inline int mul(int x,int y) {return (long long)x*y%mod;}
inline int power(int a,int b,int rs=1) {for(;b;b>>=1,a=mul(a,a)) if(b&1) rs=mul(rs,a); return rs;}namespace FFT {int w[N*8],pos[N*8],k;int A[N*8],B[N*8];inline void init(int nn) {for(k=1;k<=nn;k<<=1);for(int i=1;i<k;i++) pos[i]=(i&1) ? ((pos[i>>1]>>1)^(k>>1)) : (pos[i>>1]>>1);memset(A,0,sizeof(int)*k);memset(B,0,sizeof(int)*k);}inline void dft(int *a) {for(int i=1;i<k;i++) if(pos[i]>i) swap(a[pos[i]],a[i]);for(int bl=1;bl<k;bl<<=1) {int tl=bl<<1, wn=power(3,(mod-1)/tl);w[0]=1; for(int i=1;i<bl;i++) w[i]=mul(w[i-1],wn);for(int bg=0;bg<k;bg+=tl)for(int j=0;j<bl;j++) {int &t1=a[bg+j], &t2=a[bg+j+bl], t=mul(t2,w[j]);t2=dec(t1,t); t1=add(t1,t);}}}inline void func() {dft(A); dft(B);for(int i=0;i<k;i++) B[i]=mul(A[i],B[i]);dft(B); reverse(B+1,B+k);const int inv=power(k,mod-2);for(int i=0;i<k;i++) B[i]=mul(B[i],inv);}
}
int n,xc[N],yc[N],val[N];
struct poly {vector <int> a;inline int deg() const {return a.size()-1;}poly(int d=0,int t=0) {a.resize(d+1); a[d]=t;}inline int& operator [](const int &b) {return a[b];}inline const int& operator [](const int &b) const {return a[b];}friend inline poly operator *(const poly &a,const poly &b) {FFT::init(a.deg()+b.deg());for(int i=0;i<=a.deg();i++) FFT::A[i]=a[i];for(int i=0;i<=b.deg();i++) FFT::B[i]=b[i];FFT::func(); poly c(a.deg()+b.deg(),0);for(int i=0;i<=c.deg();i++) c[i]=FFT::B[i];return c;}friend inline poly operator *(const poly &a,const int &b) {poly c=a; for(int i=0;i<=c.deg();i++) c.a[i]=mul(c.a[i],b);return c;  }friend inline poly operator -(const poly &a,const poly &b) {poly c(max(a.deg(),b.deg()));for(int i=0;i<=a.deg();i++) c[i]=a[i];for(int i=0;i<=b.deg();i++) c[i]=dec(c[i],b[i]);return c;}  friend inline poly operator +(const poly &a,const poly &b) {poly c(max(a.deg(),b.deg()));for(int i=0;i<=a.deg();i++) c[i]=a[i];for(int i=0;i<=b.deg();i++) c[i]=add(c[i],b[i]);return c;}inline poly calc_inv(poly f,int len) {if(len==1) return poly(0,power(f[0],mod-2));poly f0=calc_inv(f.extend(len>>1),len>>1);return (f0*2-(f0*f0).extend(len-1)*f).extend(len-1);}inline poly calc_inverse(int nn) {FFT::init(nn);return calc_inv(*this,FFT::k).extend(nn);}inline poly rev() {poly c=*this; reverse(c.a.begin(),c.a.end()); return c;}friend inline poly operator %(const poly &a,const poly &b) {if(a.deg()<b.deg()) return a;poly A=a, B=b; int res=a.deg()-b.deg();A=A.rev().extend(res); B=B.rev().extend(res);poly C=B.calc_inverse(res);poly D=(A*C).extend(res).rev(); return (a-b*D).extend(b.deg()-1);}inline void dg() {if(!deg()) a[0]=0;else {for(int i=0;i<deg();i++) a[i]=mul(i+1,a[i+1]);a.pop_back();}}inline poly extend(int nn) {poly c=*this;c.a.resize(nn+1);return c;}inline int getval(int x) {int res=0, w=1;for(int i=0;i<=deg();i++)res=add(res,mul(w,a[i])), w=mul(w,x);return res;}inline void pt() {for(int i=0;i<=deg();i++) cout<<a[i]<<' '; cout<<"\n\n";}
} f[N<<2],g;
inline void build(int k,int l,int r) {if(l==r) {f[k]=poly(1,1);f[k][0]=dec(0,xc[l]);return;} int mid=(l+r)>>1;build(k<<1,l,mid);build(k<<1|1,mid+1,r);f[k]=f[k<<1]*f[k<<1|1];
}
inline void getval(poly F,int k,int l,int r,int* v) {if(l==r) {v[l]=F.getval(l); return;}int mid=(l+r)>>1;getval(F%f[k<<1],k<<1,l,mid,v);getval(F%f[k<<1|1],k<<1|1,mid+1,r,v);
}
inline pair <poly,poly> solve(int l,int r) {if(l==r) {poly a(0,val[l]);poly b(1,1);b[0]=dec(0,xc[l]);return pair <poly,poly> (a,b);} int mid=(l+r)>>1;pair <poly,poly> L=solve(l,mid), R=solve(mid+1,r);return pair <poly,poly> (L.first*R.second+L.second*R.first,L.second*R.second);
}
int main() {n=rd();for(int i=1;i<=n;i++) xc[i]=rd(), yc[i]=rd();build(1,1,n);g=f[1]; g.dg();getval(g,1,1,n,val);for(int i=1;i<=n;i++) val[i]=mul(yc[i],power(val[i],mod-2));g=solve(1,n).first;n=rd();for(int i=1;i<=n;i++) xc[i]=rd();build(1,1,n);getval(g,1,1,n,val);for(int i=1;i<=n;i++) cout<<val[i]<<' ';
}

多项式多点求值与快速插值相关推荐

@总结 - 4@ 多项式的多点求值与快速插值
目录 @0 - 参考资料@ @1 - 多点求值@ @理论推导@ @参考代码@ @例题与应用@ @2 - 快速插值@ @理论推导@ @(不建议参考的)代码@ @例题与应用@(暂无) @0 - 参考资料@ ...
洛谷P5282 【模板】快速阶乘算法（多项式多点求值+MTT）
题面传送门前置芝士 \(MTT\),多项式多点求值题解这题法老当初好像讲过--而且他还说这种题目如果模数已经给定可以直接分段打表艹过去以下是题解我们设 \[F(x)=\prod_{i=0} ...
[2021CCPC 威海G] Shinyruo and KFC (下降幂多项式乘法+下降幂转普通幂+多项式多点求值)
题意给定 n n n 个正整数 a 1 , a 2 , ⋯ , a n a_1,a_2,\cdots,a_n a1,a2,⋯,an,并给定正整数 m m m,对于每个 k ∈ [ 1 , m ...
洛谷P5050 【模板】多项式多点求值
传送门人傻常数大.jpg 因为求逆的时候没清零结果调了几个小时-- 前置芝士多项式除法,多项式求逆什么?你不会?左转你谷模板区,包教包会题解首先我们要知道一个结论\[f(x_0)\equiv ...
多项式的求逆、取模和多点求值学习小记
最近学习了多项式的求逆.取模和多点求值,这些方法能够解决很多多项式问题. 这三个操作是环环相扣的,很有趣,学完后不妨记录一下. 多项式求逆给出一个次数界为 nnn 的多项式 A(x)A(x)A(x) ...
多项式的各类计算(多项式的逆/开根/对数/exp/带余除法/多点求值)
预备知识:FFT/NTT 多项式的逆给定一个多项式 F(x)F(x)F(x),请求出一个多项式 G(x)G(x)G(x),满足 F(x)∗G(x)≡1(modxn)F(x)*G(x) \equiv ...
matlab二元多项式求值,matlab多项式代入求值
Matlab 多项式运算与方程求根 ? Matlab多项式运算无论是在线性代数中,还是信号处理.自动控制等理论中,多项式运算都有着十分重要的地位,因此,MATLAB 为多项式的操作提供了相应的函数库 ...
matlab多项式的求值,多项式求值的MATLAB实现
公茂果老师的课件中,给出了四种多项式求值的算法,下面给出代码示例: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %E-mail: [email protected] ...
nyoj 301递推求值 (矩阵+快速幂)
题目:http://acm.nyist.net/JudgeOnline/problem.php?pid=301 题意:给你一个递推公式:f(n)=a*f(n-2)+b*f(n-1)+c 并且告诉你a, ...

多项式多点求值与快速插值

多项式多点求值与快速插值相关推荐

最新文章

热门文章