【学习笔记】拉格朗日插值

仰望星空的蚂蚁

已于 2023-06-04 15:37:14 修改

阅读量618

点赞数

分类专栏：多项式文章标签：学习 c++ 算法

于 2022-09-30 23:18:16 首次发布

本文链接：https://blog.csdn.net/cqbzlydd/article/details/126995959

版权

多项式专栏收录该内容

9 篇文章 0 订阅

订阅专栏

[NOIP2020] 微信步数

首先把问题转化为将每一天中 还存活的起点数量 计入贡献。

然后考虑单独一维 $j$ 的情况。

注意到从第 $2$ 轮开始，每次死的起点数是固定的。设最后 $i$ 步新死的起点数是 $f_j$ （事实上 $f_j$ 和进行的轮数无关，只和维度 $j$ 有关），那么事实上这是个一次函数。

首先第 $1$ 轮单独算。然后第 $T + 2$ 轮贡献 $\prod_{j=1}^k(w_j-r_j+l_j-T\times v-f_j)$ 是 $k$ 次多项式，求其前缀和为 $k + 1$ 次多项式。

基于以上观察，我们可以拉插求解，复杂度 $O(nk^2)$ 。

总结：问题的关键在于想到周期性的规律，然后大力多项式。（ $y l y$ 说这是最基础的题了。。。）

#include<bits/stdc++.h>
#define ll long long
#define fi first
#define se second
using namespace std;
const int mod=1e9+7;
int n,K;
ll w[15],v[15],v2[15],f[15];
int c[500005],d[500005];
ll l[15],r[15],l2[15],r2[15],inv[15];
ll res;
ll fpow(ll x,ll y){
	ll z(1);for(;y;y>>=1){
		if(y&1)z=z*x%mod;x=x*x%mod;
	}return z;
}
ll LG(int y) {
	if(y<=K+1) {
		ll res=0;
		for(int i=0;i<=y;i++) {
			ll mul=1;
			for(int j=1;j<=K;j++){
				mul=mul*(w[j]-r[j]+l[j]-i*abs(v[j])-f[j])%mod;
			}res=(res+mul)%mod;
		}return res;
	}
	ll res=0,tmp=0;
	for(int i=0;i<=K+1;i++) {
		ll mul=1;
		for(int j=1;j<=K;j++) {
			mul=mul*(w[j]-r[j]+l[j]-i*abs(v[j])-f[j])%mod;
		}tmp=(tmp+mul)%mod;mul=tmp;
		for(int j=0;j<=K+1;j++) {
			if(i!=j) {
				if(i>=j)mul=mul*(y-j)%mod*inv[i-j]%mod;
				else mul=-mul*(y-j)%mod*inv[j-i]%mod;
			}
		}res=(res+mul)%mod;
	}return res;
} 
int main() {
	ios::sync_with_stdio(false);
	cin.tie(0),cout.tie(0);
	cin>>n>>K;for(int i=1;i<=K+1;i++)inv[i]=fpow(i,mod-2);
	for(int i=1;i<=K;i++)cin>>w[i];
	for(int i=1;i<=n;i++) {
		cin>>c[i]>>d[i],v[c[i]]+=d[i];
		l[c[i]]=min(l[c[i]],v[c[i]]),r[c[i]]=max(r[c[i]],v[c[i]]);
	}
	int ok=0;for(int i=1;i<=K;i++)if(v[i])ok=1;
	if(!ok) {cout<< -1;return 0;}
	res=1;for(int i=1;i<=K;i++)res=res*w[i]%mod;
	for(int i=1;i<=n;i++) {
        v2[c[i]]+=d[i];
		l2[c[i]]=min(l2[c[i]],v2[c[i]]),r2[c[i]]=max(r2[c[i]],v2[c[i]]);
		int ok=1;ll mx=0x3f3f3f3f;ll mul=1;for(int j=1;j<=K;j++){
            if(r2[j]-l2[j]>=w[j]){ok=0;break;}mul=mul*(w[j]-r2[j]+l2[j])%mod;
        }if(!ok)continue;res=(res+mul)%mod;
        for(int j=1;j<=K;j++) {
            if(!v[j])continue;
            if(v[j]>0)f[j]=max(0ll,r2[j]+v[j]-r[j]);
            else f[j]=max(0ll,-l2[j]-v[j]+l[j]);
        }
		for(int j=1;j<=K;j++){
            if(r[j]-l[j]+f[j]>=w[j]){ok=0;break;}
            if(v[j])mx=min(mx,(w[j]-r[j]+l[j]-f[j])/abs(v[j]));
        }if(!ok)continue;res=(res+LG(mx))%mod;
	} cout<<(res+mod)%mod;
}

