$\sum_{i = 0}^{n-1} a^i f(i)$
(math/sum_of_exp_times_poly.hpp)

• View this file on GitHub
• Last update: 2024-07-27 23:45:45+09:00
• Include: #include "math/sum_of_exp_times_poly.hpp"

説明

$K$ 次多項式 $f(n)$ について、 $f(0), f(1), \dots, f(k)$ が与えられる。 $\sum_{i = 0}^{n-1} a^i f(i)$ を $O(K + \log{mod})$ で求める。

sum_of_exp_times_poly(f, a, n)

$K$ 次多項式 $f(n)$ の $f(0), f(1), \dots, f(k)$ と $a, n$ を与えて $\sum_{i = 0}^{n-1} a^i f(i)$ を求める。

sum_of_exp_times_poly_limit(f, a)

$K$ 次多項式 $f(n)$ の $f(0), f(1), \dots, f(k)$ と $a$ を与えて $\sum_{i = 0}^{\infty} a^i f(i)$ を求める。ここで、 $a$ は $-1 < a < 1$ であるとする（$a \neq 1\pmod p$ である）。

sum_of_exp2(r, d, n)

$f(n) = n^k$ について求める。つまり、 $\sum_{i = 0}^{n-1} r^i i^d$ を求める。

sum_of_exp2_limit(r, d, n)

$f(n) = n^k$ について求める。つまり、 $\sum_{i = 0}^{\infty} r^i i^d$ を求める。

Depends on

• Binomial Coefficient (math/binomial.hpp)
• Lagrange Interpolation (math/lagrange_interpolation.hpp)
• Linear Sieve (math/linear_sieve.hpp)
• modint/base.hpp
• template/int_alias.hpp

Verified with

• test/math/Sum_of_Exponential_Times_Polynomial.test.cpp
• test/math/Sum_of_Exponential_Times_Polynomial_Limit.test.cpp

Code

#pragma once

#include <cassert>
#include <vector>

#include "../math/binomial.hpp"
#include "../math/lagrange_interpolation.hpp"
#include "../math/linear_sieve.hpp"
#include "../modint/base.hpp"
#include "../template/int_alias.hpp"

namespace ebi {

template <Modint mint>
mint sum_of_exp_times_poly(const std::vector<mint> &f, mint a, i64 n) {
    if (n == 0) return 0;
    if (a == 0) return f[0];
    if (a == 1) {
        std::vector<mint> g(f.size() + 1, 0);
        for (int i = 1; i < (int)g.size(); i++) {
            g[i] = g[i - 1] + f[i - 1];
        }
        return lagrange_interpolation(g, n);
    }
    int k = (int)f.size() - 1;
    Binomial<mint> binom(k + 1);
    std::vector<mint> g(k + 1, 0);
    {
        mint pow_a = 1;
        for (int i = 0; i < k + 1; i++) {
            g[i] = f[i] * pow_a;
            pow_a *= a;
        }
        for (int i = 0; i < k; i++) {
            g[i + 1] += g[i];
        }
    }
    mint c = 0;
    {
        mint pow_neg_a = 1;
        for (int i = 0; i < k + 1; i++) {
            c += binom.c(k + 1, i) * g[k - i] * pow_neg_a;
            pow_neg_a *= -a;
        }
    }
    c /= (1 - a).pow(k + 1);

    {
        mint inv_a_pow = 1, inv_a = a.inv();
        for (int i = 0; i < k + 1; i++) {
            g[i] = (-c + g[i]) * inv_a_pow;
            inv_a_pow *= inv_a;
        }
    }
    mint tn = lagrange_interpolation(g, n - 1);
    return tn * a.pow(n - 1) + c;
}

template <Modint mint>
mint sum_of_exp_times_poly_limit(const std::vector<mint> &f, mint a) {
    assert(a != 1);
    if (a == 0) return f[0];
    int k = (int)f.size() - 1;
    Binomial<mint> binom(k + 1);
    std::vector<mint> g(k + 1, 0);
    {
        mint pow_a = 1;
        for (int i = 0; i < k + 1; i++) {
            g[i] = f[i] * pow_a;
            pow_a *= a;
        }
        for (int i = 0; i < k; i++) {
            g[i + 1] += g[i];
        }
    }
    mint c = 0;
    {
        mint pow_neg_a = 1;
        for (int i = 0; i < k + 1; i++) {
            c += binom.c(k + 1, i) * g[k - i] * pow_neg_a;
            pow_neg_a *= -a;
        }
    }
    c /= (1 - a).pow(k + 1);
    return c;
}

template <Modint mint> mint sum_of_exp2(mint r, int d, i64 n) {
    linear_sieve sieve(d);
    auto f = sieve.pow_table<mint>(d, d);
    return sum_of_exp_times_poly(f, r, n);
}

template <Modint mint> mint sum_of_exp2_limit(mint r, int d) {
    linear_sieve sieve(d);
    auto f = sieve.pow_table<mint>(d, d);
    return sum_of_exp_times_poly_limit(f, r);
}

}  // namespace ebi

#line 2 "math/sum_of_exp_times_poly.hpp"

#include <cassert>
#include <vector>

#line 2 "math/binomial.hpp"

#include <bit>
#line 5 "math/binomial.hpp"
#include <cstdint>
#include <iostream>
#include <ranges>
#line 9 "math/binomial.hpp"

