无所畏惧的求和题解

本题是本人目前出的题中难度最高的题。

加强版紫色是肯定可以的。

题目链接：无所畏惧的求和 - 洛谷

尽情享受吧！

这道题其实做法有很多：

待定系数法 + 矩阵求解推代数公式（不够优秀）
组合数学 + 待定系数法推组合公式（受限于 $k$ 的范围）
在模意义下的特殊做法（受限于模数范围）
拉格朗日插值法（正解）
第一类斯特林数（时间复杂度可能有点问题，为 $O(k^2)$
伯努利数（奇奇怪怪）
……

那么这里讲解前四种方式

无所畏惧的求和题解

代数公式

自然数幂方求和公式？在高等教育出版社出版的《数学手册》中有这么一些公式：

采用万能的数学归纳法可以一一证明上述公式。此处不提。

但是观察上述公式，可以发现一个特征：自然数 $k$ 次幂求和公式是关于 $n$ 的 $k + 1$ 次有理多项式。

也就是

\sum_{i = 1}^{n} i^{k} = \sum_{j = 1}^{k + 1} a_{j} n^{j}

$\sum_{i=1}^n i^k = \sum_{j=1}^{k+1} a_j n^{j}$

知道上述结果之后，可采用待定系数法，也就是写出 $n = 1, 2, 3, \dots, k+1$ 的 $k + 1$ 个代数式，利用矩阵求解即可。

举个例子，对于 $k = 3$ 的情况：

\begin{aligned} \sum_{i = 1}^{1} i^{3} & = a_{1} + a_{2} + a_{3} + a_{4} & = 1 \\ \sum_{i = 1}^{2} i^{3} & = 2 a_{1} + 4 a_{2} + 8 a_{3} + 16 a_{4} & = 9 \\ \sum_{i = 1}^{3} i^{3} & = 3 a_{1} + 9 a_{2} + 27 a_{3} + 81 a_{4} & = 36 \\ \sum_{i = 1}^{4} i^{3} & = 4 a_{1} + 16 a_{2} + 64 a_{3} + 256 a_{4} & = 100 \end{aligned}

$\begin{aligned} \sum_{i=1}^1 i^3 &= a_1 + a_2 + a_3 + a_4 &= 1 \\ \sum_{i=1}^2 i^3 &= 2a_1 + 4a_2 + 8a_3 + 16a_4 &= 9 \\ \sum_{i=1}^3 i^3 &= 3a_1 + 9a_2 + 27a_3 + 81a_4 &= 36 \\ \sum_{i=1}^4 i^3 &= 4a_1 + 16a_2 + 64a_3 + 256a_4 &= 100 \\ \end{aligned}$

那么借此求出每一项的系数即可对于每一个询问在 $O(k)$ 的复杂度内完成计算。

在获取系数的部分：

vector<int> mat[K];
void get_coefficient(int k, vector<int> &cctk) {
    int h = k + 1, w = k + 2;
    cctk.assign(h + 1, 0);
    for (int i(1); i <= h; ++i) mat[i].assign(w + 1, 0);

    int sum(0);
    for (int n(1); n <= h; ++n) {
        (sum += qpow(n, k)) %= mod;

        for (int x(n), j(1); j < w; ++j, (x *= n) %= mod) {
            mat[n][j] = x;
        } mat[n][w] = sum;
    }

    // Guass 消元    
    for (int i(1); i <= h; ++i) {
        int inv = qpow(mat[i][i], mod - 2);
        for (int j(i); j <= w; ++j) 
            (mat[i][j] *= inv) %= mod;

        for (int j = 1; j <= h; ++j) {
            if (j ^ i) {
                inv = mat[j][i];
                for (int c(i); c <= w; ++c) {
                    (mat[j][c] -= mat[i][c] * inv) %= mod;
                    if (mat[j][c] < 0) mat[j][c] += mod;
                }
            }
        }
    }

    for (int i = 1; i <= h; ++i) {
        cctk[i] = mat[h + 1 - i][w];
    }
}

在最终求答案的部分：

