园龄：5年1个月粉丝：14 关注：5

拉格朗日插值小记

关于拉格朗日插值：我只会最简单的形式喵。

就是给 $n$ 个点值，就能在 $O(n^2)$ 的时间复杂度内求出当 $x=a$ 的时候的值。

其标准形式是： $\displaystyle \sum_{i=1}^n y_i \prod _{j\not =i}^n \frac{a-x_j}{x_i-x_j}$

然后它更高深的应用我暂时不会喵。

但是也许能用在一些小清新 DP 式的优化。

主要在于如果发现式子有一个消不掉的巨大数字，并且能够确定整个结果是和它有关的多项式，就能先求出小值域下的结果，然后再插值回去。

CF995F Cowmpany Cowmpensation

题意： $n$ 个点的树，给每个点分配权值 $[1,d]$ ，父亲的权值必须大于儿子，求方案数。

$f_{i,j}$ 表示整棵树的权值在 $1\sim j$ 的方案数。

$\displaystyle f_{p,i}=f_{p,i-1}+\prod_{v\in p} f_{v,i}$

所有叶子的 $f_{p,i}=i$

感性理解一下，当在叶子上时， $f_{p,i}$ 是一个关于 $i$ 的一次函数，也就是多项式的项数为 $2$ 。而考虑叶子们的父亲，就是把儿子们乘起来，变成了次数为叶子数加一的多项式。也就是说，到根的时候，最多是一个关于 $i$ ，项数为 $n$ 次的多项式。这样可以 $n^2$ 求出一部分点值，再拉格朗日插值把 $d$ 带进去。

#include <bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
typedef long long ll;
ll ksm(ll x,ll y)
{
	ll ret=1;
	while(y)
	{
		if(y&1) ret=ret*x%mod;
		x=x*x%mod; y>>=1;
	}
	return ret;
}
int n,D;
const int maxn=6005;
int to[maxn],nxt[maxn],head[maxn],num;
void add(int x,int y){num++;to[num]=y;nxt[num]=head[x];head[x]=num;}
ll f[3002][3002],a[maxn],b[maxn];
void dfs(int p,int fa)
{
	for(int i=head[p];i;i=nxt[i])
	{
		if(to[i]==fa) continue;
		dfs(to[i],p); 
	}
	for(int j=1;j<=n;j++)
	{
		ll ret=1;
		for(int i=head[p];i;i=nxt[i])
			if(to[i]!=fa) ret=ret*f[to[i]][j]%mod;
		f[p][j]=(f[p][j-1]+ret)%mod;
	}
}
ll LA(int n,int D,ll *a,ll *b)
{
	ll ret=0;
	for(int i=0;i<=n;i++) 
	{
		ll x1=1,x2=1;
		for(int j=0;j<=n;j++)
		{
			if(i==j) continue;
			x1=x1*(D-a[j])%mod; x2=x2*(a[i]-a[j])%mod;
		}
		ret=(ret+x1*ksm(x2,mod-2)%mod*b[i]%mod)%mod;
	}
	return (ret%mod+mod)%mod;
}
int main()
{
	scanf("%d%d",&n,&D);
	for(int i=2;i<=n;i++)
	{
		int f; scanf("%d",&f);
		add(f,i);
	}
	dfs(1,0);
	for(int i=0;i<=n;i++) a[i]=i,b[i]=f[1][i];
	printf("%lld\n",LA(n,D,a,b));
}

[省选联考 2022] 填树

首先题面善良地告诉你被染色的点是一条链。

那么很容易写出第一问的转移方程。

枚举最小值为 $L$ ，当前情况下的答案为 $A_{L}$ ，设 $f_{p}$ 表示以 $p$ 为根，所有符合条件的链的数量， $\displaystyle f_{i}=(b_i-a_i+1)\sum_{v \in son} f_v$ 。

$a_i,b_i$ 表示在整棵树上最小值为 $L$ 的情况下， $i$ 节点能到达的值域。 $a_i=\max\{l_i,L\},b_i=\min\{r_i,L+K-1\}$ ，如果 $a_i>b_i$ 则 $f_i=0$ 。

然后在每个点上合并儿子的时候算一下当前链和前面已经合并的 $f_i$ 匹配一下，统计到总答案里。

那么 $A_L=f_{1}-A_{L-1}$ 。这样减一下是因为 $f_{i}$ 可能统计重复，就是没有最小值取到 $L$ 的情况。

