玄学小记.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("");
    }
}