yukicoder No.665 Bernoulli Bernoulli - アルゴリズマーの備忘録

「3月中に毎日1問解いて黄色になろう！」キャンペーン5問目。結局3日に一問ペースだけど気にしない。

問題概要

1^k+2^k+3^k+...+n^k(mod 10⁹+7)を求めよ。
N≦10¹⁶、k≦10⁴

解説

何か他の人の解法はみんなベルヌーイ数とか使っているみたいだけど、自分は高校数学で解けることにすぐ気付いてしまったのでググることがなかったです。（多分やってることは同じ）

当然ながら1からnまで順番に足すことは出来ないので、Σxⁱでiを更新できないか考えます。

高校数学では2乗和の公式や3乗和の公式が出てきます。これの求め方は、例えば3乗和では、
$n^4-(n-1)^4=4n^3-6n^2+4n-1$
$(n-1)^4-(n-2)^4=4(n-1)^3-6(n-1)^2+4(n-1)-1$
$(n-2)^4-(n-3)^4=4(n-2)^3-6(n-2)^2+4(n-2)-1$
$\vdots$
$1^4-0^4=4\times 1^3-6\times 1^2+4\times 1-1$
これらを全部足し合わせると、左辺が消えまくって、
$\displaystyle{ n^4=4\sum_{i=1}^ni^3-6\sum_{i=1}^ni^2+4\sum_{i=1}^ni-n }$
$\displaystyle{ \sum_{i=1}^ni^3=\frac{n^4+6\sum_{i=1}^ni^2-4\sum_{i=1}^ni+n}{4} }$ となります。

この式変形を一般のkに対して行うと、
$\displaystyle{ \sum_{i=1}^ni^k=\frac{ n^{k+1}+\sum_{i=0}^{k-1}\left((-1)^{k+i+1}\cdot_{k+1}C_i\sum_{x=1}^nx^i\right) }{k+1} }$ となります。

この各成分について見てみると、「Σxⁱを記録した配列」「k+1までの10⁹+7を法とした逆元」「k+1までのi!とその逆元」が分かれば、O(k)で計算できることが分かります。全体でO(k²)ですね。

私はΣxⁱの配列に最初i=1まで入れて（多分modの辺りで）WAを連発したので、人力でやらせるのはダメです。

https://yukicoder.me/submissions/244064

// 1^x+2^x+...+n^x
ll dp[10010];
// xの逆元
ll inv[10010];
// x!
ll fac[10010];
// x!の逆元
ll invfac[10010];

ll modpow(ll x,ll m) {
    x %= mod;
    if (x) {
        ll ret = 1;
        while (m) {
            if (m & 1)   ret = (ret*x) % mod;
            x = (x*x) % mod;
            m >>= 1;
        }
        return ret;
    }
    return 0;
}

ll nCm(int n, int m) {
    return (fac[n] * ((invfac[n - m] * invfac[m]) % mod)) % mod;
}

int main() {
    ll n, k;    cin >> n >> k;
    inv[0] = inv[1] = fac[0] = fac[1] = invfac[0] = invfac[1] = 1;
    SREP(i, 2, k + 1) {
        inv[i] = modpow(i,mod-2);
        fac[i] = (i*fac[i - 1]) % mod;
        invfac[i] = (inv[i] * invfac[i - 1]) % mod;
    }
    dp[0] = n%mod;
    SREP(i, 1, k) {
        ll nxt = modpow(n, i + 1);
        rep(j, i) {
            nxt += mod + ((i + j + 1) & 1 ? -1 : 1)*((nCm(i + 1, j)*dp[j]) % mod);
            nxt %= mod;
        }
        dp[i] = (nxt*inv[i + 1]) % mod;
    }
    cout << dp[k] << endl;
    return 0;
}

感想

バグ埋め込んで1日が解けているのでダメですね…