然后我们猜测 $A_{L}$ 可能是关于 $L$ 的，项数为 $n+1$ 的多项式。因为可以发现最终答案就是两条链 $b_i-a_i+1$ 乘到一起然后在 $lca$ 处加起来的，而 $b_i-a_i+1$ 是关于 $L$ 的一次函数，假设 $b_i-a_i+1=len_i$ ， $len_i$ 与 $L$ 有如下关系：

$len_i=\begin{cases}0 &\text{if } L+K-1<l_i&\text{or } L>r_i\\L+K-l_i &\text{if } L+K-1\geq l_i &\text{and } l_i \leq L\leq r_i\\ K &\text{if } L\geq l_i &\text{and } L+K-1 \leq r_i\\r_i-L+1 &\text{if }l_i\leq L\leq r_i&\text{and } L+K-1\geq r_i \end{cases}$

上面就是为了向你展示 $L$ 与每个点乘进去的权值的一次函数关系。

所以 $A_L$ 是关于 $L$ 的 $n+1$ 项的多项式，不过每一项系数可能随 $L$ 的大小改变。

但是如果我们把它想象成函数图像，就会明白能够使它发生转折的只有大约 $O(n)$ 个点，剩下的时候都是连续的函数。那么我们就可以考虑拉格朗日插值：我们找到所有可能使函数转折横坐标，将它排序之后一段一段地考虑当前的函数。假如当前这段的函数图像覆盖的范围不太大（小于等于 $n$ ），那么我们完全可以把这段里所有的函数值全求出来加和；而如果覆盖的范围很大，我们求出这一段的前 $n$ 个点值，然后做一个前缀和，再用结尾的横坐标做拉格朗日插值，得到的数字相当于这一段的和。对于每一段的函数都求出来所有的和，那么总的和就是答案啦。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=1030;
const int mod=1e9+7;
ll ksm(ll x,ll y)
{
	ll ret=1;
	while(y)
	{
		if(y&1) ret=ret*x%mod;
		x=x*x%mod; y>>=1;
	}
	return ret;
}
int n,K;
int to[maxn],nxt[maxn],head[maxn],num;
void add(int x,int y){num++; nxt[num]=head[x]; to[num]=y;head[x]=num;}
int l[maxn],r[maxn];int L,R;
ll f[maxn],g[maxn];ll tmpf,tmpg,ansf,ansg;
void dfs(int p,int fa)
{
	ll ll=max(l[p],L),rr=min(r[p],R);
	if(ll<=rr) f[p]=(rr-ll+1),g[p]=((ll+rr)*(rr-ll+1)/2)%mod;
	else f[p]=0,g[p]=0;
	int tmp1=ll<=rr?rr-ll+1:0;
	int tmp2=ll<=rr?((ll+rr)*(rr-ll+1)/2)%mod:0;
	tmpf=(tmpf+f[p])%mod;
	tmpg=(tmpg+g[p])%mod;
	for(int i=head[p];i;i=nxt[i])
	{
		if(to[i]==fa) continue;
		dfs(to[i],p);
		int v=to[i];
		tmpf=(f[p]*f[v]%mod+tmpf)%mod;
		tmpg=(f[p]*g[v]%mod+g[p]*f[v]%mod+tmpg)%mod;
		f[p]=(f[p]+f[v]*tmp1%mod)%mod;
		g[p]=(g[p]+g[v]*tmp1%mod+f[v]*tmp2%mod)%mod;
	}
}
void calc(int w)
{
	//printf("%d %d\n",L,R);
	tmpf=0; tmpg=0;
	dfs(1,0);
	ansf=(ansf+w*tmpf%mod)%mod; ansg=(ansg+w*tmpg%mod)%mod;
}
ll la(int n,ll val,ll a[],ll b[])
{
	ll ret=0;
	for(int i=0;i<=n;i++)
	{
		ll xx1=1,xx2=1;
		for(int j=0;j<=n;j++)
		{
			if(i==j) continue;
			xx1=xx1*(val-a[j])%mod;
			xx2=xx2*(a[i]-a[j])%mod;
		}
		ret=(ret+xx1*ksm(xx2,mod-2)%mod*b[i]%mod)%mod;
	}
	return ret;
}
int pos[maxn],m;
ll a[maxn],b[maxn],c[maxn];
signed main()
{
//	freopen("tree.in","r",stdin);
//	freopen("tree.out","w",stdout);
	scanf("%d%d",&n,&K);
	for(int i=1;i<=n;i++)
	{
		scanf("%d%d",&l[i],&r[i]);
		pos[++m]=l[i],pos[++m]=max(0,l[i]-K),pos[++m]=r[i],pos[++m]=max(r[i]-K,0);
		pos[0]=max(pos[0],r[i]+1);
	}
	sort(pos+0,pos+m+1);m=unique(pos,pos+m+1)-pos;
	for(int i=1;i<n;i++)
	{
		int x,y; scanf("%d%d",&x,&y);
		add(x,y); add(y,x);
	}
	for(int i=0;i<m;i++)
	{
		int j;
		for(j=0;j<n+2;j++)
		{
			if(pos[i]+j==pos[i+1]) break;
			L=pos[i]+j,R=pos[i]+j+K;
			calc(1); L++; calc(-1);
			a[j]=L-1,b[j]=ansf,c[j]=ansg;
		}
		if(pos[i]+j<pos[i+1])
		{
			ansf=la(j-1,pos[i+1]-1,a,b); 
			ansg=la(j-1,pos[i+1]-1,a,c);
		}
	}
	printf("%lld\n%lld\n",(ansf%mod+mod)%mod,(ansg%mod+mod)%mod);
}

