bzoj 3453: tyvj 1858 XLkxc

题目描述：

求

$$\sum_{i=0}^{n} \sum_{j=1}^{a+id} \sum_{l=1}^{j}l^k$$

输入格式

多组数据。第一行一个整数$T$,表示数据组数；

接下来$T$行， 每行$4$个整数$k,a,n,d$,含义如公式所示。

输出格式

输出一行，一个整数，表示答案， 输出结果请对$1234567891$取模。

样例

input

4
1 51 49 1
1 52 56 1
1 53 53 1
1 51 56 1

output

4128450
5672754
5222196
5480360

限制与约定

存在$30\%$的数据，有$1 \leqslant k \leqslant 150$；
存在$40\%$的数据，有$a+nd \leqslant 10^6$；
对于$100\%$的数据，有$1 \leqslant T \leqslant 4$, $1 \leqslant k \leqslant 3000$, $1  \leqslant a,n,d < 1234567891$。

搬运工提示

您可能需要的数学知识点：

• 伯努利数与自然数幂和
• 伯努利数-百度百科

思路：
这题我一开始想用国王奇遇记的方法递推然而推了半天头昏脑胀，啥也没推出来。。。感动的落泪
在shallwe大爷的良心提示链接里提示我们用伯努利数化成多项式做。。。感觉和TA爷作业题的我那个做法有点像，下面来xjb推一发：
我们考虑将伯努利数写成多项式的形式，记$S(n,i) = \sum_{j = 1}^{n}j^{i}$ , 设$A(i,j)(i <= k + 1,j <= i + 1)$ 为$S(n,i)$ 简化为多项式后第$j$ 项的系数，根据相关知识它是一个$i + 1$ 次的多项式,而且根据上述链接我们可以在$O(k^2)$ 的时间预处理，而题目等价于求：

$$\sum_{i = 0}^{n}\sum_{j = 1}^{a + i * d}\sum_{l = 1}^{j}l^{k}$$
将内部化为多项式的形式：
$$=\sum_{i = 0}^{n}\sum_{j = 1}^{a + i * d}\sum_{l = 1}^{k + 1}A(k,l) * (j + 1)^l$$
变化一下求和顺序：
$$=\sum_{i = 0}^{n}\sum_{l = 1}^{k + 1}A(k,l) * \sum_{j = 1}^{a + 1 + i * d}j ^ {l} - (n + 1) *\sum_{l = 1}^{k + 1}A(k,l)$$
后面那个可以单次 $O(k)$ 求出，我们先不管它
再一次将内部的伯努利数化为多项式：
$$=\sum_{l = 1}^{k + 1} A(k,l) * \sum_{t = 1}^{l + 1}A(l,t)\sum_{i = 0}^{n}(a + 2 + i * d) ^ t$$
将其展开后又一次将伯努利数化为多项式最后得到：
$$\sum_{l = 1}^{k + 1} A(k,l) * \sum_{t = 1}^{l + 1}A(l,t) * F(t)$$
其中 $$F(t) = \sum_{j = 0}^{t} * \binom{t}{j} * d^j * (a + 2)^{t - j} * S(j,n)$$
这样我们完成了化简，预处理的就是组合数，多项式系数，查询的时候先算S，然后是F，最后计算关于A的贡献，以上每一步都是 $O(k^2)$ 的。
然而考场上头昏脑胀把边界写残，睡一觉改了一个地方就AC了，伤心到变形啊。。。

代码：

#include<iostream>
#include<cstring>
#include<string>
#include<cstdio>
#define N 140
using namespace std;
typedef long long LL;
const LL P = 1234567891LL;
LL k,a,n,d,K,C[N][N],kk[6],aa[6],nnn[6],ddd[6],B[N],A[N][N],inv[N],F[N],f[N],ans;
inline void inc(LL &x,LL y)
{
    x += y;
    if (x >= P) x -= P;
}

inline LL quick(LL x,LL y)
{
    LL res = 1,base = x;
    for (;y;y >>= 1)
    {
        if (y & 1) res = res * base % P;
        base = base * base % P;
    }
    return res;
}

inline void init()
{
    for (int i = 0;i <= K + 4; ++i)
        C[i][0] = 1;
    for (int i = 1;i <= K + 4; ++i)
        for (int j = 1;j <= i; ++j)
        {
            C[i][j] = C[i - 1][j] + C[i - 1][j - 1];
            if (C[i][j] >= P) C[i][j] -= P;
        } 
    inv[0] = inv[1] = 1;
    for (int i = 2;i <= K + 4; ++i)
        inv[i] = inv[P % i] * 1LL * (P - P / i) % P;
    B[0] = 1;
    for (int i = 1;i <= K + 3; ++i)
    {
        B[i] = 0;
        for (int j = 0;j < i; ++j)
            inc(B[i],B[j] * C[i + 1][j] % P);
        B[i] = (P - inv[i + 1]) * B[i] % P;
    }
    for (int i = 1;i <= K + 3; ++i)
        for (int j = 1;j <= i + 1; ++j)
            A[i][j] = inv[i + 1] * C[i + 1][j] % P * B[i + 1 - j] % P;
}

inline void DO_IT()
{
     LL nn = 1;
    //cout<<A[0][1]<<endl;
    f[0] = n + 1;
    for (int i = 1;i <= k + 2; ++i)
    {
        f[i] = 0;
        nn = n + 1;
        for (int s = 1;s <= i + 1; ++s)
        {
            inc(f[i],A[i][s] * nn % P);
            nn = nn * (n + 1) % P;
        }
    }
    //cout<<f[1]<<endl;
    LL dd,inva = quick(a + 2,P - 2),aa = a + 2,nowa;
    for (int i = 1;i <= k + 2; ++i)
    {
        F[i] = 0;
        dd = 1;
        nowa = aa;
        for (int s = 0;s <= i; ++s)
        {
            inc(F[i],dd * nowa % P * C[i][s] % P * f[s] % P);
            dd = dd * d % P;
            nowa = nowa * inva % P;
        }
        aa = aa * (a + 2) % P;
    }
     ans = 0;
     for (int t = 1;t <= k + 1; ++t)
        for (int p = 1;p <= t + 1; ++p)
            inc(ans,A[k][t] * A[t][p] % P * F[p] % P);
    LL sum = 0;
    for (int l = 1;l <= k + 1; ++l) inc(sum,A[k][l]);
    sum = sum * (n + 1) % P;
    inc(ans,P - sum);
}

int main()
{
    int T;
    scanf("%d",&T);
    K = 0; 
    for (int o = 1;o <= T;++o)
    {
        scanf("%lld%lld%lld%lld",&kk[o],&aa[o],&nnn[o],&ddd[o]);
        K = max(K,kk[o]);
    }
    init();
    for (int o = 1;o <= T; ++o)
    {
        k = kk[o]; a = aa[o]; n = nnn[o]; d = ddd[o];
        DO_IT();
        printf("%lld\n",ans);
    }
    return 0;
}