标签： 矩阵快速幂

$$ \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;
}

$$ \text{构造一个一元二次方程，使得}x_1=\frac{b+\sqrt{d}}{2},x_2=\frac{b-\sqrt{d}}{2} \\ \text{则原式}=\left( \frac{b+\sqrt{d}}{2} \right) ^n+\left( \frac{b-\sqrt{d}}{2} \right) ^n-\left( \frac{b-\sqrt{d}}{2} \right) ^n=x_{1}^{n}+x_{2}^{n}-x_{2}^{n} \\ \text{设}f\left( n \right) =x_{1}^{n}+x_{2}^{n} \\ f\left( n \right) =\left( x_1+x_2 \right) \left( x_{1}^{n-1}+x_{2}^{n-1} \right) -x_1x_2\left( x_{1}^{n-2}+x_{2}^{n-2} \right) \\ =\left( x_1+x_2 \right) f\left( n-1 \right) -x_1x_2f\left( n-2 \right) \\ x_1+x_2=b \\ x_1x_2=\frac{b^2-d}{4}\left( \text{韦达定理} \right) \\ f\left( n \right) =bf\left( n-1 \right) -\frac{b^2-d}{4}f\left( n-2 \right) \\ \text{这是一个常系数线性递推。} \\ \text{可以构造矩阵} \\ \left[ \begin{matrix} f\left( n \right)& f\left( n+1 \right)\\ \end{matrix} \right] \left[ \begin{matrix} 0& -\frac{b^2-d}{4}\\ 1& b\\ \end{matrix} \right] =\left[ \begin{matrix} f\left( n+1 \right)& f\left( n+2 \right)\\ \end{matrix} \right] \\ \text{边界条件:} \\ f\left( 0 \right) =x_{1}^{0}+x_{2}^{0}=2 \\ f\left( 1 \right) =x_1+x_2=b \\ \text{考虑}\left( \frac{b-\sqrt{d}}{2} \right) ^n\text{如何处理} \\ \text{设}g\left( n \right) =\left( \frac{b-\sqrt{d}}{2} \right) ^n \\ n=\text{0时，}g\left( n \right) =1 \\ \because b\le \sqrt{d} \\ \text{考虑}b\ne \sqrt{d}\text{的情况} \\ \therefore \frac{b-\sqrt{d}}{2}<0 \\ \because b+1>\sqrt{d} \\ \therefore b-\sqrt{d}>-1 \\ \therefore -\frac{1}{2}<\frac{b-\sqrt{d}}{2}<0 \\ n\text{为奇数时，}-\frac{1}{2}<g\left( n \right) <0 \\ n\text{为偶数时}0<g\left( n \right) <\frac{1}{2} \\ $$

#include <iostream>
#include <cstring>
#include <cinttypes>
using int_t = __int128;

using std::cin;
using std::cout;
using std::endl;

const int_t mod = 7528443412579576937ll;
const int_t inv4 = 5646332559434682703ll;
const int_t inv2 = 3764221706289788469ll;

struct Matrix
{
    int_t data[21][21];
    int_t n;
    Matrix(int_t n)
    {
        this->n = n;
        memset(data, 0, sizeof(data));
    }
    int_t &operator()(int_t r, int_t c)
    {
        return data[r][c];
    }
    int_t *operator[](int_t row)
    {
        return data[row];
    }
    Matrix operator*(Matrix another)
    {
        Matrix result(n);
        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[i][j] = (result[i][j] + another[k][j] * (*this)[i][k] % mod) % mod;
                }
            }
        }
        return result;
    }
    Matrix operator^(int_t index)
    {
        Matrix base = *this;
        Matrix result(n);
        for (int_t i = 1; i <= n; i++)
            result[i][i] = 1;
        while (index)
        {
            if (index & 1)
                result = result * base;
            base = base * base;
            index >>= 1;
        }
        return result;
    }
};

int main()
{
    int64_t b, d, n;
    cin >> b >> d >> n;
    if (n == 0)
    {
        cout << 1 << endl;
        return 0;
    }
    Matrix init(2);
    init(1, 1) = 2;
    init(1, 2) = b % mod;
    Matrix trans(2);
    trans(1, 2) = (mod - ((b * b % mod - d + mod) % mod) * inv4 % mod);
    trans(2, 2) = b;
    trans(2, 1) = 1;
    int64_t result = (init * (trans ^ n))(1, 1);
    if ((int_t)b * b != (int_t)d && n % 2 == 0)
    {
        result = (result - 1 + mod) % mod;
    }
    cout << result << endl;
    return 0;
}