test/lc_kth_term_of_linearly_recurrent_sequence.cpp

View this file on GitHub
Last update: 2025-02-05 03:54:30+09:00
Problem: https://judge.yosupo.jp/problem/kth_term_of_linearly_recurrent_sequence

Depends on

Code

// competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/kth_term_of_linearly_recurrent_sequence

#include "../template.hpp"

#include "../formal_power_series.hpp"

int main() {
    long long d, k;
    cin >> d >> k;
    vector<mint> a(d), c(d);
    cin >> a >> c;
    cout << FPS::kth_term(a, c, k) << endl;
}

#line 1 "test/lc_kth_term_of_linearly_recurrent_sequence.cpp"
// competitive-verifier: PROBLEM https://judge.yosupo.jp/problem/kth_term_of_linearly_recurrent_sequence

#line 1 "template.hpp"
#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <climits>
#include <cmath>
#include <complex>
#include <cstring>
#include <deque>
#include <functional>
#include <iomanip>
#include <iostream>
#include <map>
#include <memory>
#include <numeric>
#include <optional>
#include <queue>
#include <random>
#include <set>
#include <stack>
#include <string>
#include <tuple>
#include <unordered_map>
#include <unordered_set>
#include <vector>

#if __cplusplus >= 202002L
#include <bit>
#include <ranges>
#define TYPE(n) remove_cvref_t<decltype(n)>
#else
#define countl_zero __builtin_clzll
#define TYPE(n) remove_cv_t<remove_reference_t<decltype(n)>>
#endif

using namespace std;
using lint = long long;
using P = pair<lint, lint>;
using Pii = pair<int, int>;
using ull = unsigned long long;

struct FastIO {
    FastIO() {
        ios::sync_with_stdio(false);
        cin.tie(nullptr);
        cout << fixed << setprecision(15);
    }
} aaaAAAaaaAAA;

#define rep(i, n) for (TYPE(n) i = 0; i < (n); ++i)
#define repe(i, n) for (TYPE(n) i = 0; i <= (n); ++i)
#define rep1(i, n) for (TYPE(n) i = 1; i < (n); ++i)
#define rep1e(i, n) for (TYPE(n) i = 1; i <= (n); ++i)
#define repn(i, a, b) for (TYPE(a) i = (a); i < (b); ++i)
#define repne(i, a, b) for (TYPE(a) i = (a); i <= (b); ++i)
#define rrep(i, n) for (TYPE(n) i = (n); i >= 0; --i)
#define all(vec) begin(vec), end(vec)
#define rall(vec) rbegin(vec), rend(vec)

constexpr long long Mod = /** 1000'000'007LL /*/ 998244353LL /**/;
constexpr long long Inf = 2'000'000'000'000'000'010LL;
constexpr int IntInf = 1000'000'010;
constexpr double Pi = 3.141592653589793238;
constexpr double InvPi = 0.318309886183790671;

const int di[] = {0, -1, 0, 1, 0};
const int dj[] = {1, 0, -1, 0, 0};
pair<int, int> adj(int i, int j, int k) { return {i + di[k], j + dj[k]}; }
bool in(int i, int j, int h, int w) { return 0 <= i && i < h && 0 <= j && j < w; }

#if __has_include(<atcoder/all>)
#include <atcoder/all>

using namespace atcoder;
using mint = static_modint<Mod>;

template <int MOD>
inline istream &operator>>(istream &is, static_modint<MOD> &rhs) {
    long long tmp;
    is >> tmp;
    rhs = tmp;
    return is;
}
template <int ID>
inline istream &operator>>(istream &is, dynamic_modint<ID> &rhs) {
    long long tmp;
    is >> tmp;
    rhs = tmp;
    return is;
}
template <int MOD>
inline ostream &operator<<(ostream &os, const static_modint<MOD> &rhs) {
    return os << rhs.val();
}
template <int ID>
inline ostream &operator<<(ostream &os, const dynamic_modint<ID> &rhs) {
    return os << rhs.val();
}
// [0, n]
template <typename T>
auto enumerate_fact(int n) {
    vector<T> fact(n + 1);
    fact[0] = 1;
    for (int i = 1; i <= n; ++i) fact[i] = i * fact[i - 1];
    return fact;
}
// [0, n]
template <int MOD>
auto enumerate_inv(int n) {
    vector<static_modint<MOD>> inv(n + 1);
    inv[1] = 1;
    for (int i = 2; i <= n; ++i) inv[i] = MOD - MOD / i * inv[MOD % i];
    return inv;
}
// [0, n]
template <typename T>
auto enumerate_factinv(int n, vector<T> inv) {
    vector<T> fact_inv(n + 1);
    fact_inv[0] = 1;
    for (int i = 1; i <= n; ++i) fact_inv[i] = fact_inv[i - 1] * inv[i];
    return fact_inv;
}
// [0, n]
template <int MOD>
auto enumerate_factinv(int n) {
    return enumerate_factinv(n, enumerate_inv<MOD>(n));
}

