拉格朗日插值

插值真惨

众所周知$k+1$个点可以确定一个$k$次多项式，那么插值就是通过点值还原多项式的过程。

设给出的$k+1$个点分别是$(x_0,y_0),(x_1,y_1),...,(x_k,y_k)$，那么xjb构造一下：

设函数$f_i(x)=\frac{\prod\limits_{j\neq i}(x-x_j)}{\prod\limits_{j\neq i}(x_i-x_j)}\times y_i$

显然这个函数当$x=x_i$时值为$y_i$，$x=x_j(0\leq j\leq k且j\neq i)$时值为0。

求出所有的$f_i$，然后把它们加起来，发现得到的多项式必过给出的这些点，也就是要求的结果了。。。

即要求的多项式$F(x)=\sum\limits_{i=0}^{k}f_i(x)$

易知时间复杂度是$O(k^2)$的，看起来很暴力，但是拉格朗日插值就是这么暴力。。。

如果要动态加点的话可以稍微优化一下：

注意到每个$f_i$中$\prod\limits_{j\neq i}(x-x_j)$重复计算了，可以简化：

设$p(x)=\prod\limits_{i=0}^{k}(x-x_i)$，$q_i=\frac{y_i}{\prod\limits_{j\neq i}(x_i-x_j)}$

则$F(x)=p(x)\times\sum\limits_{i=0}^{k}\frac{q_i}{(x-x_i)}$

这样子原本的复杂度不变，但是新加点时只用重新算一遍$q_i$，复杂度为$O(k)$

这个东西貌似叫重心法？

一个小技巧：

在实际做题的时候，有时不需要根据题目给出的点值插值，而是自己带值插值，此时一般选择$(0,F(0)),...,(k,F(k))$，这时$F$的表达式可以写成：

$F(x)=\sum\limits_{i=0}^{k}\frac{F(i)x(x-1)\cdots(x-k)\times(-1)^{k-i}}{i!(k-i)!}$

就可以预处理阶乘什么的然后$O(k)$做了，但是对膜数有要求

代码（裸插值）：

 #include<algorithm>

 #include<iostream>

 #include<cstring>

 #include<cstdio>

 #include<cmath>

 #include<queue>

 #define inf 2147483647

 #define eps 1e-9

 #define mod 998244353

 using namespace std;

 typedef long long ll;

 int n,k,x[],y[];

 int fastpow(int x,int y){

     int ret=;

     for(;y;y>>=,x=(ll)x*x%mod){

         if(y&)ret=(ll)ret*x%mod;

     }

     return ret;

 }

 int lagrange(){

     int ret=,t1,t2;

     for(int i=;i<n;i++){

         t1=,t2=;

         for(int j=;j<n;j++){

             if(i!=j){

                 t1=(ll)t1*(k-x[j])%mod;

                 t2=(ll)t2*(x[i]-x[j])%mod;

             }

         }

         ret=(ret+(ll)t1*y[i]%mod*fastpow(t2,mod-)%mod)%mod;

     }

     return ret;

 }

 int main(){

     scanf("%d%d",&n,&k);

     for(int i=;i<n;i++){

         scanf("%d%d",&x[i],&y[i]);

     }

     printf("%d",(lagrange()+mod)%mod);

     return ;

 }

题目：

【BZOJ2655】Calc

DP+拉格朗日插值

 #include<algorithm>

 #include<iostream>

 #include<cstring>

 #include<cstdio>

 #include<cmath>

 #include<queue>

 #define inf 2147483647

 #define eps 1e-9

 using namespace std;

 typedef long long ll;

 ll n,m,p,t1,t2,ans=,f[][];

 ll fastpow(ll x,ll y){

     ll ret=;

     for(;y;y>>=,x=(ll)x*x%p){

         if(y&)ret=(ll)ret*x%p;

     }

     return ret;

 }

 int main(){

     scanf("%lld%lld%lld",&m,&n,&p);

     f[][]=;

     for(ll i=;i<=n*;i++){

         for(ll j=;j<=n;j++){

             if(!j)f[i][j]=f[i-][j];

             else f[i][j]=(f[i-][j]+(ll)f[i-][j-]*i%p*j%p)%p;

         }

     }

     if(m<=n*)return printf("%lld",f[m][n]),;

     for(ll i=;i<=n*;i++){

         t1=t2=;

         for(ll j=;j<=n*;j++){

             if(i==j)continue;

             t1=(ll)t1*(m-j)%p;

             t2=(ll)t2*(i-j)%p;

         }

         ans=(ans+f[i][n]*t1%p*fastpow(t2,p-)%p)%p;

     }

     ans=(ans+p)%p;

     printf("%lld",ans);

     return ;

 }

