#979079 (C++23) No.1559 Next Rational

• 問題: No.1559 Next Rational
• ユーザー: 👑 rin204
• 提出日時: 2024-05-05 16:38:41
• 言語: C++23 (gcc 12.3.0 + boost 1.83.0)
• 結果: AC
• 実行時間: 2 ms / 2,000 ms
• コード長: 10,447 bytes
• コンパイル時間: 3,873 ms
• コンパイル使用メモリ: 260,300 KB
• 実行使用メモリ: 6,820 KB
• 最終ジャッジ日時: 2024-11-27 15:44:14
• 合計ジャッジ時間: 4,460 ms
• ジャッジサーバーID （参考情報）: judge1 / judge2

入力	結果	実行時間; 実行使用メモリ
testcase_00	AC	2 ms; 6,820 KB
testcase_01	AC	2 ms; 6,820 KB
testcase_02	AC	2 ms; 6,816 KB
testcase_03	AC	2 ms; 6,820 KB
testcase_04	AC	2 ms; 6,816 KB
testcase_05	AC	2 ms; 6,816 KB
testcase_06	AC	2 ms; 6,820 KB
testcase_07	AC	2 ms; 6,816 KB
testcase_08	AC	2 ms; 6,820 KB
testcase_09	AC	2 ms; 6,816 KB
testcase_10	AC	2 ms; 6,820 KB
testcase_11	AC	2 ms; 6,820 KB
testcase_12	AC	2 ms; 6,816 KB
testcase_13	AC	2 ms; 6,820 KB
testcase_14	AC	2 ms; 6,816 KB
testcase_15	AC	2 ms; 6,820 KB
testcase_16	AC	2 ms; 6,816 KB
testcase_17	AC	2 ms; 6,816 KB

ソースコード

#define PROBLEM "https://yukicoder.me/problems/no/1559"

#include <bits/stdc++.h>
using namespace std;

template <int MOD>
struct Modint {
    int x;
    Modint() : x(0) {}
    Modint(int64_t y) {
        if (y >= 0)
            x = y % MOD;
        else
            x = (y % MOD + MOD) % MOD;
    }

    Modint &operator+=(const Modint &p) {
        x += p.x;
        if (x >= MOD) x -= MOD;
        return *this;
    }

    Modint &operator-=(const Modint &p) {
        x -= p.x;
        if (x < 0) x += MOD;
        return *this;
    }

    Modint &operator*=(const Modint &p) {
        x = int(1LL * x * p.x % MOD);
        return *this;
    }

    Modint &operator/=(const Modint &p) {
        *this *= p.inverse();
        return *this;
    }

    Modint &operator%=(const Modint &p) {
        assert(p.x == 0);
        return *this;
    }

    Modint operator-() const {
        return Modint(-x);
    }

    Modint &operator++() {
        x++;
        if (x == MOD) x = 0;
        return *this;
    }

    Modint &operator--() {
        if (x == 0) x = MOD;
        x--;
        return *this;
    }

    Modint operator++(int) {
        Modint result = *this;
        ++*this;
        return result;
    }

    Modint operator--(int) {
        Modint result = *this;
        --*this;
        return result;
    }

    friend Modint operator+(const Modint &lhs, const Modint &rhs) {
        return Modint(lhs) += rhs;
    }

    friend Modint operator-(const Modint &lhs, const Modint &rhs) {
        return Modint(lhs) -= rhs;
    }

    friend Modint operator*(const Modint &lhs, const Modint &rhs) {
        return Modint(lhs) *= rhs;
    }

    friend Modint operator/(const Modint &lhs, const Modint &rhs) {
        return Modint(lhs) /= rhs;
    }

    friend Modint operator%(const Modint &lhs, const Modint &rhs) {
        assert(rhs.x == 0);
        return Modint(lhs);
    }

    bool operator==(const Modint &p) const {
        return x == p.x;
    }

    bool operator!=(const Modint &p) const {
        return x != p.x;
    }

    bool operator<(const Modint &rhs) const {
        return x < rhs.x;
    }

    bool operator<=(const Modint &rhs) const {
        return x <= rhs.x;
    }

    bool operator>(const Modint &rhs) const {
        return x > rhs.x;
    }

    bool operator>=(const Modint &rhs) const {
        return x >= rhs.x;
    }

    Modint inverse() const {
        int a = x, b = MOD, u = 1, v = 0, t;
        while (b > 0) {
            t = a / b;
            a -= t * b;
            u -= t * v;
            std::swap(a, b);
            std::swap(u, v);
        }
        return Modint(u);
    }