template <int MOD>
struct Binomial {
    using Fp = static_modint<MOD>;
    vector<Fp> fact, inv, fact_inv;

    explicit Binomial() = default;

    // [0, n]
    void build(int n) {
        fact = enumerate_fact<Fp>(n);
        inv = enumerate_inv<MOD>(n);
        fact_inv = enumerate_factinv<Fp>(n, inv);
    }

    Fp comb(int n, int r) const {
        if (n < 0 || r < 0 || n < r) return 0;
        if (r == 0 || r == n) return 1;
        return fact[n] * fact_inv[n - r] * fact_inv[r];
    }

    Fp perm(int n, int r) const {
        if (n < 0 || r < 0 || n < r) return 0;
        return fact[n] * fact_inv[n - r];
    }

    Fp multi(int n, int r) const {
        if (n == 0 && r == 0) return 1;
        if (n < 0 || r < 0) return 0;
        return r == 0 ? 1 : comb(n + r - 1, r);
    }
};
Binomial<Mod> binomial;
inline mint fact(int n) { return binomial.fact[n]; }
inline mint comb(int n, int r) { return binomial.comb(n, r); }
inline mint perm(int n, int r) { return binomial.perm(n, r); }
inline mint multi(int n, int r) { return binomial.multi(n, r); }

mint lagrange_interpolation(const vector<mint> &y, mint t) {
    const int n = (int)y.size();

    mint res = 0;

    auto inv = enumerate_inv<Mod>(n - 1), fact_inv = enumerate_factinv(n - 1, inv);

    vector<mint> prod2(n);
    prod2.back() = 1;
    for (int i = n - 1; i > 0; --i) {
        prod2[i - 1] = (t - i) * prod2[i];
    }

    mint prod1 = 1;
    for (int i = 0; i < n; ++i) {
        mint a = y[i];
        a *= fact_inv[i] * fact_inv[n - 1 - i];
        if ((n - 1 - i) & 1) a = -a;

        res += a * prod1 * prod2[i];

        prod1 *= (t - i);
    }

    return res;
}

template <typename T>
lint inversion_number(const vector<T> vec) {
    vector<T> v = vec;
    sort(v.begin(), v.end());
    v.erase(unique(v.begin(), v.end()), v.end());

    const int n = vec.size();

    lint res = 0;

    fenwick_tree<int> b(n);
    for (int i = 0; i < n; ++i) {
        const int j = lower_bound(v.begin(), v.end(), vec[i]) - v.begin();
        res += b.sum(j + 1, n);
        b.add(j, 1);
    }

    return res;
}
#endif

// top = max
template <typename T>
using prique = priority_queue<T>;
// top = min
template <typename T>
using prique_inv = priority_queue<T, vector<T>, greater<T>>;

template <typename T, typename U>
inline istream &operator>>(istream &is, pair<T, U> &rhs) {
    return is >> rhs.first >> rhs.second;
}
template <typename T>
inline istream &operator>>(istream &is, complex<T> &c) {
    T real, imag;
    is >> real >> imag;
    c.real(real);
    c.imag(imag);
    return is;
}
template <typename T, typename U>
inline ostream &operator<<(ostream &os, const pair<T, U> &rhs) {
    return os << "{" << rhs.first << ", " << rhs.second << "}";
}

#if __cplusplus >= 202002L
template <class T>
concept Container = requires(const T &v) {
    v.begin();
    v.end();
} && !is_same_v<T, string>;
template <Container T>
inline istream &operator>>(istream &is, T &v) {
    for (auto &e : v) is >> e;
    return is;
}
template <Container T>
inline ostream &operator<<(ostream &os, const T &v) {
    for (auto itr = v.begin(), end_itr = v.end(); itr != end_itr;) {
        os << *itr;
        if (++itr != end_itr) os << " ";
    }
    return os;
}
#else
template <typename T>
inline istream &operator>>(istream &is, vector<T> &v) {
    for (auto &e : v) is >> e;
    return is;
}
template <typename T>
inline ostream &operator<<(ostream &os, const vector<T> &v) {
    for (auto itr = v.begin(), end_itr = v.end(); itr != end_itr;) {
        os << *itr;
        if (++itr != end_itr) os << " ";
    }
    return os;
}
#endif

template <typename T, typename U>
inline bool chmin(T &a, const U b) {
    return a > b ? a = b, true : false;
}
template <typename T, typename U>
inline bool chmax(T &a, const U b) {
    return a < b ? a = b, true : false;
}
template <typename T, typename U, class Pr>
inline int getid(const vector<T> &v, const U &value, Pr pred) {
    return lower_bound(v.begin(), v.end(), value, pred) - v.begin();
}
template <typename T, typename U>
inline int getid(const vector<T> &v, const U &value) {
    return getid(v, value, less<>{});
}