【BZOJ4559】【JLOI2016】成绩比较

组合数DP+拉格朗日插值

 #include<algorithm>

 #include<iostream>

 #include<cstring>

 #include<cstdio>

 #include<cmath>

 #include<queue>

 #define inf 2147483647

 #define eps 1e-9

 #define mod 1000000007

 using namespace std;

 typedef long long ll;

 int n,m,k,tmp,nw,u[],r[],g[],C[][],f[][];

 int fastpow(int x,int y){

     int ret=;

     for(;y;y>>=,x=(ll)x*x%mod){

         if(y&)ret=(ll)ret*x%mod;

     }

     return ret;

 }

 int main(){

     scanf("%d%d%d",&n,&m,&k);

     for(int i=;i<=m;i++){

         scanf("%d",&u[i]);

     }

     for(int i=;i<=m;i++){

         scanf("%d",&r[i]);

     }

     C[][]=;

     for(int i=;i<=n+;i++){

         C[i][]=;

         for(int j=;j<=i;j++){

             C[i][j]=(C[i-][j]+C[i-][j-])%mod;

         }

     }

     f[][n-]=;

     for(int i=;i<=m;i++){

         g[]=u[i];

         for(int j=;j<=n;j++){

             g[j]=(fastpow(u[i]+,j+)-+mod)%mod;

             for(int kk=;kk<j;kk++){

                 g[j]=(g[j]-(ll)C[j+][kk]*g[kk]%mod+mod)%mod;

             }

             g[j]=(ll)g[j]*fastpow(j+,mod-)%mod;

         }

         tmp=;

         nw=fastpow(u[i],r[i]-);

         int inv=fastpow(u[i],mod-);

         for(int j=;j<r[i];j++){

             tmp=(ll)(tmp+(ll)((j&)?-:)*C[r[i]-][j]*nw%mod*g[n-r[i]+j]%mod+mod)%mod;

             nw=(ll)nw*inv%mod;

         }

         for(int j=k;j<=n;j++){

             for(int kk=j;kk<=n;kk++){

                 if(n-r[i]-j>=){

                     f[i][j]=(f[i][j]+(ll)f[i-][kk]*C[kk][j]%mod*C[n-kk-][n-r[i]-j]%mod*tmp%mod)%mod;

                 }

             }

         }

     }

     printf("%d",(f[m][k]+mod)%mod);

     return ;

 }

【BZOJ3453】XLkxc

差分+拉格朗日插值

 #include<algorithm>

 #include<iostream>

 #include<cstring>

 #include<cstdio>

 #include<cmath>

 #include<queue>

 #define inf 2147483647

 #define eps 1e-9

 #define mod 1234567891

 using namespace std;

 typedef long long ll;

 ll t,k,a,n,d,f[],g[],inv[],suf[],pre[];

 ll fastpow(ll x,ll y){

     ll ret=;

     for(;y;y>>=,x=x*x%mod){

         if(y&)ret=ret*x%mod;

     }

     return ret;

 }

 void _(){

     inv[]=inv[]=;

     for(int i=;i<=;i++)inv[i]=(mod-mod/i)*inv[mod%i]%mod;

     for(int i=;i<=;i++)inv[i]=inv[i-]*inv[i]%mod;

 }

 ll lagrange(ll *s,ll x,ll n){

     ll ret=,tmp=;

     pre[]=suf[n+]=;

     for(int i=;i<=n+;i++)pre[i]=pre[i-]*(x-i+mod)%mod;

     for(int i=n+;i;i--)suf[i]=suf[i+]*(x-i+mod)%mod;

     for(int i=;i<=n+;i++){

         tmp=s[i]*pre[i-]%mod*suf[i+]%mod*inv[i-]%mod*inv[n-i+]%mod;

         if((n-i)%==)tmp=-tmp+mod;

         ret=(ret+tmp)%mod;

     }

     return ret;

 }

 int main(){

     _();

     scanf("%lld",&t);

     while(t--){

         scanf("%lld%lld%lld%lld",&k,&a,&n,&d);

         for(int i=;i<=k+;i++)g[i]=fastpow(i,k);

         for(int i=;i<=k+;i++)g[i]=(g[i-]+g[i])%mod;

         for(int i=;i<=k+;i++)g[i]=(g[i-]+g[i])%mod;

         f[]=lagrange(g,a,k+);

         for(int i=;i<=k+;i++)f[i]=(f[i-]+lagrange(g,(a+i*d)%mod,k+))%mod;

         printf("%lld\n",lagrange(f,n,k+));

     }

     return ;

 }

