library

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

---

Project maintained by tko919 Hosted on GitHub Pages — Theme by mattgraham

Pisano Period
(Math/pisano.hpp)

• View this file on GitHub
• Last update: 2025-05-25 16:11:40+09:00
• Include: #include "Math/pisano.hpp"

Depends on

• Dynamic Modint (Math/dynamic.hpp)
• Fast Division (Math/fastdiv.hpp)
• Matrix (Math/matrix.hpp)
• Miller-Rabin (Math/miller.hpp)
• Pollard-Rho (Math/pollard.hpp)
• Random (Utility/random.hpp)

Code

#pragma once
#include "Math/dynamic.hpp"
#include "Math/matrix.hpp"
#include "Math/pollard.hpp"

int Pisano(int p) {
    if (p == 2)
        return 3;
    if (p == 5)
        return 20;
    vector<ll> cand;
    if (p % 5 == 1 or p % 5 == 4)
        cand = EnumDivisors(p - 1);
    else
        cand = EnumDivisors(p * 2 + 2);
    auto check = [&](int d, int p) -> bool {
        Fp::set_mod(p);
        Matrix<Fp> mat(2, 2);
        mat[0][1] = mat[1][0] = mat[1][1] = 1;
        mat = mat.pow(d);
        return mat[0][1] == 0 and mat[1][1] == 1;
    };
    for (auto &d : cand)
        if (check(d, p)) {
            return d;
        }
    assert(0);
}

/**
 * @brief Pisano Period
 */

#line 2 "Math/fastdiv.hpp"

struct FastDiv {
    using u64 = uint64_t;
    using u128 = __uint128_t;
    constexpr FastDiv() : m(), s(), x() {}
    constexpr FastDiv(int _m)
        : m(_m), s(__lg(m - 1)), x(((u128(1) << (s + 64)) + m - 1) / m) {}
    constexpr int get() {
        return m;
    }
    constexpr friend u64 operator/(u64 n, const FastDiv &d) {
        return (u128(n) * d.x >> d.s) >> 64;
    }
    constexpr friend int operator%(u64 n, const FastDiv &d) {
        return n - n / d * d.m;
    }
    constexpr pair<u64, int> divmod(u64 n) const {
        u64 q = n / (*this);
        return {q, n - q * m};
    }
    int m, s;
    u64 x;
};

struct FastDiv64 {
    using u64 = uint64_t;
    using u128 = __uint128_t;
    u128 mod, mh, ml;
    explicit FastDiv64(u64 mod = 1) : mod(mod) {
        u128 m = u128(-1) / mod;
        if (m * mod + mod == u128(0))
            ++m;
        mh = m >> 64;
        ml = m & u64(-1);
    }
    u64 umod() const {
        return mod;
    }
    u64 modulo(u128 x) {
        u128 z = (x & u64(-1)) * ml;
        z = (x & u64(-1)) * mh + (x >> 64) * ml + (z >> 64);
        z = (x >> 64) * mh + (z >> 64);
        x -= z * mod;
        return x < mod ? x : x - mod;
    }
    u64 mul(u64 a, u64 b) {
        return modulo(u128(a) * b);
    }
};

/**
 * @brief Fast Division
 */
#line 3 "Math/dynamic.hpp"

struct Fp {
    using u64 = uint64_t;
    uint v;
    static uint get_mod() {
        return _getmod();
    }
    static void set_mod(uint _m) {
        bar = FastDiv(_m);
    }
    Fp inv() const {
        int tmp, a = v, b = get_mod(), x = 1, y = 0;
        while (b) {
            tmp = a / b, a -= tmp * b;
            swap(a, b);
            x -= tmp * y;
            swap(x, y);
        }
        if (x < 0) {
            x += get_mod();
        }
        return x;
    }
    Fp() : v(0) {}
    Fp(ll x) {
        v = x % get_mod();
        if (v < 0)
            v += get_mod();
    }
    Fp operator-() const {
        return Fp() - *this;
    }
    Fp pow(ll t) {
        assert(t >= 0);
        Fp res = 1, b = *this;
        while (t) {
            if (t & 1)
                res *= b;
            b *= b;
            t >>= 1;
        }
        return res;
    }
    Fp &operator+=(const Fp &x) {
        v += x.v;
        if (v >= get_mod())
            v -= get_mod();
        return *this;
    }
    Fp &operator-=(const Fp &x) {
        v += get_mod() - x.v;
        if (v >= get_mod())
            v -= get_mod();
        return *this;
    }
    Fp &operator*=(const Fp &x) {
        v = (u64(v) * x.v) % bar;
        return *this;
    }
    Fp &operator/=(const Fp &x) {
        (*this) *= x.inv();
        return *this;
    }
    Fp operator+(const Fp &x) const {
        return Fp(*this) += x;
    }
    Fp operator-(const Fp &x) const {
        return Fp(*this) -= x;
    }
    Fp operator*(const Fp &x) const {
        return Fp(*this) *= x;
    }
    Fp operator/(const Fp &x) const {
        return Fp(*this) /= x;
    }
    bool operator==(const Fp &x) const {
        return v == x.v;
    }
    bool operator!=(const Fp &x) const {
        return v != x.v;
    }
    friend istream &operator>>(istream &is, Fp &x) {
        return is >> x.v;
    }
    friend ostream &operator<<(ostream &os, const Fp &x) {
        return os << x.v;
    }

