玄学小记.2 ~ bzoj 4162: shlw loves matrix II

我们需要求\(g(P)\)，其中\(g\)是一个只有一项的多项式

暴力是\(n ^ 3 \log k\)的，显然过不去

怎么办？

 

特征方程优化矩阵快速幂～

考虑方程\(|xI - P| = 0\)，把det展开后可以得到一个方程\(f(x) = 0\)，这个方程称为\(P\)的特征方程，\(f\)称为\(P\)的特征多项式

基于某个厉害的定理，有\(f(P) = 0\)，也就是说矩阵\(P\)是特征方程的一个解

于是有\(g(P) = (g \mod f)(P)\)，由于\(f\)的次数只有\(n\)，而且多项式取模的复杂度很容易做到\(O(n ^ 2)\)

因此只要算出特征多项式就可以\(O(n ^ 4 + n ^2 \log k)\)计算

 

那么怎么求出特征多项式呢？

 

牛顿恒等式～

对于多项式\(F(x) = \sum_{i = 0}^{n} c_{n-i} x^i\)

若有方程\(F(x) = 0\)，它显然有\(n\)个根 \(x_1, x_2 ... x_n\)，令\(s_{k} = \sum_{i = 1}^{n} x_{i}^k\)

则有：

当\(1 \leq k \leq n\)时，有\(\sum_{i = 0}^{k - 1} c_{i} * s_{k - i} +c_{k} * k = 0\) 

当\(n < k\)时，有\(\sum_{i = 0}^{n} c_{i} * s_{k - i} = 0\)

 

这玩意对于矩阵同样成立～

基于另一个厉害的定理，特征方程\(f(x) = 0\)的根\(\lambda_{i}\)满足：\(s_k = \sum \lambda _{i}^{k} = tr(P^k)\)

于是可以通过牛顿恒等式解出特征方程每一项的系数～

这样就可以\(O(n^4)\)解出特征多项式了了……

时间复杂度\(O(n ^ 4 + n ^2 \log k)\)

#include <bits/stdc++.h>
using namespace std;

// #define int long long
const int N = 55;
const int MOD = 1000000007;
char st[20000]; int k;
struct mat
{
    int d[N][N];
    mat()
    {
        memset(d, 0, sizeof d);
    }
} X, P[N], I, ans;
mat operator * (mat a, mat b)
{
    mat c;
    for (int i = 0; i < k; ++ i)
        for (int j = 0; j < k; ++ j)
            for (int p = 0; p < k; ++ p)
                c.d[i][j] = (1ll * a.d[i][p] * b.d[p][j] + c.d[i][j]) % MOD;
    return c;
}
mat operator * (mat a, int b)
{
    mat c;
    for (int i = 0; i < k; ++ i)
        for (int j = 0; j < k; ++ j)
            c.d[i][j] = 1ll * a.d[i][j] * b % MOD;
    return c;
}
mat operator + (mat a, mat b)
{
    mat c;
    for (int i = 0; i < k; ++ i)
        for (int j = 0; j < k; ++ j)
            c.d[i][j] = (a.d[i][j] + b.d[i][j]) % MOD;
    return c;
}
int powi(int a, int b)
{
    int c = 1;
    for (; b; b /= 2, a = 1ll * a * a % MOD)
        if (b & 1) c = 1ll * c * a % MOD;
    return c;
}
int S[N], C[N];

struct poly
{
    int d[N * 2];
    poly ()
    {
        memset(d, 0, sizeof d);
    }
} pI, pS;
poly operator * (poly a, poly b)
{
    poly c;
    for (int i = 0; i < k; ++ i)
        for (int j = 0; j < k; ++ j)
            c.d[i + j] = (c.d[i + j] + 1ll * a.d[i] * b.d[j]) % MOD;
    for (int i = k + k - 2; i >= k; -- i)
    {
        int v = c.d[i];
        for (int j = 0; j <= k; ++ j)
            c.d[i - j] = (c.d[i - j] - 1ll * v * C[j] % MOD + MOD) % MOD;
    }
    return c;
}
poly powi(poly a, char b[], int l)
{
    poly c = pI;
    for (int i = l - 1; ~i; -- i, a = a * a)
        if (b[i] == '1')
            c = c * a;
    return c;
}
signed main()
{
    scanf("%s%d", st, &k);
    int l = strlen(st);
    for (int i = 0; i < k; ++ i) I.d[i][i] = 1;
    pI.d[0] = 1;
    for (int i = 0; i < k; ++ i)
        for (int j = 0; j < k; ++ j)
            scanf("%d", &X.d[i][j]);
    P[0] = I; C[0] = 1;
    for (int i = 1; i <= k; ++ i)
    {
        P[i] = P[i - 1] * X;
        for (int j = 0; j < k; ++ j)
            S[i] = (S[i] + P[i].d[j][j]) % MOD;
        for (int j = 1; j <= i; ++ j)
            C[i] = (C[i] - 1ll * C[j - 1] * S[i - j + 1] % MOD + MOD) % MOD;
        C[i] = 1ll * C[i] * powi(i, MOD - 2) % MOD;
    }
    pS.d[1] = 1;
    pS = powi(pS, st, l);


    for (int i = 0; i < k; ++ i)
        ans = ans + P[i] * pS.d[i];
    for (int i = 0; i < k; ++ i)
    {
        for (int j = 0; j < k - 1; ++ j)
            printf("%d ", ans.d[i][j]);
        printf("%d", ans.d[i][k - 1]);
        puts("");
    }
}