線形漸化式を持つ数列の第 $n$ 項目の計算
(math/fps/NthTermFind.hpp)

• View this file on GitHub
• Last update: 2026-05-19 17:02:02+09:00
• Include: #include "math/fps/NthTermFind.hpp"

概要

todo

計算量

todo

Depends on

• Matrix (linear-algebra/Matrix.hpp)
• Bostan-Mori Algorithm (math/fps/BostanMori.hpp)
• Formal Power Series (math/fps/FormalPowerSeries.hpp)

Code

#include "../../linear-algebra/Matrix.hpp"
#include "BostanMori.hpp"

// a.size() == d*2
template <class S>
S nth_term_find(const vector<S>& a, long long n) {
  int d = int(a.size()) / 2;
  Matrix<S> A(d);
  for (int i = 0; i < d; i++) {
    for (int j = 0; j < d; j++) A[i][j] = a[d - j - 1 + i];
  }
  vector<S> coef = A.linear_equation({a.begin() + d, a.end()})[0];
  FPS<S> l = {a.begin(), a.begin() + d};
  FPS<S> q(d + 1);
  q[0] = 1;
  for (int i = 0; i < d; i++) q[i + 1] = -coef[i];
  l *= q;
  return bostan_mori<S>(l, q, n);
}

// a.size() == d*2
template <class S>
S nth_term_find_2(const vector<S>& a, long long n) {
  int d = int(a.size()) / 2;
  Matrix<S> A(d);
  for (int i = 0; i < d; i++) {
    for (int j = 0; j < d; j++) A[i][j] = a[d - j - 1 + i];
  }
  vector<S> coef = A.linear_equation({a.begin() + d, a.end()})[0];
  FPS<S> l = {a.begin(), a.begin() + d};
  FPS<S> q(d + 1);
  q[0] = 1;
  for (int i = 0; i < d; i++) q[i + 1] = -coef[i];
  l = multiply(l, q, d);
  return bostan_mori_naive<S>(l, q, n);
}

#line 1 "linear-algebra/Matrix.hpp"
template <class S>
struct Matrix {
 private:
 public:
  vector<vector<S>> A;
  Matrix() {}
  Matrix(int n, int m) : A(n, vector<S>(m)) {}
  Matrix(int n) : A(n, vector<S>(n)) {}

  inline int size() const { return A.size(); }
  inline int height() const { return A.size(); }
  inline int width() const { return A[0].size(); }
  inline const vector<S>& operator[](int h) const { return (A[h]); }
  inline vector<S>& operator[](int h) { return (A[h]); }

  Matrix& operator+=(const Matrix& B) {
    int h = height();
    int w = width();
    for (int i = 0; i < h; i++) {
      for (int j = 0; j < w; j++) {
        (*this)[i][j] += B[i][j];
      }
    }
    return (*this);
  }
  Matrix& operator-=(const Matrix& B) {
    int h = height();
    int w = width();
    for (int i = 0; i < h; i++) {
      for (int j = 0; j < w; j++) {
        (*this)[i][j] -= B[i][j];
      }
    }
    return (*this);
  }

  Matrix& operator*=(const Matrix& B) {
    int h = height();
    int w = B.width();
    int c = width();
    vector<vector<S>> C(h, vector<S>(w));
    for (int i = 0; i < h; i++) {
      for (int j = 0; j < w; j++) {
        for (int k = 0; k < c; k++) {
          C[i][j] = (C[i][j] + (*this)[i][k] * B[k][j]);
        }
      }
    }
    A = move(C);
    return (*this);
  }
  Matrix operator+(const Matrix& B) const { return (Matrix(*this) += B); }
  Matrix operator-(const Matrix& B) const { return (Matrix(*this) -= B); }
  Matrix operator*(const Matrix& B) const { return (Matrix(*this) *= B); }