void work() {
    int T;
    cin >> T;

    while (T--) {
        int n, k;
        cin >> n >> k;

        vector<int> &cctk = ccts[k];
        if (cctk.empty())
            get_coefficient(k, cctk);

        int ans(0);
        n %= mod;
        for (int i = 1; i <= k + 1; ++i) {
            (ans = ans * n % mod + cctk[i]) %= mod;
        }
        cout << ans * n % mod << '\n';
    }
}

总的复杂度为 $O(k^4 + T \cdot k)$ 。~~还不够优秀，~~普通版只能过前10个点，期望得分30，加强版无法得分。

但是简单版更快……

组合公式

这个方法相对优秀一点。标程就是用的此写法。

其实不难发现，对于 $x^k$ ，我们可以改写为：

x^{k} = \sum_{i = 1}^{k} a_{i} (\binom{x}{i})

$x^k = \sum_{i=1}^{k} a_i\binom{x}i$

那么依据某些公式推导：

\sum_{x = 1}^{n} x^{k} = \sum_{i = 2}^{k + 1} a_{i - 1} (\binom{n + 1}{i})

$\sum_{x=1}^n x^k = \sum_{i = 2}^{k + 1} a_{i-1} \binom {n+1}i$

所以，类似的，我们也可以枚举 $n = 1, 2, 3, \dots, k$ 来寻找其系数：

以 $k = 3$ 为例

\begin{aligned} \sum_{x = 1}^{1} & = a_{1} (\binom{2}{2}) + a_{2} (\binom{2}{3}) + a_{3} (\binom{2}{4}) = 1 \\ \sum_{x = 1}^{2} & = a_{1} (\binom{3}{2}) + a_{2} (\binom{3}{3}) + a_{3} (\binom{3}{4}) = 9 \\ \sum_{x = 1}^{3} & = a_{1} (\binom{4}{2}) + a_{2} (\binom{4}{3}) + a_{3} (\binom{4}{4}) = 36 \end{aligned}

$\begin{aligned} \sum_{x=1}^1 &= a_1 \binom 22 + a_2 \binom 23 + a_3 \binom 24 = 1 \\ \sum_{x=1}^2 &= a_1 \binom 32 + a_2 \binom 33 + a_3 \binom 34 = 9 \\ \sum_{x=1}^3 &= a_1 \binom 42 + a_2 \binom 43 + a_3 \binom 44 = 36 \\ \end{aligned}$

同时我们规定 $\binom nr = 0\ (n \lt r)$ 。所以上式也可以写为

\begin{aligned} \sum_{x = 1}^{1} & = a_{1} (\binom{2}{2}) & = 1 \\ \sum_{x = 1}^{2} & = a_{1} (\binom{3}{2}) + a_{2} (\binom{3}{3}) & = 9 \\ \sum_{x = 1}^{3} & = a_{1} (\binom{4}{2}) + a_{2} (\binom{4}{3}) + a_{3} (\binom{4}{4}) & = 36 \end{aligned}

$\begin{aligned} \sum_{x=1}^1 &= a_1 \binom 22 &= 1 \\ \sum_{x=1}^2 &= a_1 \binom 32 + a_2 \binom 33 &= 9 \\ \sum_{x=1}^3 &= a_1 \binom 42 + a_2 \binom 43 + a_3 \binom 44 &= 36 \\ \end{aligned}$

这就是一个下三角矩阵，每一次扫一遍即可。

代码也非常简单，常数比第一种方法小很多：

void get_coefficient(int k, int * ccts) {
    int sum = 0;
    for (int n = 1; n <= k; ++n) {
        sum += qpow(n, k, MOD);

        int rest = sum;
        // 由于我们已经知道了前 (n-1) 个系数，直接通过总和一一减去即可。
        for (int i = 1; i < n; ++i) {
            // 注意 k < MOD，所以此处不需要Lucas
            (rest -= ccts[i] * C(n + 1, i + 1) % MOD) %= MOD;
        }
        if (rest < 0) rest += MOD;
        // 明显可知，C(n, n) = 1，所以剩下的即是系数
        ccts[n] = rest;
    }
}

