传送门:序列求和$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;
}
分类: 文章

quhengyi11

喵(ฅ>ω<*ฅ)

发表评论

电子邮件地址不会被公开。 必填项已用*标注

你是机器人吗? =。= *