algorithm-dev

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

View the Project on GitHub today2098/algorithm-dev

:heavy_check_mark: きたまさ法
(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,*)$ を定義する（ただし，$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,) \rightarrow D(n+1,), D(n,) \rightarrow D(2n,)$ を計算することにより，繰り返し二乗法と同じ要領で $D(0,)$ から $D(N,)$ を求めることである．

$D(n,) \rightarrow D(n+1,)$ の計算

まず，前提として $D(n,*)$ が分かっていると仮定する． このとき，

\[\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,) \rightarrow D(2n,)$ の計算

先の方法で $D(n,), D(n+1,), \ldots, D(n+k-1,*)$ が $\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,)$ から $D(N,)$ を $\mathcal{O}(k^2 \log N)$ で求めることができる．

参考文献

  1. “漸化式”. Wikipedia. https://ja.wikipedia.org/wiki/漸化式.
  2. yosupo. “きたまさ法メモ”. HatenaBlog. https://yosupo.hatenablog.com/entry/2015/03/27/025132.
  3. smijake3. “Kitamasa法”. HatenaBlog. https://smijake3.hatenablog.com/entry/2017/01/02/024712.

問題

Verified with

Code

#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<T> m_a;  // m_a[]:=(初項数列).
    std::vector<T> m_d;  // m_d[]:=(係数数列).

    // f(n)->f(n+1). O(K).
    std::vector<T> add(const std::vector<T> &x) const {
        std::vector<T> 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<T> mul(const std::vector<T> &x) const {
        std::vector<T> y(m_k, 0);
        std::vector<T> 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<T> f(size_type n) const {
        if(n == 0) {
            std::vector<T> x(m_k, 0);
            x[0] = 1;
            return x;  // f(0).
        }
        std::vector<T> &&x = mul(f(n >> 1));
        if(n & 1ULL) x = add(x);
        return x;
    }

public:
    Kitamasa() : Kitamasa({0, 1}, {1, 1}) {}  // フィボナッチ数列．
    explicit Kitamasa(std::vector<T> a, std::vector<T> 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());
    }
    template <std::input_iterator Iter>
    explicit Kitamasa(Iter a_first, Iter a_last, Iter d_first, Iter 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();
    }

    T operator[](size_type n) const { return calc(n); }

    // a[n]を求める．O((K^2)*logN).
    T calc(size_type n) const {
        const std::vector<T> &&x = f(n);
        T 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)*logN+K*(r-l)).
    std::vector<T> calc(size_type l, size_type r) const {
        assert(l < r);
        std::vector<T> res(r - l);
        std::vector<T> &&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<T> m_a;  // m_a[]:=(初項数列).
    std::vector<T> m_d;  // m_d[]:=(係数数列).

    // f(n)->f(n+1). O(K).
    std::vector<T> add(const std::vector<T> &x) const {
        std::vector<T> 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<T> mul(const std::vector<T> &x) const {
        std::vector<T> y(m_k, 0);
        std::vector<T> 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<T> f(size_type n) const {
        if(n == 0) {
            std::vector<T> x(m_k, 0);
            x[0] = 1;
            return x;  // f(0).
        }
        std::vector<T> &&x = mul(f(n >> 1));
        if(n & 1ULL) x = add(x);
        return x;
    }

public:
    Kitamasa() : Kitamasa({0, 1}, {1, 1}) {}  // フィボナッチ数列．
    explicit Kitamasa(std::vector<T> a, std::vector<T> 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());
    }
    template <std::input_iterator Iter>
    explicit Kitamasa(Iter a_first, Iter a_last, Iter d_first, Iter 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();
    }

    T operator[](size_type n) const { return calc(n); }

    // a[n]を求める．O((K^2)*logN).
    T calc(size_type n) const {
        const std::vector<T> &&x = f(n);
        T 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)*logN+K*(r-l)).
    std::vector<T> calc(size_type l, size_type r) const {
        assert(l < r);
        std::vector<T> res(r - l);
        std::vector<T> &&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
Back to top page