  int rank() {
    Matrix B(*this);
    if (B.height() == 0 or B.width() == 0) return 0;
    int res = 0;
    int h = height();
    int w = width();
    int ch = 0;
    int cw = 0;
    while (ch < h and cw < w) {
      bool ok = false;
      for (int j = cw; j < w; j++) {
        for (int i = ch; i < h; i++) {
          if (B[i][j] != 0) {
            ok = true;
            swap(B[ch], B[i]);
            S d = B[ch][j];
            for (int j2 = j; j2 < w; j2++) {
              B[ch][j2] /= d;
            }
            for (int i2 = 0; i2 < h; i2++) {
              if (B[i2][j] != 0 and i2 != ch) {
                S m = B[i2][j];
                for (int j2 = j; j2 < w; j2++) {
                  B[i2][j2] -= B[ch][j2] * m;
                }
              }
            }
            res++;
            ch++;
            cw = j + 1;
            break;
          }
        }
        if (ok) break;
      }
      if (!ok) break;
    }
    return res;
  }

  S determinant() {
    Matrix B(*this);
    if (B.height() == 0 or B.width() == 0) return 0;
    assert(B.height() == B.width());
    int h = height();
    int w = width();
    int ch = 0;
    int cw = 0;
    S div = 1;
    bool neg = false;
    while (ch < h and cw < w) {
      bool ok = false;
      for (int j = cw; j < w; j++) {
        for (int i = ch; i < h; i++) {
          if (B[i][j] != 0) {
            ok = true;
            if (ch != i) neg = !neg;
            swap(B[ch], B[i]);
            S d = B[ch][j];
            div /= d;
            for (int j2 = j; j2 < w; j2++) {
              B[ch][j2] /= d;
            }
            for (int i2 = 0; i2 < h; i2++) {
              if (B[i2][j] != 0 and i2 != ch) {
                S m = B[i2][j];
                for (int j2 = j; j2 < w; j2++) {
                  B[i2][j2] -= B[ch][j2] * m;
                }
              }
            }
            ch++;
            cw = j + 1;
            break;
          }
        }
        if (ok) {
          break;
        } else {
          return S(0);
        }
      }
      if (!ok) break;
    }
    S res = (neg ? -B[0][0] : B[0][0]) / div;
    for (int i = 1; i < h; i++) {
      res = res * B[i][i];
    }
    return res;
  }

  pair<bool, Matrix<S>> inverse() {
    int h = height();
    int w = width();
    assert(h == w);
    Matrix<S> B(h, w * 2);
    for (int i = 0; i < h; i++) {
      for (int j = 0; j < w; j++) {
        B[i][j] = (*this)[i][j];
      }
    }
    for (int i = 0; i < h; i++) {
      B[i][i + w] = 1;
    }
    w *= 2;
    int rnk = 0;

    int ch = 0;
    int cw = 0;
    while (ch < h and cw < h) {
      bool ok = false;
      for (int j = cw; j < h; j++) {
        for (int i = ch; i < h; i++) {
          if (B[i][j] != 0) {
            ok = true;
            swap(B[ch], B[i]);
            S d = B[ch][j];
            for (int j2 = j; j2 < w; j2++) {
              B[ch][j2] /= d;
            }
            for (int i2 = 0; i2 < h; i2++) {
              if (B[i2][j] != 0 and i2 != ch) {
                S m = B[i2][j];
                for (int j2 = j; j2 < w; j2++) {
                  B[i2][j2] -= B[ch][j2] * m;
                }
              }
            }
            rnk++;
            ch++;
            cw = j + 1;
            break;
          }
        }
        if (ok) break;
      }
      if (!ok) break;
    }
    Matrix<S> res(h);
    if (rnk == h) {
      for (int i = 0; i < h; i++) {
        for (int j = 0; j < h; j++) {
          res[i][j] = B[i][j + h];
        }
      }
      return {true, res};
    } else {
      return {false, res};
    }
  }

