标签: 特征方程

  • BJOI2019 勘破神机

    强行把两道不相关的题合到一起?不过还好,毕竟可以出到任意二阶递推数列(三阶的其实也可以,只要有人会解三次方程

    $$ \text{考虑第一部分:} \\ \sum_{1\le i\le n}{G\left( i,k \right)}=\sum_{1\le i\le n}{\left( \begin{array}{c} F_i\\ k\\ \end{array} \right) =\frac{1}{k!}\sum_{1\le i\le n}{F_i^{\underline{k}}}},\text{其中}F_i\text{为斐波那契数列,}F_0=\text{1,}F_1=1 \\ =\frac{1}{k!}\sum_{1\le i\le n}{\prod_{0\le j<k}{\left( F_i-i \right)}} \\ \text{构造}k\text{次多项式}H\left( x \right) =\prod_{0\le j<k}{\left( x-j \right)}=\sum_{0\le i\le k}{h_ix^i} \\ \text{则有}\sum_{1\le i\le n}{G\left( i,k \right)}=\frac{1}{k!}\sum_{1\le i\le n}{\sum_{0\le j\le k}{h_jF_{i}^{j}}} \\ =\frac{1}{k!}\sum_{0\le j\le k}{h_j\sum_{1\le i\le n}{F_{i}^{j}}} \\ \text{设}f\left( n,k \right) =\sum_{1\le i\le n}{F_{i}^{k}} \\ \text{由}F_i=\frac{5+\sqrt{5}}{10}\left( \frac{1+\sqrt{5}}{2} \right) ^n+\frac{5-\sqrt{5}}{10}\left( \frac{1-\sqrt{5}}{2} \right) ^n,\text{设}a=\frac{5+\sqrt{5}}{10},b=\frac{5-\sqrt{5}}{10} \\ x_1=\frac{1+\sqrt{5}}{2},x_2=\frac{1-\sqrt{5}}{2} \\ \text{则}F_i=ax_{1}^{n}+bx_{2}^{n} \\ \text{构造模意义下的数域}a+b\sqrt{5} \\ \text{则有}f\left( n,k \right) =\sum_{1\le i\le n}{F_{i}^{k}}=\sum_{1\le i\le n}{\left( ax_{1}^{i}+bx_{2}^{i} \right) ^k}=\sum_{1\le i\le n}{\sum_{0\le j\le k}{\left( \begin{array}{c} k\\ j\\ \end{array} \right) a^jx_{1}^{ij}b^{k-j}x_{2}^{i\left( k-j \right)}}} \\ =\sum_{0\le j\le k}{\left( \begin{array}{c} k\\ j\\ \end{array} \right) a^jb^{k-j}\sum_{1\le i\le n}{x_{1}^{ij}x_{2}^{i\left( k-j \right)}}} \\ =\sum_{0\le j\le k}{\left( \begin{array}{c} k\\ j\\ \end{array} \right) a^jb^{k-j}\sum_{1\le i\le n}{\left( x_{1}^{j}x_{2}^{k-j} \right) ^i}}=\sum_{0\le j\le k}{\left( \begin{array}{c} k\\ j\\ \end{array} \right) a^jb^{k-j}\frac{\left( x_{1}^{j}x_{2}^{k-j} \right) ^{n+1}-1}{x_{1}^{j}x_{2}^{k-j}-1}} \\ \text{第二个求和号等比数列求和即可。} \\ \text{第一部分复杂度}O\left( k^2\log k \right) \\ \text{考虑第二部分。} \\ \text{设}f\left( n \right) \text{表示长度为}n\text{的答案,打表可知} \\ f\left( 2n \right) =4f\left( 2n-2 \right) -f\left( 2n-4 \right) \\ \text{只考虑取偶数的}2n,\text{则}f\left( n \right) =4f\left( n-1 \right) -f\left( n-2 \right) \text{。} \\ \text{计算特征根}x_1=2+\sqrt{3},x_2=2-\sqrt{3} \\ \text{设}f\left( n \right) =ax_{1}^{n}+bx_{2}^{n},\text{由}f\left( 0 \right) =\text{1,}f\left( 1 \right) =\text{3得} \\ a=\frac{3+\sqrt{3}}{6},b=\frac{3-\sqrt{3}}{6} \\ \text{仍然像第一部分一样处理即可。} \\ $$

    #include <assert.h>
    #include <algorithm>
    #include <cstring>
    #include <iostream>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t mod = 998244353;
    const int_t LARGE = 1000;
    
    template <class T>
    T power(T base, int_t index) {
        base = (base % mod + mod) % mod;
        T result = 1;
        while (index) {
            if (index & 1) result = result * base % mod;
            index >>= 1;
            base = base * base % mod;
        }
        return result;
    }
    template <int_t I>
    struct Z {
        int_t a = 0, b = 0;
        Z(int_t a = 0, int_t b = 0) : a(a), b(b) {}
        Z<I> operator+(const Z<I>& rhs) const {
            return Z<I>{(a + rhs.a) % mod, (b + rhs.b) % mod};
        }
        Z<I> operator-(const Z<I>& rhs) const {
            return Z<I>{(a - rhs.a + mod) % mod, (b - rhs.b + mod) % mod};
        }
        Z<I> operator*(const Z<I>& rhs) const {
            return Z<I>{(a * (rhs.a % mod) % mod + I * b % mod * rhs.b % mod) % mod,
                        (a * rhs.b % mod + b * rhs.a % mod) % mod};
        }
        Z<I> inv() const {
            int_t under = power(
                (a * a % mod - I * b % mod * b % mod + I * mod) % mod, mod - 2);
            return Z<I>{a * under % mod, (mod - b) * under % mod};
        }
        Z<I> operator/(const Z<I>& rhs) const { return (*this) * rhs.inv(); }
        Z<I> operator%(int_t x) const { return Z<I>{a % x, b % x}; }
    };
    int_t fact[LARGE + 1], inv[LARGE + 1], factInv[LARGE + 1];
    int_t C(int_t n, int_t m) {
        return fact[n] * factInv[m] % mod * factInv[n - m] % mod;
    }
    template <int_t I>
    Z<I> sumQ(Z<I> x, int_t n) {
        if (x.a == 1 && x.b == 0) return n;
        return (power(x, n) - 1) / (x - 1);
    }
    //求斐波那契数列的k次幂的前缀和
    int_t sumF(int_t n, int_t k) {
        if (k == 0) return n + 1;
        using Z = ::Z<5>;
        const Z a = Z(5, 1) / 10;
        const Z b = Z(5, mod - 1) / 10;
        const Z x1 = Z(1, 1) / 2, x2 = Z(1, mod - 1) / 2;
        Z result = 0;
        for (int_t i = 0; i <= k; i++) {
            Z x = power(x1, i) * power(x2, k - i);
            result =
                result + power(a, i) * C(k, i) * power(b, k - i) * sumQ(x, n + 1);
        }
        assert(result.b == 0);
        return result.a;
    }
    int_t sumG(int_t n, int_t k) {
        if (k == 0) return n + 1;
        using Z = ::Z<3>;
        const Z a = Z(3, 1) / 6;
        const Z b = Z(3, mod - 1) / 6;
        const Z x1 = Z(2, 1), x2 = Z(2, mod - 1);
        Z result = 0;
        for (int_t i = 0; i <= k; i++) {
            Z x = power(x1, i) * power(x2, k - i);
            result =
                result + power(a, i) * C(k, i) * power(b, k - i) * sumQ(x, n + 1);
        }
        assert(result.b == 0);
        return result.a;
    }
    void poly_mul(const int_t* A, int_t n, const int_t* B, int_t m, int_t* C) {
        static int_t Cx[LARGE * 2 + 1];
        std::fill(Cx, Cx + n + m + 1, 0);
        for (int_t i = 0; i <= n; i++)
            for (int_t j = 0; j <= m; j++) {
                Cx[i + j] = (Cx[i + j] + A[i] * B[j] % mod) % mod;
            }
        std::copy(Cx, Cx + n + m + 1, C);
    }
    int_t poly[LARGE + 1];
    template <class Func>
    int_t process(int_t n, int_t k, Func F) {
        int_t result = 0;
        for (int_t i = 0; i <= k; i++) {
            result = (result + poly[i] * F(n, i) % mod) % mod;
        }
        result = result * factInv[k] % mod;
        return result;
    }
    int main() {
        fact[0] = fact[1] = factInv[0] = factInv[1] = inv[1] = 1;
        for (int_t i = 2; i <= LARGE; i++) {
            inv[i] = (mod - mod / i) * inv[mod % i] % mod;
            fact[i] = fact[i - 1] * i % mod;
            factInv[i] = factInv[i - 1] * inv[i] % mod;
        }
        int_t T, m, left, right, k;
        cin >> T >> m >> left >> right >> k;
        poly[0] = 1;
        for (int_t i = 0; i < k; i++) {
            int_t curr[] = {(mod - i) % mod, 1};
            poly_mul(poly, k, curr, 1, poly);
        }
        if (m == 2) {
            int_t result =
                (process(right, k, sumF) - process(left - 1, k, sumF) + mod) % mod *
                power(right - left + 1, mod - 2) % mod;
            cout << result << endl;
        } else {
            int_t result = (process(right / 2, k, sumG) -
                            process((left - 1) / 2, k, sumG) + mod) %
                           mod * power(right - left + 1, mod - 2) % mod;
            cout << result << endl;
        }
        return 0;
    }

     

  • 洛谷5110 块速递推

    $$ \text{特征根}x_1=\frac{233}{2}+\frac{13\sqrt{337}}{2},x_2=-\frac{13\sqrt{337}}{2}+\frac{233}{2} \\ \text{通项}a_n=ax_{1}^{n}+bx_{2}^{n} \\ \text{联立方程}\begin{cases} 0=a+b\\ 1=x_1a+x_2b\\ \end{cases} \\ \text{解出}a=\frac{\sqrt{337}}{4381},b=-\frac{\sqrt{337}}{4381} \\ \text{根据欧拉准则,}\left( 337 \right) ^{\frac{p-1}{2}}=\text{1,故求解二次剩余}216284168^2\equiv 337\left( mod\,\,10^9+7 \right) \\ \text{故在模}10^9+\text{7意义下,}x_1=\text{905847205,}x_2=94153035 \\ a=766769301 \\ \text{故}a_n=766769301\left( 905847205^n-94153035^n \right) \\ \text{对于一次询问,令}x=km+p\left( m\text{为块大小} \right) \text{,则有}x_{1}^{km+p}=x_{1}^{km}\times x_{1}^{p},\text{其中}p<m \\ \text{按照大小为}m\text{分块即可。} \\ \\ $$

    #include <cstdio>
    #include <iostream>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    namespace Mker {
    unsigned long long SA, SB, SC;
    void init() { scanf("%llu%llu%llu", &SA, &SB, &SC); }
    unsigned long long rand() {
        SA ^= SA << 32, SA ^= SA >> 13, SA ^= SA << 1;
        unsigned long long t = SA;
        SA = SB, SB = SC, SC ^= t ^ SA;
        return SC;
    }
    }  // namespace Mker
    const int_t mod = 1e9 + 7;
    const int_t phi = mod - 1;
    const int_t a = 766769301;
    const int_t x1 = 905847205;
    const int_t x2 = 94153035;
    const int_t BLOCK_SIZE = 65536;
    int T;
    int_t x1s[mod / BLOCK_SIZE + 1], x2s[mod / BLOCK_SIZE + 1];
    int_t x1b[BLOCK_SIZE], x2b[BLOCK_SIZE];
    int_t power(int_t base, int_t index) {
        int_t result = 1;
        while (index) {
            if (index & 1) result = result * base % mod;
            index >>= 1;
            base = base * base % mod;
        }
        return result;
    }
    int main() {
        cin >> T;
    
        x1s[0] = x2s[0] = 1;
        int_t prex1 = power(x1, BLOCK_SIZE);
        int_t prex2 = power(x2, BLOCK_SIZE);
        for (int_t i = 1; i <= mod / BLOCK_SIZE; i++) {
            x1s[i] = x1s[i - 1] * prex1 % mod;
            x2s[i] = x2s[i - 1] * prex2 % mod;
        }
        x1b[0] = x2b[0] = 1;
        for (int_t i = 1; i < BLOCK_SIZE; i++) {
            x1b[i] = x1b[i - 1] * x1 % mod;
            x2b[i] = x2b[i - 1] * x2 % mod;
        }
        int_t result = 0;
        Mker::init();
        while (T--) {
            // int_t curr = 0;
            int x = (unsigned long long int)Mker::rand() % phi;
            int_t curr = ((x1s[x >> 16] * x1b[x & ((1 << 16) - 1)] -
                           x2s[x >> 16] * x2b[x & ((1 << 16) - 1)]) %
                              mod +
                          mod) %
                         mod * a % mod;
            result ^= curr;
        }
        cout << result << endl;
        return 0;
    }

     

  • 洛谷2012 拯救世界2

    $$ \text{考虑对于每一种基因搞出来一个指数生成函,最后把他们乘起来。} \\ \text{注意到}\frac{\sin \left( ix \right)}{i}=\frac{\sinh \left( x \right)}{i}=\sum_{i\ge 0}{\frac{x^{2i+1}}{\left( 2i+1 \right) !}},\cos \left( ix \right) =\cosh \left( x \right) =\sum_{i\ge 0}{\frac{x^{2i}}{\left( 2i \right) !}} \\ \text{所以基因序列的生成函数为}F\left( x \right) =\left( \frac{\sinh \left( x \right)}{i}\cosh \left( x \right) e^x \right) ^4=\sinh \left( x \right) ^4\cosh \left( x \right) ^4e^{4x} \\ \text{大力麦克劳林展开发现} \\ F\left( x \right) =\frac{24}{\text{4!}}x^4+\frac{480}{\text{5!}}x^5+\frac{107520x^6}{\text{6!}}+\frac{1419264x^7}{\text{7!}}+\frac{18063360x^8}{\text{8!}}+\frac{225116160x^9}{\text{9!}}+…. \\ \text{所以答案即为}F^{\left( n \right)}\left( 0 \right) ,\text{同时注意到,对于}n<\text{4,答案为}0. \\ \text{经过}sympy\text{打表发现,答案极有可能是一个线性递推式,故写程序验证,发现} \\ f\left( n \right) =20f\left( n-1 \right) -80f\left( n-2 \right) -320f\left( n-3 \right) +1536f\left( n-4 \right) ,\text{同时} \\ f\left( 4 \right) =\text{24,}f\left( 5 \right) =\text{480,}f\left( 6 \right) =\text{7680,}f\left( 7 \right) =\text{107520,对于}k\le \text{3,}f\left( k \right) =0 \\ \text{然后直接矩乘即可通过本题}\left( \text{然而我一开始读入写挂导致在结尾没有回车或者空格的点}TLE \right) \\ \text{如果追求更高的速度,可以算一下这个递推的特征根,}x^4=20x^3-80x^2-320x+\text{1536,解得} \\ x=\left\{ -\text{4,4,8,}12 \right\} \\ \text{之后大力解方程求出来系数,故} \\ f\left( n-4 \right) =\left( -4 \right) ^n+6\times 4^n-64\times 8^n+81\times 12^n \\ \text{直接整数快速幂即可。} \\ $$

    from sympy import *
    x = symbols("x")
    f = ((sinh(x)*cosh(x)*exp(x))**4).diff().diff().diff().diff()
    coefs = []
    for i in range(20):
        coefs.append(f.subs([(x, 0)]))
        f = f.diff()
    print(coefs)
    
    eqs = []
    xs = list([symbols("x%d" % i) for i in range(8)])
    A = 4
    for i in range(A+4):
        a = 0
        for j in range(i, i+A):
            a += xs[j-i]*coefs[j]
        eqs.append(Eq(a, coefs[i+A]))
    print(eqs)
    print(solve(eqs))
    #pragma GCC target("avx2,sse")
    #include <assert.h>
    #include <cctype>
    #include <cstdio>
    #include <cstring>
    #include <iostream>
    using int_t = unsigned long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int mod = 1e9;
    const int phi = 400000000;
    template <class T>
    void read(T& x) {
        x = 0;
        char chr = getchar();
        assert(chr != '-');
        int counter = 0;
        while (isdigit(chr) == false) {
            int curr = getchar();
            if (curr == -1) return;
            chr = curr;
            counter++;
            assert(counter <= 10);
        }
        while (isdigit(chr)) {
            x = x * 10 + (chr - '0');
            chr = getchar();
        }
    }
    template <class T>
    void write(T x) {
        assert(x >= 0);
        if (x > 9) write(x / 10);
        putchar('0' + x % 10);
    }
    int main() {
        int counter = 0;
        while (true) {
            int_t x;
            read(x);
            if (x == 0) break;
            counter += 1;
            assert(counter <= 200000);
            if (x > phi) x = x % phi + phi;
            if (x < 4) {
                write(0);
                putchar('\n');
                continue;
            }
            x -= 4;
            int t2 = 1;
            int t3 = 1;
            int b3 = 3, b2 = 2;
            int index = x;
            while (index) {
                t2 = (int_t)t2 * ((index & 1) ? b2 : 1) % mod;
                t3 = (int_t)t3 * ((index & 1) ? b3 : 1) % mod;
                b3 = (int_t)b3 * b3 % mod;
                b2 = (int_t)b2 * b2 % mod;
                index >>= 1;
            }
            int t4 = (int_t)t2 * t2 % mod;
            int result = ((int_t)((x & 1) ? (mod - 1) : 1) + (int_t)81 * t3 + 6 -
                          (int_t)64 * t2 % mod + mod) %
                         mod * t4 % mod;
            write(result);
            putchar('\n');
        }
        return 0;
    }