$\sum_{i}r^i poly(i)$
(fps/sum-of-exp-poly.hpp)

• View this file on GitHub
• Last update: 2025-10-25 18:30:13+09:00
• Include: #include "fps/sum-of-exp-poly.hpp"

多項式 $f$ について以下の値を計算する．

• $\sum_{i=0}^{n-1}r^if(i)$

• $\sum_{i=0}^{\infty}r^if(i),\quad

  r

  \lt 1$

  • 有理数として $

    r

    \lt 1$ を満たすものが modint で与えられる

$d=\deg f$ とし，$f(0),\dots,f(d)$ が与えられるものとする．

仕組み

$r=0$ は自明．

$r=1$ のとき単に Lagrange 補間すればよい．

以下 $r\neq 0,1$ とする．

数列 ${r^if(i)}$ は $\dfrac{p(x)}{(1-rx)^{d+1}}$ の形の母関数をもつ．

従って $a_i=\sum_{i=0}^{n-1}r^if(i)$ とすれば ${a_i}$ の母関数は \(\frac{xp(x)}{(1-rx)^{d+1}(1-x)}=\frac{q(x)}{(1-rx)^{d+1}}+\frac{c}{1-x}\) と部分分数分解できる． $q(x)$ は高々 $d$ 次であるから，両辺に $(1-rx)^{d+1}$ を掛けて $x^{d+1}$ の係数を考えれば次のように $c$ が求められる． \(c=(1-r)^{-(d+1)}\sum_{i=0}^{d+1}\binom{d+1}{i}(-r)^{d+1-i}a_i\)

また $\frac{q(x)}{(1-rx)^{d+1}}$ は ${a_i-c}$ の母関数であるが，この母関数から得られる数列は $r$ の冪と多項式の和の形であるから Lagrange 補間で求められる．

無限和のとき ${a_i-c}$ は $0$ に収束するので求める値は単に $c$ である．

Depends on

• Interpolate (fps/interpolate.hpp)
• 階乗, 二項係数 (modint/factorial.hpp)

Verified with

• verify/fps/LC_sum_of_exponential_times_polynomial.test.cpp
• verify/fps/LC_sum_of_exponential_times_polynomial_limit.test.cpp

Code

#pragma once
#include "fps/interpolate.hpp"
#include "modint/factorial.hpp"

// sum_{i=0}^{infty}r^i*poly(i)
// f[i]=poly(i)
template <class mint>
mint SumOfExpPolyLimit(mint r, vector<mint>& f) {
  if (r == 0) return f[0];
  assert(r != 1);
  int k = f.size();
  vector<mint> g(k + 1, 0);
  mint prod = 1;
  for (int i = 0; i < k; i++) {
    g[i + 1] = g[i] + f[i] * prod;
    prod *= r;
  }
  using fact = Factorial<mint>;
  mint c = 0;
  prod = 1;
  for (int i = 0; i <= k; i++) {
    c += fact::binom(k, i) * prod * g[k - i];
    prod *= -r;
  }
  c /= (1 - r).pow(k);
  return c;
}
// sum_{i=0}^{n-1}r^i*poly(i)
// f[i]=poly(i)
template <class mint>
mint SumOfExpPoly(long long n, mint r, vector<mint>& f) {
  if (n <= 0) return 0;
  if (r == 0) return f[0];
  int k = f.size();
  vector<mint> g(k + 1, 0);
  mint prod = 1;
  for (int i = 0; i < k; i++) {
    g[i + 1] = g[i] + f[i] * prod;
    prod *= r;
  }
  if (r == 1) return Interpolate(g, mint(n));
  mint c = 0;
  prod = 1;
  using fact = Factorial<mint>;
  for (int i = 0; i <= k; i++) {
    c += fact::binom(k, i) * prod * g[k - i];
    prod *= -r;
  }
  c /= (1 - r).pow(k);
  for (int i = 0; i <= k; i++) g[i] -= c;
  mint ir = r.inv();
  prod = 1;
  for (int i = 1; i <= k; i++) g[i] *= (prod *= ir);
  return Interpolate(g, mint(n)) * r.pow(n) + c;
}
/**
 * @brief $\sum_{i}r^i poly(i)$
 * @docs docs/fps/sum-of-exp-poly.md
 */

#line 2 "modint/factorial.hpp"