核心部分也非常简单，只是模数较小，需要用到 $Lucas$ 定理。

int ccts[K][K]; // 用于保存系数
void work() {
    int T, n, k;
    cin >> T;

    while (T--) {
        cin >> n >> k;

        int * cctk = ccts[k];
        if (!cctk[1]) // 其实不难发现，a1 一定为 1，所以借此判断即可
            get_coefficient(k, cctk);

        int ans = 0;
        for (int i = 1; i <= k; ++i)
            (ans += cctk[i] * Lucas(n + 1, i + 1)) %= MOD;
        cout << ans << '\n';
    }
}

总的复杂度为 $O(k^3 + T \cdot \log_pn)$ ，由于 $p = 111121$ ，所以我们可以认为复杂度为 $O(2k^3 + 2T)$ 。普通版期望得分100分，加强版受限于 $k$ 的大小，期望得分30分。

本题特殊做法

由于数据原因需要取模……所以有了这么一个做法。

做法来源于 fanghaoyu

在模意义下有如下恒等式：

x^{k} = (x + p)^{k} \mod p

$x^k = (x+p)^k \mod p$

意味着

\sum_{x = 1}^{n} x^{k} \equiv ⌊ \frac{n}{p} ⌋ \sum_{x = 1}^{p} x^{k} + \sum_{x = 1}^{n % p} x^{k} (\mod p)

$\sum_{x=1}^n x^k \equiv \lfloor \frac np \rfloor \sum_{x=1}^p x^k + \sum_{x = 1}^{n \% p} x^k \pmod p$

那么，很明显，我们可以预处理出 $\lfloor \frac np \rfloor \sum_{x=1}^p x^k$

然后对于每一次求解求后一项即可，每一次查询时间为 $O(p \log k)$ 。

当然，也可以采用空间换时间的方法，保存后一项，使得查询时间复杂度为 $O(1)$ 。

而预处理需要 $O(p)$ ，所以总的复杂度大概为 $O(p k\log k + T p\log k)$ 或者 $O(p k\log k + T)$

而空间复杂度……

其实这个方法特别好卡，当初我为了提示使用 Lucas 定理，专门设置的小质数，也就成了这个做法的来源，所以只需要小小的修改一下模数，改大一点，这就熄火了 QwQ。

为了尊重这个方法的想出者，所以贴出未修改的源代码。

#include <bits/stdc++.h>
using namespace std;
const int N = 1e5 + 5,MOD = 8887;
typedef long long ll;
ll n,k,a[1001][MOD + 5];
inline ll ksm(ll base,ll pts)
{
    ll ret = 1;
    for(;pts > 0;pts >>= 1,base = base * base % MOD)
        if(pts & 1)
            ret = ret * base % MOD;
    return ret;
}
int main()
{
    ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
    int T;
    cin>>T;
    for(int t = 1;t <= 1000;t++)
        for(int i = 1;i <= MOD;i++)
            a[t][i] = a[t][i - 1] + ksm(i,t),a[t][i] %= MOD;
    while(T--)
    {
        cin>>n>>k;
        ll ans = 0;
        ans = ans + a[k][MOD] * (n / MOD) % MOD;
        ans = ans + a[k][n % MOD];
        cout<<ans<<endl;
    }
    return 0;
}

拉格朗日插值定理

在初中学过，两点确定一线（一次多项式），三点可以确定一个抛物线（二次多项式）。扩展出来， $n + 1$ 个点可以确定 $n$ 次多项式。

在上文中提及，自然数 $k$ 次幂求和公式是关于 $n$ 的 $k + 1$ 次有理多项式。

那么我们只需要确定 $k + 2$ 个点，就可以求出这个多项式。加之我们只需要求在 $n$ 处的值，所以想到使用拉格朗日插值定理求解。

定理

首先，对于一个 $k$ 次多项式，我们有了 $k + 1$ 个点 $(x_i, y_i)$ 。

定义：