P3270 [JLOI2016]成绩比较

纪念一下我自己切出来的题 qwq

我们看到了恰好有 $k$ 个人被碾压，考虑容斥。就是设 $g(i)$ 为钦定 $i$ 个人被碾压，剩下的随意的方案数。

答案就是 $\displaystyle \sum_{i=k}^n (-1)^{i-k} \dbinom i k\dbinom {n}{i} g(i)$

然后就是考虑这个 $g(i)$ 怎么求。

$\displaystyle g(i)=\prod_{j=1}^m \dbinom {n-1-i}{R_j-1} \sum_{k=1}^{U_j} k^{n-R_j} (U_j - k)^{R_j-1}$

这个意思就是先从不被碾压的人选出 $R_j-1$ 个在 D 神前面，然后枚举 D 神的分数，算方案数，然后不同科目乘起来就行。

观察到后面的和式很像数列 $k$ 次方和，它自然是一个 $n$ 次多项式 qwq

然后后面的拉插 $m$ 次预处理一下就行。

/*
　　　　　／＞　 フ   
　　　　　| 　_　 _|	
　 　　　／`ミ _x 彡	
　　 　 /　　　 　 |	
　　　 /　 ヽ　　 /	
　／￣|　　 |　|　|	
　| (￣ヽ＿_ヽ_)_)	
　＼二つ
*/
#include <bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
typedef long long ll;
ll ksm(ll x,ll y)
{
	ll ret=1;
	while(y)
	{
		if(y&1) ret=ret*x%mod;
		x=x*x%mod; y>>=1;
	}
	return ret;
}
int n,m,K;
int U[103],R[103];ll b[104],a[104],d[104],fac[103],inv[192],ans;
ll LA(ll n,ll m,ll *a,ll *b)
{
	ll ret=0;
	for(int i=1;i<=m;i++) 
	{
		ll x1=1,x2=1;
		for(int j=1;j<=m;j++)
		{
			if(i==j) continue;
			x1=x1*(n-a[j])%mod; x2=x2*(a[i]-a[j])%mod;
		}
		ret=(ret+x1*ksm(x2,mod-2)%mod*b[i]%mod)%mod;
	}
	return ret;
}
ll C(ll n,ll m)
{
	if(n<0||m<0||n<m) return 0;
	return fac[n]*inv[m]%mod*inv[n-m]%mod;
}
int main()
{
	//freopen("p.in","r",stdin);
	scanf("%d%d%d",&n,&m,&K);
	for(int i=1;i<=m;i++) scanf("%d",&U[i]);
	for(int i=1;i<=m;i++) scanf("%d",&R[i]);
	int up=0x3f3f3f3f;
	for(int i=1;i<=m;i++)
	{
		ll tmp=0;
		for(int j=1;j<=n+1;j++)
		{
			tmp=(tmp+ksm(j,n-R[i])*ksm(U[i]-j,R[i]-1)%mod)%mod;
			b[j]=tmp; a[j]=j;
		}
		d[i]=LA(U[i],n+1,a,b); 
		//fprintf(stderr,"%lld\n",d[i]);
		up=min(up,n-R[i]);
	} 
	fac[0]=1;inv[0]=1;	
	for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod,inv[i]=ksm(fac[i],mod-2);
	for(int i=K;i<=up;i++)
	{
		ll tmp=((i-K)&1?(mod-1):1)*C(n-1,i)%mod*C(i,K)%mod;
		//fprintf(stderr,"%lld\n",tmp);
		for(int j=1;j<=m;j++)
		{
			tmp=tmp*C(n-i-1,R[j]-1)%mod*d[j]%mod;
		}
		//fprintf(stderr,"%lld\n*******\n",tmp);
		ans=(ans+tmp)%mod;
	} 
	printf("%lld\n",(ans%mod+mod)%mod);
}

