拉格朗日插值法 学习笔记
拉格朗日插值法学习笔记
拉格朗日插值
对于一个$n$次函数$f(x)$，根据小学数学可知，平面上$n+1$就可以确定$f(x)$。那么，假如我们已知这$n+1$个点，我们该如何求到$f(x)$呢？
拉格朗日机智地想到，假如$n$个点为$(x_i,y_i),1\le i\le n$，那么，可以得到:
$$f(x)=\sum_{i=1}^{n} y_i \prod_{j\not= i} \frac{x-x_j}{x_i-x_j}$$
很显然，这个式子是正确的。于是，我们就可以做出【模板】拉格朗日插值了。
至于拉格朗日插值2直接卷积就好了。
拉格朗日插值法优化dp
对于一个dp式子，如果我们可以证明它是一个$g(n)$次的多项式，那我们就可以用拉格朗日插值法把它插出来。其中$g(n)$是关于$n$的多项式。
口胡可能不能体会，还是看例题吧。
$\text {The 1st}$
题目大意
给定$n,k$，对于$a_{1,2,…,n}$满足:
$$(\forall i,j,i\not=j \rightarrow a_i \not= a_j )\wedge (\forall i,a_i\le k)$$
的$\prod_{i=1}^{n} a_i$。
思路
我们考虑dp。为了方便，我们先钦定$a_{1,2,…,n}$是个严格单调上升序列，最后再乘上$n!$即可。我们就可以设$dp[i][j]$表示，前面已经选了$i$个数，最大值不超过$j$的贡献。可以得到转移方程:
$$dp[i][j]=dp[i-1][j-1]\times j + dp[i][j-1]$$
最后的答案即为$dp[n][k]\times n!$。
显然，这是一个$\Theta(nk)$的方法。
这个时候，我们想起一个东西。如果$\delta f(x)$是一个$n$次函数，那么，$f(x)$为$n+1$次函数。
其中，$\delta f(x)=f(x+1)-f(x)$（其实就是有限微积分里面那个）。
于是，我们观察到:
$$dp[i][j]-dp[i][j-1]=dp[i-1][j-1]\times j$$
下面我用$dp_i(j)$表示$dp[i]$所对应的多项式。（这叫做先猜后证么？雾）
我们设$g(i)$为$dp_i(j)$的次数，那么根据上面的式子我们就可以得到:
$$g(i)-1=g(i-1)+1$$
$$\Rightarrow g(i)=g(i-1)+2$$
又因为$g(0)=0$，所以，$g(n)=2n$。
于是，我们只需要确定$dp_n(i)$中的$2n+1$个点即可。
于是，我们就可以$\Theta(n^2)$求出$dp[n][1,2,…,2n+1]$然后$\Theta(n^2)$求到$dp_n(i)$，在把$k$代入函数中求值即可。
总时间复杂度$\Theta(n^2\log n)$。（费马小定理求逆元）
$\text {Code}$
#include <bits/stdc++.h>
using namespace std;
#define Int register int
#define ll long long
#define MAXN 1005
template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c == '-') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
int k,n,p;
int y[MAXN],dp[MAXN][MAXN];
int quick_pow (int a,int b,int c){
int res = 1;
while (b){
if (b & 1) res = 1ll * res * a % c;
a = 1ll * a * a % c;
b >>= 1;
}
return res;
}
int inv (int x,int mod){
return quick_pow (x,mod - 2,mod);
}
signed main(){
read (k,n,p);int up = n << 1 | 1;
for (Int i = 0;i <= up;++ i) dp[0][i] = 1;
for (Int i = 1;i <= n;++ i)
for (Int j = 1;j <= up;++ j)
dp[i][j] = (1ll * dp[i - 1][j - 1] * j % p + dp[i][j - 1]) % p;
for (Int i = 0;i <= up;++ i) y[i] = dp[n][i];
int ans = 0;
for (Int i = 0;i <= up;++ i){
int now = 1;
for (Int j = 0;j <= up;++ j)
if (i ^ j) now = 1ll * now * inv (i - j,p) % p * (k - j) % p;
now = (now + p) % p;
ans = (ans + 1ll * now * y[i] % p) % p;
}
for (Int i = 1;i <= n;++ i) ans = 1ll * ans * i % p;
write (ans),putchar ('\n');
return 0;
}
$\text {The 2ed}$
题目大意
给定一棵$n$个点的树，给定$d$。现在给树上每个点赋值$[1,d]$中的某个数，并满足子节点的权值不大于父节点的权值。问合法方案数。
思路
有了上一道题的经验这道题应该好做多了。
我们可以设$dp[i][j]$表示以$i$为根的子树内最大值不大于$j$的合法方案数。
可以得到转移式:
$$dp[u][i]=dp[u][i-1]+\prod_{v\in u\to to} dp[v][i]$$
答案就是$dp[1][d]$。
我们发现这是个前缀和形式。于是可以得到:
$$\delta dp_u(i)=\prod_{v} dp[v][i]$$
我们设$g(u)$为$dp_u(i)$的次数，可以得到:
$$g(u)=\sum g(v)+1$$
我们发现$g(u)$其实就是以$u$为根的子树大小。于是，我们只需要$n+1$个点就可以确定$dp_1(i)$，进而求到$dp[1][d]$。
所以我们就可以用拉格朗日插值法在$\Theta(n^2)$内求到答案。
$\text {Code}$
#include <bits/stdc++.h>
using namespace std;
#define Int register int
#define mod 1000000007
#define ll long long
#define MAXN 3005
template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c == '-') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
int quick_pow (int a,int b,int c){
int res = 1;
while (b){
if (b & 1) res = 1ll * res * a % c;
a = 1ll * a * a % c;
b >>= 1;
}
return res;
}
struct edge{
int v,nxt;
}e[MAXN];
int top = 1;
int head[MAXN];
void Add_Edge (int u,int v){
e[++ top] = edge {v,head[u]};
head[u] = top;
}
int n,d;
int g[MAXN][MAXN],sum[MAXN][MAXN];
void dfs (int u){
for (Int i = 1;i <= n + 1;++ i) sum[u][i] = 1;
for (Int i = head[u];i;i = e[i].nxt){
int v = e[i].v;
dfs (v);
for (Int j = 1;j <= n + 1;++ j) sum[u][j] = 1ll * sum[u][j] * g[v][j] % mod;
}
for (Int i = 1;i <= n + 1;++ i) g[u][i] = (g[u][i - 1] + sum[u][i]) % mod;
}
signed main(){
read (n,d);
for (Int i = 2,fa;i <= n;++ i) read (fa),Add_Edge (fa,i);
dfs (1);
int ans = 0;
for (Int i = 1;i <= n + 1;++ i){
int up = 1,down = 1;
for (Int j = 1;j <= n + 1;++ j)
if (i ^ j) up = 1ll * up * (d < j ? d + mod - j : d - j) % mod,down = 1ll * down * (i < j ? i + mod - j : i - j) % mod;
ans = (ans + 1ll * up * quick_pow (down,mod - 2,mod) % mod * g[1][i] % mod) % mod;
}
write (ans),putchar ('\n');
return 0;
}