f_{i} (x_{j}) = {\begin{cases} 1 & i = j \\ 0 & i \neq j \end{cases}

$f_i(x_j) = \begin{cases} 1 & i = j \\ 0 & i \ne j \end{cases}$

容易构造出

f_{i} (x) = \frac{\prod_{i \neq j} (x - x_{j})}{\prod_{i \neq j} (x_{i} - x_{j})}

$f_i(x) = \cfrac {\prod_{i \ne j} (x - x_j)}{\prod_{i \ne j}(x_i - x_j)}$

因此最终我们得到：

f (x) = \sum_{i = 1}^{k + 1} y_{i} f i (x)

$f(x) = \sum_{i=1}^{k+1} y_i fi(x)$

那么考虑为什么？

由于我们对于每一个 $y_if_i(x)$ ，可以保证在 $x_i$ 点处，取值为 $y_i$ ，其余的点上取值为 $0$ 。

那么 $f(x) = \sum_{i=1}^{k+1} y_i fi(x)$ 可以一一穿过这三个点。

也就是说，我们可以借此求得 $f(x)$ 了。

更严谨的证明此处略过……

连续值求解

很明显，对于每一个分母，是两个阶乘的积，所以预处理阶乘。

对于分母，预处理前缀积和后缀积即可。

核心代码如下：

typedef long long lint;
lint fac[K], rfac[K];

inline lint qpow(lint a, int x) {
    lint r = 1;
    for (; x; x >>= 1, a = a * a % mod)
        if (x & 1) r = r * a % mod;
    return r;
}

inline void prepare() {
    fac[0] = rfac[0] = 1;
    // 预处理下方阶乘
    for (int i = 1; i < K; ++i) {
        fac[i] = fac[i - 1] * i % mod;
        // rfac[i] = rfac[i - 1] * (-i) % mod;
        rfac[i] = rfac[i - 1] * (mod - i) % mod;
    }
}

// prefix, suffix
lint pfx[K], sfx[K];
int calc(int n, int k) {
    pfx[0] = sfx[k + 3] = 1;
    for (int i = 1; i <= k + 2; ++i)
        pfx[i] = pfx[i - 1] * (n - i) % mod;
    for (int i = k + 2; i; --i)
        sfx[i] = sfx[i + 1] * (n - i) % mod;

    lint y = 0, ans = 0, frup, frdown;
    for (int i = 1; i <= k + 2; ++i) {
        y = (y + qpow(i, k)) % mod;
        frup = pfx[i - 1] * sfx[i + 1] % mod;
        frdown = fac[i - 1] * rfac[k + 2 - i] % mod;
        ans = (ans + y * frup % mod * qpow(frdown, mod - 2) % mod) % mod;
    }

    return (ans + mod) % mod;
}

时间复杂度： $O(\sum k)$ ，期望得分 100。

posted @ 2023-03-31 17:42 jeefy 阅读(87) 评论(3) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

相关博文：

· 算法学习笔记(40): 具体数学

· 算法学习笔记(16): 组合数学基础

· LOJ3075 「2019 集训队互测 Day 3」组合数求和

· 刷题笔记（2023.9.21）

· 组合数学初步

阅读排行：
· TypeScript + Deepseek 打造卜卦网站：技术与玄学的结合
· 阿里巴巴 QwQ-32B真的超越了 DeepSeek R-1吗？
· 如何调用 DeepSeek 的自然语言处理 API 接口并集成到在线客服系统
· 【译】Visual Studio 中新的强大生产力特性
· 2025年我用 Compose 写了一个 Todo App

公告

昵称： jeefy
园龄： 2年2个月
粉丝： 32
关注： 13

+加关注

2025年3月

日

一

二

三

四

五

六

jeefy

信仰，梦想，希望！

无所畏惧的求和题解

无所畏惧的求和题解

代数公式

组合公式

本题特殊做法

拉格朗日插值定理

定理

连续值求解

公告

最新随笔

我的标签

合集 (7)

随笔分类 (105)

阅读排行榜

推荐排行榜