    Modint pow(int64_t k) const {
        Modint ret(1);
        Modint y(x);
        while (k > 0) {
            if (k & 1) ret *= y;
            y *= y;
            k >>= 1;
        }
        return ret;
    }

    friend std::ostream &operator<<(std::ostream &os, const Modint &p) {
        return os << p.x;
    }

    friend std::istream &operator>>(std::istream &is, Modint &p) {
        int64_t y;
        is >> y;
        p = Modint<MOD>(y);
        return (is);
    }

    static int get_mod() {
        return MOD;
    }
};

struct Arbitrary_Modint {
    int x;
    static int MOD;

    static void set_mod(int mod) {
        MOD = mod;
    }

    Arbitrary_Modint() : x(0) {}
    Arbitrary_Modint(int64_t y) {
        if (y >= 0)
            x = y % MOD;
        else
            x = (y % MOD + MOD) % MOD;
    }

    Arbitrary_Modint &operator+=(const Arbitrary_Modint &p) {
        x += p.x;
        if (x >= MOD) x -= MOD;
        return *this;
    }

    Arbitrary_Modint &operator-=(const Arbitrary_Modint &p) {
        x -= p.x;
        if (x < 0) x += MOD;
        return *this;
    }

    Arbitrary_Modint &operator*=(const Arbitrary_Modint &p) {
        x = int(1LL * x * p.x % MOD);
        return *this;
    }

    Arbitrary_Modint &operator/=(const Arbitrary_Modint &p) {
        *this *= p.inverse();
        return *this;
    }

    Arbitrary_Modint &operator%=(const Arbitrary_Modint &p) {
        assert(p.x == 0);
        return *this;
    }

    Arbitrary_Modint operator-() const {
        return Arbitrary_Modint(-x);
    }

    Arbitrary_Modint &operator++() {
        x++;
        if (x == MOD) x = 0;
        return *this;
    }

    Arbitrary_Modint &operator--() {
        if (x == 0) x = MOD;
        x--;
        return *this;
    }

    Arbitrary_Modint operator++(int) {
        Arbitrary_Modint result = *this;
        ++*this;
        return result;
    }

    Arbitrary_Modint operator--(int) {
        Arbitrary_Modint result = *this;
        --*this;
        return result;
    }

    friend Arbitrary_Modint operator+(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) {
        return Arbitrary_Modint(lhs) += rhs;
    }

    friend Arbitrary_Modint operator-(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) {
        return Arbitrary_Modint(lhs) -= rhs;
    }

    friend Arbitrary_Modint operator*(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) {
        return Arbitrary_Modint(lhs) *= rhs;
    }

    friend Arbitrary_Modint operator/(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) {
        return Arbitrary_Modint(lhs) /= rhs;
    }

    friend Arbitrary_Modint operator%(const Arbitrary_Modint &lhs, const Arbitrary_Modint &rhs) {
        assert(rhs.x == 0);
        return Arbitrary_Modint(lhs);
    }

    bool operator==(const Arbitrary_Modint &p) const {
        return x == p.x;
    }

    bool operator!=(const Arbitrary_Modint &p) const {
        return x != p.x;
    }

    bool operator<(const Arbitrary_Modint &rhs) {
        return x < rhs.x;
    }

    bool operator<=(const Arbitrary_Modint &rhs) {
        return x <= rhs.x;
    }

    bool operator>(const Arbitrary_Modint &rhs) {
        return x > rhs.x;
    }

    bool operator>=(const Arbitrary_Modint &rhs) {
        return x >= rhs.x;
    }

    Arbitrary_Modint inverse() const {
        int a = x, b = MOD, u = 1, v = 0, t;
        while (b > 0) {
            t = a / b;
            a -= t * b;
            u -= t * v;
            std::swap(a, b);
            std::swap(u, v);
        }
        return Arbitrary_Modint(u);
    }

    Arbitrary_Modint pow(int64_t k) const {
        Arbitrary_Modint ret(1);
        Arbitrary_Modint y(x);
        while (k > 0) {
            if (k & 1) ret *= y;
            y *= y;
            k >>= 1;
        }
        return ret;
    }

    friend std::ostream &operator<<(std::ostream &os, const Arbitrary_Modint &p) {
        return os << p.x;
    }

