矩阵迭代加速

一、矩阵快速幂

矩阵快速幂和普通快速幂一样，只需要创建矩阵类型，对乘法进行重载

二、矩阵加速

差分方程问题

1. 斐波那契数列：

   $$a_x = \begin{split}\begin{cases} 1 & & x \in \{1,2\} \\ a_{x-1} + a_{x - 2} & & x \ge 3 \end{cases}\end{split}$$

   快速求第 $ n $ 项，可以使用通项公式，也可以使用矩阵快速幂

   $$\begin{split}a_x &= 1 \times a_{x-1} + 1\times a_{x-2} \\ a_{x-1} &= 1 \times a_{x-1} + 0 \times a_{x-2} \end{split}$$

   因此可以得到

   $$\begin{bmatrix} a_x \\ a_{x-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \times \begin{bmatrix} a_{x-1} \\ a_{x-2} \end{bmatrix}$$

   基本上这个系数矩阵的第一行表示的是差分关系，然后其他的都是每行只有一个 $ 1 $ ，是一个 $ I $ 矩阵

2. 类斐波那契数列：

   $$a_x = \begin{split}\begin{cases} 1 & & x \in \{1,2,3\} \\ a_{x-1} + a_{x - 3} & & x \ge 4 \end{cases}\end{split}$$

$$\begin{split}a_x &= 1 \times a_{x-1} + 0\times a_{x-2} + 1 \times a_{x-3} \\ a_{x-1} &= 1 \times a_{x-1} + 0 \times a_{x-2} + 0 \times a_{x-3}\\ a_{x-2} &= 0 \times a_{x-1} + 1 \times a_{x-2} + 0 \times a_{x-3} \end{split}$$ $$\begin{bmatrix} a_x \\ a_{x-1} \\ a_{x-2} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix} \times \begin{bmatrix} a_{x-1} \\ a_{x-2} \\ a_{x-3} \\ \end{bmatrix}$$

三、例题

1. P3216 [HNOI2011]数学作业

发现数据范围非常大， $ O(n) $ 根本解决不了，需要 $ \log(n) $ 的做法。

$ f_n = f_{n-1} \times 10^{k} + n $ ，因此很容易可以得到这个矩阵乘法，然后使用加速幂

$$\left[\begin{array}{c} f_{n} \\ n \\ 1 \end{array}\right]=\left[\begin{array}{ccc} 10^{k} & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{array}\right] \times\left[\begin{array}{c} f_{n-1} \\ n-1 \\ 1 \end{array}\right]$$

但是发现，指数 k 的值为 $ \left\lfloor\log_{10}i\right\rfloor+1 $ （也就是 $ len(i) $ 的值是会变化的，没法直接矩乘。由于 $ len(i) $ 的值一共只有 $ 18 $ 种左右，所以我们直接枚举 然后分 $ len(i) $ 块矩乘后合并。

#include <bits/stdc++.h>
#define ll long long
namespace FAST_IO
{
    template <class T>
    inline void read(T &x)
    {
        T flag = 1;
        x = 0;
        char ch = getchar();
        while (ch < '0' || ch > '9')
        {
            if (ch == '-')
                flag = -1;
            ch = getchar();
        }
        while (ch >= '0' && ch <= '9')
        {
            x = (x << 3) + (x << 1) + (ch ^ 48), ch = getchar();
        }
        x *= flag;
    }

    template <class T, class... _T>
    inline void read(T &x, _T &...y)
    {
        return read(x), read(y...);
    }

} // namespace FAST_IO
using namespace std;
using namespace FAST_IO;
ll mod = 1e9 + 7;
const int INF = 0x3f3f3f3f;
const ll INF_LL = 0x3f3f3f3f3f3f3f3f;
const double eqs = 1e-5;
const int maxn = 10 + 10;
const int maxm = 1e5 + 10;
ll t, n, m, k;