template <typename T>
T gcd(const vector<T> &vec) {
    T res = vec.front();
    for (T e : vec) {
        res = gcd(res, e);
        if (res == 1) return 1;
    }
    return res;
}
template <typename T>
T gcd(initializer_list<T> init) {
    auto first = init.begin(), last = init.end();
    T res = *(first++);
    for (auto itr = first; itr != last; ++itr) {
        res = gcd(res, *itr);
        if (res == 1) return 1;
    }
    return res;
}
template <typename T>
T lcm(const vector<T> &vec) {
    T res = vec.front();
    for (T e : vec) res = lcm(res, e);
    return res;
}
template <typename T>
T lcm(initializer_list<T> init) {
    auto first = init.begin(), last = init.end();
    T res = *(first++);
    for (auto itr = first; itr != last; ++itr) {
        res = lcm(res, *itr);
    }
    return res;
}

inline void YesNo(bool b) { cout << (b ? "Yes\n" : "No\n"); }
inline void YESNO(bool b) { cout << (b ? "YES\n" : "NO\n"); }
inline void takaao(bool b) { cout << (b ? "Takahashi\n" : "Aoki\n"); }
inline void aotaka(bool b) { cout << (b ? "Aoki\n" : "Takahashi\n"); }

// [l, r]
template <typename T>
T rand(T l, T r) {
    static mt19937 mt(random_device{}());
    if constexpr (is_integral_v<T>) {
        uniform_int_distribution<T> dist(l, r);
        return dist(mt);
    } else if constexpr (is_floating_point_v<T>) {
        uniform_real_distribution<T> dist(l, r);
        return dist(mt);
    }
}

bool is_prime_naive(lint x) {
    for (lint i = 2; i * i <= x; ++i) {
        if (x % i == 0) return false;
    }
    return 1 < x;
}

vector<lint> divisors(lint n) {
    vector<lint> f, l;
    for (lint x = 1; x * x <= n; ++x) {
        if (n % x == 0) {
            f.push_back(x);
            if (x * x != n) l.push_back(n / x);
        }
    }
    f.reserve(f.size() + l.size());
    copy(l.rbegin(), l.rend(), back_inserter(f));
    return f;
}

lint phi(lint n) {
    lint r = n;
    for (lint i = 2; i * i <= n; ++i) {
        if (n % i == 0) {
            r -= r / i;
            while (n % i == 0) n /= i;
        }
    }
    if (n > 1) r -= r / n;
    return r;
}

lint floor_sqrt(lint n) {
    if (n == 0 || n == 1) return n;
    lint x0 = 1LL << ((65 - countl_zero(static_cast<uint64_t>(n))) >> 1);
    lint x1 = (x0 + n / x0) >> 1;
    while (x1 < x0) {
        x0 = x1;
        x1 = (x0 + n / x0) >> 1;
    }
    return x0;
}
lint ceil_sqrt(lint n) {
    const lint f = floor_sqrt(n);
    if (f * f == n) return f;
    return f + 1;
}

template <typename T>
constexpr bool is_intersect(T l1, T r1, T l2, T r2) {
    return l1 <= r2 && l2 <= r1;
}
template <typename T>
constexpr bool is_intersect2(T l1, T r1, T l2, T r2) {
    return l1 < r2 && l2 < r1;
}

lint modinv(lint a, lint m = Mod) {
    lint b = m, u = 1, v = 0;
    while (b != 0) {
        lint t = a / b;
        a -= t * b;
        swap(a, b);
        u -= t * v;
        swap(u, v);
    }
    u %= m;
    if (u < 0) u += m;
    return u;
}

template <typename T>
T modpow(T x, T n, T m = Mod) {
    if (m == 1) return 0;
    T res = 1;
    x %= m;
    if (x < 0) x += m;
    while (n > 0) {
        if (n & 1) res = res * x % m;
        x = x * x % m;
        n >>= 1;
    }
    return res;
}

template <typename T>
T intpow(T x, T n) {
    T res = 1;
    while (n > 0) {
        if (n & 1) res *= x;
        x *= x;
        n >>= 1;
    }
    return res;
}

template <typename T>
vector<T> compressed(vector<T> v) {
    sort(v.begin(), v.end());
    v.erase(unique(v.begin(), v.end()), v.end());
    return v;
}

class Sieve {
   private:
    const int max_n;
    vector<int> sieve;

   public:
    explicit Sieve(int max_n) : max_n(max_n), sieve(max_n + 1) {
        iota(sieve.begin(), sieve.end(), 0);

        for (int i = 2; i * i <= max_n; ++i) {
            if (sieve[i] < i) continue;

            for (int j = i * i; j <= max_n; j += i) {
                if (sieve[j] == j) sieve[j] = i;
            }
        }
    }

    unordered_map<int, int> calc(int x) const {
        assert(x <= max_n);
        unordered_map<int, int> res;
        while (x > 1) {
            ++res[sieve[x]];
            x /= sieve[x];
        }
        return res;
    }