【xsy1537】五颜六色的幻想乡

拆边矩阵树定理+拉格朗日插值

 #include<algorithm>

 #include<iostream>

 #include<cstring>

 #include<cstdio>

 #include<cmath>

 #include<queue>

 #define inf 2147483647

 #define eps 1e-9

 #define mod 1000000007

 using namespace std;

 typedef long long ll;

 int N,n,m,s[][],f[],g[],h[],num[],x[],y[],z[];

 int fastpow(int x,int y){

     int ret=;

     for(;y;y>>=,x=(ll)x*x%mod){

         if(y&)ret=(ll)ret*x%mod;

     }

     return ret;

 }

 int calc(int xx){

     int ret=,pw=fastpow(xx,n),tmp,inv;

     memset(s,,sizeof(s));

     for(int i=;i<=m;i++){

         tmp=(z[i]==)?xx:(z[i]==)?pw:;

         s[x[i]][y[i]]=(s[x[i]][y[i]]-tmp)%mod;

         s[y[i]][x[i]]=(s[y[i]][x[i]]-tmp)%mod;

         s[x[i]][x[i]]=(s[x[i]][x[i]]+tmp)%mod;

         s[y[i]][y[i]]=(s[y[i]][y[i]]+tmp)%mod;

     }

     for(int i=;i<n;i++){

         tmp=;

         for(int j=i;j<n;j++){

             if(s[j][i]){

                 tmp=j;

                 break;

             }

         }

         if(!tmp)return ;

         if(tmp!=i){

             ret=-ret;

             for(int j=i;j<n;j++)swap(s[i][j],s[tmp][j]);

         }

         ret=(ll)ret*s[i][i]%mod;

         inv=fastpow(s[i][i],mod-);

         for(int j=i+;j<n;j++){

             if(s[j][i]){

                 int t=(ll)s[j][i]*inv%mod;

                 for(int k=i;k<n;k++){

                     s[j][k]=(s[j][k]-(ll)s[i][k]*t%mod+mod)%mod;

                 }

             }

         }

     }

     return ret;

 }

 void gao(){

     f[]=;

     for(int i=;i<=N;i++){

         for(int j=i;j>=;j--)f[j+]=f[j];

         f[]=;

         for(int j=;j<=i;j++)f[j]=(f[j]-(ll)f[j+]*i%mod+mod)%mod;

     }

     for(int i=;i<=N;i++){

         int tmp=;

         h[N+]=;

         for(int j=N;j>=;j--){

             if(i!=j)tmp=(ll)tmp*(i-j)%mod;

             h[j]=(f[j+]+(ll)h[j+]*i)%mod;

         }

         tmp=(ll)fastpow(tmp,mod-)*num[i]%mod;

         for(int j=;j<=N;j++){

             g[j]=(g[j]+(ll)tmp*h[j])%mod;

         }

     }

 }

 int main(){

     scanf("%d%d",&n,&m);

     for(int i=;i<=m;i++){

         scanf("%d%d%d",&x[i],&y[i],&z[i]);

     }

     N=n*n-n;

     for(int i=;i<=N;i++){

         num[i]=calc(i);

     }

     gao();

     for(int i=;i<=n;i++){

         for(int j=;j<n-i;j++){

             printf("%d\n",(g[i+j*n]+mod)%mod);

         }

     }

     return ;

 }

【BZOJ4162】shlw loves matrix II