  Matrix<S> linear_equation(vector<S> b) {
    Matrix A(*this);

    int rnk = 0;
    assert(A.height() == b.size());
    int h = height();
    int w = width();
    int ch = 0;
    int cw = 0;
    vector<int> pivot_row(w, -1);
    while (ch < h and cw < w) {
      bool ok = false;
      for (int j = cw; j < w; j++) {
        for (int i = ch; i < h; i++) {
          if (A[i][j] != 0) {
            ok = true;
            swap(A[ch], A[i]);
            swap(b[ch], b[i]);
            S d = A[ch][j];
            for (int j2 = j; j2 < w; j2++) {
              A[ch][j2] /= d;
            }
            b[ch] /= d;
            for (int i2 = 0; i2 < h; i2++) {
              S m = A[i2][j];
              if (A[i2][j] != 0 and i2 != ch) {
                for (int j2 = j; j2 < w; j2++) {
                  A[i2][j2] -= A[ch][j2] * m;
                }
              }
              if (i2 != ch) b[i2] -= b[ch] * m;
            }
            pivot_row[j] = ch;
            rnk++;
            ch++;
            cw = j + 1;
            break;
          }
        }
        if (ok) break;
      }
      if (!ok) break;
    }

    for (int i = rnk; i < h; i++) {
      if (b[i] != 0) return Matrix<S>(0);
    }
    Matrix<S> sol(w - rnk + 1, w);
    int idx = 1;
    for (int j = 0; j < w; j++) {
      if (pivot_row[j] != -1) {
        sol[0][j] = b[pivot_row[j]];
      } else {
        sol[idx][j] = 1;
        for (int i = 0; i < w; i++) {
          if (pivot_row[i] != -1) {
            sol[idx][i] = -A[pivot_row[i]][j];
          }
        }
        idx++;
      }
    }
    return sol;
  }
};
#line 1 "math/fps/FormalPowerSeries.hpp"
template <class S>
struct FPS;

template <class S>
struct SFPS : vector<pair<int, S>> {
  using vector<pair<int, S>>::vector;
  using vector<pair<int, S>>::operator=;

  FPS<S> log(int deg);
  FPS<S> exp(int deg);
  FPS<S> pow(long long m, int deg);
};

template <class S>
struct FPS : vector<S> {
  using vector<S>::vector;
  using vector<S>::operator=;

  FPS<S> inv() const {
    int n = int((*this).size());
    FPS<S> res = {(*this)[0].inv()};
    while (int(res.size()) < n) {
      int m = int(res.size());
      // f = f[0, 2m)
      FPS<S> f((*this).begin(), (*this).begin() + min(n, m << 1));
      FPS<S> inv_f(res);
      f.resize(m << 1);
      internal::butterfly(f);
      inv_f.resize(m << 1);
      internal::butterfly(inv_f);
      // f = f*g
      for (int i = 0; i < m << 1; i++) f[i] *= inv_f[i];
      internal::butterfly_inv(f);

      // f = f[m, 2m)
      f.erase(f.begin(), f.begin() + m);
      f.resize(m << 1);
      // f = f*g
      internal::butterfly(f);
      for (int i = 0; i < m << 1; i++) f[i] *= inv_f[i];
      internal::butterfly_inv(f);
      S m2i = S(m << 1).inv();
      m2i *= -m2i;
      for (int i = 0; i < m << 1; i++) f[i] *= m2i;
      res.insert(res.end(), f.begin(), f.begin() + m);
    }
    return {res.begin(), res.begin() + n};
  }

  FPS<S> exp() const {
    int n = int((*this).size());
    FPS<S> res = {1};
    assert((*this)[0] == 0);
    for (int siz = 1; siz < n; siz <<= 1) {
      FPS<S> f(siz << 1);
      f[0] = 1;
      res.resize(siz << 1);
      FPS<S> lg = res.log();
      for (int i = 0; i < siz << 1; i++) f[i] -= lg[i];
      for (int i = 0; i < min(siz << 1, n); i++) f[i] += (*this)[i];
      res *= f;
    }
    return {res.begin(), res.begin() + n};
  }

  FPS<S> log() const {
    FPS<S> res = *this;
    res.log_inplace();
    return res;
  }

