拉格朗日插值&&快速插值
拉格朗日插值
插值真惨
众所周知k+1k+1个点可以确定一个kk次多项式,那么插值就是通过点值还原多项式的过程。
设给出的k+1k+1个点分别是(x0,y0),(x1,y1),...,(xk,yk)(x0,y0),(x1,y1),...,(xk,yk),那么xjb构造一下:
设函数fi(x)=∏j≠i(x−xj)∏j≠i(xi−xj)×yifi(x)=∏j≠i(x−xj)∏j≠i(xi−xj)×yi
显然这个函数当x=xix=xi时值为yiyi,x=xj(0≤j≤k且j≠i)x=xj(0≤j≤k且j≠i)时值为0。
求出所有的fifi,然后把它们加起来,发现得到的多项式必过给出的这些点,也就是要求的结果了。。。
即要求的多项式F(x)=∑i=0kfi(x)F(x)=∑i=0kfi(x)
易知时间复杂度是O(k2)O(k2)的,看起来很暴力,但是拉格朗日插值就是这么暴力。。。
如果要动态加点的话可以稍微优化一下:
注意到每个fifi中∏j≠i(x−xj)∏j≠i(x−xj)重复计算了,可以简化:
设p(x)=∏i=0k(x−xi)p(x)=∏i=0k(x−xi),qi=yi∏j≠i(xi−xj)qi=yi∏j≠i(xi−xj)
则F(x)=p(x)×∑i=0kqi(x−xi)F(x)=p(x)×∑i=0kqi(x−xi)
这样子原本的复杂度不变,但是新加点时只用重新算一遍qiqi,复杂度为O(k)O(k)
这个东西貌似叫重心法?
一个小技巧:
在实际做题的时候,有时不需要根据题目给出的点值插值,而是自己带值插值,此时一般选择(0,F(0)),...,(k,F(k))(0,F(0)),...,(k,F(k)),这时FF的表达式可以写成:
F(x)=∑i=0kF(i)x(x−1)⋯(x−k)×(−1)k−ii!(k−i)!F(x)=∑i=0kF(i)x(x−1)⋯(x−k)×(−1)k−ii!(k−i)!
就可以预处理阶乘什么的然后O(k)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
先讲个东西:多点求值
多点求值
给出一个kk次多项式A(x)A(x)和nn个点x0,x1,...,xn−1x0,x1,...,xn−1,求多项式在这nn个点的值;
显然暴力是O(nk)O(nk)的,因此要优化。
考虑一个简(shen)单(xian)的构造:构造多项式Bi(x)=x−xiBi(x)=x−xi,Ci(x)=A(x)modBi(x)Ci(x)=A(x)modBi(x),易知Bi(xi)=0Bi(xi)=0,所以A(xi)=Ci(xi)A(xi)=Ci(xi);
但是这样计算BiBi和CiCi依然是O(n)O(n)的,还需要优化;
可以考虑类似多项式求逆取模那样分治倍增。
先把当前点集X={x0,x1,...,xn−1}X={x0,x1,...,xn−1}分成两半:
X0={x0,x1,...,x⌊n2⌋−1}X0={x0,x1,...,x⌊n2⌋−1}
X1={x⌊n2⌋,x⌊n2⌋+1,...,xn−1}X1={x⌊n2⌋,x⌊n2⌋+1,...,xn−1}
然后构造多项式B0,B1B0,B1类似把BB分成两半,AA同理:
B0(x)=∏i=0⌊n2⌋−1(x−xi)B0(x)=∏i=0⌊n2⌋−1(x−xi)
B1(x)=∏i=⌊n2⌋n−1(x−xi)B1(x)=∏i=⌊n2⌋n−1(x−xi)
A0(x)=A(x)modB0(x)A0(x)=A(x)modB0(x)
A1(x)=A(x)modB1(x)A1(x)=A(x)modB1(x)
容易看出,这样分治下去当x∈X0x∈X0时,A(x)=A0(x)A(x)=A0(x),否则A(x)=A1(x)A(x)=A1(x),可以每次递归处理。
每层时间复杂度是O(nlogn)O(nlogn)的,根据主定理,总的时间复杂度就是:
T(n)=2T(n2)+O(nlogn)=O(nlog2n)T(n)=2T(n2)+O(nlogn)=O(nlog2n)
快速插值
拉格朗日插值是O(n2)O(n2)的,只有特殊情况才能优化到O(n)O(n),而快速插值法可以在任意情况下做到O(nlog2n)O(nlog2n)的时间复杂度的插值。
快速插值法是基于拉格朗日插值法进行数学优化的,那么先来回顾一下拉格朗日插值公式:
F(x)=∑i=0n∏j≠i(x−xj)∏j≠i(xi−xj)×yiF(x)=∑i=0n∏j≠i(x−xj)∏j≠i(xi−xj)×yi
主要问题就是怎么快速求分子分母。
先考虑分母,设gi=∏j≠i(xi−xj)gi=∏j≠i(xi−xj),则:
gi=limx→xi∏j=0n(x−xj)x−xigi=limx→xi∏j=0n(x−xj)x−xi
=(∏j=0n(x−xj))′|x=xi=(∏j=0n(x−xj))′|x=xi
于是就可以分治求出∏j=0n(x−xj)∏j=0n(x−xj),求导后对所有xixi多点求值即可;
分子也可以类似分治求,处理好分母然后顺手合并答案就好了。
时间复杂度易证是O(nlog2n)O(nlog2n)
代码(多点求值+快速插值):
太难了,先咕着