    vector<int> enumerate_prime(int x) const {
        assert(x <= max_n);
        vector<int> primes;
        for (int i = 2; i <= x; ++i) {
            if (sieve[i] == i) {
                primes.push_back(i);
            }
        }
        return primes;
    }
};

struct UnionFind {
    int n;
    vector<int> par, rank, siz, es;  // [root(i)]
    int c;

    UnionFind() = default;

    explicit UnionFind(int _n) : n(_n), par(_n), rank(_n), siz(_n, 1), es(_n), c(_n) { iota(par.begin(), par.end(), 0); }

    int root(int x) {
        while (par[x] != x) x = par[x] = par[par[x]];
        return x;
    }

    bool same(int x, int y) { return root(x) == root(y); }

    void unite(int x, int y) {
        if (x == y) return;

        x = root(x);
        y = root(y);
        if (x == y)
            ++es[x];
        else {
            c--;
            if (rank[x] < rank[y]) {
                par[x] = y;
                siz[y] += siz[x];
                es[y] += es[x] + 1;
            } else {
                par[y] = x;
                if (rank[x] == rank[y]) ++rank[x];
                siz[x] += siz[y];
                es[x] += es[y] + 1;
            }
        }
    }

    int size(int x) { return siz[root(x)]; }

    vector<int> roots() {
        vector<int> res;
        res.reserve(c);

        for (int i = 0; i < n; ++i) {
            if (par[i] == i) {
                res.push_back(i);
            }
        }

        return res;
    }

    vector<vector<int>> groups() {
        vector<vector<int>> res(n);

        for (int i = 0; i < n; ++i)
            if (par[i] == i) res[i].reserve(siz[i]);
        for (int i = 0; i < n; ++i) res[root(i)].push_back(i);

        res.erase(remove_if(res.begin(), res.end(), [](const vector<int> &v) { return v.empty(); }), res.end());

        return res;
    }
};

template <typename T>
class BinaryIndexedTree {
   private:
    int n;
    vector<T> dat;

   public:
    BinaryIndexedTree() = default;

    explicit BinaryIndexedTree(const int size) : n(size), dat(size + 1) {}

    explicit BinaryIndexedTree(const vector<T> &vec) : n(vec.size()), dat(n + 1) {
        for (int i = 0; i < n; ++i) dat[i + 1] = vec[i];
        for (int i = 1; i <= n; ++i) {
            const int j = i + (i & -i);
            if (j <= n) dat[j] += dat[i];
        }
    }

    // 0-indexed
    void add(const int a, const T v) {
        for (int x = a + 1; x <= n; x += (x & -x)) dat[x] += v;
    }

    // 0-indexed
    void set(const int a, const T v) { add(a, v - get(a)); }

    void reset() { fill(dat.begin(), dat.end(), T(0)); }

    // [0, a)
    T sum(const int a) const {
        T res = 0;
        for (int x = a; x > 0; x -= (x & -x)) res += dat[x];
        return res;
    }

    // [a, b)
    T sum(const int a, const int b) const { return sum(b) - sum(a); }

    T get(const int i) const { return sum(i, i + 1); }

    // min i s.t. sum(i) >= w
    int lower_bound(T w) const {
        int x = 0, r = 1;
        while (r < n) r <<= 1;
        for (int i = r; i > 0; i >>= 1) {
            if (x + i <= n && dat[x + i] < w) {
                w -= dat[x + i];
                x += i;
            }
        }
        return x + 1;
    }
    // min i s.t. sum(i) > w
    int upper_bound(T w) const {
        int x = 0, r = 1;
        while (r < n) r <<= 1;
        for (int i = r; i > 0; i >>= 1) {
            if (x + i <= n && dat[x + i] <= w) {
                w -= dat[x + i];
                x += i;
            }
        }
        return x + 1;
    }
};

constexpr int msb(long long x) { return 63 - countl_zero(static_cast<uint64_t>(x)); }
#line 4 "test/lc_kth_term_of_linearly_recurrent_sequence.cpp"

#line 1 "formal_power_series.hpp"
#include <atcoder/convolution>

template <int MOD>
struct FormalPowerSeries : public vector<atcoder::static_modint<MOD>> {
    using Fp = atcoder::static_modint<MOD>;

    using vector<Fp>::vector;
    using vector<Fp>::operator=;
    using F = FormalPowerSeries;

    FormalPowerSeries(const vector<Fp> &vec) { *this = vec; }

    void shrink() {
        while (!this->empty() && this->back() == 0) this->pop_back();
    }

    F operator+() const noexcept { return *this; }
    F operator-() const noexcept {
        F res(*this);

        for (auto &&e : res) e = -e;

        return res;
    }

    F operator*(const Fp &k) const noexcept { return F(*this) *= k; }
    F operator/(const Fp &k) const { return F(*this) /= k; }

    F operator+(const F &g) const noexcept { return F(*this) += g; }
    F operator-(const F &g) const noexcept { return F(*this) -= g; }