其实这是道常系数齐次线性递推模板题哒

 #include<algorithm>

 #include<iostream>

 #include<cstring>

 #include<cstdio>

 #include<cmath>

 #include<queue>

 #define inf 2147483647

 #define eps 1e-9

 #define mod 1000000007

 using namespace std;

 typedef long long ll;

 typedef double db;

 struct lambda{

     int id,x;

     lambda(){}

     lambda(int _id,int _x){

         id=_id,x=_x;

     }

 }p[];

 int n,a[][],sq[][],tmp[][],dt[][],x[],pw[],t[],f[],g[],ans[][];

 char s[];

 int fastpow(int x,int y){

     int ret=;

     for(;y;y>>=,x=(ll)x*x%mod){

         if(y&)ret=(ll)ret*x%mod;

     }

     return ret;

 }

 int det(){

     int ret=,nw,tmp;

     for(int i=;i<=n;i++){

         nw=i;

         for(int j=i;j<=n;j++){

             if(dt[j][i]){

                 nw=j;

                 break;

             }

         }

         if(nw!=i){

             for(int j=i;j<=n;j++){

                 swap(dt[i][j],dt[nw][j]);

             }

             ret=-ret;

         }

         for(int j=i+;j<=n;j++){

             if(dt[j][i]){

                 tmp=(ll)dt[j][i]*fastpow(dt[i][i],mod-)%mod;

                 for(int k=i;k<=n;k++){

                     dt[j][k]=(dt[j][k]-(ll)dt[i][k]*tmp%mod)%mod;

                 }

             }

         }

         ret=(ll)ret*dt[i][i]%mod;

     }

     return (ret+mod)%mod;

 }

 void mul(int *s,int *x,int *y){

     for(int i=;i<=n*;i++)t[i]=;

     for(int i=;i<n;i++){

         for(int j=;j<n;j++){

             //add(t[i+j],(ll)x[i]*y[j]%mod);

             t[i+j]=(t[i+j]+(ll)x[i]*y[j]%mod)%mod;

         }

     }

     int mo=fastpow(f[n],mod-);

     for(int i=n*-;i>=n;i--){

         for(int j=n-;j>=;j--){

             //add(t[i+j-n],mod-(ll)f[j]*t[i]%mod*mo%mod);

             t[i+j-n]=(t[i+j-n]-(ll)f[j]*t[i]%mod*mo%mod)%mod;

         }

     }

     for(int i=;i<n;i++){

         s[i]=t[i];

     }

 }

 int main(){

     scanf("%s%d",s,&n);

     for(int i=;i<=n;i++){

         for(int j=;j<=n;j++){

             scanf("%d",&a[i][j]);

         }

     }

     for(int i=;i<=n;i++){

         memcpy(dt,a,sizeof(a));

         for(int j=;j<=n;j++)dt[j][j]-=i;

         p[i]=lambda(i,det());

     }

     for(int i=;i<=n;i++){

         for(int j=;j<=n;j++)g[j]=;

         g[]=p[i].x;

         for(int j=;j<=n;j++){

             if(i!=j){

                 g[]=(ll)g[]*fastpow(p[j].id-p[i].id,mod-)%mod;

             }

         }

         for(int j=;j<=n;j++){

             if(i!=j){

                 for(int k=n;k;k--){

                     g[k]=((ll)g[k]*p[j].id%mod-g[k-])%mod;

                 }

                 g[]=(ll)g[]*p[j].id%mod;

             }

         }

         for(int j=;j<=n;j++){

             f[j]=(f[j]+g[j])%mod;

         }

     }

     x[]=pw[]=;

     for(int i=strlen(s)-;i>=;i--){

         if(s[i]=='')mul(pw,pw,x);

         mul(x,x,x);

     }

     for(int i=;i<=n;i++){

         sq[i][i]=;

     }

     for(int i=;i<n;i++){

         for(int j=;j<=n;j++){

             for(int k=;k<=n;k++){

                 //add(ans[j][k],(ll)sq[j][k]*pw[i]%mod);

                 ans[j][k]=(ans[j][k]+(ll)sq[j][k]*pw[i]%mod)%mod;

             }

         }

         memset(tmp,,sizeof(tmp));

         for(int j=;j<=n;j++){

             for(int k=;k<=n;k++){

                 for(int l=;l<=n;l++){

                     //add(tmp[j][k],(ll)sq[j][l]*a[l][k]%mod);

                     tmp[j][k]=(tmp[j][k]+(ll)sq[j][l]*a[l][k]%mod)%mod;

                 }

             }

         }

         memcpy(sq,tmp,sizeof(tmp));

     }

     for(int i=;i<=n;i++){

         for(int j=;j<=n;j++){

             printf("%d ",(ans[i][j]%mod+mod)%mod);

         }

         puts("");

     }

     return ;

 }

