P3263 题解

首先，我们可以发现 \(x_1 = \dfrac{b + \sqrt{d}}{2}\) 和 \(x_2 = \dfrac{b - \sqrt{d}}{2}\) 是二次方程 \(x^2 - bx + \dfrac{b^2 - d}{4} = 0\) 的根。

则根据韦达定理：

\[\begin{cases} x_1 + x_2 = b \\ x_1 x_2 = \dfrac{b^2 - d}{4} \end{cases} \]

有了这个，就可以递推算 \(x_1^n + x_2^n\)：

\[\begin{aligned} x_1^n + x_2^n &= (x_1 + x_2)(x_1^{n - 1} + x_2^{n - 1}) - x_1 x_2(x_1^{n - 2} + x_2^{n - 2}) \\ &= b(x_1^{n - 1} + x_2^{n - 1}) - \dfrac{(b^2 - d)(x_1^{n - 2} + x_2^{n - 2})}{4} \end{aligned} \]

设 \(f_n = x_1^n + x_2^n\)，则：

\[\begin{cases} f_0 = 2 \\ f_1 = b \\ f_n = b \cdot f_{n - 1} - \dfrac{b^2 - d}{4} \cdot f_{n - 2} \end{cases} \]

于是我们可以使用矩阵快速幂计算 \(f_n\)。

最后考虑如何通过 \(f_n\) 计算答案，可以发现：

• 若 \(2 \nmid n\)，\(-1 < x_2^n \le 0\)。
• 否则：\(0 \le x_2^n < 1\)。

所以：

• 若 \(n = 0\)：答案为 \(1\)。
• 若 \(b^2 \neq d\) 且 \(2 \mid n\)：答案为 \(f_n - 1\)。
• 否则：答案为 \(f_n\)。

时间复杂度：\(O(\log n)\)。

#include <algorithm>
#include <iostream>
#include <ostream>
#include <vector>

using namespace std;
using i64 = long long;
using i128 = __int128_t;

constexpr i64 MOD = 7'528'443'412'579'576'937;

inline i128 mul(i128 a, i128 b)
{
    return a * b % MOD;
}

inline i128 add(i128 a, i128 b)
{
    return (a + b) % MOD;
}

ostream& operator<<(ostream& os, i128 x)
{
    if (x == 0)
        return os << x;

    string s;

    while (x > 0) {
        s += x % 10 + '0';
        x /= 10;
    }

    reverse(s.begin(), s.end());
    return os << s;
}

struct Matrix {
    vector<vector<i128>> data;

    Matrix(int n, int m)
        : data(n, vector(m, (i128)0))
    {
    }

    Matrix(const vector<vector<i128>>& val)
        : data(val)
    {
    }

    auto& operator[](int idx)
    {
        return data[idx];
    }

    const auto& operator[](int idx) const
    {
        return data[idx];
    }

    friend Matrix operator*(const Matrix& base, const Matrix& other)
    {
        Matrix res(base.data.size(), other[0].size());

        for (int i = 0; i < (int)base.data.size(); ++i) {
            for (int k = 0; k < (int)base[i].size(); ++k) {
                auto curr = base[i][k];

                for (int j = 0; j < (int)other[k].size(); ++j) {
                    res[i][j] = add(res[i][j], mul(curr, other[k][j]));
                }
            }
        }

        return res;
    }

    inline Matrix& operator*=(const Matrix& other) // left mul
    {
        return *this = *this * other;
    }

    Matrix operator^(i64 n) const
    {
        Matrix res(data.size(), data.size()), base(*this);

        for (int i = 0; i < (int)data.size(); ++i)
            res[i][i] = 1;

        for (; n > 0; n >>= 1, base *= base) {
            if (n & 1)
                res *= base;
        }

        return res;
    }

    inline Matrix& operator^=(i64 n)
    {
        return *this = *this ^ n;
    }
};

int main()
{
    i64 b, d, n;
    cin >> b >> d >> n;

    if (n == 0) {
        cout << "1\n";
        return 0;
    }

    Matrix F(1, 2), G(2, 2);

    F[0] = { b, 2 };
    G.data = {
        { b, 1 },
        { (d - b * b) / 4, 0 },
    };

    F *= G ^ (n - 1);
    cout << F[0][0] - (int)(b * b != d && (n & 1) == 0) << '\n';

    return 0;
}