    F operator<<(const int d) const noexcept { return F(*this) <<= d; }
    F operator>>(const int d) const noexcept { return F(*this) >>= d; }

    F operator*(const F &g) const { return F(*this) *= g; }
    F operator/(const F &g) const { return F(*this) /= g; }
    F operator%(const F &g) const { return F(*this) %= g; }

    F &operator*=(const Fp &k) noexcept {
        for (auto &&e : *this) e *= k;
        return *this;
    }
    F &operator/=(const Fp &k) {
        assert(k != 0);
        *this *= k.inv();
        return *this;
    }

    F &operator+=(const F &g) noexcept {
        const int n = this->size(), m = g.size();
        this->resize(max(n, m));

        for (int i = 0; i < m; ++i) (*this)[i] += g[i];

        return *this;
    }
    F &operator-=(const F &g) noexcept {
        const int n = this->size(), m = g.size();
        this->resize(max(n, m));

        for (int i = 0; i < m; ++i) (*this)[i] -= g[i];

        return *this;
    }

    F &operator<<=(const int d) {
        const int n = this->size();

        this->insert(this->begin(), d, Fp(0));

        return *this;
    }
    F &operator>>=(const int d) {
        const int n = this->size();

        if (n <= d)
            this->clear();
        else
            this->erase(this->begin(), this->begin() + d);

        return *this;
    }

    F &operator*=(const F &g) {
        const auto f = atcoder::convolution(std::move(*this), g);
        return *this = F(f.begin(), f.end());
    }
    F &operator/=(const F &g) {
        if (this->size() < g.size()) {
            this->clear();
            return *this;
        }
        const int n = this->size() - g.size() + 1;
        return *this = (rev().pre(n) * g.rev().inv(n)).pre(n).rev(n);
    }
    F &operator%=(const F &g) {
        *this -= *this / g * g;

        this->shrink();

        return *this;
    }

    bool zero() const noexcept {
        bool res = true;
        for (const auto &e : *this) {
            res &= (e.val() == 0);
        }
        return res;
    }

    Fp eval(const Fp &x) const noexcept {
        Fp res = this->back();

        for (auto itr = ++this->rbegin(), itr_rend = this->rend(); itr != itr_rend; ++itr) {
            res *= x;
            res += *itr;
        }

        return res;
    }

    F pre(int d) const { return F(this->begin(), this->begin() + min((int)this->size(), d)); }
    F rev(int d = -1) const {
        F res(*this);
        if (d != -1) res.resize(d, Fp(0));
        reverse(res.begin(), res.end());
        return res;
    }

    F inv(int d = -1) const {
        int n = this->size();
        assert(n != 0 && this->front() != Fp(0));
        if (d == -1) d = n;
        assert(d > 0);

        F res = {1 / this->front()};
        res.reserve(2 * d);

        int m = res.size();
        while (m < d) {
            F f(this->begin(), this->begin() + min(n, 2 * m));
            F r(res);

            f.resize(2 * m);
            r.resize(2 * m);
            atcoder::internal::butterfly(f);
            atcoder::internal::butterfly(r);
            for (int i = 0; i < 2 * m; ++i) f[i] *= r[i];
            atcoder::internal::butterfly_inv(f);

            f.erase(f.begin(), f.begin() + m);
            f.resize(2 * m);
            atcoder::internal::butterfly(f);
            for (int i = 0; i < 2 * m; ++i) f[i] *= r[i];
            atcoder::internal::butterfly_inv(f);

            Fp iz = Fp(1) / (2 * m);
            iz *= -iz;

            for (int i = 0; i < m; ++i) f[i] *= iz;

            res.insert(res.end(), f.begin(), f.begin() + m);

            m <<= 1;
        }

        return {res.begin(), res.begin() + d};
    }

    F &multiply_inplace(const F &g, int d = -1) {
        if (d == -1) d = this->size();
        assert(d >= 0);

        *this = atcoder::convolution(move(*this), g);
        this->resize(d);

        return *this;
    }
    F multiply(const F &g, int d = -1) const { return F(*this).multiply_inplace(g, d); }

    F &diff_inplace() noexcept {
        const int n = this->size();

        for (int i = 1; i < n; ++i) (*this)[i - 1] = (*this)[i] * i;
        this->back() = 0;

        return *this;
    }
    F diff() const noexcept { return F(*this).diff_inplace(); }

    F &integral_inplace() noexcept {
        constexpr int p = Fp::mod();
        const int n = this->size();

        vector<Fp> inv(n);
        inv[1] = 1;
        for (int i = 2; i < n; ++i) inv[i] = -inv[p % i] * (p / i);

        for (int i = n - 1; i > 0; --i) (*this)[i] = (*this)[i - 1] * inv[i];

        this->front() = 0;

        return *this;
    }
    F integral() const noexcept { return F(*this).integral_inplace(); }