template <class mint>
struct Factorial {
  static void reserve(int n) {
    inv(n);
    fact(n);
    fact_inv(n);
  }
  static mint inv(int n) {
    static long long mod = mint::get_mod();
    static vector<mint> buf({0, 1});
    assert(n != 0);
    if (mod != mint::get_mod()) {
      mod = mint::get_mod();
      buf = vector<mint>({0, 1});
    }
    while ((int)buf.capacity() <= n) buf.reserve(buf.capacity() * 2);
    while ((int)buf.size() <= n) {
      long long k = buf.size(), q = (mod + k - 1) / k;
      buf.push_back(q * buf[k * q - mod]);
    }
    return buf[n];
  }
  static mint fact(int n) {
    static long long mod = mint::get_mod();
    static vector<mint> buf({1, 1});
    assert(n >= 0);
    if (mod != mint::get_mod()) {
      mod = mint::get_mod();
      buf = vector<mint>({1, 1});
    }
    while ((int)buf.capacity() <= n) buf.reserve(buf.capacity() * 2);
    while ((int)buf.size() <= n) {
      long long k = buf.size();
      buf.push_back(buf.back() * k);
    }
    return buf[n];
  }
  static mint fact_inv(int n) {
    static long long mod = mint::get_mod();
    static vector<mint> buf({1, 1});
    assert(n >= 0);
    if (mod != mint::get_mod()) {
      mod = mint::get_mod();
      buf = vector<mint>({1, 1});
    }
    if ((int)buf.size() <= n) inv(n);
    while ((int)buf.capacity() <= n) buf.reserve(buf.capacity() * 2);
    while ((int)buf.size() <= n) {
      long long k = buf.size();
      buf.push_back(buf.back() * inv(k));
    }
    return buf[n];
  }
  static mint binom(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact(n) * fact_inv(r) * fact_inv(n - r);
  }
  static mint binom_naive(int n, int r) {
    if (r < 0 || r > n) return 0;
    mint res = fact_inv(r);
    for (int i = 0; i < r; i++) res *= n - i;
    return res;
  }
  static mint multinom(const vector<int>& r) {
    int n = 0;
    for (auto& x : r) {
      if (x < 0) return 0;
      n += x;
    }
    mint res = fact(n);
    for (auto& x : r) res *= fact_inv(x);
    return res;
  }
  static mint P(int n, int r) {
    if (r < 0 || r > n) return 0;
    return fact(n) * fact_inv(n - r);
  }
  // partition n items to r groups (allow empty group)
  static mint H(int n, int r) {
    if (n < 0 || r < 0) return 0;
    return r == 0 ? 1 : binom(n + r - 1, r);
  }
};
/**
 * @brief 階乗, 二項係数
 */
#line 3 "fps/interpolate.hpp"

// f(0),f(1),...,f(n-1) -> f(x)
template <class mint>
mint Interpolate(const vector<mint>& f, mint x) {
  int n = f.size();
  vector<mint> l(n, 1), r(n, 1);
  for (int i = 0; i + 1 < n; i++) l[i + 1] = l[i] * (x - i);
  for (int i = n - 1; i > 0; i--) r[i - 1] = r[i] * (x - i);
  using fact = Factorial<mint>;
  mint s = 0;
  for (int i = 0; i < n; i++) {
    mint v = f[i] * l[i] * r[i] * fact::fact_inv(i) * fact::fact_inv(n - 1 - i);
    if ((n - i) & 1)
      s += v;
    else
      s -= v;
  }
  return s;
}
/**
 * @brief Interpolate
 * @docs docs/fps/interpolate.md
 */
#line 4 "fps/sum-of-exp-poly.hpp"

// sum_{i=0}^{infty}r^i*poly(i)
// f[i]=poly(i)
template <class mint>
mint SumOfExpPolyLimit(mint r, vector<mint>& f) {
  if (r == 0) return f[0];
  assert(r != 1);
  int k = f.size();
  vector<mint> g(k + 1, 0);
  mint prod = 1;
  for (int i = 0; i < k; i++) {
    g[i + 1] = g[i] + f[i] * prod;
    prod *= r;
  }
  using fact = Factorial<mint>;
  mint c = 0;
  prod = 1;
  for (int i = 0; i <= k; i++) {
    c += fact::binom(k, i) * prod * g[k - i];
    prod *= -r;
  }
  c /= (1 - r).pow(k);
  return c;
}
// sum_{i=0}^{n-1}r^i*poly(i)
// f[i]=poly(i)
template <class mint>
mint SumOfExpPoly(long long n, mint r, vector<mint>& f) {
  if (n <= 0) return 0;
  if (r == 0) return f[0];
  int k = f.size();
  vector<mint> g(k + 1, 0);
  mint prod = 1;
  for (int i = 0; i < k; i++) {
    g[i + 1] = g[i] + f[i] * prod;
    prod *= r;
  }
  if (r == 1) return Interpolate(g, mint(n));
  mint c = 0;
  prod = 1;
  using fact = Factorial<mint>;
  for (int i = 0; i <= k; i++) {
    c += fact::binom(k, i) * prod * g[k - i];
    prod *= -r;
  }
  c /= (1 - r).pow(k);
  for (int i = 0; i <= k; i++) g[i] -= c;
  mint ir = r.inv();
  prod = 1;
  for (int i = 1; i <= k; i++) g[i] *= (prod *= ir);
  return Interpolate(g, mint(n)) * r.pow(n) + c;
}
/**
 * @brief $\sum_{i}r^i poly(i)$
 * @docs docs/fps/sum-of-exp-poly.md
 */

Back to top page