写在前面
昨天 ABC 的 F 题,结论推出来了,猜到是拉格朗日,奈何我只会板子,不会分析次数;
赛后知道正解后感觉大受震撼。我还是太 naive 了 /kk
前置知识
-
组合数学
-
拉格朗日插值法。
Description
给你 (n, m, k),定义 (f(n, m)) 为:
求 (f(n,m) mod 10^9+7)
数据范围:
(0 le n le 10^{18})
(0 le m le 30)
(1 le k le 2.5 imes 10^{6})
Solution
先对第一个样例搞一下。
我们用下面这个表格来算一下 (f(3,4)) 最终递归调用 (f(i,0)) 多少次
结合 (f(n, m) = f(n-1,m) + f(n, m-1) (n > 0, m > 0)) 的递归式
发现这个玩意很像平面上 ((0,0)) 到 ((n,m)) 的方案数。
而我们求的 (f(i,0)) 就对应着 ((i, 1)) 到 ((n, m)) 的方案数。
根据组合数学,可以算出 (f(i,0)) 调用了 (C_{n+m-i-1}^{m-1}) 次
又因为 (f(i,0) = i^k)
然后我们就能够用一个式子表示出答案:
根据 @GuidingStar 的说法
前缀和是 (1) 次, (inom{n+m-i-1}{m-1}) 是 (m-1) 次,(i^k) 是 (k) 次, 一共 (k+m) 次。
吾以为 之所以 (inom{n+m-i-1}{m-1}) 是 (m-1) 次,是因为这个式子约分之后分子与 (i) 相关的一共有 (m-1) 项,至于下面的 ((m-1)!) 与 (i) 无关。
我们设 (F(n) = sum_{i = 0}^{n} inom{n+m-i-1}{m-1} i^k)
这是一个 (k+m) 次多项式,我们只要找到 (k+m+1) 个点就能确定他的式子。
这里我们取 (x = 0 sim k + m) 这一段连续的区间,用暴力刷表的方法求出 ((x, F(x)))
然后利用拉格朗日插值法求解即可。
用到 (x) 取值连续时的优化,如果不会可以参考这篇博客。
总复杂度为 (O((k+m) log k)),可以通过。
Code
/*
Work by: Suzt_ilymics
Problem:
Knowledge: 拉格朗日插值法
Time: O((k+m)log mod)
*/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<queue>
#define int long long
#define orz cout<<"lkp AK IOI!"<<endl
using namespace std;
const int MAXN = 5e6+5;
const int INF = 1e9+7;
const int mod = 1e9+7;
int n, m, k, Ans = 0;
int S[MAXN], fac[MAXN], inv[MAXN], pre[MAXN], suf[MAXN];
int ans[MAXN];
int read(){
int s = 0, f = 0;
char ch = getchar();
while(!isdigit(ch)) f |= (ch == '-'), ch = getchar();
while(isdigit(ch)) s = (s << 1) + (s << 3) + ch - '0' , ch = getchar();
return f ? -s : s;
}
int Pow(int x, int p, int mod) {
int res = 1;
while(p) {
if(p & 1) res = res * x % mod;
x = x * x % mod;
p >>= 1;
}
return res;
}
signed main()
{
n = read(), m =read(), k = read();
if(!m) {
printf("%lld
", Pow(n % mod, k, mod));
return 0;
}
if(!n) {
printf("0
");
return 0;
}
ans[0] = 0;
for(int i = 1; i <= k + m; ++i) S[i] = Pow(i, k, mod);
for(int i = 1; i <= m; ++i) {
for(int j = 1; j <= k + m; ++j) {
ans[j] = (ans[j - 1] + S[j]) % mod;
}
for(int j = 1; j <= k + m; ++j) {
S[j] = ans[j];
}
}
if(n <= k + m) {
printf("%lld
", ans[n]);
return 0;
}
// for(int i = 1; i <= k + 2; ++i) inv[i] = Pow(fac[i], mod - 2, mod);
pre[0] = n % mod, suf[k + m + 1] = 1, fac[0] = 1;
for(int i = 1; i <= k + m; ++i) pre[i] = pre[i - 1] * ((n - i) % mod) % mod;
for(int i = k + m; i >= 0; --i) suf[i] = suf[i + 1] * ((n - i) % mod) % mod;
for(int i = 1; i <= k + m; ++i) fac[i] = fac[i - 1] * i % mod;
for(int i = 0; i <= k + m; ++i) {
int a = (i == 0 ? 1 : pre[i - 1]) * suf[i + 1] % mod;
int b = fac[i] * ((k + m - i) & 1 ? -1 : 1) * fac[k + m - i] % mod;
Ans = Ans + S[i] * a % mod * Pow(b, mod - 2, mod) % mod;
Ans %= mod;
}
printf("%lld
", (Ans + mod) % mod);
return 0;
}