拉格朗日插值
插值真惨
众所周知$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)$做了，但是对膜数有要求
代码（裸插值）：
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 #define mod 998244353
10 using namespace std; 11 typedef long long ll; 12 int n,k,x[2001],y[2001]; 13 int fastpow(int x,int y){ 14 int ret=1; 15 for(;y;y>>=1,x=(ll)x*x%mod){ 16 if(y&1)ret=(ll)ret*x%mod; 17 } 18 return ret; 19 } 20 int lagrange(){ 21 int ret=0,t1,t2; 22 for(int i=0;i<n;i++){ 23 t1=1,t2=1; 24 for(int j=0;j<n;j++){ 25 if(i!=j){ 26 t1=(ll)t1*(k-x[j])%mod; 27 t2=(ll)t2*(x[i]-x[j])%mod; 28 } 29 } 30 ret=(ret+(ll)t1*y[i]%mod*fastpow(t2,mod-2)%mod)%mod; 31 } 32 return ret; 33 } 34 int main(){ 35 scanf("%d%d",&n,&k); 36 for(int i=0;i<n;i++){ 37 scanf("%d%d",&x[i],&y[i]); 38 } 39 printf("%d",(lagrange()+mod)%mod); 40 return 0; 41 }
题目：
【BZOJ2655】Calc
DP+拉格朗日插值
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 using namespace std; 10 typedef long long ll; 11 ll n,m,p,t1,t2,ans=0,f[1001][1001]; 12 ll fastpow(ll x,ll y){ 13 ll ret=1; 14 for(;y;y>>=1,x=(ll)x*x%p){ 15 if(y&1)ret=(ll)ret*x%p; 16 } 17 return ret; 18 } 19 int main(){ 20 scanf("%lld%lld%lld",&m,&n,&p); 21 f[0][0]=1; 22 for(ll i=1;i<=n*2;i++){ 23 for(ll j=0;j<=n;j++){ 24 if(!j)f[i][j]=f[i-1][j]; 25 else f[i][j]=(f[i-1][j]+(ll)f[i-1][j-1]*i%p*j%p)%p; 26 } 27 } 28 if(m<=n*2)return printf("%lld",f[m][n]),0; 29 for(ll i=0;i<=n*2;i++){ 30 t1=t2=1; 31 for(ll j=0;j<=n*2;j++){ 32 if(i==j)continue; 33 t1=(ll)t1*(m-j)%p; 34 t2=(ll)t2*(i-j)%p; 35 } 36 ans=(ans+f[i][n]*t1%p*fastpow(t2,p-2)%p)%p; 37 } 38 ans=(ans+p)%p; 39 printf("%lld",ans); 40 return 0; 41 }
【BZOJ4559】【JLOI2016】成绩比较
组合数DP+拉格朗日插值
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 #define mod 1000000007
10 using namespace std; 11 typedef long long ll; 12 int n,m,k,tmp,nw,u[201],r[201],g[201],C[201][201],f[201][201]; 13 int fastpow(int x,int y){ 14 int ret=1; 15 for(;y;y>>=1,x=(ll)x*x%mod){ 16 if(y&1)ret=(ll)ret*x%mod; 17 } 18 return ret; 19 } 20 int main(){ 21 scanf("%d%d%d",&n,&m,&k); 22 for(int i=1;i<=m;i++){ 23 scanf("%d",&u[i]); 24 } 25 for(int i=1;i<=m;i++){ 26 scanf("%d",&r[i]); 27 } 28 C[0][0]=1; 29 for(int i=1;i<=n+1;i++){ 30 C[i][0]=1; 31 for(int j=1;j<=i;j++){ 32 C[i][j]=(C[i-1][j]+C[i-1][j-1])%mod; 33 } 34 } 35 f[0][n-1]=1; 36 for(int i=1;i<=m;i++){ 37 g[0]=u[i]; 38 for(int j=1;j<=n;j++){ 39 g[j]=(fastpow(u[i]+1,j+1)-1+mod)%mod; 40 for(int kk=0;kk<j;kk++){ 41 g[j]=(g[j]-(ll)C[j+1][kk]*g[kk]%mod+mod)%mod; 42 } 43 g[j]=(ll)g[j]*fastpow(j+1,mod-2)%mod; 44 } 45 tmp=0; 46 nw=fastpow(u[i],r[i]-1); 47 int inv=fastpow(u[i],mod-2); 48 for(int j=0;j<r[i];j++){ 49 tmp=(ll)(tmp+(ll)((j&1)?-1:1)*C[r[i]-1][j]*nw%mod*g[n-r[i]+j]%mod+mod)%mod; 50 nw=(ll)nw*inv%mod; 51 } 52 for(int j=k;j<=n;j++){ 53 for(int kk=j;kk<=n;kk++){ 54 if(n-r[i]-j>=0){ 55 f[i][j]=(f[i][j]+(ll)f[i-1][kk]*C[kk][j]%mod*C[n-kk-1][n-r[i]-j]%mod*tmp%mod)%mod; 56 } 57 } 58 } 59 } 60 printf("%d",(f[m][k]+mod)%mod); 61 return 0; 62 }
【BZOJ3453】XLkxc
差分+拉格朗日插值
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 #define mod 1234567891
10 using namespace std; 11 typedef long long ll; 12 ll t,k,a,n,d,f[210],g[210],inv[210],suf[210],pre[210]; 13 ll fastpow(ll x,ll y){ 14 ll ret=1; 15 for(;y;y>>=1,x=x*x%mod){ 16 if(y&1)ret=ret*x%mod; 17 } 18 return ret; 19 } 20 void _(){ 21 inv[0]=inv[1]=1; 22 for(int i=2;i<=200;i++)inv[i]=(mod-mod/i)*inv[mod%i]%mod; 23 for(int i=2;i<=200;i++)inv[i]=inv[i-1]*inv[i]%mod; 24 } 25 ll lagrange(ll *s,ll x,ll n){ 26 ll ret=0,tmp=0; 27 pre[0]=suf[n+2]=1; 28 for(int i=1;i<=n+1;i++)pre[i]=pre[i-1]*(x-i+mod)%mod; 29 for(int i=n+1;i;i--)suf[i]=suf[i+1]*(x-i+mod)%mod; 30 for(int i=1;i<=n+1;i++){ 31 tmp=s[i]*pre[i-1]%mod*suf[i+1]%mod*inv[i-1]%mod*inv[n-i+1]%mod; 32 if((n-i)%2==0)tmp=-tmp+mod; 33 ret=(ret+tmp)%mod; 34 } 35 return ret; 36 } 37 int main(){ 38 _(); 39 scanf("%lld",&t); 40 while(t--){ 41 scanf("%lld%lld%lld%lld",&k,&a,&n,&d); 42 for(int i=0;i<=k+3;i++)g[i]=fastpow(i,k); 43 for(int i=1;i<=k+3;i++)g[i]=(g[i-1]+g[i])%mod; 44 for(int i=1;i<=k+3;i++)g[i]=(g[i-1]+g[i])%mod; 45 f[0]=lagrange(g,a,k+2); 46 for(int i=1;i<=k+5;i++)f[i]=(f[i-1]+lagrange(g,(a+i*d)%mod,k+2))%mod; 47 printf("%lld\n",lagrange(f,n,k+4)); 48 } 49 return 0; 50 }
【xsy1537】五颜六色的幻想乡
拆边矩阵树定理+拉格朗日插值
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 #define mod 1000000007
10 using namespace std; 11 typedef long long ll; 12 int N,n,m,s[51][51],f[3001],g[3001],h[3001],num[3001],x[3001],y[3001],z[3001]; 13 int fastpow(int x,int y){ 14 int ret=1; 15 for(;y;y>>=1,x=(ll)x*x%mod){ 16 if(y&1)ret=(ll)ret*x%mod; 17 } 18 return ret; 19 } 20 int calc(int xx){ 21 int ret=1,pw=fastpow(xx,n),tmp,inv; 22 memset(s,0,sizeof(s)); 23 for(int i=1;i<=m;i++){ 24 tmp=(z[i]==1)?xx:(z[i]==2)?pw:1; 25 s[x[i]][y[i]]=(s[x[i]][y[i]]-tmp)%mod; 26 s[y[i]][x[i]]=(s[y[i]][x[i]]-tmp)%mod; 27 s[x[i]][x[i]]=(s[x[i]][x[i]]+tmp)%mod; 28 s[y[i]][y[i]]=(s[y[i]][y[i]]+tmp)%mod; 29 } 30 for(int i=1;i<n;i++){ 31 tmp=0; 32 for(int j=i;j<n;j++){ 33 if(s[j][i]){ 34 tmp=j; 35 break; 36 } 37 } 38 if(!tmp)return 0; 39 if(tmp!=i){ 40 ret=-ret; 41 for(int j=i;j<n;j++)swap(s[i][j],s[tmp][j]); 42 } 43 ret=(ll)ret*s[i][i]%mod; 44 inv=fastpow(s[i][i],mod-2); 45 for(int j=i+1;j<n;j++){ 46 if(s[j][i]){ 47 int t=(ll)s[j][i]*inv%mod; 48 for(int k=i;k<n;k++){ 49 s[j][k]=(s[j][k]-(ll)s[i][k]*t%mod+mod)%mod; 50 } 51 } 52 } 53 } 54 return ret; 55 } 56 void gao(){ 57 f[0]=1; 58 for(int i=0;i<=N;i++){ 59 for(int j=i;j>=0;j--)f[j+1]=f[j]; 60 f[0]=0; 61 for(int j=0;j<=i;j++)f[j]=(f[j]-(ll)f[j+1]*i%mod+mod)%mod; 62 } 63 for(int i=0;i<=N;i++){ 64 int tmp=1; 65 h[N+1]=0; 66 for(int j=N;j>=0;j--){ 67 if(i!=j)tmp=(ll)tmp*(i-j)%mod; 68 h[j]=(f[j+1]+(ll)h[j+1]*i)%mod; 69 } 70 tmp=(ll)fastpow(tmp,mod-2)*num[i]%mod; 71 for(int j=0;j<=N;j++){ 72 g[j]=(g[j]+(ll)tmp*h[j])%mod; 73 } 74 } 75 } 76 int main(){ 77 scanf("%d%d",&n,&m); 78 for(int i=1;i<=m;i++){ 79 scanf("%d%d%d",&x[i],&y[i],&z[i]); 80 } 81 N=n*n-n; 82 for(int i=0;i<=N;i++){ 83 num[i]=calc(i); 84 } 85 gao(); 86 for(int i=0;i<=n;i++){ 87 for(int j=0;j<n-i;j++){ 88 printf("%d\n",(g[i+j*n]+mod)%mod); 89 } 90 } 91 return 0; 92 }
【BZOJ4162】shlw loves matrix II
其实这是道常系数齐次线性递推模板题哒
1 #include<algorithm>
2 #include<iostream>
3 #include<cstring>
4 #include<cstdio>
5 #include<cmath>
6 #include<queue>
7 #define inf 2147483647
8 #define eps 1e-9
9 #define mod 1000000007
10 using namespace std; 11 typedef long long ll; 12 typedef double db; 13 struct lambda{ 14 int id,x; 15 lambda(){} 16 lambda(int _id,int _x){ 17 id=_id,x=_x; 18 } 19 }p[55]; 20 int n,a[55][55],sq[55][55],tmp[55][55],dt[55][55],x[110],pw[110],t[110],f[55],g[55],ans[55][55]; 21 char s[10001]; 22 int fastpow(int x,int y){ 23 int ret=1; 24 for(;y;y>>=1,x=(ll)x*x%mod){ 25 if(y&1)ret=(ll)ret*x%mod; 26 } 27 return ret; 28 } 29 int det(){ 30 int ret=1,nw,tmp; 31 for(int i=1;i<=n;i++){ 32 nw=i; 33 for(int j=i;j<=n;j++){ 34 if(dt[j][i]){ 35 nw=j; 36 break; 37 } 38 } 39 if(nw!=i){ 40 for(int j=i;j<=n;j++){ 41 swap(dt[i][j],dt[nw][j]); 42 } 43 ret=-ret; 44 } 45 for(int j=i+1;j<=n;j++){ 46 if(dt[j][i]){ 47 tmp=(ll)dt[j][i]*fastpow(dt[i][i],mod-2)%mod; 48 for(int k=i;k<=n;k++){ 49 dt[j][k]=(dt[j][k]-(ll)dt[i][k]*tmp%mod)%mod; 50 } 51 } 52 } 53 ret=(ll)ret*dt[i][i]%mod; 54 } 55 return (ret+mod)%mod; 56 } 57 void mul(int *s,int *x,int *y){ 58 for(int i=0;i<=n*2;i++)t[i]=0; 59 for(int i=0;i<n;i++){ 60 for(int j=0;j<n;j++){ 61 //add(t[i+j],(ll)x[i]*y[j]%mod);
62 t[i+j]=(t[i+j]+(ll)x[i]*y[j]%mod)%mod; 63 } 64 } 65 int mo=fastpow(f[n],mod-2); 66 for(int i=n*2-2;i>=n;i--){ 67 for(int j=n-1;j>=0;j--){ 68 //add(t[i+j-n],mod-(ll)f[j]*t[i]%mod*mo%mod);
69 t[i+j-n]=(t[i+j-n]-(ll)f[j]*t[i]%mod*mo%mod)%mod; 70 } 71 } 72 for(int i=0;i<n;i++){ 73 s[i]=t[i]; 74 } 75 } 76 int main(){ 77 scanf("%s%d",s,&n); 78 for(int i=1;i<=n;i++){ 79 for(int j=1;j<=n;j++){ 80 scanf("%d",&a[i][j]); 81 } 82 } 83 for(int i=0;i<=n;i++){ 84 memcpy(dt,a,sizeof(a)); 85 for(int j=1;j<=n;j++)dt[j][j]-=i; 86 p[i]=lambda(i,det()); 87 } 88 for(int i=0;i<=n;i++){ 89 for(int j=0;j<=n;j++)g[j]=0; 90 g[0]=p[i].x; 91 for(int j=0;j<=n;j++){ 92 if(i!=j){ 93 g[0]=(ll)g[0]*fastpow(p[j].id-p[i].id,mod-2)%mod; 94 } 95 } 96 for(int j=0;j<=n;j++){ 97 if(i!=j){ 98 for(int k=n;k;k--){ 99 g[k]=((ll)g[k]*p[j].id%mod-g[k-1])%mod; 100 } 101 g[0]=(ll)g[0]*p[j].id%mod; 102 } 103 } 104 for(int j=0;j<=n;j++){ 105 f[j]=(f[j]+g[j])%mod; 106 } 107 } 108 x[1]=pw[0]=1; 109 for(int i=strlen(s)-1;i>=0;i--){ 110 if(s[i]=='1')mul(pw,pw,x); 111 mul(x,x,x); 112 } 113 for(int i=1;i<=n;i++){ 114 sq[i][i]=1; 115 } 116 for(int i=0;i<n;i++){ 117 for(int j=1;j<=n;j++){ 118 for(int k=1;k<=n;k++){ 119 //add(ans[j][k],(ll)sq[j][k]*pw[i]%mod);
120 ans[j][k]=(ans[j][k]+(ll)sq[j][k]*pw[i]%mod)%mod; 121 } 122 } 123 memset(tmp,0,sizeof(tmp)); 124 for(int j=1;j<=n;j++){ 125 for(int k=1;k<=n;k++){ 126 for(int l=1;l<=n;l++){ 127 //add(tmp[j][k],(ll)sq[j][l]*a[l][k]%mod);
128 tmp[j][k]=(tmp[j][k]+(ll)sq[j][l]*a[l][k]%mod)%mod; 129 } 130 } 131 } 132 memcpy(sq,tmp,sizeof(tmp)); 133 } 134 for(int i=1;i<=n;i++){ 135 for(int j=1;j<=n;j++){ 136 printf("%d ",(ans[i][j]%mod+mod)%mod); 137 } 138 puts(""); 139 } 140 return 0; 141 }
剩下道【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)$
代码（多点求值+快速插值）：
太难了，先咕着