JLOI2015 有意义的字符串

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