    friend std::istream &operator>>(std::istream &is, Arbitrary_Modint &p) {
        int64_t y;
        is >> y;
        p = Arbitrary_Modint(y);
        return (is);
    }

    static int get_mod() {
        return MOD;
    }
};
int Arbitrary_Modint::MOD = 998244353;

using modint9 = Modint<998244353>;
using modint1 = Modint<1000000007>;
using modint  = Arbitrary_Modint;

template <typename T>
T modinv(T a, T MOD) {
    T b = MOD;
    T u = 1;
    T v = 0;
    while (b > 0) {
        T t = a / b;
        a -= t * b;
        u -= t * v;
        std::swap(a, b);
        std::swap(u, v);
    }
    if (a != 1) return -1;
    if (u < 0) u += MOD;
    return u;
}

template <typename T>
std::vector<T> berlekampMessy(const std::vector<T> &A) {
    int n = A.size();
    std::vector<T> B(1, -1);
    std::vector<T> C(1, -1);
    T y = 1;
    for (int j = 1; j <= n; j++) {
        int l = C.size();
        int m = B.size();
        T x   = 0;
        for (int i = 0; i < l; i++) {
            x += C[i] * A[j - l + i];
        }
        B.push_back(0);
        m++;
        if (x == 0) continue;
        T freq = x / y;
        if (l < m) {
            std::vector<T> D(m - l, T(0));
            D.insert(D.end(), C.begin(), C.end());
            for (int i = 0; i < m; i++) {
                D[m - 1 - i] -= freq * B[m - 1 - i];
            }
            std::swap(B, C);
            std::swap(C, D);
            y = x;
        } else {
            for (int i = 0; i < m; i++) {
                C[l - 1 - i] -= freq * B[m - 1 - i];
            }
        }
    }

    std::reverse(C.begin(), C.end());
    for (auto &c : C) c = -c;
    return C;
}

/*
return x^n mod f
引数: f_reversed: f(x)の係数(逆順, f_reversed[0] = 1)
*/
template <typename T>
std::vector<T> monomial_mod_polynomial(long long n, const std::vector<T> &f_reversed) {
    assert(!f_reversed.empty() and f_reversed[0] == 1);
    int K = f_reversed.size() - 1;
    if (!K) return {};
    int D = 64 - __builtin_clzll(n);
    std::vector<T> ret(K, 0);
    ret[0]         = 1;
    auto self_conv = [](std::vector<T> &x) -> std::vector<T> {
        int d = x.size();
        std::vector<T> ret(2 * d - 1);
        for (int i = 0; i < d; i++) {
            ret[2 * i] += x[i] * x[i];
            for (int j = 0; j < i; j++) ret[i + j] += 2 * x[i] * x[j];
        }
        return ret;
    };
    for (int d = D - 1; d >= 0; d--) {
        ret = self_conv(ret);
        for (int i = 2 * K - 2; i >= K; i--) {
            for (int j = 1; j <= K; j++) {
                ret[i - j] -= ret[i] * f_reversed[j];
            }
        }
        ret.resize(K);
        if ((n >> d) & 1) {
            std::vector<T> c(K);
            c[0] = -ret[K - 1] * f_reversed[K];
            for (int i = 1; i < K; i++) {
                c[i] = ret[i - 1] - ret[K - 1] * f_reversed[K - i];
            }
            ret = c;
        }
    }
    return ret;
}

template <typename T>
T guess_kth_term(const std::vector<T> &A, long long k, bool debug = false) {
    assert(k >= 0);
    if (k < int(A.size())) return A[k];

    const auto F = berlekampMessy(A);
    if (debug) {
        std::cerr << "F.size() = " << F.size() << std::endl;
    }
    const auto G = monomial_mod_polynomial<T>(k, F);
    T ret        = 0;
    for (size_t i = 0; i < G.size(); i++) ret += G[i] * A[i];
    return ret;
}
using mint = modint1;

void solve() {
    long long n, a, b, k;
    cin >> n >> a >> b >> k;
    const int N = 1000;
    vector<mint> A(N);
    A[0] = a;
    A[1] = b;
    for (int i = 2; i < N; i++) {
        A[i] = (A[i - 1] * A[i - 1] + k) / A[i - 2];
    }
    auto ans = guess_kth_term(A, n - 1, true);
    cout << ans << endl;
}

int main() {
    cin.tie(0)->sync_with_stdio(0);
    cout << fixed << setprecision(12);
    int t;
    t = 1;
    // cin >> t;
    while (t--) solve();
    return 0;
}