This documentation is automatically generated by online-judge-tools/verification-helper
#include "algorithm/Math/Algebra/kitamasa.hpp"
初めの $k$ 項 $\lbrace a_0, a_1, \ldots, a_{k-1} \rbrace$ と $k$ 階線形漸化式
\[a_n = d_0 \cdot a_{n-k} + d_1 \cdot a_{n-(k-1)} + \cdots + d_{k-1} \cdot a_{n-1}\]によって定まる数列の任意の項 $a_N$ を $\mathcal{O}(k^2 \log N)$ で求める.
説明のため,
\[a_n = D(n,0) \cdot a_0 + D(n,1) \cdot a_1 + \cdots + D(n,k-1) \cdot a_{k-1} = \sum_{i=0}^{k-1} D(n,i) \cdot a_{i}\]となる $D(n,\ast)$ を定義する(ただし,$n \geq 0$ ).
このとき,
\[\begin{equation} D(0,i) = \begin{cases} 1 &\text{if \ $i=0$ ,} \\ 0 &\text{otherwise} \end{cases} \notag \end{equation}\]である.
きたまさ法の方針は,$D(n,\ast) \rightarrow D(n+1,\ast)$ と $D(n,\ast) \rightarrow D(2n,\ast)$ を計算することにより,繰り返し二乗法と同じ要領で $D(0,\ast)$ から $D(N,\ast)$ を求めることである.
まず,前提として $D(n,\ast)$ が分かっていると仮定する. このとき,
\[\begin{align} a_{n+1} &= D(n,0) \cdot a_1 + D(n,1) \cdot a_2 + \cdots + D(n,k-2) \cdot a_{k-1} + D(n,k-1) \cdot a_k \notag\\ &= D(n,0) \cdot a_1 + D(n,1) \cdot a_2 + \cdots + D(n,k-2) \cdot a_{k-1} + D(n,k-1) \cdot (d_0 \cdot a_0 + d_1 \cdot a_1 + \cdots + d_{k-1} \cdot a_{k-1}) \notag\\ &= D(n,k-1) \cdot d_0 \cdot a_0 + (D(n,0) + D(n,k-1) \cdot d_1) \cdot a_1 + (D(n,1) + D(n,k-1) \cdot d_2) \cdot a_2 + \cdots + (D(n,k-2) + D(n,k-1) \cdot d_{k-1}) \cdot a_{k-1} \notag\\ \end{align}\]となるから,
\[\begin{equation} D(n+1,i) = \begin{cases} D(n,k-1) \cdot d_0 &\text{if \ $i=0$ ,} \\ D(n,i-1) + D(n,k-1) \cdot d_i &\text{otherwise} \end{cases} \notag \end{equation}\]と求まる. これは $\mathcal{O}(k)$ で計算できる.
先の方法で $D(n,\ast), D(n+1,\ast), \ldots, D(n+k-1,\ast)$ が $\mathcal{O}(k^2)$ で列挙できるとする. このとき,
\[\begin{align} a_{2n} &= D(n,0) \cdot a_n + D(n,1) \cdot a_{n+1} + \cdots + D(n,k-1) \cdot a_{n+k-1} \notag\\ &= D(n,0) \cdot (D(n,0) \cdot a_0 + D(n,1) \cdot a_1 + \cdots + D(n,k-1) \cdot a_{k-1}) \notag\\ &\quad + D(n,1) \cdot (D(n+1,0) \cdot a_0 + D(n+1,1) \cdot a_1 + \cdots + D(n+1,k-1) \cdot a_{k-1}) \notag\\ &\quad + \cdots \notag\\ &\quad + D(n,k-1) \cdot (D(n+k-1,0) \cdot a_0 + D(n+k-1,1) \cdot a_1 + \cdots + D(n+k-1,k-1) \cdot a_{k-1}) \notag\\ &= (D(n,0) \cdot D(n,0) + D(n,1) \cdot D(n+1,0) + \cdots + D(n,k-1) \cdot D(n+k-1,0)) \cdot a_0 \notag\\ &\quad + (D(n,0) \cdot D(n,1) + D(n,1) \cdot D(n+1,1) + \cdots + D(n,k-1) \cdot D(n+k-1,1)) \cdot a_1 \notag\\ &\quad + \cdots \notag\\ &\quad + (D(n,0) \cdot D(n,k-1) + D(n,1) \cdot D(n+1,k-1) + \cdots + D(n,k-1) \cdot D(n+k-1,k-1)) \cdot a_{k-1} \notag\\ \end{align}\]となるから,
\[D(2n,i) = \sum_{j=0}^{k-1} D(n,j) \cdot D(n+j,i)\]と求まる. これは $\mathcal{O}(k^2)$ で計算できる.
よって,繰り返し二乗法と同じ要領で $D(0,\ast)$ から $D(N,\ast)$ を $\mathcal{O}(k^2 \log N)$ で求めることができる.
#ifndef ALGORITHM_KITAMASA_HPP
#define ALGORITHM_KITAMASA_HPP 1
#include <cassert>
#include <iterator>
#include <utility>
#include <vector>
namespace algorithm {
// きたまさ法.
template <typename T>
class Kitamasa {
public:
using value_type = T;
using size_type = std::size_t;
private:
size_type m_k; // m_k:=(階数).
std::vector<value_type> m_a; // m_a[]:=(初項数列).
std::vector<value_type> m_d; // m_d[]:=(係数数列).
// f(n)->f(n+1). O(K).
std::vector<value_type> add(const std::vector<value_type> &x) const {
std::vector<value_type> y(m_k);
y[0] = x[m_k - 1] * m_d[0];
for(size_type i = 1; i < m_k; ++i) y[i] = x[i - 1] + x[m_k - 1] * m_d[i];
return y;
}
// f(n)->f(2*n). O(K^2).
std::vector<value_type> mul(const std::vector<value_type> &x) const {
std::vector<value_type> y(m_k, 0);
std::vector<value_type> t = x;
for(size_type i = 0; i < m_k; ++i) {
for(size_type j = 0; j < m_k; ++j) y[j] += x[i] * t[j];
if(i == m_k - 1) break;
t = add(t);
}
return y;
}
// f(n)を返す.O((K^2)*logN).
std::vector<value_type> f(size_type n) const {
if(n == 0) {
std::vector<value_type> x(m_k, 0);
x[0] = 1;
return x; // f(0).
}
std::vector<value_type> &&x = mul(f(n >> 1));
if(n & 1ULL) x = add(x);
return x;
}
public:
Kitamasa() : Kitamasa({0, 1}, {1, 1}) {} // フィボナッチ数列.
template <std::input_iterator InputIter>
explicit Kitamasa(InputIter a_first, InputIter a_last, InputIter d_first, InputIter d_last) : m_a(a_first, a_last), m_d(d_first, d_last) {
assert(m_a.size() >= 1);
assert(m_a.size() == m_d.size());
m_k = m_a.size();
}
explicit Kitamasa(std::vector<value_type> a, std::vector<value_type> d) : m_k(a.size()), m_a(std::move(a)), m_d(std::move(d)) {
assert(m_a.size() >= 1);
assert(m_a.size() == m_d.size());
}
value_type operator[](size_type n) const { return calc(n); }
// a[n]を求める.O((K^2)*log N).
value_type calc(size_type n) const {
std::vector<value_type> &&x = f(n);
value_type res = 0;
for(size_type i = 0; i < m_k; ++i) res += m_a[i] * x[i];
return res;
}
// a[l:r]を求める.O((K^2)*log(N) + K*(r-l)).
std::vector<value_type> calc(size_type l, size_type r) const {
assert(l < r);
std::vector<value_type> res(r - l);
std::vector<value_type> &&x = f(l);
for(size_type i = l; i < r; ++i) {
for(size_type j = 0; j < m_k; ++j) res[i - l] += m_a[j] * x[j];
if(i == r - 1) break;
x = add(x);
}
return res;
}
};
} // namespace algorithm
#endif
#line 1 "algorithm/Math/Algebra/kitamasa.hpp"
#include <cassert>
#include <iterator>
#include <utility>
#include <vector>
namespace algorithm {
// きたまさ法.
template <typename T>
class Kitamasa {
public:
using value_type = T;
using size_type = std::size_t;
private:
size_type m_k; // m_k:=(階数).
std::vector<value_type> m_a; // m_a[]:=(初項数列).
std::vector<value_type> m_d; // m_d[]:=(係数数列).
// f(n)->f(n+1). O(K).
std::vector<value_type> add(const std::vector<value_type> &x) const {
std::vector<value_type> y(m_k);
y[0] = x[m_k - 1] * m_d[0];
for(size_type i = 1; i < m_k; ++i) y[i] = x[i - 1] + x[m_k - 1] * m_d[i];
return y;
}
// f(n)->f(2*n). O(K^2).
std::vector<value_type> mul(const std::vector<value_type> &x) const {
std::vector<value_type> y(m_k, 0);
std::vector<value_type> t = x;
for(size_type i = 0; i < m_k; ++i) {
for(size_type j = 0; j < m_k; ++j) y[j] += x[i] * t[j];
if(i == m_k - 1) break;
t = add(t);
}
return y;
}
// f(n)を返す.O((K^2)*logN).
std::vector<value_type> f(size_type n) const {
if(n == 0) {
std::vector<value_type> x(m_k, 0);
x[0] = 1;
return x; // f(0).
}
std::vector<value_type> &&x = mul(f(n >> 1));
if(n & 1ULL) x = add(x);
return x;
}
public:
Kitamasa() : Kitamasa({0, 1}, {1, 1}) {} // フィボナッチ数列.
template <std::input_iterator InputIter>
explicit Kitamasa(InputIter a_first, InputIter a_last, InputIter d_first, InputIter d_last) : m_a(a_first, a_last), m_d(d_first, d_last) {
assert(m_a.size() >= 1);
assert(m_a.size() == m_d.size());
m_k = m_a.size();
}
explicit Kitamasa(std::vector<value_type> a, std::vector<value_type> d) : m_k(a.size()), m_a(std::move(a)), m_d(std::move(d)) {
assert(m_a.size() >= 1);
assert(m_a.size() == m_d.size());
}
value_type operator[](size_type n) const { return calc(n); }
// a[n]を求める.O((K^2)*log N).
value_type calc(size_type n) const {
std::vector<value_type> &&x = f(n);
value_type res = 0;
for(size_type i = 0; i < m_k; ++i) res += m_a[i] * x[i];
return res;
}
// a[l:r]を求める.O((K^2)*log(N) + K*(r-l)).
std::vector<value_type> calc(size_type l, size_type r) const {
assert(l < r);
std::vector<value_type> res(r - l);
std::vector<value_type> &&x = f(l);
for(size_type i = l; i < r; ++i) {
for(size_type j = 0; j < m_k; ++j) res[i - l] += m_a[j] * x[j];
if(i == r - 1) break;
x = add(x);
}
return res;
}
};
} // namespace algorithm