[APIO2016] 划艇

拉插做法比较复杂。~~滚回去学多项式基础了~~

考虑简单 $d p$ 。然后是关键的观察：

$1.1$ 从 $[0, m]$ 中选 $n$ 个数，满足非零数单调递增的方案数是 $\binom{n+m}{n}$

证明：考虑枚举 $0$ 的个数，答案是 $\sum_{i=0}^n\binom{n}{i}\binom{m}{n-i}$ 。容易发现它等价于从 $n + m$ 个数中选取 $n$ 个数，也就是 $\binom{n+m}{n}$ 。

如果钦定第 $n$ 个数非零的话，答案是 $\binom{n+m}{n}-\binom{n+m-1}{n-1}$ 。

设 $f_{i,j}$ 表示第 $i$ 个学校派出的划艇数目在第 $j$ 段的方案数。显然维护前缀和转移。复杂度 $O(n^3)$ 。

其实做到这一步可以想到拉插，形式化的说，每一段 $a_i,a_{i+1})$ 对应一个 $n$ 次多项式。

怎么观察到这一步的呢。注意看每一段里面具体的 $j$ 只影响组合数的值，前面公共的那一坨可以看成常数，而组合数又是与 $j$ 有关的多项式，于是可以拉插。自然一段内的前缀和也是 $n + 1$ 次多项式。

具体实现看代码。

~~我信仰什么，我便实现哪种方法。~~

#include<bits/stdc++.h>
#define ll long long
#define fi first
#define se second
using namespace std;
const int mod=1e9+7;
int n,a[505],b[505],lsh[1005],cnt;
ll f[505][1005],inv[505],res;
ll fpow(ll x,ll y) {
 ll z(1);
 for(;y;y>>=1) {
  if(y&1)z=z*x%mod;
  x=x*x%mod;
 }return z;
}
int main() {
 ios::sync_with_stdio(false);
 cin>>n;for(int i=1;i<=n;i++) inv[i]=fpow(i,mod-2);
 for(int i=1;i<=n;i++) {
  cin>>a[i]>>b[i],b[i]++;
  lsh[++cnt]=a[i],lsh[++cnt]=b[i];
 }sort(lsh+1,lsh+1+cnt),cnt=unique(lsh+1,lsh+1+cnt)-lsh-1;
 for(int i=1;i<=n;i++) {
  a[i]=lower_bound(lsh+1,lsh+1+cnt,a[i])-lsh;
  b[i]=lower_bound(lsh+1,lsh+1+cnt,b[i])-lsh-1;
 }
 for(int i=0;i<=cnt;i++) f[0][i]=1;
 for(int i=1;i<=n;i++) {
  for(int j=a[i];j<=b[i];j++) {
   ll L=1,mul=lsh[j+1]-lsh[j];
   for(int k=i-1;k>=0;k--) {
    f[i][j]=(f[i][j]+mul*f[k][j-1]%mod)%mod;
    if(a[k]<=j&&j<=b[k]) {
     mul=mul*(lsh[j+1]-lsh[j]+L)%mod*inv[L+1]%mod,L++;
    }
   }  }
  for(int j=a[i];j<=cnt;j++) {
   f[i][j]=(f[i][j]+f[i][j-1])%mod;
  }
  res=(res+f[i][cnt])%mod; }cout<<res;
}

大常数拉插。~~好逊啊~~