P7116 [NOIP2020] 微信步数

这题卡常，卡着卡着，拉插就消失了。

首先，如果走完 $n$ 步之后没有坐标发生变化且存在期间没有走出界的起点，那么就一定走不出去。

我们考虑每个时间产生的贡献。显然就是这个时间点还剩下的点的个数。

我们把维度分开考虑，对于每一个时间点 $j$ ，假如维度 $i$ 能保持在场地内的数量为 $a_{i,j}$ ，那么当前时间点能产生的贡献为 $\prod_{i=1}^k a_{i,j}$ 。

那么我们考虑计算每一个维度向左或向右走的最远距离 $l_{i,j},r_{i,j}$ ，所以 $a_{i,j}=w_i-r_{i,j}+l_{i,j}$ 。我们一步一步枚举，如果某一步 $a_{i,j}\leq 0$ ，那么就结束，否则统计所有的乘积，就是答案。期望得分 30pts。

这样很慢的原因是每一轮都要一步一步枚举。而可以发现，从第二轮开始，每个维度增加或减少的坐标都相同。所以我们枚举步数 $x$ ，计算出还能走到这一步多少次。每一次就是一轮过去了。而这个的贡献就是 $\displaystyle \sum_{i=0}^{T} \prod_{j=1}^{k} (c_j-i\times d_j-f_{x,j})$ 。其中 $c_j$ 表示这一维度第一轮过后还剩多少人， $d_j$ 表示这一维度每轮要离开多少人， $f_{i,j}$ 表示再走 $i$ 步 $j$ 维度消失多少个。然后可以发现里面的乘积形式的东西是个与 $i$ 有关的 $k+1$ 次多项式。然后就可以拉插了。

但是问题就是拉插写不好就是 $O(nm^2\log)$ 的太慢。其实我们发现 $T$ 不会超过 $\max w_i$ ，所以我们预处理 $\sum i^k$ ，然后把里面的多项式暴力乘在一起，然后直接算就行。

这样就没拉插什么事了 qwq

P4463 [集训队互测 2012] calc

这题先咕着，因为我不会证明它是个多项式...

P5469 [NOI2019] 机器人

一开始看错题自闭了好久

一个挺巧妙的区间 DP：由于区间内最高点左边的点不可能到右边来，所以设 $f[i][j][k]$ 表示 $i\sim j$ 区间内最大值为 $k$ 的方案数。转移就是枚举可以的起点 $mid$ , $f[i][j][k]=\sum f[i][mid-1][k]\times f[mid+1][r][k-1]$ , $mid$ 必须是一个合法的转移点，也就是 $abs(mid-l-r+mid)\leq 2$

然后发现 $mid$ 合法点不太多，~~打表发现不超过3000个~~ 所以此时复杂度为 $O(3000mx)$ ，mx 为 $A_i,B_i$ 的最大值。

我们观察 $A_i=1,B_i=1e9$ ，然后我们发现这个式子中 $f[l][l][x]=1=x^0$ ，所以最后的式子一定是关于 $x$ 的不超过 $n$ 次多项式。那么我们可以直接拉格朗日插值。

然后我们考虑还像填树那样，把所有分界点都找出来，这样 $f[i][i][x]$ 保证是 $x^0$ ，然后每次都求出 $n+1$ 个点值，拉插求前缀和即可。

然后是痛苦卡常时间。

upd on 2023/05/26

补一下关于如何拉格朗日插值系数。

