传送门：序列求和 $1$
求：$\sum_{i = 1}^n i ^ k\quad(1 <= N <= 10^{18}, 1 <= K <= 2000)$
我才不会告诉你们那天晚上我去学伯努利数的时候弃坑了
感觉伯努利数的写法构造性太强了 $QwQ$而且似乎也只能写自然数幂和（其实听说它的深刻含义是指数型生成函数？不过我不太熟那东西），因此我后来就去学了拉格朗日插值法（毕竟普适性强一点）
其实一句话概括拉格朗日插值法的话它就是一个无脑构造（类似中国剩余定理的构造），其实很简单了，不过在这之前我要先抛出来一个定理：
一个 $n$次多项式可以通过 $n+1$个不同的点唯一确定这个多项式。
咋证呢？
其实反证一下就行了，假设有两个不同的 $n$次多项式在这 $n+1$个点的取值都相同，那么这两个多项式相减一下，这 $n+1$个点就是新多项式的所有根。而新多项式是一个 $n$次多项式，由代数基本定理可以知道 $n$次多项式是不可能有 $n+1$个根的，因此原命题得证。
拉格朗日插值法的思路是通过 $n+1$个点构造一个 $n$次多项式 $A(x)$，也就是说当你确定一个多项式的次数 $n$以及它在 $n+1$个点的取值后，你就可以计算出这个多项式啦~
设 $n+1$个点分别为 $(x_1,y_1),(x_2,y_2)\cdots(x_n,y_n),(x_{n+1},y_{n+1})$
则我们可以暴力构造一个多项式，它长这个样子
$$A(x)=\sum_{i=1}^{n+1} y_i \ell_{i}(x)$$
其中
$$\ell_{i}(x_j)=[i==j]$$
显然把这 $n+1$个点带进去后，点所对应的值都是等于多项式值的。
那么问题是我们怎样构造这个 $\ell_i(x)$？
其实超暴力的啦~
我们令
$$\ell_i(x)=\prod_{j\ne i} \frac{x-x_j}{x_i-x_j}$$
是不是就符合上面那个构造了呢？（读者不妨自己带几个值进去试一试）
嗯那拉格朗日插值法的总式就是
$$A(x)=\sum_{i=1}^{n+1} y_i\prod_{j\ne i} \frac{x-x_j}{x_i-x_j}$$
话说回来我们好像忘记了一件重要的事 $QAQ$
我们要求的是自然数幂和呀~
可是 $k$次的自然数幂和的多项式次数是多少呢？
$emmmmmmmmm$
由常识可知，$k=0$时结果是 $n$，$k=1$时结果是 $\frac{n(n-1)} 2$，那我们不妨猜想 $k$次的自然数幂和的多项式次数是 $k+1$吧。
证明：懒得写了怎么办，口糊吧。其实用数学归纳法随便搞搞 (二项式展开一下)，很容易弄出来的，具体的话可以参考这道数学题的第一问 (抱 $sammy$酱大腿 $QwQ$)
嗯那自然数幂求和就解决啦~
等等别跑，你还没说代码实现怎么搞呢！
唔~
首先我们明确一下我们实际上要求的就是这个式子的值：
$$\sum_{i=1}^{k+1} y_i\prod_{j\ne i} \frac{n-x_j}{x_i-x_j}$$
你需要先手动处理 $k+2$个点值（我处理的是 $0~k+1$），时间复杂度 $\Theta(k)$
然后你可以发现，其实每次的分母就是 $i!(n-i)!$，分子就是 $\frac{\prod_{j=0}^{k+1}(n-x_j)}{n-x_i}$
因此我们可以在 $\Theta(klogq)$的时间求出这个式子的值（其中 $q$是模数，$logq$是取逆元）
对于这道题来说是 $\Theta(TKlogq)$。然后成功被卡掉。
好吧其实有更快一点的算法可以去掉 $logq$，我们预处理一个 $1$乘到 $k+1$的数组，再预处理一个从 $k+1$乘到 $1$的数组，就可以 $\Theta(1)$求分子啦~
因此现在的复杂度是 $\Theta(TK)$
顺便说一下这是一道双倍经验，序列求和 $4$也可以用一样的代码 A 掉。
代码：
#include<bits/stdc++.h>
#define fo(i, a, b) for (int i = (a); i <= (b); ++i)
#define fd(i, a, b) for (int i = (a); i >= (b); --i)
#define mod 1000000007
#define ll long long
#define K 50005
int T, sgn;
ll fac[K + 3], y[K + 3], inv[K + 3], ans[K], p[K], q[K];
struct ANS{
ll n;
int k, id;
friend inline bool operator < (const ANS &x, const ANS &y)
{
return x.k < y.k;
}
}t[K];
inline ll pow (ll x, int p = mod - 2)
{
ll sum = 1;
while (p)
{
if (p & 1) sum = sum * x % mod;
x = x * x % mod, p >>= 1;
}
return sum;
}
inline ll solve (ll n, int k)
{
ll ans = 0, num = n;
++k;
// fo (i, 1, k)
// num = num * (n - i) % mod;
sgn = k & 1 ? -1 : 1;
p[0] = n;
fo (i, 1, k) p[i] = p[i - 1] * (n - i) % mod;
q[k + 1] = 1;
fd (i, k, 1) q[i] = q[i + 1] * (n - i) % mod;
fo (i, 1, k)
{
sgn = -sgn;
ans = (ans + y[i] * fac[i] % mod * fac[k - i] % mod * p[i - 1] % mod * q[i + 1] * sgn) % mod;
}
return (ans + mod) % mod;
}
int main ()
{
scanf("%d", &T);
fac[0] = 1;
inv[1] = 1;
fo (i, 2, K) inv[i] = (mod - mod / i) * inv[mod % i] % mod;
fo (i, 1, K) fac[i] = fac[i - 1] * inv[i] % mod;
y[0] = 0;
fo (i, 1, T)
{
scanf("%lld %d", &t[i].n, &t[i].k);
t[i].id = i;
}
std::sort(t + 1, t + T + 1);
fo (i, 1, T)
{
t[i].n = t[i].n % mod;
if (t[i - 1].k < t[i].k)
fo (j, 1, t[i].k + 1) y[j] = (y[j - 1] + pow(1LL * j, t[i].k)) % mod;
if (t[i].k + 1 >= t[i].n)
ans[t[i].id] = y[t[i].n];
else
ans[t[i].id] = solve(t[i].n, t[i].k);
}
fo (i, 1, T)
printf("%lld\n", ans[i]);
return 0;
}
0 条评论