多项式及其求和

        本文思路来源：杜教—-多项式及求和。待更新。

引入

        本文主要讨论如何求解如下的求和式：

$$\sum_{i=1}^ni^d$$

解法1

        规定：

$$S_d(n)=\sum_{i=0}^{n-1}i^d$$

        构造向量：

$$V_n=\begin{bmatrix} S_d(n) & n^0 &\cdots & n^d \end{bmatrix}$$

        注意到：

$$S_d(n)=S_d(n-1)+n^d$$

        以及根据二项式定理：

$$(n+1)^i=\sum_{j=0}^iC_i^jn^j$$

        于是我们可以预处理组合数，然后构造矩阵$A$满足：

$$V_{n+1}=V_nA$$

        $A$是一个$(d+2)*(d+2)$的矩阵，容易得到：

$$V_{n+1}=V_1A^n$$

        利用矩阵快速幂可以在$O(d^3logn)$的时间复杂度下解决这个问题。

#include <bits/stdc++.h>
using namespace std;
const int mod = 1e9 + 7;
typedef long long ll;
struct Matrix
{
    int op[51][51], n;
    void operator*=(const Matrix& s)
    {
        int tmp[51][51];
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) tmp[i][j] = 0;
        }
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                for (int k = 1; k <= n; k++) tmp[i][j] = (tmp[i][j] + 1ll * op[i][k] * s.op[k][j] % mod) % mod;
            }
        }
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) op[i][j] = tmp[i][j];
        }
    }
} sta, ans;
int C[51][51];
int main()
{
    int d;
    ll n;
    scanf("%d%lld", &d, &n);
    C[0][0] = C[1][0] = C[1][1] = 1;
    for (int i = 2; i <= 50; i++) {
        C[i][0] = C[i][i] = 1;
        for (int j = 1; j < i; j++) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod;
    }
    ans.n = sta.n = d + 2, sta.op[1][1] = 1, sta.op[d + 2][1] = 1;
    for (int i = 2; i <= d + 2; i++) {
        for (int j = 2; j <= i; j++) sta.op[j][i] = C[i - 2][j - 2];
    }
    for (int i = 1; i <= d + 2; i++) ans.op[i][i] = 1;
    while (n) {
        if (n & 1) ans *= sta;
        sta *= sta, n >>= 1;
    }
    int tot = 0;
    for (int i = 2; i <= d + 2; i++) tot = (tot + ans.op[i][1]) % mod;
    printf("%d", tot);
    return 0;
}

解法2

        注意到：

$$(j+1)^{d+1}-j^{d+1}=\sum_{i=0}^dC_{d+1}^ij^i$$

        对其两边求和：

$$\sum_{j=0}^{n-1}((j+1)^{d+1}-j^{d+1})=\sum_{j=1}^{n-1}\sum_{i=0}^dC_{d+1}^ij^i$$

        即为：

$$n^{d+1}=\sum_{i=0}^dC_{d+1}^iS_i(n)$$

        将$S_d(n)$移出来：

$$(d+1)S_d(n)=n^{d+1}-\sum_{i=0}^{d-1}C_{d+1}^iS_i(n)$$

        之后就可以递推出$S_d(n)$，时间复杂度$O(d^2+dlogd)$。这种方法需要$2$到$d+1$总是与模数互质。
        本方法在$d=0$时需要特判掉，否则会导致答案多一。这是因为我们在定义$S_d(n)$时，本身就会导致在$d=0$时出现$0^0$的特殊情况。

#include <bits/stdc++.h>
using namespace std;
const int mod = 1e9 + 7;
typedef long long ll;
int C[1005][1005], S[1005];
inline int qpow(ll x, int y)
{
    int ans = 1, sta = x;
    while (y) {
        if (y & 1) ans = 1ll * ans * sta % mod;
        sta = 1ll * sta * sta % mod, y >>= 1;
    }
    return ans;
}
inline int inv(int x)
{
    return qpow(x, mod - 2);
}
int main()
{
    int d;
    ll n;
    scanf("%d%lld", &d, &n);
    C[0][0] = C[1][0] = C[1][1] = 1;
    for (int i = 2; i <= 1000; i++) {
        C[i][0] = C[i][i] = 1;
        for (int j = 1; j < i; j++) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod;
    }
    ++n, S[0] = n % mod;
    for (int i = 1; i <= d; i++) {
        S[i] = 0;
        for (int j = 0; j < i; j++) S[i] = (S[i] + 1ll * S[j] * C[i + 1][j] % mod) % mod;
        S[i] = (qpow(n, i + 1) - S[i]) % mod;
        S[i] = 1ll * S[i] * inv(i + 1) % mod;
    }
    if (d == 0) {  //very important!!
        printf("%d", S[0] - 1);
    } else {
        printf("%d", (S[d] + mod) % mod);
    }
    return 0;
}

