#line 1 "test/yosupo-judge/kth_term_of_linearly_recurrent_sequence/main.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/kth_term_of_linearly_recurrent_sequence"
#include <iostream>
#line 2 "Mylib/Convolution/ntt_convolution.cpp"
#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>
#line 4 "Mylib/Number/Mint/mint.cpp"
namespace haar_lib {
template <int32_t M>
class modint {
uint32_t val_;
public:
constexpr static auto mod() { return M; }
constexpr modint() : val_(0) {}
constexpr modint(int64_t n) {
if (n >= M)
val_ = n % M;
else if (n < 0)
val_ = n % M + M;
else
val_ = n;
}
constexpr auto &operator=(const modint &a) {
val_ = a.val_;
return *this;
}
constexpr auto &operator+=(const modint &a) {
if (val_ + a.val_ >= M)
val_ = (uint64_t) val_ + a.val_ - M;
else
val_ += a.val_;
return *this;
}
constexpr auto &operator-=(const modint &a) {
if (val_ < a.val_) val_ += M;
val_ -= a.val_;
return *this;
}
constexpr auto &operator*=(const modint &a) {
val_ = (uint64_t) val_ * a.val_ % M;
return *this;
}
constexpr auto &operator/=(const modint &a) {
val_ = (uint64_t) val_ * a.inv().val_ % M;
return *this;
}
constexpr auto operator+(const modint &a) const { return modint(*this) += a; }
constexpr auto operator-(const modint &a) const { return modint(*this) -= a; }
constexpr auto operator*(const modint &a) const { return modint(*this) *= a; }
constexpr auto operator/(const modint &a) const { return modint(*this) /= a; }
constexpr bool operator==(const modint &a) const { return val_ == a.val_; }
constexpr bool operator!=(const modint &a) const { return val_ != a.val_; }
constexpr auto &operator++() {
*this += 1;
return *this;
}
constexpr auto &operator--() {
*this -= 1;
return *this;
}
constexpr auto operator++(int) {
auto t = *this;
*this += 1;
return t;
}
constexpr auto operator--(int) {
auto t = *this;
*this -= 1;
return t;
}
constexpr static modint pow(int64_t n, int64_t p) {
if (p < 0) return pow(n, -p).inv();
int64_t ret = 1, e = n % M;
for (; p; (e *= e) %= M, p >>= 1)
if (p & 1) (ret *= e) %= M;
return ret;
}
constexpr static modint inv(int64_t a) {
int64_t b = M, u = 1, v = 0;
while (b) {
int64_t t = a / b;
a -= t * b;
std::swap(a, b);
u -= t * v;
std::swap(u, v);
}
u %= M;
if (u < 0) u += M;
return u;
}
constexpr static auto frac(int64_t a, int64_t b) { return modint(a) / modint(b); }
constexpr auto pow(int64_t p) const { return pow(val_, p); }
constexpr auto inv() const { return inv(val_); }
friend constexpr auto operator-(const modint &a) { return modint(M - a.val_); }
friend constexpr auto operator+(int64_t a, const modint &b) { return modint(a) + b; }
friend constexpr auto operator-(int64_t a, const modint &b) { return modint(a) - b; }
friend constexpr auto operator*(int64_t a, const modint &b) { return modint(a) * b; }
friend constexpr auto operator/(int64_t a, const modint &b) { return modint(a) / b; }
friend std::istream &operator>>(std::istream &s, modint &a) {
s >> a.val_;
return s;
}
friend std::ostream &operator<<(std::ostream &s, const modint &a) {
s << a.val_;
return s;
}
template <int N>
static auto div() {
static auto value = inv(N);
return value;
}
explicit operator int32_t() const noexcept { return val_; }
explicit operator int64_t() const noexcept { return val_; }
};
} // namespace haar_lib
#line 7 "Mylib/Convolution/ntt_convolution.cpp"
namespace haar_lib {
template <typename T, int PRIM_ROOT, int MAX_SIZE>
class number_theoretic_transform {
public:
using value_type = T;
constexpr static int primitive_root = PRIM_ROOT;
constexpr static int max_size = MAX_SIZE;
private:
const int MAX_POWER_;
std::vector<T> BASE_, INV_BASE_;
public:
number_theoretic_transform() : MAX_POWER_(__builtin_ctz(MAX_SIZE)),
BASE_(MAX_POWER_ + 1),
INV_BASE_(MAX_POWER_ + 1) {
static_assert((MAX_SIZE & (MAX_SIZE - 1)) == 0, "MAX_SIZE must be power of 2.");
static_assert((T::mod() - 1) % MAX_SIZE == 0);
T t = T::pow(PRIM_ROOT, (T::mod() - 1) >> (MAX_POWER_ + 2));
T s = t.inv();
for (int i = MAX_POWER_; --i >= 0;) {
t *= t;
s *= s;
BASE_[i] = -t;
INV_BASE_[i] = -s;
}
}
void run(std::vector<T> &f, bool INVERSE = false) const {
const int n = f.size();
assert((n & (n - 1)) == 0 and n <= MAX_SIZE); // データ数は2の冪乗個
if (INVERSE) {
for (int b = 1; b < n; b <<= 1) {
T w = 1;
for (int j = 0, k = 1; j < n; j += 2 * b, ++k) {
for (int i = 0; i < b; ++i) {
const auto s = f[i + j];
const auto t = f[i + j + b];
f[i + j] = s + t;
f[i + j + b] = (s - t) * w;
}
w *= INV_BASE_[__builtin_ctz(k)];
}
}
const T t = T::inv(n);
for (auto &x : f) x *= t;
} else {
for (int b = n >> 1; b; b >>= 1) {
T w = 1;
for (int j = 0, k = 1; j < n; j += 2 * b, ++k) {
for (int i = 0; i < b; ++i) {
const auto s = f[i + j];
const auto t = f[i + j + b] * w;
f[i + j] = s + t;
f[i + j + b] = s - t;
}
w *= BASE_[__builtin_ctz(k)];
}
}
}
}
template <typename U>
std::vector<T> convolve(std::vector<U> f, std::vector<U> g, bool is_same = false) const {
const int m = f.size() + g.size() - 1;
int n = 1;
while (n < m) n *= 2;
std::vector<T> f2(n);
for (int i = 0; i < (int) f.size(); ++i) f2[i] = (int64_t) f[i];
run(f2);
if (is_same) {
for (int i = 0; i < n; ++i) f2[i] *= f2[i];
run(f2, true);
} else {
std::vector<T> g2(n);
for (int i = 0; i < (int) g.size(); ++i) g2[i] = (int64_t) g[i];
run(g2);
for (int i = 0; i < n; ++i) f2[i] *= g2[i];
run(f2, true);
}
return f2;
}
template <typename U>
std::vector<T> operator()(std::vector<U> f, std::vector<U> g, bool is_same = false) const {
return convolve(f, g, is_same);
}
};
template <typename T>
std::vector<T> convolve_general_mod(std::vector<T> f, std::vector<T> g) {
static constexpr int M1 = 167772161, P1 = 3;
static constexpr int M2 = 469762049, P2 = 3;
static constexpr int M3 = 1224736769, P3 = 3;
auto res1 = number_theoretic_transform<modint<M1>, P1, 1 << 20>().convolve(f, g);
auto res2 = number_theoretic_transform<modint<M2>, P2, 1 << 20>().convolve(f, g);
auto res3 = number_theoretic_transform<modint<M3>, P3, 1 << 20>().convolve(f, g);
const int n = res1.size();
std::vector<T> ret(n);
const int64_t M12 = (int64_t) modint<M2>::inv(M1);
const int64_t M13 = (int64_t) modint<M3>::inv(M1);
const int64_t M23 = (int64_t) modint<M3>::inv(M2);
for (int i = 0; i < n; ++i) {
const int64_t r[3] = {(int64_t) res1[i], (int64_t) res2[i], (int64_t) res3[i]};
const int64_t t0 = r[0] % M1;
const int64_t t1 = (r[1] - t0 + M2) * M12 % M2;
const int64_t t2 = ((r[2] - t0 + M3) * M13 % M3 - t1 + M3) * M23 % M3;
ret[i] = T(t0) + T(t1) * M1 + T(t2) * M1 * M2;
}
return ret;
}
} // namespace haar_lib
#line 4 "Mylib/IO/input_vector.cpp"
namespace haar_lib {
template <typename T>
std::vector<T> input_vector(int N) {
std::vector<T> ret(N);
for (int i = 0; i < N; ++i) std::cin >> ret[i];
return ret;
}
template <typename T>
std::vector<std::vector<T>> input_vector(int N, int M) {
std::vector<std::vector<T>> ret(N);
for (int i = 0; i < N; ++i) ret[i] = input_vector<T>(M);
return ret;
}
} // namespace haar_lib
#line 3 "Mylib/Math/linearly_recurrent_sequence.cpp"
#include <cstdint>
#line 5 "Mylib/Math/linearly_recurrent_sequence.cpp"
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 3 "Mylib/Number/Mod/mod_pow.cpp"
namespace haar_lib {
constexpr int64_t mod_pow(int64_t n, int64_t p, int64_t m) {
int64_t ret = 1;
while (p > 0) {
if (p & 1) (ret *= n) %= m;
(n *= n) %= m;
p >>= 1;
}
return ret;
}
} // namespace haar_lib
#line 3 "Mylib/Number/Prime/primitive_root.cpp"
namespace haar_lib {
constexpr int primitive_root(int p) {
int pf[30] = {};
int k = 0;
{
int n = p - 1;
for (int64_t i = 2; i * i <= p; ++i) {
if (n % i == 0) {
pf[k++] = i;
while (n % i == 0) n /= i;
}
}
if (n != 1)
pf[k++] = n;
}
for (int g = 2; g <= p; ++g) {
bool ok = true;
for (int i = 0; i < k; ++i) {
if (mod_pow(g, (p - 1) / pf[i], p) == 1) {
ok = false;
break;
}
}
if (not ok) continue;
return g;
}
return -1;
}
} // namespace haar_lib
#line 9 "test/yosupo-judge/kth_term_of_linearly_recurrent_sequence/main.test.cpp"
namespace hl = haar_lib;
constexpr int mod = 998244353;
constexpr int prim_root = hl::primitive_root(mod);
using mint = hl::modint<mod>;
using NTT = hl::number_theoretic_transform<mint, prim_root, 1 << 21>;
const static auto ntt = NTT();
int main() {
std::cin.tie(0);
std::ios::sync_with_stdio(false);
int d;
std::cin >> d;
int64_t k;
std::cin >> k;
auto a = hl::input_vector<mint>(d);
auto c = hl::input_vector<mint>(d);
std::reverse(c.begin(), c.end());
auto ans = hl::linearly_recurrent_sequence<mint, ntt>(a, c, k);
std::cout << ans << "\n";
return 0;
}
test/yosupo-judge/kth_term_of_linearly_recurrent_sequence/main.test.cpp
Input vector
(Mylib/IO/input_vector.cpp)