  FPS<S> pow(__int128_t m) const {
    __int128_t n = int((*this).size());
    if (m == 0) {
      FPS<S> res(n);
      if (n) res[0] = 1;
      return res;
    }
    // 定数項を1にする
    for (__int128_t i = 0; i < n; i++) {
      if ((*this)[i] != 0) {
        S coef = (*this)[i];
        FPS<S> f((*this).begin() + i, (*this).end());
        if (coef != 1) {
          for (int j = 0; j < n - i; j++) f[j] /= coef;
        }
        f.log_inplace();
        for (int j = 0; j < n - i; j++) f[j] *= m;
        f.exp_inplace();
        coef = coef.pow(m);
        for (int j = 0; j < n - i; j++) f[j] *= coef;
        FPS<S> res(min(__int128_t(m) * i, n), 0);
        if (res.size() < n) res.insert(res.end(), f.begin(), f.begin() + min(__int128_t(n), n - res.size()));
        return res;
      }
      if (__int128_t(i + 1) * m >= n) return FPS<S>(n, 0);
    }
    return FPS<S>(n, 0);
  }

  FPS<S> differentiate() const {
    int n = int((*this).size());
    FPS<S> res(n);
    for (int i = 0; i < n - 1; i++) res[i] = (*this)[i + 1] * S(i + 1);
    res[n - 1] = 0;
    return res;
  }

  FPS<S> integrate() const {
    int n = int((*this).size());
    vector<S> iv(n);
    iv[1] = 1;
    for (int i = 2; i < n; i++) iv[i] = iv[S::mod() % i] * (-(S::mod() / i));
    FPS<S> res(n);
    res[0] = 0;
    for (int i = 0; i < n - 1; i++) res[i + 1] = (*this)[i] * iv[i + 1];
    return res;
  }

  void integrate_inplace() {
    int n = int((*this).size());
    static vector<S> inv{0, 1};
    if (int(inv.size()) < n) {
      int old_siz = inv.size();
      inv.resize(n);
      int mod = S::mod();
      for (int i = old_siz; i < n; i++) inv[i] = -inv[mod % i] * (mod / i);
    }
    for (int i = n - 1; i > 0; i--) (*this)[i] = (*this)[i - 1] * inv[i];
    (*this)[0] = 0;
  }

  void differentiate_inplace() {
    int n = int((*this).size());
    if (n == 0) return;
    for (int i = 0; i < n - 1; i++) {
      (*this)[i] = (*this)[i + 1] * S(i + 1);
    }
    (*this)[n - 1] = 0;
  }

  void inv_inplace() {
    *this = this->inv();
  }
  void exp_inplace() {
    *this = this->exp();
  }

  void log_inplace() {
    assert(!this->empty() and (*this)[0] == 1);
    FPS<S> inv_f = this->inv();
    this->differentiate_inplace();
    *this *= inv_f;
    this->integrate_inplace();
  }

  void pow_inplace(__int128_t m) {
    *this = this->pow(m);
  }

  FPS<S>& operator*=(const FPS<S>& g) {
    int n = int((*this).size());
    *this = convolution(*this, g);
    (*this).resize(n);
    return *this;
  }

  FPS<S>& operator/=(FPS<S>& g) {
    int n = int((*this).size());
    *this = convolution(*this, g.inv());
    (*this).resize(n);
    return *this;
  }

  FPS<S>& operator<<=(int k) {
    int n = int((*this).size());
    if (k >= n) {
      (*this).assign(n, 0);
      return *this;
    }
    for (int i = n - 1; i >= k; i--) (*this)[i] = (*this)[i - k];
    for (int i = 0; i < k; i++) (*this)[i] = 0;
    return *this;
  }

  FPS<S>& operator>>=(int k) {
    int n = int((*this).size());
    if (k >= n) {
      (*this).assign(n, 0);
      return *this;
    }
    for (int i = 0; i < n - k; i++) (*this)[i] = (*this)[i + k];
    for (int i = n - k; i < n; i++) (*this)[i] = 0;
    return *this;
  }