剩下道【NOI2012】骑行川藏，神仙题，咕咕咕

快速插值

本部分参考了yww的博客：Orzyww

先讲个东西：多点求值

多点求值

给出一个$k$次多项式$A(x)$和$n$个点$x_0,x_1,...,x_{n-1}$，求多项式在这$n$个点的值；

显然暴力是$O(nk)$的，因此要优化。

考虑一个简（shen）单（xian）的构造：构造多项式$B_i(x)=x-x_i$，$C_i(x)=A(x) \mod B_i(x)$，易知$B_i(x_i)=0$，所以$A(x_i)=C_i(x_i)$；

但是这样计算$B_i$和$C_i$依然是$O(n)$的，还需要优化；

可以考虑类似多项式求逆取模那样分治倍增。

先把当前点集$X=\{x_0,x_1,...,x_{n-1}\}$分成两半：

$X_0=\{x_0,x_1,...,x_{\lfloor\frac{n}{2}\rfloor -1}\}$

$X_1=\{x_{\lfloor\frac{n}{2}\rfloor},x_{\lfloor\frac{n}{2}\rfloor+1},...,x_{n-1}\}$

然后构造多项式$B_0,B_1$类似把$B$分成两半，$A$同理：

$B_0(x)=\prod\limits_{i=0}^{\lfloor\frac{n}{2}\rfloor-1}(x-x_i)$

$B_1(x)=\prod\limits_{i=\lfloor\frac{n}{2}\rfloor}^{n-1}(x-x_i)$

$A_0(x)=A(x) \mod B_0(x)$

$A_1(x)=A(x) \mod B_1(x)$

容易看出，这样分治下去当$x∈X_0$时，$A(x)=A_0(x)$，否则$A(x)=A_1(x)$，可以每次递归处理。

每层时间复杂度是$O(nlogn)$的，根据主定理，总的时间复杂度就是：

$T(n)=2T(\frac{n}{2})+O(nlogn)=O(nlog^2n)$

快速插值

拉格朗日插值是$O(n^2)$的，只有特殊情况才能优化到$O(n)$，而快速插值法可以在任意情况下做到$O(nlog^2n)$的时间复杂度的插值。

快速插值法是基于拉格朗日插值法进行数学优化的，那么先来回顾一下拉格朗日插值公式：

$F(x)=\sum\limits_{i=0}^{n}\frac{\prod\limits_{j\neq i}(x-x_j)}{\prod\limits_{j\neq i}(x_i-x_j)}\times y_i$

主要问题就是怎么快速求分子分母。

先考虑分母，设$g_i=\prod\limits_{j\neq i}(x_i-x_j)$，则：

$g_i=\lim\limits_{x \to x_i}\frac{\prod\limits_{j=0}^{n}(x-x_j)}{x-x_i}$

$=(\prod\limits_{j=0}^{n}(x-x_j))'|_{x=x_i}$

于是就可以分治求出$\prod\limits_{j=0}^{n}(x-x_j)$，求导后对所有$x_i$多点求值即可；

分子也可以类似分治求，处理好分母然后顺手合并答案就好了。

时间复杂度易证是$O(nlog^2n)$

代码（多点求值+快速插值）：

太难了，先咕着

巴特西

拉格朗日插值&&快速插值

拉格朗日插值

一个小技巧：

代码（裸插值）：

题目：

【BZOJ2655】Calc

【BZOJ4559】【JLOI2016】成绩比较

【BZOJ3453】XLkxc

【xsy1537】五颜六色的幻想乡

【BZOJ4162】shlw loves matrix II

快速插值

多点求值

快速插值

代码（多点求值+快速插值）：

最新文章

热门文章