// 使用的时候注意一下矩阵乘法自带mod，矩阵的len，type等
struct Matrix
{
    using type = ll;
    static const int len = maxn;
    type a[len][len];
    int m, n; // 矩阵的行列
    int dim;  // 矩阵的维度
    Matrix(int n, int m)
    {
        this->m = m, this->n = n;
        clean();
    }
    Matrix(int n)
    {
        this->m = this->n = n;
        clean();
    }
    Matrix(int n, char ch)
    {
        this->m = this->n = n;
        if (ch == 'i' || ch == 'I')
        {
            for (int i = 1; i <= n; i++)
                for (int j = 1; j <= m; j++)
                    a[i][j] = (i == j);
        }
    }

    void clean()
    {
        for (int i = 1; i <= n; i++)
            for (int j = 1; j <= m; j++)
                a[i][j] = 0;
    }

    Matrix operator=(const Matrix &a)
    {

        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= m; j++)
            {
                this->a[i][j] = a.a[i][j];
            }
        }
        return *this;
    }

    Matrix operator%(const ll &mod) const
    {
        Matrix ret(n, m);
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= m; j++)
            {
                ret.a[i][j] = a[i][j] % mod;
            }
        }
        return ret;
    }

    Matrix operator%=(ll mod)
    {
        *this = *this % mod;
        return *this;
    }

    Matrix operator*(const Matrix &a) const
    {
        Matrix ret(n, a.m);
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= a.m; j++)
            {
                for (int k = 1; k <= m; k++)
                {
                    ret.a[i][j] = (ret.a[i][j] + this->a[i][k] * a.a[k][j] % mod) % mod;
                }
            }
        }
        return ret;
    }
    Matrix operator*=(const Matrix &a)
    {
        *this = (*this) * a;
        return *this;
    }

    type *operator[](const int x)
    {
        return a[x];
    }

    Matrix qpow(ll k)
    {
        Matrix ans(n, 'i'), tmp = *this;
        while (k)
        {
            if (k & 1)
                ans = ans * tmp;
            tmp = tmp * tmp;
            k >>= 1;
        }
        return ans;
    }
};

int lenof(ll n)
{
    int i;
    for (i = 0; n; i++)
        n /= 10;
    return i;
}
ll pow10(int x)
{
    ll ans = 1;
    for (int i = 1; i <= x; i++)
        ans *= 10;
    return ans;
}
int main()
{
// #define COMP_DATA
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    read(n, m);
    Matrix a(3, 1);
    mod = m;
    int len = lenof(n);
    a[1][1] = a[2][1] = a[3][1] = 1;
    for (int i = 0; i < len - 1; i++)
    {
        Matrix b(3, 3);
        b[1][1] = pow10(i + 1) % mod;
        b[1][2] = b[1][3] = b[2][2] = b[2][3] = b[3][3] = 1;
        m = pow10(i + 1) - pow10(i);
        a = (b.qpow(m - (i == 0 ? 1 : 0))) * a;
    }
    Matrix b(3, 3);
    b[1][1] = pow10(len) % mod;
    b[1][2] = b[1][3] = b[2][2] = b[2][3] = b[3][3] = 1;
    m = n - pow10(len - 1) + 1;
    a = (b.qpow(m - (len - 1 == 0 ? 1 : 0))) * a;
    printf("%lld\n", a[1][1]);
    return 0;
}

2. P2886 [USACO07NOV]Cow Relays G

还是使用矩阵乘法，矩阵乘法重载一下乘号，变成 $ \min $ 这种形式， $ a^k $ 表示恰好经过 $ k $ 条边最短路。

注意求的是最短路，然后需要初始化为 $ \inf $ ，注意这个里面 $ ans $ ，需要把对角线置为 $ 0 $ ， $ ans $ 作为单位矩阵，需要保持不变

无向图得到的是一个对称阵，运算一下发现，刚好对角线为 $ 0 $ 的时候才满足