    F &log_inplace(int d = -1) {
        const int n = this->size();
        assert(this->front() == 1);
        if (d == -1) d = n;
        assert(d >= 0);

        this->resize(d);

        const F f_inv = this->inv();

        this->diff_inplace().multiply_inplace(f_inv).integral_inplace();

        return *this;
    }
    F log(int d = -1) const { return F(*this).log_inplace(d); }

    F &exp_inplace(int d = -1) {
        const int n = this->size();
        assert(this->front() == 0);
        if (d == -1) d = n;
        assert(d >= 0);

        F g = {1}, g_fft;

        this->resize(d);
        this->front() = 1;

        const F h_diff = this->diff();
        for (int m = 1; m < d; m <<= 1) {
            F f_fft(this->begin(), this->begin() + m);
            f_fft.resize(2 * m);
            atcoder::internal::butterfly(f_fft);
            if (m > 1) {
                F _f(m);
                for (int i = 0; i < m; ++i) _f[i] = f_fft[i] * g_fft[i];
                atcoder::internal::butterfly_inv(_f);
                _f.erase(_f.begin(), _f.begin() + m / 2);
                _f.resize(m);
                atcoder::internal::butterfly(_f);
                for (int i = 0; i < m; ++i) _f[i] *= g_fft[i];
                atcoder::internal::butterfly_inv(_f);
                _f.resize(m / 2);

                _f /= Fp(-m) * m;
                g.insert(g.end(), _f.begin(), _f.begin() + m / 2);
            }

            F t(this->begin(), this->begin() + m);
            t.diff_inplace();
            {
                F r(h_diff.begin(), h_diff.begin() + (m - 1));

                r.resize(m);
                atcoder::internal::butterfly(r);
                for (int i = 0; i < m; ++i) r[i] *= f_fft[i];
                atcoder::internal::butterfly_inv(r);
                r /= -m;

                t += r;
                t.insert(t.begin(), t.back());
                t.pop_back();
            }

            if (2 * m < d || m == 1) {
                t.resize(2 * m);
                atcoder::internal::butterfly(t);
                g_fft = g;
                g_fft.resize(2 * m);
                atcoder::internal::butterfly(g_fft);
                for (int i = 0; i < 2 * m; ++i) t[i] *= g_fft[i];
                atcoder::internal::butterfly_inv(t);
                t.resize(m);
                t /= 2 * m;
            } else {
                F g1(g.begin() + m / 2, g.end());
                F s1(t.begin() + m / 2, t.end());
                t.resize(m / 2);
                g1.resize(m);
                atcoder::internal::butterfly(g1);
                t.resize(m);
                atcoder::internal::butterfly(t);
                s1.resize(m);
                atcoder::internal::butterfly(s1);
                for (int i = 0; i < m; ++i) s1[i] = g_fft[i] * s1[i] + g1[i] * t[i];
                for (int i = 0; i < m; ++i) t[i] *= g_fft[i];
                atcoder::internal::butterfly_inv(t);
                atcoder::internal::butterfly_inv(s1);
                for (int i = 0; i < m / 2; ++i) t[i + m / 2] += s1[i];
                t /= m;
            }

            F v(this->begin() + m, this->begin() + min(d, 2 * m));
            v.resize(m);
            t.insert(t.begin(), m - 1, 0);
            t.push_back(0);
            t.integral_inplace();
            for (int i = 0; i < m; ++i) v[i] -= t[m + i];

            v.resize(2 * m);
            atcoder::internal::butterfly(v);
            for (int i = 0; i < 2 * m; ++i) v[i] *= f_fft[i];
            atcoder::internal::butterfly_inv(v);
            v.resize(m);
            v /= 2 * m;

            for (int i = 0; i < min(d - m, m); ++i) (*this)[m + i] = v[i];
        }

        return *this;
    }
    F exp(int d = -1) const { return F(*this).exp_inplace(d); }

    F &pow_inplace(long long m, long long d = -1) {
        const long long n = this->size();
        if (d == -1) d = n;
        assert(d > 0);

        if (m == 0) {
            F res(d);
            res[0] = 1;
            return *this = res;
        }
        if (zero()) {
            return *this = F(d);
        }

        long long k = 0;
        while (k < n && (*this)[k] == 0) ++k;
        if (k >= (d + m - 1) / m) return *this = F(d);

        const Fp c_inv = (*this)[k].inv();
        const Fp c_pow = (*this)[k].pow(m);

        this->erase(this->begin(), this->begin() + k);
        *this *= c_inv;
        this->log_inplace(d);
        *this *= m;
        this->exp_inplace(d);
        *this *= c_pow;

        this->insert(this->begin(), m * k, 0);
        this->resize(d);

        return *this;
    }
    F pow(long long m, int d = -1) const { return F(*this).pow_inplace(m, d); }

