線形漸化的数列の第$k$項 (src/Math/bostan_mori.hpp)

View this file on GitHub
View document part on GitHub
Last update: 2024-10-14 11:50:44+09:00
Include: #include "src/Math/bostan_mori.hpp"

名前	概要	計算量
`div_at(p,q,k)`	$\displaystyle \lbrack x^k\rbrack\frac{p(x)}{q(x)}$ の値を返す. $q(x)$の定数項は $1$ か逆元が存在しないとだめ．	$O(d^2\log k)$ ただし $d = \max(\deg(p), \deg(q))$
`linear_recurrence(c,a,k)`	線形漸化式 $\displaystyle b_n = \begin{cases}a_n & n \lt d \newline c_0 b_{n-1} + c_1 b_{n-2} + … + c_{d-1} b_{n-d} & n \geq d \end{cases}$ を満たす数列 $(b_n)$ の第 $k$ 項 $b_k$ の値を返す.	$O(d^2\log k)$

参考

線形漸化的数列のN項目の計算

Verified with

Code

#pragma once
#include <vector>
#include <cassert>
#include <cstdint>
template <class K> K div_at(std::vector<K> p, std::vector<K> q, uint64_t k) {
 int n= p.size() - 1, m= q.size() - 1;
 for (assert(q[0] != K(0));; --n)
  if (n < 0 || p[n] != K()) break;
 for (;; --m)
  if (m < 0 || q[m] != K()) break;
 const int l= std::max(n, m) + 1;
 p.resize(l), q.resize(l);
 for (std::vector<K> np; k > (uint64_t)m; q.swap(p), p.swap(np), k>>= 1) {
  np.assign(l, K());
  if (k & 1) {
   for (int i= 0; i < l; i+= 2)
    for (int j= 1; j < l; j+= 2) np[(i + j) >> 1]+= p[j] * q[i] - p[i] * q[j];
  } else {
   for (int i= 0; i < l; i+= 2)
    for (int j= 0; j < l; j+= 2) np[(i + j) >> 1]+= p[i] * q[j];
   for (int i= 1; i < l; i+= 2)
    for (int j= 1; j < l; j+= 2) np[(i + j) >> 1]-= p[i] * q[j];
  }
  p.assign(l, K());
  for (int i= 0; i < l; i+= 2)
   for (int j= 0; j < i; j+= 2) p[(i + j) >> 1]+= q[i] * q[j];
  for (int i= 1; i < l; i+= 2)
   for (int j= 1; j < i; j+= 2) p[(i + j) >> 1]-= q[i] * q[j];
  for (int i= l; i--;) p[i]+= p[i];
  for (int i= 0; i < l; i+= 2) p[i]+= q[i] * q[i];
  for (int i= 1; i < l; i+= 2) p[i]-= q[i] * q[i];
 }
 K iv= K(1) / q[0];
 for (unsigned j= 0; j <= k; p[j++]*= iv)
  for (int i= j; i; --i) p[j]-= p[j - i] * q[i];
 return p[k];
}
// a[n] = c[0] * a[n-1] + c[1] * a[n-2] + ... + c[d-1] * a[n-d]
// return a[k]
template <class K> K linear_recurrence(std::vector<K> c, const std::vector<K> &a, uint64_t k) {
 if (k < a.size()) return a[k];
 const size_t d= c.size();
 assert(d <= a.size());
 for (auto &x: c) x= -x;
 std::vector<K> p(d);
 c.insert(c.begin(), K(1));
 for (int i= d; i--;)
  for (int j= i; j >= 0; --j) p[i]+= c[j] * a[i - j];
 return div_at<K>(p, c, k);
}

#line 2 "src/Math/bostan_mori.hpp"
#include <vector>
#include <cassert>
#include <cstdint>
template <class K> K div_at(std::vector<K> p, std::vector<K> q, uint64_t k) {
 int n= p.size() - 1, m= q.size() - 1;
 for (assert(q[0] != K(0));; --n)
  if (n < 0 || p[n] != K()) break;
 for (;; --m)
  if (m < 0 || q[m] != K()) break;
 const int l= std::max(n, m) + 1;
 p.resize(l), q.resize(l);
 for (std::vector<K> np; k > (uint64_t)m; q.swap(p), p.swap(np), k>>= 1) {
  np.assign(l, K());
  if (k & 1) {
   for (int i= 0; i < l; i+= 2)
    for (int j= 1; j < l; j+= 2) np[(i + j) >> 1]+= p[j] * q[i] - p[i] * q[j];
  } else {
   for (int i= 0; i < l; i+= 2)
    for (int j= 0; j < l; j+= 2) np[(i + j) >> 1]+= p[i] * q[j];
   for (int i= 1; i < l; i+= 2)
    for (int j= 1; j < l; j+= 2) np[(i + j) >> 1]-= p[i] * q[j];
  }
  p.assign(l, K());
  for (int i= 0; i < l; i+= 2)
   for (int j= 0; j < i; j+= 2) p[(i + j) >> 1]+= q[i] * q[j];
  for (int i= 1; i < l; i+= 2)
   for (int j= 1; j < i; j+= 2) p[(i + j) >> 1]-= q[i] * q[j];
  for (int i= l; i--;) p[i]+= p[i];
  for (int i= 0; i < l; i+= 2) p[i]+= q[i] * q[i];
  for (int i= 1; i < l; i+= 2) p[i]-= q[i] * q[i];
 }
 K iv= K(1) / q[0];
 for (unsigned j= 0; j <= k; p[j++]*= iv)
  for (int i= j; i; --i) p[j]-= p[j - i] * q[i];
 return p[k];
}
// a[n] = c[0] * a[n-1] + c[1] * a[n-2] + ... + c[d-1] * a[n-d]
// return a[k]
template <class K> K linear_recurrence(std::vector<K> c, const std::vector<K> &a, uint64_t k) {
 if (k < a.size()) return a[k];
 const size_t d= c.size();
 assert(d <= a.size());
 for (auto &x: c) x= -x;
 std::vector<K> p(d);
 c.insert(c.begin(), K(1));
 for (int i= d; i--;)
  for (int j= i; j >= 0; --j) p[i]+= c[j] * a[i - j];
 return div_at<K>(p, c, k);
}