  private:
    static FastDiv bar;
    static uint _getmod() {
        return bar.get();
    }
};
FastDiv Fp::bar(998244353);

void rd(Fp &x) {
    fastio::rd(x.v);
}
void wt(Fp x) {
    fastio::wt(x.v);
}

/**
 * @brief Dynamic Modint
 */
#line 2 "Math/matrix.hpp"

template <class T> struct Matrix {
    int h, w;
    vector<vector<T>> val;
    T det;
    Matrix() {}
    Matrix(int n) : h(n), w(n), val(vector<vector<T>>(n, vector<T>(n))) {}
    Matrix(int n, int m)
        : h(n), w(m), val(vector<vector<T>>(n, vector<T>(m))) {}
    vector<T> &operator[](const int i) {
        return val[i];
    }
    Matrix &operator+=(const Matrix &m) {
        assert(h == m.h and w == m.w);
        rep(i, 0, h) rep(j, 0, w) val[i][j] += m.val[i][j];
        return *this;
    }
    Matrix &operator-=(const Matrix &m) {
        assert(h == m.h and w == m.w);
        rep(i, 0, h) rep(j, 0, w) val[i][j] -= m.val[i][j];
        return *this;
    }
    Matrix &operator*=(const Matrix &m) {
        assert(w == m.h);
        Matrix<T> res(h, m.w);
        rep(i, 0, h) rep(j, 0, m.w) rep(k, 0, w) res.val[i][j] +=
            val[i][k] * m.val[k][j];
        *this = res;
        return *this;
    }
    Matrix operator+(const Matrix &m) const {
        return Matrix(*this) += m;
    }
    Matrix operator-(const Matrix &m) const {
        return Matrix(*this) -= m;
    }
    Matrix operator*(const Matrix &m) const {
        return Matrix(*this) *= m;
    }
    Matrix pow(ll k) {
        Matrix<T> res(h, h), c = *this;
        rep(i, 0, h) res.val[i][i] = 1;
        while (k) {
            if (k & 1)
                res *= c;
            c *= c;
            k >>= 1;
        }
        return res;
    }
    vector<int> gauss(int c = -1) {
        det = 1;
        if (val.empty())
            return {};
        if (c == -1)
            c = w;
        int cur = 0;
        vector<int> res;
        rep(i, 0, c) {
            if (cur == h)
                break;
            rep(j, cur, h) if (val[j][i] != 0) {
                swap(val[cur], val[j]);
                if (cur != j)
                    det *= -1;
                break;
            }
            det *= val[cur][i];
            if (val[cur][i] == 0)
                continue;
            rep(j, 0, h) if (j != cur) {
                T z = val[j][i] / val[cur][i];
                rep(k, i, w) val[j][k] -= val[cur][k] * z;
            }
            res.push_back(i);
            cur++;
        }
        return res;
    }
    Matrix inv() {
        assert(h == w);
        Matrix base(h, h * 2), res(h, h);
        rep(i, 0, h) rep(j, 0, h) base[i][j] = val[i][j];
        rep(i, 0, h) base[i][h + i] = 1;
        base.gauss(h);
        det = base.det;
        rep(i, 0, h) rep(j, 0, h) res[i][j] = base[i][h + j] / base[i][i];
        return res;
    }
    bool operator==(const Matrix &m) {
        assert(h == m.h and w == m.w);
        rep(i, 0, h) rep(j, 0, w) if (val[i][j] != m.val[i][j]) return false;
        return true;
    }
    bool operator!=(const Matrix &m) {
        assert(h == m.h and w == m.w);
        rep(i, 0, h) rep(j, 0, w) if (val[i][j] == m.val[i][j]) return false;
        return true;
    }
    friend istream &operator>>(istream &is, Matrix &m) {
        rep(i, 0, m.h) rep(j, 0, m.w) is >> m[i][j];
        return is;
    }
    friend ostream &operator<<(ostream &os, Matrix &m) {
        rep(i, 0, m.h) {
            rep(j, 0, m.w) os << m[i][j]
                              << (j == m.w - 1 and i != m.h - 1 ? '\n' : ' ');
        }
        return os;
    }
};

