拉格朗日插值&&快速插值
拉格朗日插值
插值真惨
众所周知$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)$
代码（多点求值+快速插值）：
太难了，先咕着
