ilovegzd

首先你需要知道 矩阵运算。

过程

让我们看看数据范围：$0<b^2\leq d<(b+1)^2\leq 10^{18},n\leq 10^{18}$，并且 $b\bmod 2=1,d\bmod 4=1$。

看到这个数据范围我们肯定要一个 $O(\log n)$ 左右的算法。

然后我们对两边取个根号，得到 $0<b\leq\sqrt{d}<b+1\leq 10^9$。

沿用 P5136 的解法，如果我们设一个数列

$$F_n=(\dfrac{b+\sqrt{d}}{2})^n+(\dfrac{b-\sqrt{d}}{2})^n $$

我们能在 $O(\log n)$ 的时间复杂度内使用矩阵快速幂计算 $F_n$。

计算 $(\frac{b-\sqrt{d}}{2})^n$

$$\begin{aligned} b&\leq\sqrt{d}<b+1\\ 0&\leq\sqrt{d}-b<1\\ -1&< b-\sqrt{d}\leq 0\\ -1<-\frac{1}{2}&<\frac{b-\sqrt{d}}{2}\leq 0 \end{aligned} $$

于是

$$(\frac{b-\sqrt{d}}{2})^n \left\{\begin{aligned} \leq 0,n\bmod 2=1\\ \geq 0,n\bmod 2=0 \end{aligned}\right. $$

$O(\log n)$ 计算 $F_n$

第二波式子：

$$\begin{aligned} F_n&=(\dfrac{b+\sqrt{d}}{2})^n+(\dfrac{b-\sqrt{d}}{2})^n\\ &=\dfrac{b^2+2b\sqrt{d}+d}{4}\times (\dfrac{b+\sqrt{d}}{2})^{n-2}+\dfrac{b^2-2b\sqrt{d}+d}{4}\times (\dfrac{b-\sqrt{d}}{2})^{n-2}\\ &=\dfrac{2b^2+2b\sqrt{d}+d-b^2}{4}\times (\dfrac{b+\sqrt{d}}{2})^{n-2}+\dfrac{2b^2-2b\sqrt{d}+d-b^2}{4}\times (\dfrac{b-\sqrt{d}}{2})^{n-2}\\ &=\dfrac{2b^2+2b\sqrt{d}}{4}\times (\dfrac{b+\sqrt{d}}{2})^{n-2}+\dfrac{d-b^2}{4}\times (\dfrac{b+\sqrt{d}}{2})^{n-2}+\dfrac{2b^2-2b\sqrt{d}}{4}\times (\dfrac{b-\sqrt{d}}{2})^{n-2}+\dfrac{d-b^2}{4}\times (\dfrac{b-\sqrt{d}}{2})^{n-2}\\ &=\dfrac{b(b+\sqrt{d})}{2}\times (\dfrac{b+\sqrt{d}}{2})^{n-2}+\dfrac{d-b^2}{4}\times (\dfrac{b+\sqrt{d}}{2})^{n-2}+\dfrac{b(b-\sqrt{d})}{2}\times (\dfrac{b-\sqrt{d}}{2})^{n-2}+\dfrac{d-b^2}{4}\times (\dfrac{b-\sqrt{d}}{2})^{n-2}\\ &=b((\dfrac{b+\sqrt{d}}{2})^{n-1}+(\dfrac{b-\sqrt{d}}{2})^{n-1})+\dfrac{d-b^2}{4}((\dfrac{b+\sqrt{d}}{2})^{n-2}+(\dfrac{b-\sqrt{d}}{2})^{n-2})\\ &=bF_{n-1}+\dfrac{d-b^2}{4}F_{n-2}\\ \end{aligned} $$

由于 $b\bmod 2=1,d\bmod 4=1$，得到 $\frac{d-b^2}{4}\in\mathbb{Z}$。

递推矩阵即为

$$\left\lbrack\begin{matrix} b&1\\ \frac{d-b^2}{4}&0\\ \end{matrix}\right\rbrack $$

初始矩阵即为

$$\left\lbrack\begin{matrix} F_2&F_1\\ 0&0\\ \end{matrix}\right\rbrack=\left\lbrack\begin{matrix} d+b^2&b\\ 0&0\\ \end{matrix}\right\rbrack $$

然后就可以矩阵快速幂了。

代码：

mat base, res;
constexpr mint one(1), zero(0);
mint ans;
unsigned long long n = 0;
mint b, d;
inline void solve() {
    cin >> b >> d >> n;
    if (n == 0) cout << '1' << '\n';
    else if (n == 1) cout << mint((unsigned long long)(b.val() + sqrt(d.val())) >> 1) << '\n';
    else {
        base.mat[1][1] = zero;
        base.mat[0][0] = b;
        base.mat[0][1] = one;
        base.mat[1][0] = (d - b * b).val() >> 2;
        ans = !(n & 1);
        ans = -ans;
        n -= 2;
        res.mat[0][0] = (b * b + d).val() >> 1;
        res.mat[1][0] = res.mat[1][1] = zero;
        res.mat[0][1] = b;
        for (; n; base *= base, n >>= 1) if (n & 1) res *= base;
        ans += res.mat[0][0];
        cout << ans << '\n';
    }
}
int main() {
    base.mat.resize(2);
    base.mat[0].resize(2);
    base.mat[1].resize(2);
    res.mat.resize(2);
    res.mat[0].resize(2);
    res.mat[1].resize(2);
    solve();
}