公式：$f(x)=sum_{i=1}^{n} y_{i} prod_{i eq j} frac{x-x_{j}}{x_{i}-x_{j}}$.
这个式子正常算的话是 $O(n^2)$ 的，如果遇到 $x$ 是连续的情况可以优化到 $O(n log n)$.
但是有些时候我们只知道 $f(x)$ 在 $x=k$ 时的点值是不够的，有时必须求出这个多项式每一位系数.
多项式快速插值可以做到 $O(n log^2 n)$，但是快速插值非常非常难写，用处并不多.
相比之下，有一种简易的写法可以在 $O(n^2)$ 的时间复杂度内通过 $n$ 个不同的点来还原一个 $n-1$ 次多项式.
插值公式中 $prod_{j eq i} (x-x_{j})$ 是比较难求的，其他地方由于都是基于整数的运算，所以比较简单.
先令 $f_{i,j}$ 表示考虑前 $i$ 个点 $(x,y)$，$x^j$ 前的系数.
那么有转移：$f_{i,j}=f_{i-1,j-1}+f_{i-1,j} imes (-x_{i})$ 即分别表示当前位的贡献为 $x^1 / -x_{i}$.
求出这个后，我们枚举 $i$，然后想 $O(n)$ 计算 $h(x)=prod_{i eq j} (x-x_{j})$.
令 $k1[i],k2[i]$ 分别表示 $h(x)$ 的 $x^i$ 前的系数，强制让第 $i$ 位贡献 $x^1$ 时 $x^i$ 前的系数.
由于有 $k2$ 这个强制贡献的状态，转移就比较简单：
$k1[i] leftarrow k2[i+1]$
$f_{n,i}=k1[i] imes (-x_{i}) +k2[i] Rightarrow k2[i]=f_{n,i}+k1[i] imes (x_{i})$.
算出 $k2$ 后把 $y_{i}$ 及插值公式中分母的贡献乘上然后累加到答案数组中即可.
应用：
求 $sum_{i=1}^{n} i^k$.
这是一个关于 $n$ 的 $k+1$ 次多项式.
所以可以取 $k+2$ 个点带进去，然后用拉格朗日插值法来求值.
具体，$f(k)=sum_{i=1}^{n} y_{i} prod_{}^{i eq j}frac{k-x_{j}}{x_{i}-x_{j}}$
由于点可以做到取 $x$ 连续的，所以提前预处理前缀/后缀积极可以做到 $O(n log n)$.
code:
#include <cstdio>
#include <vector>
#include <cstring>
#include <algorithm>
#define N 1000009
#define ll long long
#define mod 1000000007
#define setIO(s) freopen(s".in","r",stdin)
using namespace std;
int f[N];
int ifac[N],fac[N],pre[N],suf[N],inv[N],n,K;
int qpow(int x,int y) {
int tmp=1;
for(;y;y>>=1,x=(ll)x*x%mod)
if(y&1) tmp=(ll)tmp*x%mod;
return tmp;
}
int INV(int x) { return qpow(x,mod-2); }
void init() {
ifac[0]=fac[0]=inv[1]=1;
for(int i=2;i<N;++i) inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod;
inv[0]=1;
for(int i=1;i<N;++i) {
fac[i]=(ll)fac[i-1]*i%mod;
ifac[i]=(ll)ifac[i-1]*inv[i]%mod;
}
pre[0]=1,suf[n+1]=1;
for(int i=1;i<=n;++i) pre[i]=(ll)pre[i-1]*(K-i+mod)%mod;
for(int i=n;i>=1;--i) suf[i]=(ll)suf[i+1]*(K-i+mod)%mod;
}
int sol() {
int ans=0;
for(int i=1;i<=n;++i) {
int a1=(ll)ifac[i-1]*ifac[n-i]%mod;
if((n-i)&1) a1=(ll)a1*(mod-1)%mod;
int a2=(ll)pre[i-1]*suf[i+1]%mod;
(ans+=(ll)f[i]*a1%mod*a2%mod)%=mod;
}
return ans;
}
int main() {
// setIO("input");
scanf("%d%d",&K,&n),n+=2;
init();
for(int i=1;i<=n;++i) {
f[i]=(ll)(f[i-1]+qpow(i,n-2))%mod;
}
printf("%d
",sol());
return 0;
}
还原多项式系数
#include <cstdio>
#include <cstring>
#include <algorithm>
#define N 2008
#define ll long long
#define mod 998244353
#define setIO(s) freopen(s".in","r",stdin)
using namespace std;
int f[N][N],k1[N],k2[N],s[N],n;
struct point {
int x,y;
point(int x=0,int y=0):x(x),y(y){}
}a[N];
int ADD(int x,int y) {
return (ll)(x+y)%mod;
}
int DEC(int x,int y) {
return (ll)(x-y+mod)%mod;
}
int MUL(int x,int y) {
return (ll)x*y%mod;
}
int qpow(int x,int y) {
int tmp=1;
for(;y;y>>=1,x=MUL(x,x))
if(y&1) {
tmp=MUL(tmp,x);
}
return tmp;
}
int get_inv(int x) { return qpow(x,mod-2); }
void init() {
f[0][0]=1;
for(int i=1;i<=n;++i) {
for(int j=1;j<=i;++j) {
f[i][j]=ADD(f[i-1][j-1],MUL(mod-a[i].x,f[i-1][j]));
}
f[i][0]=MUL(f[i-1][0],mod-a[i].x);
}
}
int main() {
// setIO("input");
int X;
scanf("%d%d",&n,&X);
for(int i=1;i<=n;++i) {
scanf("%d%d",&a[i].x,&a[i].y);
}
init();
for(int i=1;i<=n;++i) {
for(int j=0;j<=n;++j) k2[j]=f[n][j];
for(int j=n-1;j>=0;--j) {
k1[j]=k2[j+1];
k2[j]=ADD(k2[j],MUL(k1[j],a[i].x));
}
int inv=1;
for(int j=1;j<=n;++j)
if(i!=j) {
inv=(ll)inv*(a[i].x-a[j].x+mod)%mod;
}
inv=get_inv(inv);
for(int j=0;j<=n-1;++j) {
(s[j]+=(ll)inv*a[i].y%mod*k1[j]%mod)%=mod;
}
}
int ans=0;
for(int i=n-1;i>=0;--i) {
ans=(ll)((ll)ans*X%mod+s[i])%mod;
}
printf("%d
",ans);
return 0;
}
例题
CF917D Stranger Trees
给你一颗树，求 $n$ 个点有多少个生成树满足该生成树与给定树有 $k$ 条边是重合的.
题解：
先对完全图构建矩阵，然后将原树上的边 $(x,y)$ 在矩阵中的边权标记成 $x^1$，其余边权为 $1$.
矩阵树定理求的是所有生成树边权乘积之和，那么要是可以对含 $x$ 的矩阵求行列式的话可以直接得出答案.
但是复杂度太高，而且难写（写不了）
所以用 $n$ 个不同的整数来替换那个 $x^1$，然后跑出来 $n$ 个结果，用拉格朗日插值还原出多项式的系数即可.
#include <cstdio>
#include <vector>
#include <cstring>
#include <algorithm>
#define N 103
#define ll long long
#define mod 1000000007
#define setIO(s) freopen(s".in","r",stdin)
using namespace std;
int n;
int A[N],B[N];
int f[N][N],k1[N],k2[N],ans[N];
int deg[N][N],con[N][N],a[N][N];
struct point {
int x,y;
point(int x=0,int y=0):x(x),y(y){}
}p[N];
int qpow(int x,int y) {
int tmp=1;
for(;y;y>>=1,x=(ll)x*x%mod)
if(y&1) {
tmp=(ll)tmp*x%mod;
}
return tmp;
}
int get_inv(int x) {
return qpow(x,mod-2);
}
int ADD(int x,int y) {
return (ll)(x+y)%mod;
}
int DEC(int x,int y) {
return (ll)(x-y+mod)%mod;
}
int MUL(int x,int y) {
return (ll)x*y%mod;
}
int gauss() {
int ans=1;
for(int i=1;i<n;++i) {
for(int j=i+1;j<n;++j) {
while(a[j][i]) {
int t=a[i][i]/a[j][i];
for(int k=i;k<n;++k) {
a[i][k]=DEC(a[i][k],MUL(t,a[j][k]));
}
swap(a[j],a[i]);
ans=(ll)ans*(mod-1)%mod;
}
}
if(!a[i][i]) {
return 0;
}
}
for(int i=1;i<n;++i) {
ans=(ll)ans*a[i][i]%mod;
}
return ans;
}
int cal(int val) {
for(int i=1;i<=n;++i) {
for(int j=1;j<=n;++j) {
a[i][j]=mod-1;
}
}
for(int i=1;i<=n;++i) {
a[i][i]=n-1;
}
for(int i=1;i<n;++i) {
int x=A[i],y=B[i];
a[x][x]=(ll)(DEC(a[x][x],1)+val)%mod;
a[y][y]=(ll)(DEC(a[y][y],1)+val)%mod;
a[x][y]=(ll)(a[x][y]+1-val+mod)%mod;
a[y][x]=(ll)(a[y][x]+1-val+mod)%mod;
}
return gauss();
}
void init() {
f[0][0]=1;
for(int i=1;i<=n;++i) {
for(int j=1;j<=i;++j)
f[i][j]=ADD(f[i-1][j-1],MUL(f[i-1][j],mod-p[i].x));
f[i][0]=(ll)f[i-1][0]*(mod-p[i].x)%mod;
}
}
int main() {
// setIO("input");
scanf("%d",&n);
int x,y,z;
for(int i=1;i<n;++i) {
scanf("%d%d",&A[i],&B[i]);
}
for(int i=1;i<=n;++i) {
p[i].x=i;
p[i].y=cal(i);
}
init();
for(int i=1;i<=n;++i) {
int inv=1;
for(int j=1;j<=n;++j) {
if(i!=j) inv=(ll)inv*(p[i].x-p[j].x+mod)%mod;
}
inv=get_inv(inv);
for(int j=0;j<=n;++j) {
k2[j]=f[n][j];
}
for(int j=n-1;j>=0;--j) {
k1[j]=k2[j+1];
k2[j]=ADD(k2[j],MUL(p[i].x,k1[j]));
}
for(int j=0;j<=n-1;++j) {
ans[j]=ADD(ans[j],(ll)k1[j]*inv%mod*p[i].y%mod);
}
}
for(int i=0;i<n;++i) {
printf("%d ",ans[i]);
}
return 0;
}
LuoguP4463 [集训队互测2012] calc
朴素的 DP 非常好列：$f[i][j]$ 表示选了 $i$ 个数，且值域为 $[1,j]$ 的总价值和.
那么有 $f[i][j]=f[i-1][j-1] imes j+f[i][j-1]$，直接算的话复杂度是 $O(nD)$ 的.
但是我们可以猜测这是一个关于 $j$ 的 $g_{i}$ 次多项式.
有一个结论：对于 $n$ 次多项式 $h(x)$，满足 $h(x)-h(x-1)$ 是 $n-1$ 次多项式.
那么有 $f[i][j]-f[i][j-1]=f[i-1][j-1] imes j$.
将 $g$ 带入，有 $g_{i}-1=g_{i-1}+1$.
即 $g_{i}=g_{i-1}+2$，说明这是一个关于 $j$ 的 $2 imes i$ 次多项式.
那么我们就求出 $f[n][1...2n+1]$ 后将值带入，然后拉格朗日插值来插一下就行了.
code:
#include <cstdio>
#include <cstring>
#include <algorithm>
#define N 2002
#define ll long long
#define setIO(s) freopen(s".in","r",stdin)
using namespace std;
int D,n,mod,tot,f[N][N],fac[N];
void init() {
fac[0]=1;
for(int i=1;i<N;++i) {
fac[i]=(ll)fac[i-1]*i%mod;
}
}
struct point {
int x,y;
point(int x=0,int y=0):x(x),y(y){}
}a[N];
int qpow(int x,int y) {
int tmp=1;
for(;y;y>>=1,x=(ll)x*x%mod)
if(y&1) tmp=(ll)tmp*x%mod;
return tmp;
}
int get_inv(int x) {
return qpow(x,mod-2);
}
int calc() {
int ans=0;
for(int i=1;i<=tot;++i) {
int inv=1,up=1;
for(int j=1;j<=tot;++j) {
if(i==j) continue;
up=(ll)up*(D-a[j].x+mod)%mod;
inv=(ll)inv*(a[i].x-a[j].x+mod)%mod;
}
inv=get_inv(inv);
(ans+=(ll)a[i].y*up%mod*inv%mod)%=mod;
}
return ans;
}
int main() {
// setIO("input");
scanf("%d%d%d",&D,&n,&mod);
init();
for(int i=0;i<=2*n+1;++i) f[0][i]=1;
for(int i=1;i<=n;++i) {
for(int j=1;j<=2*n+1;++j) {
f[i][j]=(ll)(f[i][j-1]+(ll)f[i-1][j-1]*j%mod)%mod;
}
}
for(int i=1;i<=2*n+1;++i) {
a[++tot]=point(i,f[n][i]);
}
printf("%d
",(ll)calc()*fac[n]%mod);
return 0;
}