  FPS<S>& operator*=(const SFPS<S>& g) {
    int n = (*this).size();
    int k = int(g.size());
    auto [d, c] = g.front();
    int start = 0;
    if (d == 0) {
      start = 1;
    } else {
      c = 0;
    }
    for (int i = n - 1; i >= 0; i--) {
      (*this)[i] *= c;
      for (int j = start; j < k; j++) {
        const auto& [d_, c_] = g[j];
        if (d_ > i) break;
        (*this)[i] += (*this)[i - d_] * c_;
      }
    }
    return *this;
  }

  FPS<S>& operator/=(const SFPS<S>& g) {
    int n = (*this).size();
    int k = int(g.size());
    auto [d, c] = g.front();
    assert(d == 0 and c != 0);
    S inv = c.inv();
    for (int i = 0; i < n; i++) {
      for (int j = 1; j < k; j++) {
        const auto& [d, c] = g[j];
        if (d > i) break;
        (*this)[i] -= (*this)[i - d] * c;
      }
      (*this)[i] *= inv;
    }
    return *this;
  }

  FPS<S> operator<<(int k) const { return FPS<S>(*this) <<= k; }
  FPS<S> operator>>(int k) const { return FPS<S>(*this) >>= k; }
};

template <class S>
FPS<S> SFPS<S>::log(int deg) {
  FPS<S> f(deg);
  assert((*this)[0].first == 0 and (*this)[0].second == 1 and (*this).back().first < deg);
  int k = (*this).size();
  for (int i = 0; i < k; i++) {
    const auto& [d, c] = (*this)[i];
    f[d] = c;
  }
  f.differentiate_inplace();
  f /= (*this);
  f.integrate_inplace();
  return f;
}

template <class S>
FPS<S> SFPS<S>::exp(int deg) {
  SFPS df = (*this);
  int k = (*this).size();
  vector<S> inv(deg);
  inv[1] = 1;
  for (int i = 2; i < deg; i++) inv[i] = inv[S::mod() % i] * (-(S::mod() / i));

  // df = f'
  for (int i = 0; i < k; i++) {
    const auto& [d, c] = df[i];
    df[i] = {d - 1, d * c};
  }

  // F = exp(f)
  // F' = f'F
  FPS<S> F(deg);
  F[0] = 1;
  for (int i = 0; i < deg - 1; i++) {
    S conv_sum = 0;
    for (int j = 0; j < k; j++) {
      const auto& [d, c] = df[j];
      if (d > i) break;
      conv_sum += c * F[i - d];
    }
    F[i + 1] = conv_sum * inv[i + 1];
  }
  return F;
}

template <class S>
FPS<S> SFPS<S>::pow(long long m, int deg) {
  if (m == 0) {
    FPS<S> res(deg);
    if (deg) res[0] = 1;
    return res;
  }
  vector<S> inv(deg);
  inv[1] = 1;
  for (int i = 2; i < deg; i++) inv[i] = inv[S::mod() % i] * (-(S::mod() / i));

  int k = (*this).size();
  // F = f ^ m
  FPS<S> F(deg);
  // F' = m(f^(n-1))f'
  // fF' = mFf'

  // 定数項を1にする
  for (int i = 0; i < k; i++) {
    const auto& [d, c] = (*this)[i];
    if (c != 0) {
      SFPS f((*this).begin() + i, (*this).end());
      for (int j = 0; j < f.size(); j++) {
        f[j].first -= d;
        f[j].second /= c;
      }

      FPS<S> F(deg);
      F[0] = 1;
      for (int j = 0; j < deg - 1; j++) {
        S dF_j = 0;
        for (int l = 0; l < f.size(); l++) {
          const auto& [d_, c_] = f[l];
          if (d_ == 0) continue;
          if (d_ - 1 > j) break;
          dF_j += c_ * F[j - d_ + 1] * (S(m) * d_ - (j - d_ + 1));
        }
        F[j + 1] = dF_j * inv[j + 1];
      }
      S coef_pw = S(c).pow(m);
      for (int j = 0; j < deg; j++) F[j] *= coef_pw;

      FPS<S> res(min(__int128_t(m) * d, __int128_t(deg)), 0);
      if (res.size() < deg) res.insert(res.end(), F.begin(), F.begin() + min(deg, deg - int(res.size())));
      return res;
    }

    if (__int128_t(d + 1) * m >= deg) return FPS<S>(deg, 0);
  }

  return FPS<S>(deg, 0);
}

