cplib-cpp

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

View the Project on GitHub hitonanode/cplib-cpp

:heavy_check_mark: Sum of exponential times polynomial ($\sum_{i=0}^{N - 1} r^i f(i)$)
(formal_power_series/sum_of_exponential_times_polynomial.hpp)

一般に $f(x)$ を高々 $d$ 次の多項式として，$\displaystyle \sum_{i=0}^{N - 1} r^i f(i)$ を求める．先頭 $d + 1$ 項の値が既知であれば計算量は $O(d)$．関数は分かるが具体的な値が未知の場合は $O\left(d \left(\log d\right)^2\right)$ で先頭 $d + 1$ 項の値の multipoint evaluation を行う必要があると思われる．

導出の説明

min_25 さん によると，求めたい値は $C + g(x) r^x$, $g(x)$ はなんらかの多項式，という形になる．したがって，$r = 1$ の場合はそのまま Lagrange 補間を行えばよく，$r \neq 1$ の場合は $C$ や $r^x$ の分を打ち消した $g(x)$ に対して Lagrange 補間を行えばよい．$C$ の値は sum_of_sum_of_exponential_times_polynomial_limit() 関数によって $O(d)$ で求められる．

使用方法

先頭の $d + 1$ 項を計算して sum_of_exponential_times_polynomial() 関数に与える．

using mint = ModInt<998244353>;
vector<mint> xs(d + 1);
for (int i = 0; i <= d; i++) xs[i] = i;
vector<mint> f = some_function();
vector<mint> initial_terms = MultipointEvaluation<mint>(xs).evaluate_polynomial(f);

mint sum = sum_of_exponential_times_polynomial<mint>(r, initial_terms, 12345678910111213LL);

単項式の場合，先頭 $d + 1$ 項の列挙が線形篩により $O(d)$ で可能である．使用方法は以下の通り．

mint r;
int d;
vector<mint> initial_terms = Sieve(d).enumerate_kth_pows<mint>(d, d);
cout << sum_of_exponential_times_polynomial<mint>(r, initial_terms, 12345678910111213LL) << '\n';

リンク

Depends on

Required by

Verified with

Code

#pragma once
#include "lagrange_interpolation.hpp"
#include "sum_of_exponential_times_polynomial_limit.hpp"
#include <cassert>
#include <vector>

// CUT begin
// $d$ 次以下の多項式 $f(x)$ と定数 $r$ について，
// $\sum_{i=0}^{N-1} r^i f(i)$ の値を $[f(0), ..., f(d - 1), f(d)]$ の値から $O(d)$ で計算．
// Reference:
// - https://judge.yosupo.jp/problem/sum_of_exponential_times_polynomial
// - https://min-25.hatenablog.com/entry/2015/04/24/031413
template <typename MODINT>
MODINT sum_of_exponential_times_polynomial(MODINT r, const std::vector<MODINT> &init, long long N) {
    if (N <= 0) return 0;
    if (init.size() == 1) return init[0] * (r == 1 ? MODINT(N) : (1 - r.pow(N)) / (1 - r));

    std::vector<MODINT> S(init.size() + 1);
    MODINT rp = 1;
    for (int i = 0; i < int(init.size()); i++) {
        S[i + 1] = S[i] + init[i] * rp;
        rp *= r;
    }
    if (N < (long long)S.size()) return S[N];

    if (r == 1) return interpolate_iota<MODINT>(S, N);

    const MODINT Sinf = sum_of_exponential_times_polynomial_limit<MODINT>(r, init);
    MODINT rinv = r.inv(), rinvp = 1;
    for (int i = 0; i < int(S.size()); i++) {
        S[i] = (S[i] - Sinf) * rinvp;
        rinvp *= rinv;
    }
    return interpolate_iota<MODINT>(S, N) * r.pow(N) + Sinf;
};
#line 2 "formal_power_series/lagrange_interpolation.hpp"
#include <vector>

