Skip to main content

__debug's Home Keep it simple, stupid

计算自然数幂和 II

经典问题. 新学习了几个做法.

(很久之前也写过一篇: 传送门)

问题

给定 \(n\), \(k\), 求 \(\sum_{i = 0}^{n}i^k\). 其中 \(n\) 非常大.

对质数 \(P\) 取模. (\(P\) 不是质数的时候可以参考以前写的那篇, 有一个 \(O(k^2)\) 的解法)

做法

#1. 多项式差分

(差分的相关定义可以参考 组合数学 或是 具体数学, 或者也可以看看 这篇博客)

考虑 \(f(n) = n^k\) 这一函数的差分表的第 \(0\) 条对角线 \(c_0 ... c_k\).

我们有

\begin{align*} \sum_{i = 0}^{n} f(i) &= \sum_{i = 0}^{n} \sum_{j = 0}^{k} c_j \binom{i}{j} \\ &= \sum_{j = 0}^{k} \sum_{i = 0}^{n} c_j \binom{i}{j} \\ &= \sum_{j = 0}^{k} c_j \binom{n + 1}{j + 1} \end{align*}

如果模数为质数, 就可以 \(O(k^2)\) 计算了.

不过 \(P\) 不是质数的时候不是很会做…

static long long D[MAXK][MAXK];

for (long long i = 0; i <= K; ++i) {
    D[0][i] = fpm(i, K);
}
for (long long i = 1; i <= K; ++i) {
    for (long long j = 0; j <= K - i; ++j) {
	D[i][j] = mod(D[i - 1][j + 1] - D[i - 1][j]);
    }
}

long long ans = 0;
long long clast = 1;
for (long long i = 0; i <= K; ++i) {
    clast = clast * (mod(N - i + 1) % MOD) % MOD * inv(i + 1) % MOD;
    modr(ans += D[i][0] * clast % MOD);
}

return ans;

#2. 拉格朗日插值

由上面那个做法我们可以知道 \(\sum_{i = 0}^{n} i^k\) 是一个关于 \(n\) 的 \(k + 1\) 次多项式. 拉格朗日插值一下不就好了吗…

我们把 \(n = 0 ... k + 1\) 的函数值算出来, 由于 \(n\) 的取值很特殊, 可以 \(O(k)\) 或 \(O(k \log k)\) 算出 \(n_0\) 时的函数值.

至于之前算函数值, 可以简单地用快速幂做到 \(O(k \log k)\), 也可以用线性筛做到 \(O(k)\).

下面的是 \(O(n \log n)\) 的.

static long long val[MAXK];

val[0] = 0;
for (long long i = 1; i <= K + 1; ++i) {
    val[i] = mod(val[i - 1] + fpm(i, K));
}

if (N <= K + 1)
    return val[N];

static long long pre[MAXK], suf[MAXK];

pre[0] = N % MOD;
for (long long i = 1; i <= K + 1; ++i) {
    pre[i] = pre[i - 1] * ((N - i) % MOD) % MOD;
}
suf[K + 2] = 1;
for (long long i = K + 1; i >= 1; --i) {
    suf[i] = suf[i + 1] * ((N - i) % MOD) % MOD;
}

long long ans = 0;
for (long long i = 0; i <= K + 1; ++i) {
    long long t = (i ? pre[i - 1] : 1) * suf[i + 1] % MOD * ifac[i] % MOD * ifac[K + 1 - i] % MOD * val[i] % MOD;
    if ((K + 1 - i) & 1)
	modr(ans -= t);
    else
	modr(ans += t);
}

return ans;

#3. Bernoulli 数

不想填了…

Comments

Comments powered by Disqus