    // a.size() = n, c.size() = n + m - 1
    // res[j] = sum a[i] * c[i + j] (0 <= j < m)
    static vector<Fp> middle_product(vector<Fp> a, vector<Fp> c, bool c_reversed = false, bool b_reversed = false) {
        int n = a.size(), m = c.size() + 1 - n;
        if (m <= 0) return {};
        if (min(n, m) <= 60) return middle_product_naive(a, c, c_reversed, b_reversed);
        return middle_product_fft(a, c, c_reversed, b_reversed);
    }

    template <typename U>
    vector<Fp> eval(const vector<U> &x) const {
        const int n = this->size();
        const int m1 = x.size();

        if (m1 == 1) return {eval(x[0])};

        int m = 1;
        while (m < m1) m <<= 1;

        vector t(2 * m, vector<Fp>{1});
        for (int i = m; i < m + m1; ++i) {
            t[i].resize(2);
            t[i][0] = -x[i - m];
            t[i][1] = 1;
        }
        for (int i = m - 1; i >= 1; --i) t[i] = atcoder::convolution(t[i << 1], t[i << 1 | 1]);

        F t1 = F(t[1]).rev().inv(n);

        vector<Fp> f(*this);
        f.resize(n + m1 - 1);

        vector<Fp> a = middle_product(t1, f, false, false);

        vector b(2 * m, vector<Fp>{});
        b[1] = a;

        for (int i = 1; i < m; ++i) {
            b[i << 1 | 1] = middle_product(t[i << 1], b[i], true, true);
            b[i << 1] = middle_product(t[i << 1 | 1], b[i], true, true);
        }

        vector<Fp> res(m1);
        for (int i = m; i < m + m1; ++i) res[i - m] = b[i][0];

        return res;
    }

    friend F operator*(const Fp &k, const F &f) noexcept { return f * k; }

    static F prod(vector<F> fs) {
        const int n = fs.size();
        if (n == 0) return {1};

        for (int i = n - 1; i > 0; --i) {
            fs[i / 2] = atcoder::convolution(fs[i / 2], fs[i]);
        }
        return fs[0];
    }

    // prod (a_i x + b_i)
    template <typename U, typename V>
    static F prod(const vector<U> &a, const vector<V> &b) {
        const size_t n = a.size();

        vector<F> fs(n, F{0, 0});
        for (size_t i = 0; i < n; ++i) {
            fs[i][0] = b[i];
            fs[i][1] = a[i];
        }
        return prod(fs);
    }

    // pre: O(sqrt(MOD) log^2 MOD)
    // O(sqrt(MOD))
    static Fp factorial(long long n) {
        if (n >= MOD) return 0;
        static constexpr int v = 1 << 15;  // v * v >= MOD
        static vector<Fp> ps(v + 1);       // ps[i] = prod_1^{iv} j

        static bool init = false;
        if (!init) {
            init = true;
            vector<int> k(v);
            iota(k.begin(), k.end(), 1);
            const F f = prod(vector<int>(v, 1), k);  // (x+1)(x+2)...(x+v)
            vector<int> xs(v);
            for (int i = 0; i < v; ++i) xs[i] = v * i;
            vector<Fp> es = f.eval(xs);
            ps[0] = 1;
            for (int i = 0; i < v; ++i) ps[i + 1] = ps[i] * es[i];
        }
        const int m = min<int>(v, (n + 1) / v);

        Fp p = ps[m];
        for (int i = m * v + 1; i <= n; ++i) p *= i;
        return p;
    }

    // dft.size() == 2 * n, dft[0 : n] = DFT(f)
    // dft <- DFT(f + [0] * n)
    // time complexity: FFT(n)
    static void fft_doubling(vector<Fp> &dft, vector<Fp> f, const Fp r_2n) {
        const int n = dft.size() >> 1;
        Fp rp = 1;
        for (auto &e : f) {
            e *= rp;
            rp *= r_2n;
        }
        atcoder::internal::butterfly(f);
        copy(f.begin(), f.end(), dft.begin() + n);
    }

    // dft.size() == 2 * n
    // dft <- DFT(IDFT(dft[0 : n]) + [0] * n)
    // time complexity: 2 FFT(n)
    static void fft_doubling(vector<Fp> &dft, const Fp r, const Fp n_inv) {
        const int n = dft.size() >> 1;

        vector<Fp> b(n);
        copy(dft.begin(), dft.begin() + n, b.begin());
        atcoder::internal::butterfly_inv(b);
        Fp rp = 1;
        for (auto &e : b) {
            e *= rp * n_inv;
            rp *= r;
        }
        atcoder::internal::butterfly(b);
        copy(b.begin(), b.end(), dft.begin() + n);
    }