/**
 * @brief Matrix
 */
#line 2 "Math/miller.hpp"

struct m64 {
    using i64 = int64_t;
    using u64 = uint64_t;
    using u128 = __uint128_t;

    static u64 mod;
    static u64 r;
    static u64 n2;

    static u64 get_r() {
        u64 ret = mod;
        rep(_,0,5) ret *= 2 - mod * ret;
        return ret;
    }

    static void set_mod(u64 m) {
        assert(m < (1LL << 62));
        assert((m & 1) == 1);
        mod = m;
        n2 = -u128(m) % m;
        r = get_r();
        assert(r * mod == 1);
    }
    static u64 get_mod() { return mod; }

    u64 a;
    m64() : a(0) {}
    m64(const int64_t &b) : a(reduce((u128(b) + mod) * n2)){};

    static u64 reduce(const u128 &b) {
        return (b + u128(u64(b) * u64(-r)) * mod) >> 64;
    }
    u64 get() const {
        u64 ret = reduce(a);
        return ret >= mod ? ret - mod : ret;
    }
    m64 &operator*=(const m64 &b) {
        a = reduce(u128(a) * b.a);
        return *this;
    }
    m64 operator*(const m64 &b) const { return m64(*this) *= b; }
    bool operator==(const m64 &b) const {
        return (a >= mod ? a - mod : a) == (b.a >= mod ? b.a - mod : b.a);
    }
    bool operator!=(const m64 &b) const {
        return (a >= mod ? a - mod : a) != (b.a >= mod ? b.a - mod : b.a);
    }
    m64 pow(u128 n) const {
        m64 ret(1), mul(*this);
        while (n > 0) {
        if (n & 1) ret *= mul;
        mul *= mul;
        n >>= 1;
        }
        return ret;
    }
};
typename m64::u64 m64::mod, m64::r, m64::n2;

bool Miller(ll n){
    if(n<2 or (n&1)==0)return (n==2);
    m64::set_mod(n);
    ll d=n-1; while((d&1)==0)d>>=1;
    vector<ll> seeds;
    if(n<(1<<30))seeds={2, 7, 61};
    else seeds={2, 325, 9375, 28178, 450775, 9780504};
    for(auto& x:seeds){
        if(n<=x)break;
        ll t=d;
        m64 y=m64(x).pow(t);
        while(t!=n-1 and y!=1 and y!=n-1){
            y*=y;
            t<<=1;
        }
        if(y!=n-1 and (t&1)==0)return 0;
    } return 1;
}

/**
 * @brief Miller-Rabin
 */
#line 2 "Utility/random.hpp"

