Kth term of linearly recurrent sequence
(Mylib/Math/linearly_recurrent_sequence.cpp)

• View this file on GitHub
• Last update: 2021-04-23 23:44:44+09:00
• Link: View error logs on GitHub Actions

Operations

$a_i \equiv \sum_{j=0}^{d - 1}c_j a_{i - d + j}$を満たす数列aの$a_0 \ldots a_{d-1}$と係数$c_0 \ldots c_{d-1}$から$a_n$を求める。

Requirements

Notes

Problems

References

• https://arxiv.org/abs/2008.08822

Verified with

• test/yosupo-judge/kth_term_of_linearly_recurrent_sequence/main.test.cpp

Code

#pragma once
#include <cassert>
#include <cstdint>
#include <vector>

namespace haar_lib {
  template <typename T, auto &convolve>
  T linearly_recurrent_sequence(const std::vector<T> &a, const std::vector<T> &c, int64_t k) {
    assert(a.size() == c.size());

    const int d = a.size();

    std::vector<T> Q(d + 1);
    Q[0] = 1;
    for (int i = 0; i < d; ++i) Q[d - i] = -c[i];

    std::vector<T> P = convolve(a, Q);
    P.resize(d);

    while (k > 0) {
      auto q = Q;
      for (size_t i = 1; i < q.size(); i += 2) q[i] = -q[i];
      auto U = convolve(P, q);
      auto A = convolve(Q, q);

      if (k % 2 == 0) {
        for (int i = 0; i < d; ++i) P[i] = i * 2 < (int) U.size() ? U[i * 2] : 0;
      } else {
        for (int i = 0; i < d; ++i) P[i] = i * 2 + 1 < (int) U.size() ? U[i * 2 + 1] : 0;
      }

      for (int i = 0; i <= d; ++i) Q[i] = i * 2 < (int) A.size() ? A[i * 2] : 0;

      k >>= 1;
    }

    return P[0];
  }
}  // namespace haar_lib

#line 2 "Mylib/Math/linearly_recurrent_sequence.cpp"
#include <cassert>
#include <cstdint>
#include <vector>

namespace haar_lib {
  template <typename T, auto &convolve>
  T linearly_recurrent_sequence(const std::vector<T> &a, const std::vector<T> &c, int64_t k) {
    assert(a.size() == c.size());

    const int d = a.size();

    std::vector<T> Q(d + 1);
    Q[0] = 1;
    for (int i = 0; i < d; ++i) Q[d - i] = -c[i];

    std::vector<T> P = convolve(a, Q);
    P.resize(d);

    while (k > 0) {
      auto q = Q;
      for (size_t i = 1; i < q.size(); i += 2) q[i] = -q[i];
      auto U = convolve(P, q);
      auto A = convolve(Q, q);

      if (k % 2 == 0) {
        for (int i = 0; i < d; ++i) P[i] = i * 2 < (int) U.size() ? U[i * 2] : 0;
      } else {
        for (int i = 0; i < d; ++i) P[i] = i * 2 + 1 < (int) U.size() ? U[i * 2 + 1] : 0;
      }

      for (int i = 0; i <= d; ++i) Q[i] = i * 2 < (int) A.size() ? A[i * 2] : 0;

      k >>= 1;
    }

    return P[0];
  }
}  // namespace haar_lib

Back to top page