#include<bits/stdc++.h>
#define ll long long
#define fi first
#define se second
using namespace std;
const int mod=1e9+7;
int n,a[505],b[505],lsh[1005],cnt;
ll fac[505],inv[505],res;
ll g[1005][505],val[1005],sum[1005];
ll iv[1005][505];
ll fpow(ll x,ll y) {
 ll z(1);
 for(;y;y>>=1) {
  if(y&1)z=z*x%mod;
  x=x*x%mod;
 }return z;
}
void add(ll &x,ll y) {x=(x+y)%mod;}
ll LG(int x,ll y) {
 if(lsh[x+1]-lsh[x]<=n+2) {
  ll tmp=0;for(int j=0;j<lsh[x+1]-lsh[x];j++)add(tmp,g[x][j]);
  return tmp;
 }ll tmp=0,res=0,mul=1;
 for(int i=0;i<=n+1;i++)mul=mul*(y-i)%mod;
 for(int i=0;i<=n+1;i++) {
  add(tmp,g[x][i]);
  if(n+1-i&1) {
   res=(res-tmp*mul%mod*iv[x][i]%mod*inv[i]%mod*inv[n+1-i]%mod)%mod;
  }
  else {
   res=(res+tmp*mul%mod*iv[x][i]%mod*inv[i]%mod*inv[n+1-i]%mod)%mod;
  }
 }return (res+mod)%mod;
}
int main() {
// freopen("fuck.in","r",stdin);
 ios::sync_with_stdio(false);
 cin>>n;fac[0]=1;for(int i=1;i<=n+1;i++)fac[i]=fac[i-1]*i%mod;
 inv[n+1]=fpow(fac[n+1],mod-2);for(int i=n+1;i>=1;i--)inv[i-1]=inv[i]*i%mod;
 for(int i=1;i<=n;i++) {
  cin>>a[i]>>b[i],b[i]++;
  lsh[++cnt]=a[i],lsh[++cnt]=b[i];
 }sort(lsh+1,lsh+1+cnt),cnt=unique(lsh+1,lsh+1+cnt)-lsh-1;
 for(int i=1;i<=n;i++) {
  a[i]=lower_bound(lsh+1,lsh+1+cnt,a[i])-lsh;
  b[i]=lower_bound(lsh+1,lsh+1+cnt,b[i])-lsh;
 }
 for(int i=1;i<cnt;i++) {
     for(int j=0;j<=n+1;j++) {
      iv[i][j]=fpow(lsh[i+1]-lsh[i]-1-j,mod-2);
  }
 }
 g[0][0]=1;for(int i=0;i<=cnt;i++)sum[i]=1;
 for(int i=1;i<=n;i++) {
  for(int j=b[i];j<=cnt;j++) val[j]=sum[j]-sum[j-1];
  for(int j=a[i];j<b[i];j++) {
   add(g[j][0],sum[j-1]);
   for(int k=1;k<=n+1;k++) add(g[j][k],g[j][k-1]);
   val[j]=LG(j,lsh[j+1]-lsh[j]-1);
  }for(int j=a[i];j<b[i];j++)sum[j]=(sum[j-1]+val[j])%mod;
  for(int j=b[i];j<=cnt;j++)sum[j]=(sum[j-1]+val[j])%mod;
 }cout<<(sum[cnt]+mod-1)%mod;
}

[省选联考 2022] 填树

这题挺烦的。

把 $a_i$ , $a_i-K$ , $b_i$ , $b_i-K$ 塞进去离散化，很容易看出这是一个分段多项式。对每一段分别求答案，实现时还是对前缀和插值。

树 $d p$ 求链的时候在 $L C A$ 处统计答案。第一个问题的多项式次数是 $n + 1$ ，第二个问题的多项式次数是 $n + 2$ 。

复杂度 $O(n^3)$ 。

细节贼多。恶心死我了。