template <class S>
FPS<S> multiply(const FPS<S>& a, const FPS<S>& b, int d = -1) {
  int siz = int(a.size()) + int(b.size()) - 1;
  FPS<S> c(siz);
  for (int i = 0; i < int(a.size()); i++) {
    for (int j = 0; j < int(b.size()); j++) {
      if (d != -1 and i + j >= d) break;
      c[i + j] += a[i] * b[j];
    }
  }
  if (d != -1) c.resize(d);
  return c;
}
#line 2 "math/fps/BostanMori.hpp"

template <class S>
S bostan_mori(FPS<S> p, FPS<S> q, long long n) {
  auto filter = [&](const FPS<S>& p, int start) {
    FPS<S> ret((p.size() + 1) / 2);
    for (int i = 0; i * 2 + start < int(p.size()); i++) ret[i] = p[i * 2 + start];
    return ret;
  };

  while (n > 0) {
    FPS<S> q_ = q;
    for (int i = 1; i < int(q_.size()); i += 2) q_[i] = -q_[i];
    auto pq = convolution(p, q_);
    auto qq = convolution(q, q_);
    p = filter(FPS<S>(pq.begin(), pq.end()), n & 1);
    q = filter(FPS<S>(qq.begin(), qq.end()), 0);
    n >>= 1;
  }
  return p[0] / q[0];
}

template <class S>
S bostan_mori_naive(FPS<S> p, FPS<S> q, long long n) {
  auto filter = [&](const FPS<S>& p, int start) {
    FPS<S> ret((p.size() + 1) / 2);
    for (int i = 0; i * 2 + start < int(p.size()); i++) ret[i] = p[i * 2 + start];
    return ret;
  };
  while (n > 0) {
    FPS<S> q_ = q;
    for (int i = 1; i < int(q_.size()); i += 2) q_[i] = -q_[i];
    p = filter(multiply(p, q_), n & 1);
    q = filter(multiply(q, q_), 0);
    n >>= 1;
  }
  return p[0] / q[0];
}
#line 3 "math/fps/NthTermFind.hpp"

// a.size() == d*2
template <class S>
S nth_term_find(const vector<S>& a, long long n) {
  int d = int(a.size()) / 2;
  Matrix<S> A(d);
  for (int i = 0; i < d; i++) {
    for (int j = 0; j < d; j++) A[i][j] = a[d - j - 1 + i];
  }
  vector<S> coef = A.linear_equation({a.begin() + d, a.end()})[0];
  FPS<S> l = {a.begin(), a.begin() + d};
  FPS<S> q(d + 1);
  q[0] = 1;
  for (int i = 0; i < d; i++) q[i + 1] = -coef[i];
  l *= q;
  return bostan_mori<S>(l, q, n);
}

// a.size() == d*2
template <class S>
S nth_term_find_2(const vector<S>& a, long long n) {
  int d = int(a.size()) / 2;
  Matrix<S> A(d);
  for (int i = 0; i < d; i++) {
    for (int j = 0; j < d; j++) A[i][j] = a[d - j - 1 + i];
  }
  vector<S> coef = A.linear_equation({a.begin() + d, a.end()})[0];
  FPS<S> l = {a.begin(), a.begin() + d};
  FPS<S> q(d + 1);
  q[0] = 1;
  for (int i = 0; i < d; i++) q[i + 1] = -coef[i];
  l = multiply(l, q, d);
  return bostan_mori_naive<S>(l, q, n);
}

Back to top page