namespace Random {
mt19937_64 randgen(chrono::steady_clock::now().time_since_epoch().count());
using u64 = unsigned long long;
u64 get() {
    return randgen();
}
template <typename T> T get(T L) { // [0,L]

    return get() % (L + 1);
}
template <typename T> T get(T L, T R) { // [L,R]

    return get(R - L) + L;
}
double uniform() {
    return double(get(1000000000)) / 1000000000;
}
string str(int n) {
    string ret;
    rep(i, 0, n) ret += get('a', 'z');
    return ret;
}
template <typename Iter> void shuffle(Iter first, Iter last) {
    if (first == last)
        return;
    int len = 1;
    for (auto it = first + 1; it != last; it++) {
        len++;
        int j = get(0, len - 1);
        if (j != len - 1)
            iter_swap(it, first + j);
    }
}
template <typename T> vector<T> select(int n, T L, T R) { // [L,R]

    if (n * 2 >= R - L + 1) {
        vector<T> ret(R - L + 1);
        iota(ALL(ret), L);
        shuffle(ALL(ret));
        ret.resize(n);
        return ret;
    } else {
        unordered_set<T> used;
        vector<T> ret;
        while (SZ(used) < n) {
            T x = get(L, R);
            if (!used.count(x)) {
                used.insert(x);
                ret.push_back(x);
            }
        }
        return ret;
    }
}

void relabel(int n, vector<pair<int, int>> &es) {
    shuffle(ALL(es));
    vector<int> ord(n);
    iota(ALL(ord), 0);
    shuffle(ALL(ord));
    for (auto &[u, v] : es)
        u = ord[u], v = ord[v];
}
template <bool directed, bool multi, bool self>
vector<pair<int, int>> genGraph(int n, int m) {
    vector<pair<int, int>> cand, es;
    rep(u, 0, n) rep(v, 0, n) {
        if (!self and u == v)
            continue;
        if (!directed and u > v)
            continue;
        cand.push_back({u, v});
    }
    if (m == -1)
        m = get(SZ(cand));
    // chmin(m, SZ(cand));

    vector<int> ord;
    if (multi)
        rep(_, 0, m) ord.push_back(get(SZ(cand) - 1));
    else {
        ord = select(m, 0, SZ(cand) - 1);
    }
    for (auto &i : ord)
        es.push_back(cand[i]);
    relabel(n, es);
    return es;
}
vector<pair<int, int>> genTree(int n) {
    vector<pair<int, int>> es;
    rep(i, 1, n) es.push_back({get(i - 1), i});
    relabel(n, es);
    return es;
}
}; // namespace Random


/**
 * @brief Random
 */
#line 4 "Math/pollard.hpp"

vector<ll> Pollard(ll n) {
    if (n <= 1)
        return {};
    if (Miller(n))
        return {n};
    if ((n & 1) == 0) {
        vector<ll> v = Pollard(n >> 1);
        v.push_back(2);
        return v;
    }
    for (ll x = 2, y = 2, d;;) {
        ll c = Random::get(2LL, n - 1);
        do {
            x = (__int128_t(x) * x + c) % n;
            y = (__int128_t(y) * y + c) % n;
            y = (__int128_t(y) * y + c) % n;
            d = __gcd(x - y + n, n);
        } while (d == 1);
        if (d < n) {
            vector<ll> lb = Pollard(d), rb = Pollard(n / d);
            lb.insert(lb.end(), ALL(rb));
            return lb;
        }
    }
}

vector<pair<ll, int>> Pollard2(ll n) {
    auto ps = Pollard(n);
    sort(ALL(ps));
    using P = pair<ll, int>;
    vector<P> pes;
    for (auto &p : ps) {
        if (pes.empty() or pes.back().first != p) {
            pes.push_back({p, 1});
        } else {
            pes.back().second++;
        }
    }
    return pes;
}

vector<ll> EnumDivisors(ll n) {
    auto pes = Pollard2(n);
    vector<ll> ret;
    auto rec = [&](auto &rec, int id, ll d) -> void {
        if (id == SZ(pes)) {
            ret.push_back(d);
            return;
        }
        rec(rec, id + 1, d);
        rep(e, 0, pes[id].second) {
            d *= pes[id].first;
            rec(rec, id + 1, d);
        }
    };
    rec(rec, 0, 1);
    sort(ALL(ret));
    return ret;
}

/**
 * @brief Pollard-Rho
 */
#line 5 "Math/pisano.hpp"

int Pisano(int p) {
    if (p == 2)
        return 3;
    if (p == 5)
        return 20;
    vector<ll> cand;
    if (p % 5 == 1 or p % 5 == 4)
        cand = EnumDivisors(p - 1);
    else
        cand = EnumDivisors(p * 2 + 2);
    auto check = [&](int d, int p) -> bool {
        Fp::set_mod(p);
        Matrix<Fp> mat(2, 2);
        mat[0][1] = mat[1][0] = mat[1][1] = 1;
        mat = mat.pow(d);
        return mat[0][1] == 0 and mat[1][1] == 1;
    };
    for (auto &d : cand)
        if (check(d, p)) {
            return d;
        }
    assert(0);
}

/**
 * @brief Pisano Period
 */

Back to top page