#include <bits/stdc++.h>
#define ll long long
namespace FAST_IO
{
    template <class T>
    inline void read(T &x)
    {
        T flag = 1;
        x = 0;
        char ch = getchar();
        while (ch < '0' || ch > '9')
        {
            if (ch == '-')
                flag = -1;
            ch = getchar();
        }
        while (ch >= '0' && ch <= '9')
        {
            x = (x << 3) + (x << 1) + (ch ^ 48), ch = getchar();
        }
        x *= flag;
    }

    template <class T, class... _T>
    inline void read(T &x, _T &...y)
    {
        return read(x), read(y...);
    }

} // namespace FAST_IO
using namespace std;
using namespace FAST_IO;
const ll mod = 1e9 + 7;
const int INF = 0x3f3f3f3f;
const ll INF_LL = 0x3f3f3f3f3f3f3f3f;
const double eqs = 1e-5;
const int maxn = 3e2 + 10;
const int maxm = 1e5 + 10;
// 使用的时候注意一下矩阵乘法自带mod，矩阵的len，type等
struct Matrix
{
    using type = ll;
    static const int len = maxn;
    type a[len][len];
    int m, n; // 矩阵的行列
    int dim;  // 矩阵的维度
    Matrix(int n, int m)
    {
        this->m = m, this->n = n;
        clean();
    }
    Matrix(int n)
    {
        this->m = this->n = n;
        clean();
    }
    Matrix(int n, char ch)
    {
        this->m = this->n = n;
        if (ch == 'i' || ch == 'I')
        {
            for (int i = 1; i <= n; i++)
                for (int j = 1; j <= m; j++)
                    a[i][j] = (i == j);
        }
        else if (ch == 'u' || ch == 'U')
        {
            for (int i = 1; i <= n; i++)
                for (int j = 1; j <= m; j++)
                    a[i][j] = INF_LL;
        }
    }

    void clean()
    {
        for (int i = 1; i <= n; i++)
            for (int j = 1; j <= m; j++)
                a[i][j] = 0;
    }
    type *operator[](const int x)
    {
        return a[x];
    }

    Matrix operator=(const Matrix &a)
    {

        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= m; j++)
            {
                this->a[i][j] = a.a[i][j];
            }
        }
        return *this;
    }

    Matrix operator%(const ll &mod) const
    {
        Matrix ret(n, m);
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= m; j++)
            {
                ret.a[i][j] = a[i][j] % mod;
            }
        }
        return ret;
    }

    Matrix operator%=(ll mod)
    {
        *this = *this % mod;
        return *this;
    }

    Matrix operator*(const Matrix &a) const
    {
        Matrix ret(n, 'u');
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= a.m; j++)
            {
                for (int k = 1; k <= m; k++)
                {
                    // ret.a[i][j] = (ret.a[i][j] + this->a[i][k] * a.a[k][j] % mod) % mod;
                    ret.a[i][j] = min(ret.a[i][j], (this->a)[i][k] + a.a[k][j]);
                }
            }
        }
        return ret;
    }
    Matrix operator*=(const Matrix &a)
    {
        *this = (*this) * a;
        return *this;
    }

    Matrix qpow(ll k)
    {
        Matrix ans(n, 'u'), tmp = *this;
        for (int i = 1; i <= n; i++)
            ans[i][i] = 0;
        while (k)
        {
            if (k & 1)
                ans = ans * tmp;
            tmp = tmp * tmp;
            k >>= 1;
        }
        return ans;
    }
};

int n, t, s, e;
map<int, int> id;
int main()
{
// #define COMP_DATA
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    read(n, t, s, e);
    int cnt = 0;
    Matrix g(t << 1, 'u');
    for (int i = 1, u, v, w; i <= t; i++)
    {
        read(w, u, v);
        if (!id.count(u))
            id[u] = ++cnt;
        if (!id.count(v))
            id[v] = ++cnt;
        int uu = id[u], vv = id[v];
        g[uu][vv] = g[vv][uu] = min((ll)w, g[uu][vv]);
    }
    g.n = g.m = cnt;
    g = g.qpow(n);
    printf("%lld\n", g[id[s]][id[e]]);
    return 0;
}