#line 2 "modint/base.hpp"

#include <concepts>
#line 5 "modint/base.hpp"
#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 11 "math/binomial.hpp"

namespace ebi {

template <Modint mint> struct Binomial {
  private:
    static void extend(int len = -1) {
        int sz = (int)fact.size();
        if (len < 0)
            len = 2 * sz;
        else if (len <= sz)
            return;
        else
            len = std::max(2 * sz, (int)std::bit_ceil(std::uint32_t(len)));
        len = std::min(len, mint::mod());
        assert(sz <= len);
        fact.resize(len);
        inv_fact.resize(len);
        for (int i : std::views::iota(sz, len)) {
            fact[i] = fact[i - 1] * i;
        }
        inv_fact[len - 1] = fact[len - 1].inv();
        for (int i : std::views::iota(sz, len) | std::views::reverse) {
            inv_fact[i - 1] = inv_fact[i] * i;
        }
    }

  public:
    Binomial() = default;

    Binomial(int n) {
        extend(n + 1);
    }

    static mint f(int n) {
        if (n >= (int)fact.size()) [[unlikely]] {
            extend(n + 1);
        }
        return fact[n];
    }

    static mint inv_f(int n) {
        if (n >= (int)fact.size()) [[unlikely]] {
            extend(n + 1);
        }
        return inv_fact[n];
    }

    static mint c(int n, int r) {
        if (r < 0 || n < r) return 0;
        return f(n) * inv_f(r) * inv_f(n - r);
    }

    static mint neg_c(int k, int d) {
        assert(d > 0);
        return c(k + d - 1, d - 1);
    }

    static mint p(int n, int r) {
        if (r < 0 || n < r) return 0;
        return f(n) * inv_f(n - r);
    }

    static mint catalan_number(int n) {
        return c(2 * n, n) * inv(n + 1);
    }

    static mint inv(int n) {
        return inv_f(n) * f(n - 1);
    }

    static void reserve(int n) {
        extend(n + 1);
    }

  private:
    static std::vector<mint> fact, inv_fact;
};

template <Modint mint>
std::vector<mint> Binomial<mint>::fact = std::vector<mint>(2, 1);

template <Modint mint>
std::vector<mint> Binomial<mint>::inv_fact = std::vector<mint>(2, 1);

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

#line 4 "math/lagrange_interpolation.hpp"

/*
    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 "template/int_alias.hpp"

#line 4 "template/int_alias.hpp"

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 11 "math/sum_of_exp_times_poly.hpp"

namespace ebi {

template <Modint mint>
mint sum_of_exp_times_poly(const std::vector<mint> &f, mint a, i64 n) {
    if (n == 0) return 0;
    if (a == 0) return f[0];
    if (a == 1) {
        std::vector<mint> g(f.size() + 1, 0);
        for (int i = 1; i < (int)g.size(); i++) {
            g[i] = g[i - 1] + f[i - 1];
        }
        return lagrange_interpolation(g, n);
    }
    int k = (int)f.size() - 1;
    Binomial<mint> binom(k + 1);
    std::vector<mint> g(k + 1, 0);
    {
        mint pow_a = 1;
        for (int i = 0; i < k + 1; i++) {
            g[i] = f[i] * pow_a;
            pow_a *= a;
        }
        for (int i = 0; i < k; i++) {
            g[i + 1] += g[i];
        }
    }
    mint c = 0;
    {
        mint pow_neg_a = 1;
        for (int i = 0; i < k + 1; i++) {
            c += binom.c(k + 1, i) * g[k - i] * pow_neg_a;
            pow_neg_a *= -a;
        }
    }
    c /= (1 - a).pow(k + 1);

    {
        mint inv_a_pow = 1, inv_a = a.inv();
        for (int i = 0; i < k + 1; i++) {
            g[i] = (-c + g[i]) * inv_a_pow;
            inv_a_pow *= inv_a;
        }
    }
    mint tn = lagrange_interpolation(g, n - 1);
    return tn * a.pow(n - 1) + c;
}

template <Modint mint>
mint sum_of_exp_times_poly_limit(const std::vector<mint> &f, mint a) {
    assert(a != 1);
    if (a == 0) return f[0];
    int k = (int)f.size() - 1;
    Binomial<mint> binom(k + 1);
    std::vector<mint> g(k + 1, 0);
    {
        mint pow_a = 1;
        for (int i = 0; i < k + 1; i++) {
            g[i] = f[i] * pow_a;
            pow_a *= a;
        }
        for (int i = 0; i < k; i++) {
            g[i + 1] += g[i];
        }
    }
    mint c = 0;
    {
        mint pow_neg_a = 1;
        for (int i = 0; i < k + 1; i++) {
            c += binom.c(k + 1, i) * g[k - i] * pow_neg_a;
            pow_neg_a *= -a;
        }
    }
    c /= (1 - a).pow(k + 1);
    return c;
}

template <Modint mint> mint sum_of_exp2(mint r, int d, i64 n) {
    linear_sieve sieve(d);
    auto f = sieve.pow_table<mint>(d, d);
    return sum_of_exp_times_poly(f, r, n);
}

template <Modint mint> mint sum_of_exp2_limit(mint r, int d) {
    linear_sieve sieve(d);
    auto f = sieve.pow_table<mint>(d, d);
    return sum_of_exp_times_poly_limit(f, r);
}

}  // namespace ebi

Back to top page