首先考虑拉插的原理：我们知道了 $n$ 个点的点值，然后需要构造一个多项式。一种构造方案是构造 $n$ 个多项式，对于多项式 $f_i(X)$ ，当 $X=x_i$ 时等于 $y_i$ ，其余时候都等于 $0$ 。那么可以显然的构造一个东西就是 $\prod\limits_{j\not =i } (X-x_j)$ ，这个东西当 $X\not = x_i$ 时一定等于 $0$ ，但是当 $X=x_i$ 时则不一定为 $y_i$ ，而是一个别的数字。所以我们将这个数字除掉然后再乘上 $y_i$ 就行了。

所以 $f_i(X)=y_i\frac {\prod\limits_{j\not = i }(X-x_j)}{\prod \limits_{j\not = i}(x_i-x_j)}$

$f_i(X)=y_i\prod \limits_{j\not = i} \frac {X-x_j}{x_i-x_j}$

然后正常的拉插就将 $X$ 带入求和就可以了。

那么我们想要求系数怎么办。

我们可以发现 $y_i$ 和 $\prod \limits _{j\not =i } \frac {1}{x_i-x_j}$ 都是常数，可以直接求。

剩下的一部分是一个多项式，如果对每个 $i$ 都直接展开求系数是 $O(n^3)$ 的。所以我们考虑先 $O(n^2)$ 求出 $\prod (X-x_i)$ 的系数，然后对每个 $i$ 除以 $X-x_i$ ，因为需要除的多项式比较特殊，可以在的时间内求出来。所以求的系数是 $O(n^2)$ 的。

然后求出的 $n$ 个多项式对应项系数相乘就行。

#include <bits/stdc++.h>
using namespace std;
const int maxn=2e4+3;
typedef long long ll;
const int mod=998244353;
ll a[maxn],b[maxn];int n,K;
ll x[maxn],y[maxn];ll f[maxn],g[maxn],c[maxn],d[maxn];
ll ksm(ll x,ll y)
{
    ll ret=1;
    while(y)
    {
        if(y&1) ret=ret*x%mod; 
        x=x*x%mod; y>>=1;
    }
    return ret;
}
void solve()
{
    for(int i=1;i<=n;i++)
    {
        g[i]=y[i];
        for(int j=1;j<=n;j++)
        {
            if(i==j) continue;
            g[i]=g[i]*ksm((x[i]-x[j]+mod)%mod,mod-2)%mod;
        }
    }
    a[0]=1;
    for(int i=1;i<=n;i++)
    {
        for(int j=i;j>=1;j--)
            a[j]=(a[j]*(mod-x[i])%mod+a[j-1])%mod;
        a[0]=a[0]*(mod-x[i])%mod;
    }
    
    for(int i=1;i<=n;i++)
    {
        for(int j=0;j<=n;j++) c[j]=a[j];
        for(int j=n;j>=1;j--)
        {
            d[j-1]=c[j];
            c[j-1]=(c[j-1]+c[j]*x[i]%mod)%mod;
        }
    //     for(int j=0;j<=n;j++) cerr<<d[j]<<" ";
    // cerr<<endl;
        for(int j=0;j<n;j++) f[j]=(f[j]+d[j]*g[i]%mod)%mod;
    }
}
int main()
{
    //freopen("p.in","r",stdin);
    scanf("%d%d",&n,&K);
    for(int i=1;i<=n;i++)
    {
        scanf("%lld%lld",&x[i],&y[i]);
    }
    solve();
    ll tmp=1,ans=0;
    for(int i=0;i<n;i++) ans=(ans+tmp*f[i]%mod)%mod,tmp=tmp*K%mod;
    printf("%lld\n",ans); 
}

上一篇2.15模拟赛

下一篇单位根反演小记

本文作者：cc0000

本文链接：https://www.cnblogs.com/cc0000/p/17133059.html

posted @ 2023-02-18 16:59 cc0000 阅读(157) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

cc0000

一杯茶，一包烟，一道破题调一天

拉格朗日插值小记

CF995F Cowmpany Cowmpensation

[省选联考 2022] 填树

P3270 [JLOI2016]成绩比较

P7116 [NOIP2020] 微信步数

P4463 [集训队互测 2012] calc

P5469 [NOI2019] 机器人

公告

常用链接

随笔分类

随笔档案

阅读排行榜

评论排行榜

最新评论