拉格朗日插值
拉格朗日插值
为啥都开始学数学了？睡着了。呼呼呼(～ o ～)~zZ
定义
众所周知，平面上 \(n+1\) 个点可以唯一确定一个 \(n\) 次多项式。
已知 \(n+1\) 个点 \((x_i,y_i)\)，拉格朗日插值支持在 \(O(n^2)\) 的复杂度内得到一个 \(n\) 次多项式的各项系数。
在 \(x_i\) 为连续段时，可以在 \(O(n)\) 的复杂度内求另一个点值。
公式如下：
\[F(x)=\sum_{i=1}^{n} y_i \prod_{j \ne i} \frac{x-x_j}{x_i-x_j}
\]
这是 \(O(n^2)\) 求点值，直接模拟：
code
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define fi first
#define se second
const int N = 2e3+5,mod = 998244353;
inline LL qpow(LL a,int b)
{
LL res=1;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod; b>>=1;
}
return res;
}
int n;
LL k,inv[N],x[N],y[N],ans;
int main()
{
// freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
scanf("%d%lld",&n,&k);
for(int i=1;i<=n;i++) scanf("%lld%lld",&x[i],&y[i]);
for(int i=1;i<=n;i++)
{
LL t1=1,t2=1;
for(int j=1;j<=n;j++) if(j!=i)
{
t1=(t1*(k-x[j])%mod+mod)%mod; t2=(t2*(x[i]-x[j])%mod+mod)%mod;
}
ans=(ans+t1*y[i]%mod*qpow(t2,mod-2)%mod)%mod;
}
printf("%lld\n",ans);
return 0;
}
这是 \(O(n)\) 求点值，因为 \(\prod x-x_j\) 可以用前缀积和后缀积优化，\(\prod x_i-x_j\) 就是两部分阶乘：
code
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define fi first
#define se second
const int N = 1e6+5,mod = 1e9+7;
int n,k;
LL y[N],d[N],c[N],fac[N];
inline LL qpow(LL a,int b)
{
LL res=1;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod; b>>=1;
}
return res;
}
int main()
{
// freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
scanf("%d%d",&n,&k);
for(int i=1;i<=k+2;i++) y[i]=(y[i-1]+qpow(i,k))%mod;
if(n<=k+2) return printf("%lld\n",y[n]),0;
d[0]=fac[0]=fac[1]=c[k+3]=1;
for(int i=2;i<=k+2;i++) fac[i]=(mod-mod/i)*fac[mod%i]%mod;
for(int i=1;i<=k+2;i++) d[i]=d[i-1]*(n-i)%mod,fac[i]=fac[i-1]*fac[i]%mod;
for(int i=k+2;i>=1;i--) c[i]=c[i+1]*(n-i)%mod;
LL ans=0;
for(int i=1;i<=k+2;i++)
{
LL t1=c[i+1]*d[i-1]%mod,t2=fac[i-1]*((k+2-i)&1?(mod-fac[k+2-i]):(fac[k+2-i]))%mod;
ans=(ans+t1*t2%mod*y[i]%mod)%mod;
}
printf("%lld\n",ans);
return 0;
}
插多项式系数，先把常数 \(y_i\prod \frac{1}{x_i-x_j}\) 预处理，然后计算 \(\prod x-x_j\)，
递推转移，设 \(\prod x-x_j\) 的第 \(i\) 次项系数是 \(f_i\)，那么有 \(f_i=f_{i-1}-x_j f_{i}\)，注意后效性。
至于 \(j \ne i\) 的限制可以退背包做，令 \(ff_{i}\) 为退完的数组，\(ff_i=(f_i-ff_{i-1}) \times x_j^{-1}\)。
code
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define fi first
#define se second
const int N = 2e3+5,mod = 998244353;
inline LL qpow(LL a,int b)
{
LL res=1;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod; b>>=1;
}
return res;
}
int n;
LL k,inv[N],x[N],y[N],ans,g[N],f[N],p[N],ff[N];
int main()
{
// freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
scanf("%d%lld",&n,&k);
for(int i=1;i<=n;i++) scanf("%lld%lld",&x[i],&y[i]);
for(int i=1;i<=n;i++)
{
g[i]=1;
for(int j=1;j<=n;j++) if(i!=j) g[i]=g[i]*(x[i]-x[j]+mod)%mod;
g[i]=qpow(g[i],mod-2)*y[i]%mod;
}
f[0]=1;
for(int i=1;i<=n;i++)
{
for(int j=n-1;j>=1;j--) f[j]=(f[j-1]-f[j]*x[i]%mod+mod)%mod;
f[0]=(mod-f[0]*x[i]%mod)%mod;
}
for(int i=1;i<=n;i++)
{
LL in=qpow((mod-x[i])%mod,mod-2);
ff[0]=f[0]*in%mod; p[0]=(p[0]+g[i]*ff[0])%mod;
for(int j=1;j<n;j++)
{
ff[j]=(f[j]-ff[j-1]+mod)%mod*in%mod;
p[j]=(p[j]+ff[j]*g[i])%mod;
}
}
LL ans=0,tmp=1;
for(int i=0;i<=n;i++) ans=(ans+p[i]*tmp)%mod,tmp=tmp*k%mod;
printf("%lld\n",ans);
return 0;
}
例题
只做了 \(x\) 道题，除了模板好像还是模板？
自然数幂和
\(\sum i^k\) 是一个 \(k+1\) 次多项式，证明找一下递推式是不是就好了。
所以可以直接插，这个用得很多。\(code\) 在上面。
教科书般的亵渎
一开始误以为答案就是一个 \(k+1\) 次多项式，直接插就对了。是我高估拉插了 QwQ
事实证明在你都不知道式子是啥的时候，它不可能是一个确定的多项式。
考虑拆贡献，可以看做对 \(0\) 使用炎拳亵渎。这样每次的贡献就是后面所有的 \(m+1\) 次，减去后面 \(0\) 的贡献。
然后你发现要求 \(\sum_{i=1}^{x} i^k\)，然后直接拉插求这个就好了。
code
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define fi first
#define se second
const int N = 105,mod = 1e9+7;
int T,m;
LL inv[N],n,a[N],pre[N],nxt[N],y[N];
inline LL qpow(LL a,int b)
{
LL res=1;
while(b)
{
if(b&1) res=a*res%mod;
a=a*a%mod; b>>=1;
}
return res;
}
inline LL cal(LL n,int k)
{
LL res=0; k+=2;
pre[0]=nxt[k+1]=1;
for(int i=1;i<=k;i++) y[i]=(y[i-1]+qpow(i,k-2))%mod;
if(n<=k) return y[n];
n%=mod;
for(int i=1;i<=k;i++) pre[i]=pre[i-1]*(n-i)%mod;
for(int i=k;i>=1;i--) nxt[i]=nxt[i+1]*(n-i)%mod;
for(int i=1;i<=k;i++)
{
LL t1=pre[i-1]*nxt[i+1]%mod,t2=inv[i-1]*((k-i)&1?mod-inv[k-i]:inv[k-i])%mod;
res=(res+t1*t2%mod*y[i])%mod;
}
return res;
}
int main()
{
// freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
scanf("%d",&T);
inv[0]=inv[1]=1;
for(int i=2;i<N;i++) inv[i]=(mod-mod/i)*inv[mod%i]%mod;
for(int i=2;i<N;i++) inv[i]=inv[i-1]*inv[i]%mod;
while(T--)
{
scanf("%lld%d",&n,&m);
for(int i=1;i<=m;i++) scanf("%lld",&a[i]); a[++m]=0;
sort(a+1,a+1+m);
LL ans=0;
for(int i=1;i<=m;i++)
{
ans=(ans+cal(n-a[i],m))%mod;
for(int j=i+1;j<=m;j++) ans=(ans+mod-qpow(a[j]-a[i],m))%mod;
}
printf("%lld\n",ans);
}
return 0;
}
XLkxc
经过上一道题的教训，这次不敢直接猜答案就是多项式了，但是这个明显就是多项式啊！是我低估拉插了 QwQ
直接插三次就做完了，复杂度单次 \(O(n^3)\)，稍微推一下式子能有 \(O(n^2)\) 的，好像还能更优？
code
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define fi first
#define se second
const int N = 205,mod = 1234567891;
LL a,n,d;
int k,T;
LL pre[3][N],nxt[3][N],y[N],z[N],w[N],inv[N];
inline LL qpow(LL a,int b)
{
LL res=1;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod; b>>=1;
}
return res;
}
inline LL S(LL n,int k)
{
LL res=0;
k+=2; pre[0][0]=nxt[0][k+1]=1;
for(int i=1;i<=k;i++) pre[0][i]=pre[0][i-1]*(n-i)%mod,y[i]=(y[i-1]+qpow(i,k-2))%mod;
for(int i=k;i>=1;i--) nxt[0][i]=nxt[0][i+1]*(n-i)%mod;
// printf("%lld %d ",n,k);
if(n<=k) return /*printf("%lld\n",y[n]),*/y[n];
for(int i=1;i<=k;i++)
{
LL t1=pre[0][i-1]*nxt[0][i+1]%mod,t2=inv[i-1]*((k-i)&1?mod-inv[k-i]:inv[k-i])%mod;
// printf("%lld %lld\n",t1,t2);
res=(res+t1*t2%mod*y[i])%mod;
}
return /*printf("%lld\n",res),*/res;
}
inline LL F(LL n,int k)
{
LL res=0;
k+=3; pre[1][0]=nxt[1][k+1]=1;
for(int i=1;i<=k;i++) pre[1][i]=pre[1][i-1]*(n%mod-i+mod)%mod,z[i]=(z[i-1]+S(i,k-3))%mod;
for(int i=k;i>=1;i--) nxt[1][i]=nxt[1][i+1]*(n%mod-i+mod)%mod;
if(n<=k) return z[n];
for(int i=1;i<=k;i++)
{
LL t1=pre[1][i-1]*nxt[1][i+1]%mod,t2=inv[i-1]*((k-i)&1?mod-inv[k-i]:inv[k-i])%mod;
res=(res+t1*t2%mod*z[i])%mod;
}
return res;
}
inline LL E(LL n,int k)
{
LL res=0; k+=4;
pre[2][0]=nxt[2][k+1]=1; w[0]=F(a,k-4);
for(int i=1;i<=k;i++) pre[2][i]=pre[2][i-1]*(n-i)%mod,w[i]=(w[i-1]+F(a+i*d,k-4))%mod;
for(int i=k;i>=1;i--) nxt[2][i]=nxt[2][i+1]*(n-i)%mod;
if(n<=k) return w[n];
for(int i=1;i<=k;i++)
{
LL t1=pre[2][i-1]*nxt[2][i+1]%mod,t2=inv[i-1]*((k-i)&1?mod-inv[k-i]:inv[k-i])%mod;
res=(res+t1*t2%mod*w[i])%mod;
}
return res;
}
int main()
{
// freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
scanf("%d",&T);
inv[0]=inv[1]=1;
for(int i=2;i<N;i++) inv[i]=(mod-mod/i)*inv[mod%i]%mod;
for(int i=2;i<N;i++) inv[i]=inv[i-1]*inv[i]%mod;
while(T--)
{
scanf("%d%lld%lld%lld",&k,&a,&n,&d);
printf("%lld\n",E(n,k));
// for(int i=1;i<=k;i++) printf("%lld ",y[i]); putchar('\n');
}
return 0;
}
成绩比较
抽象题目，有很多个限制。
首先发现有恰好碾压 k 个人。所以考虑二项式反演。
不要怕写暴力，虽然式子很吓人但是很好想。
设 \(f_x\) 为恰好碾压 \(x\) 个人，\(g_x\) 为钦定碾压 \(x\) 个人。
有 \(f_k = \sum_{i=k}^{n-1} \binom{i}{k} (-1)^{i-k} g_i\)。
考虑计算 \(g_k\)。
\[g_k=\binom{n-1}{k} \prod_{i=1}^{m}\ \sum_{j=1}^{U_i}\ \binom{n-1-k}{R_i-1}\ j^{n-R_i}\ (U_i-j)^{R_i-1}
\]
先选 \(k\) 个人，枚举科目，枚举 B 神的分数，没被碾压有 \(\binom{n-1-k}{R_i-1}\) 个人，选 \(R_i-1\) 个比 B 神高，然后枚举分数。
化简：
\[g_k=\binom{n-1}{k} \prod_{i=1}^{m}\ \binom{n-1-k}{R_i-1} \sum_{j=1}^{U_i}\ \ j^{n-R_i}\ (U_i-j)^{R_i-1}
\]
考虑如何计算 \(\sum_{j=1}^{U_i}\ \ j^{n-R_i}\ (U_i-j)^{R_i-1}\)，然后你发现一下是 \(n\) 次的多项式，直接插。
code
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define fi first
#define se second
#define LL long long
const int N = 105,mod = 1e9+7;
int n,m,k,u[N],r[N];
LL fac[N],inv[N],h[N],g;
inline LL C (int x,int y) {return x>=y?fac[x]*inv[y]%mod*inv[x-y]%mod:0;}
LL pre[N],nxt[N],y[N];
inline LL qpow(LL a,int b)
{
LL res=1;
while(b)
{
if(b&1) res=res*a%mod;
a=a*a%mod; b>>=1;
}
return res;
}
inline LL cal(int s)
{
int k=n+1; LL res=0;
pre[0]=nxt[k+1]=1;
for(int i=1;i<=k;i++) pre[i]=pre[i-1]*(u[s]-i)%mod;
for(int i=k;i>=1;i--) nxt[i]=nxt[i+1]*(u[s]-i)%mod;
for(int i=1;i<=k;i++) y[i]=(y[i-1]+qpow(i,n-r[s])*qpow(u[s]-i,r[s]-1)%mod)%mod;
if(u[s]<=k) return y[u[s]];
for(int i=1;i<=k;i++)
{
LL t1=pre[i-1]*nxt[i+1]%mod,t2=inv[i-1]*((k-i)&1?mod-inv[k-i]:inv[k-i])%mod;
res=(res+t1*t2%mod*y[i])%mod;
}
// printf("%lld\n",res);
return res;
}
int main()
{
// freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
scanf("%d%d%d",&n,&m,&k);
fac[0]=fac[1]=inv[0]=inv[1]=1;
for(int i=2;i<=n;i++) fac[i]=fac[i-1]*i%mod,inv[i]=(mod-mod/i)*inv[mod%i]%mod;
for(int i=2;i<=n;i++) inv[i]=inv[i-1]*inv[i]%mod;
for(int i=1;i<=m;i++) scanf("%d",&u[i]);
for(int i=1;i<=m;i++) scanf("%d",&r[i]);
for(int i=1;i<=m;i++) h[i]=cal(i);
// for(int i=1;i<=m;i++) printf("%lld ",h[i]); putchar('\n');
LL ans=0;
for(int i=k;i<=n-1;i++)
{
g=C(n-1,i);
// printf("%lld$$\n",g);
for(int j=1;j<=m;j++) g=g*C(n-1-i,r[j]-1)%mod*h[j]%mod/*,printf("%lld\n",C(n-1-i,r[j]-1))*/;
ans=(ans+C(i,k)*g%mod*((i-k)&1?-1:1)+mod)%mod;
}
printf("%lld\n",ans);
return 0;
}
calc
子序列问题并不好做，但其实排个序就变成集合选数问题了，最后贡献乘一个阶乘。
dp 是朴素的 \(f_{i,j}\) 表示在 \([1,i]\) 的数中选 \(j\) 个的权值和。
有 \(f_{i,j}=f_{i-1,j-1}\times i+f_{i-1,j}\)。
发现一下这是一个 \(2i+1\) 次多项式，考虑每次乘了一个系数并且求了一遍和（感觉有点抽象）。
直接插。
code
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define fi first
#define se second
const int N = 1005;
int n,k,mod;
LL inv[N],pre[N],nxt[N],f[N][N];
LL cal(int n,int k)
{
k++; k+=::n; LL res=0;
pre[0]=nxt[k+1]=1;
for(int i=1;i<=k;i++) pre[i]=pre[i-1]*(n-i)%mod;
for(int i=k;i>=1;i--) nxt[i]=nxt[i+1]*(n-i)%mod;
f[0][0]=1;
for(int i=1;i<=k;i++){
f[i][0]=1;
for(int j=1;j<=::n;j++)
f[i][j]=(f[i-1][j-1]*i+f[i-1][j])%mod;
}
if(n<=k) return f[n][::n];
for(int i=1;i<=k;i++)
{
LL t1=pre[i-1]*nxt[i+1]%mod,t2=inv[i-1]*((k-i)&1?mod-inv[k-i]:inv[k-i])%mod;
res=(res+t1*t2%mod*f[i][::n])%mod;
}
return res;
}
int main()
{
// freopen("in.in","r",stdin);
// freopen("out.out","w",stdout);
scanf("%d%d%d",&k,&n,&mod); inv[0]=inv[1]=1;
for(int i=2;i<=n*2;i++) inv[i]=(mod-mod/i)*inv[mod%i]%mod;
for(int i=2;i<=n*2;i++) inv[i]=inv[i]*inv[i-1]%mod;
LL ans=cal(k,n);
for(int i=1;i<=n;i++) ans=ans*i%mod;
printf("%lld\n",ans);
return 0;
}
浙公网安备 33010602011771号