拉格朗日插值小记
关于拉格朗日插值：我只会最简单的形式喵。
就是给 \(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)\) 的时间内求出来。所以求 \(\prod \limits_{j\not = 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);
}
浙公网安备 33010602011771号