解法3

        解法三需要用到伯努利数。设$B_i$表示第$i$个伯努利数，定义为：

$$\frac{x}{e^x-1}=\sum_{n=1}^\infty\frac{B_n}{n!}x^n$$

        定义伯努利多项式：

$$\sum_{n=0}^\infty\frac{\beta_n(t)}{n!}x^n=\frac{x}{e^x-1}e^{tx}=(\sum_{n=1}^\infty\frac{B_n}{n!}x^n)(\sum_{n=1}^\infty\frac{t^n}{n!}x^n)$$

        对右侧卷积得到：

$$\beta_n(t)=\sum_{k=0}^nC_n^kB_{n-k}t^k$$

        作差得到：

$$\sum_{n=0}^\infty\frac{\beta_n(t+1)-\beta_n(t)}{n!}x^n=\frac{xe^{(t+1)x}}{e^x-1}-\frac{xe^{tx}}{e^x-1}=xe^{tx}=\sum_{n=0}^\infty\frac{t^n}{n!}x^{n+1}$$

        对比系数可得：

$$\beta_{n+1}(t+1)-\beta_{n+1}(t)=(n+1)t^n$$

        求和：

$$S_d(n)=\frac{\beta_{d+1}(n)-\beta_{d+1}(0)}{d+1}$$

        然后代入上文得到的公式：

$$S_d(n)=\frac{1}{d+1}\sum_{i=0}^dC_{d+1}^iB_in^{d+1-i}$$

        在这个式子中，令$n=1,d\not=0$即可得到伯努利数递推公式：

$$B_n=-\frac{1}{n+1}\sum_{k=0}^{n-1}C_{n+1}^kB_k$$

        使用以上两式配合$B_0=1$递推边界可以得到在$O(d^2)$复杂度内处理组合数和伯努利数，进而$O(dlogd)$求解。本方法也需要在$d=0$时特判。
        在需要大量求和的时候，本方法由于进行了预处理，比解法2更加高效。

#include <bits/stdc++.h>

using namespace std;
const int mod = 1e9 + 7, N = 5000000;
typedef long long ll;
int C[60][60], S[60], B[60];
inline int qpow(int x, int y)
{
    int ans = 1, sta = x;
    while (y) {
        if (y & 1) ans = 1ll * ans * sta % mod;
        sta = 1ll * sta * sta % mod, y >>= 1;
    }
    return ans;
}
inline int inv(int x)
{
    return qpow(x, mod - 2);
}

int main()
{
    int d;
    ll n;
    scanf("%d%lld", &d, &n);
    C[0][0] = C[1][0] = C[1][1] = 1;
    for (int i = 2; i < 60; i++) {
        C[i][0] = C[i][i] = 1;
        for (int j = 1; j < i; j++) C[i][j] = (C[i - 1][j - 1] + C[i - 1][j]) % mod;
    }
    B[0] = 1;
    for (int i = 1; i < 60; i++) {
        for (int j = 0; j < i; j++) B[i] = (B[i] + 1ll * C[i + 1][j] * B[j] % mod) % mod;
        B[i] = -1ll * B[i] * inv(i + 1) % mod;
    }
    int s = 0;
    for (int i = 0; i <= d; i++) {
        s = (s + 1ll * C[d + 1][i] * B[i] % mod * qpow(n + 1, d + 1 - i) % mod) % mod;
    }
    s = 1ll * s * inv(d + 1) % mod;
    if (d == 0) {
        printf("%d\n", s - 1);
    } else {
        printf("%d\n", (s + mod) % mod);
    }
    return 0;
}