Library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub ebi-fly13/Library

:heavy_check_mark: $\sum_{i = 0}^{n-1} i^k$
(math/sum_of_powers_iota.hpp)

説明

$\sum_{i = 0}^{n-1} i ^ k$ を $O(k)$ で求める。

Depends on

Verified with

Code

#pragma once

#include <cassert>

#include "../math/lagrange_interpolation.hpp"
#include "../math/linear_sieve.hpp"
#include "../modint/base.hpp"

namespace ebi {

template <Modint mint> mint sum_of_powers_iota(long long n, int k) {
    assert(n > 0 && k >= 0);
    linear_sieve sieve(k + 1);
    auto pow_table = sieve.pow_table<mint>(k + 1, k);
    for (int i = 0; i < k + 1; i++) {
        pow_table[i + 1] += pow_table[i];
    }
    return lagrange_interpolation(pow_table, n - 1);
}

}  // namespace ebi
#line 2 "math/sum_of_powers_iota.hpp"

#include <cassert>

#line 2 "math/lagrange_interpolation.hpp"

#include <vector>

/*
    reference: https://atcoder.jp/contests/abc208/editorial/2195
    verify: https://atcoder.jp/contests/abc208/tasks/abc208_f
*/

namespace ebi {

template <class mint>
mint lagrange_interpolation(const std::vector<mint> &f, long long n) {
    const int d = int(f.size()) - 1;  // Nのd次以下の多項式
    mint fact = 1;
    std::vector<mint> inv_fact(d + 1);
    for (int i = 1; i < d + 1; ++i) {
        fact *= i;
    }
    inv_fact[d] = fact.inv();
    for (int i = d; i > 0; i--) {
        inv_fact[i - 1] = inv_fact[i] * i;
    }
    std::vector<mint> l(d + 1), r(d + 1);
    l[0] = 1;
    for (int i = 0; i < d; ++i) {
        l[i + 1] = l[i] * (n - i);
    }
    r[d] = 1;
    for (int i = d; i > 0; --i) {
        r[i - 1] = r[i] * (n - i);
    }
    mint res = 0;
    for (int i = 0; i < d + 1; ++i) {
        res += mint((d - i) % 2 == 1 ? -1 : 1) * f[i] * l[i] * r[i] *
               inv_fact[i] * inv_fact[d - i];
    }
    return res;
}

}  // namespace ebi
#line 2 "math/linear_sieve.hpp"

#line 2 "modint/base.hpp"

#include <concepts>
#include <iostream>
#include <utility>

namespace ebi {

template <class T>
concept Modint = requires(T a, T b) {
    a + b;
    a - b;
    a * b;
    a / b;
    a.inv();
    a.val();
    a.pow(std::declval<long long>());
    T::mod();
};

template <Modint mint> std::istream &operator>>(std::istream &os, mint &a) {
    long long x;
    os >> x;
    a = x;
    return os;
}

template <Modint mint>
std::ostream &operator<<(std::ostream &os, const mint &a) {
    return os << a.val();
}

}  // namespace ebi
#line 2 "template/int_alias.hpp"

#include <cstdint>

namespace ebi {

using ld = long double;
using std::size_t;
using i8 = std::int8_t;
using u8 = std::uint8_t;
using i16 = std::int16_t;
using u16 = std::uint16_t;
using i32 = std::int32_t;
using u32 = std::uint32_t;
using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

}  // namespace ebi
#line 5 "math/linear_sieve.hpp"

/*
    reference: https://37zigen.com/linear-sieve/
    verify:    https://atcoder.jp/contests/abc162/submissions/25095562
*/

#line 13 "math/linear_sieve.hpp"

namespace ebi {

struct linear_sieve {
  private:
    using u64 = std::uint64_t;
    int n;
    std::vector<int> sieve;
    std::vector<int> prime;

  public:
    linear_sieve(int _n) : n(_n), sieve(std::vector<int>(_n + 1, -1)) {
        for (int i = 2; i <= n; i++) {
            if (sieve[i] < 0) {
                sieve[i] = i;
                prime.emplace_back(i);
            }
            for (auto p : prime) {
                if (u64(p) * u64(i) > u64(n) || p > sieve[i]) break;
                sieve[p * i] = p;
            }
        }
    }

    std::vector<int> prime_table() const {
        return prime;
    }

    std::vector<std::pair<int, int>> prime_power_table(int m) const {
        assert(m <= n);
        std::vector<std::pair<int, int>> table(m + 1, {1, 1});
        for (int i = 2; i <= m; i++) {
            int p = sieve[i];
            table[i] = {p, p};
            if (sieve[i / p] == p) {
                table[i] = table[i / p];
                table[i].second *= p;
            }
        }
        return table;
    }

    std::vector<std::pair<int, int>> factorize(int x) {
        assert(x <= n);
        std::vector<std::pair<int, int>> res;
        while (x > 1) {
            int p = sieve[x];
            int exp = 0;
            if (p < 0) {
                res.emplace_back(x, 1);
                break;
            }
            while (sieve[x] == p) {
                x /= p;
                exp++;
            }
            res.emplace_back(p, exp);
        }
        return res;
    }

    std::vector<int> divisors(int x) {
        assert(x <= n);
        std::vector<int> res;
        res.emplace_back(1);
        auto pf = factorize(x);
        for (auto p : pf) {
            int sz = (int)res.size();
            for (int i = 0; i < sz; i++) {
                int ret = 1;
                for (int j = 0; j < p.second; j++) {
                    ret *= p.first;
                    res.emplace_back(res[i] * ret);
                }
            }
        }
        return res;
    }

    template <class T> std::vector<T> fast_zeta(const std::vector<T> &f) {
        std::vector<T> F = f;
        int sz = f.size();
        assert(sz <= n + 1);
        for (int i = 2; i < sz; i++) {
            if (sieve[i] != i) continue;
            for (int j = (sz - 1) / i; j >= 1; j--) {
                F[j] += F[j * i];
            }
        }
        return F;
    }

    template <class T> std::vector<T> fast_mobius(const std::vector<T> &F) {
        std::vector<T> f = F;
        int sz = F.size();
        assert(sz <= n + 1);
        for (int i = 2; i < sz; i++) {
            if (sieve[i] != i) continue;
            for (int j = 1; j * i < sz; j++) {
                f[j] -= f[j * i];
            }
        }
        return f;
    }

    template <Modint mint> std::vector<mint> pow_table(int m, int k) {
        assert(m <= n && k >= 0);
        std::vector<mint> table(m + 1, 1);
        table[0] = (k == 0);
        for (int i = 2; i <= m; i++) {
            if (sieve[i] == i) {
                table[i] = mint(i).pow(k);
                continue;
            }
            table[i] = table[sieve[i]] * table[i / sieve[i]];
        }
        return table;
    }

    template <Modint mint> std::vector<mint> inv_table() {
        return pow_table(mint::mod() - 2);
    }
};

}  // namespace ebi

#line 8 "math/sum_of_powers_iota.hpp"

namespace ebi {

template <Modint mint> mint sum_of_powers_iota(long long n, int k) {
    assert(n > 0 && k >= 0);
    linear_sieve sieve(k + 1);
    auto pow_table = sieve.pow_table<mint>(k + 1, k);
    for (int i = 0; i < k + 1; i++) {
        pow_table[i + 1] += pow_table[i];
    }
    return lagrange_interpolation(pow_table, n - 1);
}

}  // namespace ebi
Back to top page