This documentation is automatically generated by online-judge-tools/verification-helper
#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})$ で求める。
$K$ 次多項式 $f(n)$ の $f(0), f(1), \dots, f(k)$ と $a, n$ を与えて $\sum_{i = 0}^{n-1} a^i f(i)$ を求める。
$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$ である)。
$f(n) = n^k$ について求める。つまり、 $\sum_{i = 0}^{n-1} r^i i^d$ を求める。
$f(n) = n^k$ について求める。つまり、 $\sum_{i = 0}^{\infty} r^i i^d$ を求める。
#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