officeyutong

标签： 常系数线性递推

  • 常系数线性递推的新做法 – 计算[x^k]P(x)/Q(x)

    $$ a_n=\sum_{1\le i\le k}{f_ia_{n-i}} \\ \left\{ a_0,a_1….a_{k-1} \right\} \text{已知} \\ \\ \text{设}F\left( x \right) \text{表示从第}k\text{项开始的该数列的生成函数} \\ \sum_{i\ge k}{a_ix^i}=\sum_{i\ge k}{\sum_{1\le j\le k}{f_ja_{i-j}x^i}} \\ =\sum_{1\le j\le k}{f_j\sum_{i\ge k}{a_{i-j}x^i}} \\ =\sum_{1\le j\le k}{f_j\sum_{i\ge k-j}{a_ix^{i+j}}} \\ =\sum_{1\le j\le k}{f_jx^j\sum_{i\ge k-j}{a_ix^i}} \\ =\sum_{1\le j\le k}{f_jx^j\left( F\left( x \right) +\sum_{k-j\le i\le k-1}{a_ix^i} \right)} \\ \\ F\left( x \right) =F\left( x \right) \sum_{1\le j\le k}{f_jx^j}+\sum_{1\le j\le k}{f_jx^j\sum_{k-j\le i\le k-1}{a_ix^i}} \\ F\left( x \right) =\frac{\sum_{1\le j\le k}{f_jx^j\sum_{k-j\le i\le k-1}{a_ix^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{1\le j\le k}{f_jx^j\sum_{k-j\le i-j\le k-1}{a_{i-j}x^{i-j}}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{k\le i\le 2k-1}{\begin{array}{c} x^i\sum_{i-k+1\le j\le k}{a_{i-j}f_j}\\ \end{array}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ F\left( x \right) +\sum_{0\le i\le k-1}{a_ix^i}=\frac{\left( 1-\sum_{1\le j\le k}{f_jx^j} \right) \left( \sum_{0\le i\le k-1}{a_ix^i} \right) +\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k}{f_j\sum_{0\le i\le k-1}{a_ix^{i+j}}}+\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k}{f_j\sum_{j\le i\le j+k-1}{a_{i-j}x^i}}+\sum_{1\le j\le k}{f_j\sum_{k\le i\le k-1+j}{a_{i-j}x^i}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}+\sum_{1\le j\le k}{f_j\left( \sum_{k\le i\le k-1+j}{a_{i-j}x^i}-\sum_{j\le i\le j+k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}+\sum_{1\le j\le k}{f_j\left( \sum_{k\le i\le k-1+j}{a_{i-j}x^i}-\sum_{k\le i\le j+k-1}{a_{i-j}x^i}-\sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \frac{\sum_{0\le i\le k-1}{a_ix^i}+\sum_{1\le j\le k}{f_j\left( -\sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k}{f_j\left( \sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le j\le k-1}{f_j\left( \sum_{j\le i\le k-1}{a_{i-j}x^i} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{\sum_{0\le i\le k-1}{a_ix^i}-\sum_{1\le i\le k-1}{x^i\sum_{1\le j\le i}{a_{i-j}f_j}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{a_0+\sum_{1\le i\le k-1}{x^ia_i}-\sum_{1\le i\le k-1}{x^i\sum_{1\le j\le i}{a_{i-j}f_j}}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ =\frac{a_0+\sum_{1\le i\le k-1}{x^i\left( a_i-\sum_{1\le j\le i}{a_{i-j}f_j} \right)}}{1-\sum_{1\le j\le k}{f_jx^j}} \\ \text{我们得到了原数列}\left\{ a_i \right\} \text{的生成函数}G\left( x \right) =\frac{P\left( x \right)}{Q\left( x \right)} \\ \text{考虑计算这个有理函数的}k\text{次项} \\ \text{令}G_k\left( x \right) =\frac{P\left( x \right)}{Q\left( x \right)}=\frac{P\left( x \right) P\left( -x \right)}{Q\left( x \right) Q\left( -x \right)}=\frac{xA\left( x^2 \right) +B\left( x^2 \right)}{C\left( x^2 \right)} \\ \text{即分母只有偶数次方项},\text{分子的奇数次方项和偶数次方项拆开} \\ \text{如果}k\text{为奇数},\text{那么}x^k\text{之可能存在于}x\frac{A\left( x^2 \right)}{C\left( x^2 \right)}\text{中}\left( \text{只有这边才有奇数次方} \right) \\ \text{同理},\text{如果}k\text{为偶数，那么}x^k\text{只能出现在}\frac{B\left( x^2 \right)}{C\left( x^2 \right)}\text{中，} \\ \text{然后根据}k\text{的奇偶性，继续递归} \\ \text{如果}k\text{是奇数},\text{那么计算}G_{\frac{k-1}{2}}\left( x \right) =\frac{A\left( x \right)}{C\left( x \right)}\left( \text{项的次数除以}2 \right) ,\text{答案为}\left[ x^{\frac{k-1}{2}} \right] G_{\frac{k-1}{2}}\left( x \right) \\ \text{如果}k\text{是偶数},\text{那么计算}G_{\frac{k}{2}}\left( x \right) =\frac{B\left( x \right)}{C\left( x \right)}\text{，答案为}\left[ x^k \right] G_{\frac{k}{2}}\left( x \right) \\ \text{然后有两种选择}:\text{递归到}k=0\text{时，计算}\frac{P\left( 0 \right)}{G\left( 0 \right)}\mathrm{mod}p\text{，此时不需要多项式求逆运算} \\ \text{但常数较大}\left( \text{递归次数多} \right) \\ \text{递归到}k<n\text{时，直接求逆计算。} \\ \\ \\ \\ \\ $$

    代码1: 递归到$k=0$时执行整数运算，较慢

    #include <cstring>
#include <iostream>
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;

const int_t LARGE = 2e6;

const int_t mod = 998244353;
const int_t g = 3;
int_t power(int_t b, int_t i) {
    int_t r = 1;
    if (i < 0)
        i = (i % (mod - 1) + mod - 1) % (mod - 1);
    while (i) {
        if (i & 1)
            r = r * b % mod;
        b = b * b % mod;
        i >>= 1;
    }
    return r;
}

void makeflip(int_t* arr, int_t size2) {
    int_t len = (1 << size2);
    arr[0] = 0;
    for (int_t i = 1; i < len; i++) {
        arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
    }
}

int_t upper2n(int_t x) {
    int_t r = 0;
    while ((1 << r) < x)
        r++;
    return r;
}
template <int_t arg = 1>
void transform(int_t* A, int_t size2, int_t* flip) {
    int_t len = (1 << size2);
    for (int_t i = 0; i < len; i++) {
        int_t r = flip[i];
        if (r > i)
            std::swap(A[i], A[r]);
    }
    for (int_t i = 2; i <= len; i *= 2) {
        int_t mr = power(g, arg * (mod - 1) / i);
        for (int_t j = 0; j < len; j += i) {
            int_t curr = 1;
            for (int_t k = 0; k < i / 2; k++) {
                int_t u = A[j + k], v = curr * A[j + k + i / 2] % mod;
                A[j + k] = (u + v) % mod;
                A[j + k + i / 2] = (u - v + mod) % mod;
                curr = curr * mr % mod;
            }
        }
    }
}
void poly_inv(const int_t* A, int_t n, int_t* result) {
    /*
    计算B(x)*A(x) = 1 mod x^n, 其中A(x)已知
    假设已知A(x)*B(x) = 1 mod x^{ceil(n/2)}
    假设C(x)*A(x) = 1 mod x^n
    (A(x)B(x)-1)^2 = A^2(x)B^2(x)-2A(x)B(x)+1= 0
    A(x)B^2(x)-2B(x)+C(x) = 0
    C(x) = 2B(x)-A(x)B^2(x)
    */
    if (n == 1) {
        result[0] = power(A[0], -1);
        return;
    }
    int_t next = n / 2 + n % 2;
    poly_inv(A, next, result);
    //次数不要选错了，应该用n次的A和B去卷
    int_t size2 = upper2n(n + 2 * next + 1);
    static int_t X[LARGE];
    static int_t Y[LARGE];
    int_t len = (1 << size2);
    memset(X + next, 0, sizeof(X[0]) * (len - n));
    memset(Y + next, 0, sizeof(Y[0]) * (len - next));
    memcpy(X, A, sizeof(A[0]) * n);
    memcpy(Y, result, sizeof(result[0]) * next);
    static int_t fliparr[LARGE];
    makeflip(fliparr, size2);
    transform<1>(X, size2, fliparr);
    transform<1>(Y, size2, fliparr);

    for (int_t i = 0; i < len; i++) {
        X[i] = (2 * Y[i] - X[i] * Y[i] % mod * Y[i] % mod + mod) % mod;
    }
    transform<-1>(X, size2, fliparr);
    const int_t inv = power(len, -1);
    for (int_t i = 0; i < n; i++)
        result[i] = X[i] * inv % mod;
}
int_t poly_divat(const int_t* F, const int_t* G, int_t n, int_t k) {
    /*
        n次多项式F和G
        计算F(x)/G(x)的k次项前系数
        考虑F(x)*G(-x)/G(x)*G(-x)，分母只有偶数次项，写为C(x^2)；分子写成xA(x^2)+B(x^2)，如果k是奇数，那么递归(A,C,n,k/2)，如果k是偶数，那么递归(B,C,n,k/2)
        到k<=n时直接计算
    */
    int_t size2 = upper2n(2 * n + 1);
    int_t len = 1 << size2;
    static int_t fliparr[LARGE];
    makeflip(fliparr, size2);
    static int_t T1[LARGE], T2[LARGE], T3[LARGE];
    for (int_t i = 0; i < len; i++) {
        if (i <= n)
            T1[i] = F[i], T2[i] = G[i];
        else
            T1[i] = T2[i] = T3[i] = 0;
    }
    const int_t inv = power(len, -1);
    while (k != 0) {
        for (int_t i = 0; i < len; i++) {
            if (i <= n) {
                T3[i] = T2[i] * (i % 2 ? (mod - 1) : 1);
            } else
                T3[i] = 0;
        }
        transform(T1, size2, fliparr);
        transform(T2, size2, fliparr);
        transform(T3, size2, fliparr);
        for (int_t i = 0; i < len; i++) {
            T1[i] = T1[i] * T3[i] % mod;
            T2[i] = T2[i] * T3[i] % mod;
        }
        transform<-1>(T1, size2, fliparr);
        transform<-1>(T2, size2, fliparr);
        for (int_t i = 0; i < len; i++) {
            if (i * 2 < len) {
                T2[i] = T2[i * 2] * inv % mod;
            } else
                T2[i] = 0;
        }
        int_t b = k % 2;
        for (int_t i = 0; i < len; i++) {
            if (i % 2 == b) {
                T1[i / 2] = T1[i] * inv % mod;
            }

            if (i > 0)  //防止把T1[0]改为0
                T1[i] = 0;
        }
        k >>= 1;
    }

    return T1[0] * power(T2[0], -1) % mod;
}
void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
    int_t size2 = upper2n(n + m + 1);
    int_t len = 1 << size2;
    static int_t T1[LARGE], T2[LARGE];
    for (int_t i = 0; i < len; i++) {
        if (i <= n)
            T1[i] = A[i];
        else
            T1[i] = 0;
        if (i <= m)
            T2[i] = B[i];
        else
            T2[i] = 0;
    }
    static int_t fliparr[LARGE];
    makeflip(fliparr, size2);
    transform(T1, size2, fliparr);
    transform(T2, size2, fliparr);
    for (int_t i = 0; i < len; i++)
        T1[i] = T1[i] * T2[i] % mod;
    transform<-1>(T1, size2, fliparr);
    int_t inv = power(len, -1);
    for (int_t i = 0; i <= n + m; i++)
        C[i] = T1[i] * inv % mod;
}
int main() {
    std::ios::sync_with_stdio(false);
    static int_t A[LARGE], F[LARGE];
    static int_t T1[LARGE], T2[LARGE];
    int_t n, k;
    cin >> n >> k;
    for (int_t i = 1; i <= k; i++) {
        cin >> F[i];
        F[i] = (F[i] % mod + mod) % mod;
    }
    F[0] = 0;
    for (int_t i = 0; i < k; i++) {
        cin >> A[i];
        A[i] = (A[i] % mod + mod) % mod;
    }
    A[k] = 0;
    poly_mul(A, k, F, k, T1);
    T1[0] = A[0];
    for (int_t i = 1; i <= k - 1; i++) {
        T1[i] = (A[i] - T1[i] + mod) % mod;
    }
    for (int_t i = k; i <= 2 * k; i++)
        T1[i] = 0;
    T2[0] = 1;
    for (int_t i = 1; i <= k; i++)
        T2[i] = (mod - F[i]) % mod;
    int_t r = poly_divat(T1, T2, k, n);
    cout << r << endl;
    return 0;
}
/*
(1+2*x)/(1+x+x^2)
2 5

1 2 0
1 1 1


ans = -2 = 998244351

*/

    代码2: 递归到$k<n$时直接求逆运算，较快，大约是代码1的一半时间

    #include <cstring>
#include <iostream>
using int_t = long long int;
using std::cin;
using std::cout;
using std::endl;

const int_t LARGE = 2e6;

const int_t mod = 998244353;
const int_t g = 3;
int_t power(int_t b, int_t i) {
    int_t r = 1;
    if (i < 0)
        i = (i % (mod - 1) + mod - 1) % (mod - 1);
    while (i) {
        if (i & 1)
            r = r * b % mod;
        b = b * b % mod;
        i >>= 1;
    }
    return r;
}

void makeflip(int_t* arr, int_t size2) {
    int_t len = (1 << size2);
    arr[0] = 0;
    for (int_t i = 1; i < len; i++) {
        arr[i] = (arr[i >> 1] >> 1) | ((i & 1) << (size2 - 1));
    }
}

int_t upper2n(int_t x) {
    int_t r = 0;
    while ((1 << r) < x)
        r++;
    return r;
}
template <int_t arg = 1>
void transform(int_t* A, int_t size2, int_t* flip) {
    int_t len = (1 << size2);
    for (int_t i = 0; i < len; i++) {
        int_t r = flip[i];
        if (r > i)
            std::swap(A[i], A[r]);
    }
    for (int_t i = 2; i <= len; i *= 2) {
        int_t mr = power(g, arg * (mod - 1) / i);
        for (int_t j = 0; j < len; j += i) {
            int_t curr = 1;
            for (int_t k = 0; k < i / 2; k++) {
                int_t u = A[j + k], v = curr * A[j + k + i / 2] % mod;
                A[j + k] = (u + v) % mod;
                A[j + k + i / 2] = (u - v + mod) % mod;
                curr = curr * mr % mod;
            }
        }
    }
}
void poly_inv(const int_t* A, int_t n, int_t* result) {
    /*
    计算B(x)*A(x) = 1 mod x^n, 其中A(x)已知
    假设已知A(x)*B(x) = 1 mod x^{ceil(n/2)}
    假设C(x)*A(x) = 1 mod x^n
    (A(x)B(x)-1)^2 = A^2(x)B^2(x)-2A(x)B(x)+1= 0
    A(x)B^2(x)-2B(x)+C(x) = 0
    C(x) = 2B(x)-A(x)B^2(x)
    */
    if (n == 1) {
        result[0] = power(A[0], -1);
        return;
    }
    int_t next = n / 2 + n % 2;
    poly_inv(A, next, result);
    //次数不要选错了，应该用n次的A和B去卷
    int_t size2 = upper2n(n + 2 * next + 1);
    static int_t X[LARGE];
    static int_t Y[LARGE];
    int_t len = (1 << size2);
    memset(X + next, 0, sizeof(X[0]) * (len - n));
    memset(Y + next, 0, sizeof(Y[0]) * (len - next));
    memcpy(X, A, sizeof(A[0]) * n);
    memcpy(Y, result, sizeof(result[0]) * next);
    static int_t fliparr[LARGE];
    makeflip(fliparr, size2);
    transform<1>(X, size2, fliparr);
    transform<1>(Y, size2, fliparr);

    for (int_t i = 0; i < len; i++) {
        X[i] = (2 * Y[i] - X[i] * Y[i] % mod * Y[i] % mod + mod) % mod;
    }
    transform<-1>(X, size2, fliparr);
    const int_t inv = power(len, -1);
    for (int_t i = 0; i < n; i++)
        result[i] = X[i] * inv % mod;
}
int_t poly_divat(const int_t* F, const int_t* G, int_t n, int_t k) {
    /*
        n次多项式F和G
        计算F(x)/G(x)的k次项前系数
        考虑F(x)*G(-x)/G(x)*G(-x)，分母只有偶数次项，写为C(x^2)；分子写成xA(x^2)+B(x^2)，如果k是奇数，那么递归(A,C,n,k/2)，如果k是偶数，那么递归(B,C,n,k/2)
        到k<=n时直接计算
    */
    int_t size2 = upper2n(2 * n + 1);
    int_t len = 1 << size2;
    static int_t fliparr[LARGE];
    makeflip(fliparr, size2);
    static int_t T1[LARGE], T2[LARGE], T3[LARGE];
    for (int_t i = 0; i < len; i++) {
        if (i <= n)
            T1[i] = F[i], T2[i] = G[i];
        else
            T1[i] = T2[i] = T3[i] = 0;
    }
    const int_t inv = power(len, -1);
    while (k >= n) {
        for (int_t i = 0; i < len; i++) {
            if (i <= n) {
                T3[i] = T2[i] * (i % 2 ? (mod - 1) : 1);
            } else
                T3[i] = 0;
        }
        transform(T1, size2, fliparr);
        transform(T2, size2, fliparr);
        transform(T3, size2, fliparr);
        for (int_t i = 0; i < len; i++) {
            T1[i] = T1[i] * T3[i] % mod;
            T2[i] = T2[i] * T3[i] % mod;
        }
        transform<-1>(T1, size2, fliparr);
        transform<-1>(T2, size2, fliparr);
        for (int_t i = 0; i < len; i++) {
            if (i * 2 < len) {
                T2[i] = T2[i * 2] * inv % mod;
            } else
                T2[i] = 0;
        }
        int_t b = k % 2;
        for (int_t i = 0; i < len; i++) {
            if (i % 2 == b) {
                T1[i / 2] = T1[i] * inv % mod;
            }

            if (i > 0)  //防止把T1[0]改为0
                T1[i] = 0;
        }
        k >>= 1;
    }

    poly_inv(T2, k + 1, T3);
    int_t result = 0;
    //计算结果的k次项
    for (int_t i = 0; i <= k; i++) {
        result = (result + T1[i] * T3[k - i] % mod) % mod;
    }
    return result;
}
void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
    int_t size2 = upper2n(n + m + 1);
    int_t len = 1 << size2;
    static int_t T1[LARGE], T2[LARGE];
    for (int_t i = 0; i < len; i++) {
        if (i <= n)
            T1[i] = A[i];
        else
            T1[i] = 0;
        if (i <= m)
            T2[i] = B[i];
        else
            T2[i] = 0;
    }
    static int_t fliparr[LARGE];
    makeflip(fliparr, size2);
    transform(T1, size2, fliparr);
    transform(T2, size2, fliparr);
    for (int_t i = 0; i < len; i++)
        T1[i] = T1[i] * T2[i] % mod;
    transform<-1>(T1, size2, fliparr);
    int_t inv = power(len, -1);
    for (int_t i = 0; i <= n + m; i++)
        C[i] = T1[i] * inv % mod;
}
int main() {
    std::ios::sync_with_stdio(false);
    static int_t A[LARGE], F[LARGE];
    static int_t T1[LARGE], T2[LARGE];
    int_t n, k;
    cin >> n >> k;
    for (int_t i = 1; i <= k; i++) {
        cin >> F[i];
        F[i] = (F[i] % mod + mod) % mod;
    }
    F[0] = 0;
    for (int_t i = 0; i < k; i++) {
        cin >> A[i];
        A[i] = (A[i] % mod + mod) % mod;
    }
    A[k] = 0;
    poly_mul(A, k, F, k, T1);
    T1[0] = A[0];
    for (int_t i = 1; i <= k - 1; i++) {
        T1[i] = (A[i] - T1[i] + mod) % mod;
    }
    for (int_t i = k; i <= 2 * k; i++)
        T1[i] = 0;
    T2[0] = 1;
    for (int_t i = 1; i <= k; i++)
        T2[i] = (mod - F[i]) % mod;
    int_t r = poly_divat(T1, T2, k, n);
    cout << r << endl;
    return 0;
}
/*
(1+2*x)/(1+x+x^2)
2 5

1 2 0
1 1 1


ans = -2 = 998244351

*/

     

  • 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;
}