    // [x^k] (p/q) (deg p < deg q)
    static Fp coeff(vector<Fp> p, vector<Fp> q, long long k) {
        static const atcoder::internal::fft_info<Fp> info;
        static const Fp inv2 = Fp::raw((Fp::mod() + 1) / 2);

        const int n = atcoder::internal::bit_ceil((unsigned int)(q.size()));

        p.resize(2 * n), q.resize(2 * n);
        atcoder::internal::butterfly(p);
        atcoder::internal::butterfly(q);

        const int w = __builtin_ctz((unsigned int)(n));
        const Fp n_inv = Fp::raw(n).inv();
        const Fp r_z = info.root[w + 1];
        const Fp ir_z = info.iroot[w + 1];

        vector<Fp> ir_p(n, 1);
        for (int i = 0; i < n; ++i) {
            Fp ir_z_p = ir_z;
            for (int j = w - 1; j >= 0; --j) {
                if (i >> j & 1) ir_p[i] *= ir_z_p;
                ir_z_p *= ir_z_p;
            }
        }

        Fp inv2_p = 1;
        while (k > 0) {
            inv2_p *= inv2;

            if (k & 1) {
                for (int i = 0; i < n; ++i) p[i] = ir_p[i] * (p[i << 1] * q[i << 1 | 1] - p[i << 1 | 1] * q[i << 1]);
            } else {
                for (int i = 0; i < n; ++i) p[i] = p[i << 1] * q[i << 1 | 1] + p[i << 1 | 1] * q[i << 1];
            }
            for (int i = 0; i < n; ++i) q[i] = q[i << 1] * q[i << 1 | 1];

            fft_doubling(p, r_z, n_inv);
            fft_doubling(q, r_z, n_inv);

            k >>= 1;
            if (k < 2 * n) break;
        }
        atcoder::internal::butterfly_inv(p);
        atcoder::internal::butterfly_inv(q);

        F f_q(q);
        f_q = f_q.inv(k + 1);

        Fp conv = 0;
        for (int i = 0; i <= k; ++i) conv += p[i] * f_q[k - i];

        return inv2_p * conv;
    }

    static Fp kth_term(const vector<Fp> &a, const vector<Fp> &c, long long k) {
        const int d = a.size();
        vector<Fp> q(d + 1);
        q[0] = 1;
        for (int i = 0; i < d; ++i) q[i + 1] = -c[i];
        vector<Fp> p = atcoder::convolution(a, q);
        p.resize(d);
        return coeff(p, q, k);
    }

    static vector<Fp> middle_product_naive(vector<Fp> a, vector<Fp> c, bool c_reversed = false, bool b_reversed = false) {
        if (c_reversed) reverse(c.begin(), c.end());
        int n = a.size(), m = c.size() + 1 - n;
        vector<Fp> b(m);
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < m; ++j) {
                b[j] += a[i] * c[i + j];
            }
        }
        if (b_reversed) reverse(b.begin(), b.end());
        return b;
    }

    static vector<Fp> middle_product_fft(vector<Fp> a, vector<Fp> c, bool c_reversed = false, bool b_reversed = false) {
        int n = a.size(), m = c.size() + 1 - n;
        if (!c_reversed) reverse(c.begin(), c.end());

        int z = atcoder::internal::bit_ceil((unsigned int)(n + m));
        a.resize(z), c.resize(z);
        atcoder::internal::butterfly(a), atcoder::internal::butterfly(c);
        for (int i = 0; i < z; ++i) a[i] *= c[i];
        atcoder::internal::butterfly_inv(a);
        a.resize(n + m - 1);
        a.erase(a.begin(), a.begin() + n - 1);
        if (!b_reversed) reverse(a.begin(), a.end());
        const Fp iz = Fp::raw(z).inv();
        for (auto &e : a) e *= iz;
        return a;
    }
};

using FPS = FormalPowerSeries<Mod>;
struct Comp {
    bool operator()(const FPS &a, const FPS &b) const { return a.size() > b.size(); }
};
#line 6 "test/lc_kth_term_of_linearly_recurrent_sequence.cpp"

int main() {
    long long d, k;
    cin >> d >> k;
    vector<mint> a(d), c(d);
    cin >> a >> c;
    cout << FPS::kth_term(a, c, k) << endl;
}

Test cases

Env	Name	Status	Elapsed	Memory
g++	example_00	AC	5 ms	4 MB
g++	max_random_00	AC	510 ms	13 MB
g++	max_random_01	AC	505 ms	13 MB
g++	max_random_02	AC	503 ms	13 MB
g++	near_65536_00	AC	250 ms	8 MB
g++	near_65536_01	AC	490 ms	12 MB
g++	near_65536_02	AC	500 ms	12 MB
g++	random_00	AC	127 ms	6 MB
g++	random_01	AC	492 ms	12 MB
g++	random_02	AC	249 ms	8 MB
g++	small_00	AC	6 ms	4 MB
g++	small_01	AC	5 ms	4 MB
g++	small_02	AC	5 ms	4 MB
g++	small_03	AC	5 ms	4 MB
g++	small_04	AC	5 ms	4 MB
g++	small_05	AC	7 ms	4 MB
g++	small_06	AC	7 ms	4 MB
g++	small_07	AC	7 ms	4 MB
g++	small_08	AC	7 ms	4 MB
g++	small_09	AC	6 ms	4 MB