#include<bits/stdc++.h>
#define ll long long
#define fi first
#define se second
#define pb push_back
using namespace std;
const int mod=1e9+7;
int n,K,ct;
ll l[805],r[805],lsh[805],sum,inv2;
ll inv[805],fac[805],dp[805],dp2[805];
vector<int>g[805];
ll fpow(ll x,ll y){
 ll z(1);for(;y;y>>=1){
  if(y&1)z=z*x%mod;x=x*x%mod;
 }return z;
}
void init(int n){
 fac[0]=1,inv2=fpow(2,mod-2);for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i%mod;
 inv[n]=fpow(fac[n],mod-2);for(int i=n;i>=1;i--)inv[i-1]=inv[i]*i%mod;
}
void dfs(int u,ll x,ll K,int topf){
 ll L=max(x,l[u]),R=min(x+K,r[u]-1);
 if(L<=R)dp[u]=R-L+1,sum=(sum+R-L+1)%mod;
 for(auto v:g[u]){
  if(v!=topf){
   dfs(v,x,K,u);
   if(L<=R){
    sum=(sum+dp[u]*dp[v]%mod)%mod;
    dp[u]=(dp[u]+(R-L+1)*dp[v]%mod)%mod;
   }
  }
 }
}
ll Sum(ll x){
 return x*(x+1)%mod*inv2%mod;
} 
void dfs2(int u,ll x,ll K,int topf){
 ll L=max(x,l[u]),R=min(x+K,r[u]-1);
 if(L<=R)dp[u]=R-L+1,dp2[u]=Sum(R)-Sum(L-1),sum=(sum+Sum(R)-Sum(L-1))%mod;
 for(auto v:g[u]){
  if(v!=topf){
   dfs2(v,x,K,u);
   if(L<=R){
    sum=(sum+dp2[u]*dp[v]%mod+dp2[v]*dp[u]%mod)%mod;
    dp[u]=(dp[u]+(R-L+1)*dp[v]%mod)%mod;
    dp2[u]=(dp2[u]+dp2[v]*(R-L+1)%mod+dp[v]*(Sum(R)-Sum(L-1))%mod)%mod;
   }
  }
 }
}
ll calc(ll x,ll K){
 memset(dp,0,sizeof dp),sum=0,dfs(1,x,K,0);
 return sum;
}
ll calc2(ll x,ll K){
 memset(dp,0,sizeof dp),memset(dp2,0,sizeof dp2),sum=0,dfs2(1,x,K,0); '<<sum<<endl;
 return sum;
}
ll solve(ll K){
 ll res=0;
 for(int i=1;i<ct;i++){
  int cnt=lsh[i+1]-lsh[i];
  if(cnt<=n+2){
      for(int j=0;j<cnt;j++)res=(res+calc(lsh[i]+j,K))%mod;
  }
  else{
   ll tmp=0,mul=1;for(int j=0;j<=n+1;j++)mul=mul*(cnt-1-j)%mod;
   for(int j=0;j<=n+1;j++){
    tmp=(tmp+calc(lsh[i]+j,K))%mod;
    if(n+1-j&1)res=(res-tmp*mul%mod*fpow(cnt-1-j,mod-2)%mod*inv[n+1-j]%mod*inv[j]%mod)%mod;
    else res=(res+tmp*mul%mod*fpow(cnt-1-j,mod-2)%mod*inv[n+1-j]%mod*inv[j]%mod)%mod;
   }
  }
 }return (res+mod)%mod;
}
ll solve2(ll K){
 ll res=0;
 for(int i=1;i<ct;i++){
  int cnt=lsh[i+1]-lsh[i];
  if(cnt<=n+3){
   for(int j=0;j<cnt;j++)res=(res+calc2(lsh[i]+j,K))%mod;
  }
  else {
   ll tmp=0,mul=1;for(int j=0;j<=n+2;j++)mul=mul*(cnt-1-j)%mod;
   for(int j=0;j<=n+2;j++){
    tmp=(tmp+calc2(lsh[i]+j,K))%mod;
    if(n+2-j&1)res=(res-tmp*mul%mod*fpow(cnt-1-j,mod-2)%mod*inv[n+2-j]%mod*inv[j]%mod)%mod;
    else res=(res+tmp*mul%mod*fpow(cnt-1-j,mod-2)%mod*inv[n+2-j]%mod*inv[j]%mod)%mod;
   }
  }
 }return (res+mod)%mod;
} 
signed main(){
 freopen("tree.in","r",stdin);
 freopen("tree.out","w",stdout);
 ios::sync_with_stdio(false);
 cin.tie(0),cout.tie(0);
 cin>>n>>K,init(2*n);for(int i=1;i<=n;i++){
  cin>>l[i]>>r[i],r[i]++,lsh[++ct]=l[i],lsh[++ct]=r[i];
  if(l[i]>K)lsh[++ct]=l[i]-K;if(r[i]>K)lsh[++ct]=r[i]-K;
 }lsh[++ct]=1;sort(lsh+1,lsh+1+ct),ct=unique(lsh+1,lsh+1+ct)-lsh-1;
 for(int i=1;i<n;i++){
  int u,v;cin>>u>>v,g[u].pb(v),g[v].pb(u);
 }
 cout<<(solve(K)-solve(K-1)+calc(1,K-1)+mod)%mod<<endl;
 cout<<((solve2(K)-solve2(K-1)+calc2(1,K-1)+mod)%mod+mod)%mod; 
}

[NOI2019] 机器人

玄学 $d p$ 优化。。。

不难看出这是一个区间 $d p$ /kk

记 $f_{l,r,M}$ 表示区间 $[l, r]$ 中每个位置高度均不大于 $M$ 的方案数。

注意到从区间最大高度出发能移动到整个区间，并且将整个区间分成了两个独立的子问题，那么 $f_{l,r,M}=\sum_{k=1}^Mf_{l,mid-1,k}\times f_{mid+1,r,k-1}[a_{mid}\le k\le b_{mid}]$ 。答案是 $f_{1,n,\max b_i}$ 。时间复杂度 $O(n^2A)$ 。