【动态规划】【拉格朗日插值优化dp】集训队互测2012 calc
题目描述
一个序列 \(a_1,a_2,\dots,a_n\) 是合法的，当且仅当：
- \(a_1,a_2,\dots,a_n\) 都是 \([1,k]\) 中的整数。
- \(a_1,a_2,\dots,a_n\) 互不相等。
一个序列的值定义为它里面所有数的乘积，即 \(a_1\times a_2\times\dots\times a_n\)。
求所有不同合法序列的值的和对 \(p\) 取模后的结果。两个序列不同当且仅当他们任意一位不同。
对于 \(100\%\) 的数据，\(k \le 10 ^ 9\)，\(n \le 500\)，\(p \le 10 ^ 9\)，保证 \(p\) 为素数，保证 \(n + 1 < k < p\)。
前置芝士：拉格朗日插值
算法描述
考虑dp，我们发现如果乱序不是很好dp，所以我们假设序列递增，最后乘上一个\(n!\)就好了。我们设\(dp_{i,j}\)表示前\(i\)个数用\([1,j]\)凑成的合法序列值之和，则有以下转移：
第\(i\)个数是\(j\)：\(dp_{i,j} += dp_{i - 1,j - 1} * j\)
第\(i\)个数不是\(j\)：\(dp_{i,j} += dp_{i,j - 1}\)
所以
\[dp_{i,j} = dp_{i - 1,j - 1} * j + dp_{i,j - 1} \]最后答案就是\(dp_{n,k} * n!\)，时间复杂度\(O(nk)\)
黑题切了
然而我们发现\(k \leq 10^9\)，并不能通过此题，那么如何加速这个式子呢？
考虑到\(k\)很大，我们试图将复杂度变得只和\(n\)有关，这个dp式子有一个性质：它只由乘法和加法组成，而且是一个没有\(min、max\)，没有多种情况的递推式。并且在\(n\)一定的情况下，\(dp_{n,j}\)只与\(j\)有关，我们不妨大胆假设\(dp_{n,j}\)是一个关于\(j\)的多项式，设它的次数是\(g(n)\)。
首先有一结论：
设多项式\(f(x)\)的次数是\(n\)，那么它的差分\(f(x) - f(x - 1)\)的次数是\(n - 1\)
所以我们可以知道，\(dp_{n,j} - dp_{n,j - 1}\)的次数是\(g(n) - 1\)
通过递推式，我们可以得出：
\[dp_{n,j} - dp_{n,j - 1} = j * dp_{n - 1,j - 1} \]因为\(dp\)是关于\(j\)的多项式，所以乘\(j\)会使次数加1，可得\(dp_{n,j} - dp_{n,j - 1}\)的次数是\(g(n - 1) + 1\)
联立得：
\[g(n) - 1 = g(n - 1) + 1 \]\[g(n) - g(n - 1) = 2 \]因为\(g(0) = 0\)（前0个数种类为1，是0次多项式，即常数），所以：
\[g(n) = 2n \]所以我们最终确定了\(dp_{n,j}\)是一个\(2n\)次多项式，要求它在\(k\)处的值，而\(k\)又很大，可以使用拉格朗日插值，\(dp\)算出\((1,dp_{n,1})\)到\((2n + 1,dp_{n,2n + 1})\)这\(2n + 1\)个连续点，然后代入插值就可以求出\(dp_{n,k}\)的值，然后乘上\(n!\)即可。这样，时间复杂度就转化成了\(O(n^2)\)，可以通过本题。
前文已经阐述过，这样优化的关键条件就是\(dp\)式子的特殊性质（其实所有\(dp\)优化都是这样），本题的特殊性质就是\(dp\)是一个只和乘法和加法有关的，没有多种情况、分类讨论的简单递推式，并且目标维中有一个很大，将其设为多项式的自变量，用拉插就可以完成计算。
Code
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1005;
ll n,k,p,dp[N][N],y[N],pre[N],suf[N],inv[N];
inline ll lag(ll len,ll X)
{
ll ret = 0; pre[0] = 1;suf[len + 1] = 1;
for(int i = 1;i <= len;i++) pre[i] = pre[i - 1] * ((X - i) % p + p) % p;
for(int i = len;i >= 1;i--) suf[i] = suf[i + 1] * ((X - i) % p + p) % p;
for(int i = 1;i <= len;i++)
{
ll res = y[i] * pre[i - 1] % p * suf[i + 1] % p * inv[i - 1] % p * inv[len - i] % p;
if((len - i) & 1) res = p - res;
ret = (ret + res) % p;
}
return ret;
}
int main()
{
cin>>k>>n>>p;
inv[0] = inv[1] = 1;
for(int i = 2;i <= N - 1;i++) inv[i] = (p - (p / i)) * inv[p % i] % p;
for(int i = 1;i <= N - 1;i++) inv[i] = inv[i] * inv[i - 1] % p;
memset(dp,0,sizeof(dp));
for(int i = 0;i <= 2 * n + 1;i++) dp[0][i] = 1;
for(int i = 1;i <= n;i++)
for(int j = 1;j <= 2 * n + 1;j++)
dp[i][j] = dp[i - 1][j - 1] * j % p + dp[i][j - 1],dp[i][j] %= p;
for(int i = 1;i <= 2 * n + 1;i++) y[i] = dp[n][i];
ll ans = lag(2 * n + 1,k);
for(int i = 1;i <= n;i++) ans = ans * i % p;
cout<<ans;
return 0;
}