洛谷4000 斐波那契数列

$$ \text{求}F\left( n \right) \,\,mod\,\,p,n\le 10^{30000000},p\le 2^{31}-1 \\ \text{斐波那契数列在模}p\text{意义下的循环节长度为}O\left( p \right) \text{级别。} \\ \text{首先将}p\text{进行质因数分解，得到}p_{1}^{k_1}p_{2}^{k_2}p_{3}^{k_3}… \\ \text{设}T\left( k \right) \text{为模}k\text{意义下的循环节长度，则}T\left( p \right) =\prod{T\left( p_i^{k_i} \right)} \\ \text{同时对于}p^k,\text{有}T\left( p^k \right) =T\left( p \right) \times p^{k-1} \\ \text{考虑如何求出模质数的循环节长度。} \\ \text{对于2,3,5直接特判。} \\ \text{对于其他质数}p\text{，如果5是模}p\text{意义下的二次剩余}\left( \text{即}5^{\frac{p-1}{2}}\equiv 1\left( mod\,\,p \right) \right) ,\text{则模}p\text{意义下的循环节长度是}p-\text{1的约数。} \\ \text{部分证明：} \\ \text{斐波那契数列的通项公式为}F_n=\frac{\sqrt{5}}{5}\left( \left( \frac{\left( 1+\sqrt{5} \right)}{2} \right) ^n-\left( \frac{\left( 1-\sqrt{5} \right)}{2} \right) ^n \right) \\ \text{设5在模}p\text{意义下的二次剩余为}x\text{，则}F_n\text{一定可以表示为}k\left( x_0^n-x_1^n \right) \text{的形式。} \\ \text{其中}x_0^n\text{和}x_1^n\text{的周期一定是}\varphi \left( p \right) \text{的约数。} \\ \text{故}F_n\text{的周期是}\varphi \left( p \right) \text{的约数。} \\ \\ \text{反之，模}p\text{意义下的循环节长度是}2\left( p+1 \right) \text{的约数。} \\ \text{枚举约数验证即可。} \\ \text{这个我不会证} \\ \\ $$

#include <climits>
#include <cstring>
#include <iostream>
#include <map>
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;
struct Matrix mat_mul(const struct Matrix& a, const struct Matrix& b,
                      int_t mod);
int_t power(int_t base, int_t index, int_t mod);
struct Matrix {
    int_t data[4][4];
    Matrix() { memset(data, 0, sizeof(data)); }
    int_t& operator()(int_t r, int_t c) { return data[r][c]; }
    const int_t N = 2;
    Matrix& operator=(const Matrix& mat2) {
        memcpy(data, mat2.data, sizeof(data));
        return *this;
    }
    Matrix power(int_t index, int_t mod) {
        Matrix result;
        Matrix base = *this;
        result(1, 1) = result(2, 2) = 1;
        while (index) {
            if (index & 1) result = mat_mul(result, base, mod);
            base = mat_mul(base, base, mod);
            index >>= 1;
        }
        return result;
    }
};

Matrix mat_mul(const Matrix& a, const Matrix& b, int_t mod) {
    const auto N = a.N;
    Matrix result;
    for (int_t i = 1; i <= N; i++) {
        for (int_t j = 1; j <= N; j++) {
            for (int_t k = 1; k <= N; k++) {
                result.data[i][j] =
                    (result.data[i][j] + a.data[i][k] * b.data[k][j] % mod) %
                    mod;
            }
        }
    }
    return result;
}
Matrix F(int_t k, int_t mod) {
    Matrix mat;
    mat(1, 1) = 0;
    mat(1, 2) = 1;
    mat(2, 1) = 1;
    mat(2, 2) = 1;
    return mat.power(k, mod);
}
//计算模p意义下的循环周期，p是质数
int_t getPeriod(int_t p) {
    if (p == 2)
        return 3;
    else if (p == 3)
        return 8;
    else if (p == 5)
        return 20;
    int_t fac;
    if (power(5, (p - 1) / 2, p) == 1) {
        fac = p - 1;
    } else {
        fac = 2 * (p + 1);
    }
    const auto check = [&](int_t x) -> bool {
        Matrix res = ::F(x + 1, p);
        return res(1, 1) == 0 && res(1, 2) == 1;
    };
    int_t result = 1e12;
    for (int_t i = 1; i * i <= fac; i++) {
        if (fac % i == 0) {
            if (i != 1 && check(i)) {
                result = std::min(result, i);
            }
            if (fac / i != i && check(fac / i)) {
                result = std::min(result, fac / i);
            }
        }
    }
    return result;
}
std::map<int_t, int_t> divide(int_t x) {
    std::map<int_t, int_t> result;
    for (int_t i = 2; i * i <= x; i++) {
        if (x % i == 0) {
            result[i] = 0;
            while (x % i == 0) {
                x /= i;
                result[i]++;
            }
        }
    }
    if (x > 1) result[x] += 1;
    return result;
}
int main() {
    static char buf[30000000 + 10];
    scanf("%s", buf);
    int_t p;
    cin >> p;
    int_t length = 1;
    for (const auto& pair : divide(p)) {
        length *= ::getPeriod(pair.first) *
                  power(pair.first, pair.second - 1, int_t(1) << 40);
    }
    int_t sum = 0;
    for (auto ptr = buf; *ptr != '\0'; ptr++) {
        sum = (sum * 10 + *ptr - '0') % length;
    }
    cout << F(sum , p)(1, 2) << endl;
    return 0;
}

int_t power(int_t base, int_t index, int_t mod) {
    int_t result = 1;
    while (index) {
        if (index & 1) result = result * base % mod;
        index >>= 1;
        base = base * base % mod;
    }
    return result;
}