// CUT begin
// Lagrange interpolation
// Input: [f(0), ..., f(N-1)] (length = N), deg(f) < N
// Output: f(x_eval)
// Complexity: O(N)
// Verified: https://atcoder.jp/contests/arc033/tasks/arc033_4
template <typename MODINT> MODINT interpolate_iota(const std::vector<MODINT> ys, MODINT x_eval) {
    const int N = ys.size();
    if (x_eval.val() < N) return ys[x_eval.val()];
    std::vector<MODINT> facinv(N);
    facinv[N - 1] = MODINT(N - 1).fac().inv();
    for (int i = N - 1; i > 0; i--) facinv[i - 1] = facinv[i] * i;
    std::vector<MODINT> numleft(N);
    MODINT numtmp = 1;
    for (int i = 0; i < N; i++) {
        numleft[i] = numtmp;
        numtmp *= x_eval - i;
    }
    numtmp = 1;
    MODINT ret = 0;
    for (int i = N - 1; i >= 0; i--) {
        MODINT tmp = ys[i] * numleft[i] * numtmp * facinv[i] * facinv[N - 1 - i];
        ret += ((N - 1 - i) & 1) ? (-tmp) : tmp;
        numtmp *= x_eval - i;
    }
    return ret;
}
#line 2 "formal_power_series/sum_of_exponential_times_polynomial_limit.hpp"
#include <cassert>
#line 4 "formal_power_series/sum_of_exponential_times_polynomial_limit.hpp"

// CUT begin
// $d$ 次以下の多項式 $f(x)$ と定数 $r$ について，
// $\sum_{i=0}^\infty r^i f(i)$ の値を $[f(0), ..., f(d - 1), f(d)]$ の値から $O(d)$ で計算．
// Requirement: r != 1
// https://judge.yosupo.jp/problem/sum_of_exponential_times_polynomial_limit
// Document: https://hitonanode.github.io/cplib-cpp/formal_power_series/sum_of_exponential_times_polynomial_limit.hpp
template <typename MODINT>
MODINT sum_of_exponential_times_polynomial_limit(MODINT r, std::vector<MODINT> init) {
    assert(r != 1);
    if (init.empty()) return 0;
    if (init.size() == 1) return init[0] / (1 - r);
    auto &bs = init;
    const int d = int(bs.size()) - 1;
    MODINT rp = 1;
    for (int i = 1; i <= d; i++) rp *= r, bs[i] = bs[i] * rp + bs[i - 1];
    MODINT ret = 0;
    rp = 1;
    for (int i = 0; i <= d; i++) {
        ret += bs[d - i] * MODINT(d + 1).nCr(i) * rp;
        rp *= -r;
    }
    return ret / MODINT(1 - r).pow(d + 1);
};
#line 6 "formal_power_series/sum_of_exponential_times_polynomial.hpp"

// CUT begin
// $d$ 次以下の多項式 $f(x)$ と定数 $r$ について，
// $\sum_{i=0}^{N-1} r^i f(i)$ の値を $[f(0), ..., f(d - 1), f(d)]$ の値から $O(d)$ で計算．
// Reference:
// - https://judge.yosupo.jp/problem/sum_of_exponential_times_polynomial
// - https://min-25.hatenablog.com/entry/2015/04/24/031413
template <typename MODINT>
MODINT sum_of_exponential_times_polynomial(MODINT r, const std::vector<MODINT> &init, long long N) {
    if (N <= 0) return 0;
    if (init.size() == 1) return init[0] * (r == 1 ? MODINT(N) : (1 - r.pow(N)) / (1 - r));

    std::vector<MODINT> S(init.size() + 1);
    MODINT rp = 1;
    for (int i = 0; i < int(init.size()); i++) {
        S[i + 1] = S[i] + init[i] * rp;
        rp *= r;
    }
    if (N < (long long)S.size()) return S[N];

    if (r == 1) return interpolate_iota<MODINT>(S, N);

    const MODINT Sinf = sum_of_exponential_times_polynomial_limit<MODINT>(r, init);
    MODINT rinv = r.inv(), rinvp = 1;
    for (int i = 0; i < int(S.size()); i++) {
        S[i] = (S[i] - Sinf) * rinvp;
        rinvp *= rinv;
    }
    return interpolate_iota<MODINT>(S, N) * r.pow(N) + Sinf;
};
Back to top page