矩乘入门学习笔记

<更新档案>


<无>
<第一次更新>
<第二次更新>
<第三次更新>


<前言>

前几天和hxc神仙说想入门一下矩乘。

神仙甩手给我三道题说这几道题会写了那就入门了。

分别是luoguP3702、luoguP5303、CF1182E.

稍微看了一下luogu评分,两黑一紫。。

稍微写了一道P5303就不会。最后终于是在看题解的情况下找到了我对矩乘理解上的错误。。。怪不得一直调不出来。

顺便说一下,我是之前看蓝书入门的,本文对小白可能并不友好,并不讲解基础内容,默认会板子。

矩阵入门参考hxc神仙的blog

最后一道是ACM的考试题。

<正文>

P5303[GXOI/GZOI2019]逼死强迫症

大佬们一看题目就知道是斐波那契的变式。

我看了好久也没看出来。

50%

这题不用矩乘加速的话,就是个递推。

(f_i)表示(2 imes i)的方格满足条件的方案数。

如何从已有状态转移?

  • 我们先考虑从i-1转移到i状态。新加入一列(1 imes 2)的空地

    可以直接放入一块竖的砖,方案数就是加上(f_{i-1})的值,因为不改变(1 imes 1)的方格的取值。

  • 我们再考虑从i-2转移。

    你也可以在这么空出的(4 imes 4)的格子中放入两个横的(1 imes 2)的砖,依旧不会改变(1 imes 1)的方格,方案数加上(f_{i-2})

你会发现上面的方案都没有包括(1 imes 1)的方格在新加入列的情况,其他情况都已囊括。

目前转移为

[f_i=f_{i-1}+f_{i-2}+x ]

其中x为我们接下来要讨论的值。

(g_i)表示(2 imes i)只用(1 imes2)方格填充的方案数。

  • 此时转移就只包括没有包括(1 imes 1)的方格在新加入列的情况,也就是上面的讨论。

    并且初值是(g_0=1,g_1=1),发现这就是个斐波那契数列特别的,(g_0=1)

  • 那么(1 imes1)方格在新加入列的情况就是(2g_{i-1}),为什么呢?

    如图,黄色部分为斐波那契数列的g方案数,在两边塞入正方形。就变成了这样:

    此时只有一种合法方案。但是左边就变成了任意摆放的g。即:


    可以证明左边的方格只能放在(1)(i-2)的列中且每列恰好可以放一个(因为奇偶性质)。
    设它在第k列,则左边方案数为(g_{k-1}),而k取1到i-2.注意右边的矩形上下都可以放,乘个2.
    所以方案数为:

    [2(g_0+g_1+g_2+g_3……g_{i-3}) ]

    g为0开始的斐波那契。
    即前缀和。由斐波那契前缀和公式可知上式

    [egin{aligned} &= 2 g_{i-1} -2\ end{aligned} ]


故总的转态转移方程为:

[f_i=f_{i-1}+f_{i-2}+ 2 g_{i-1}-2 ]

这样,你可以拿到50point(^ - ^)


100%

然后是矩乘优化递推。

我们发现(f_i)与两项f有关,一项g有关(但因为g的递推要两个位),一项常数有关。

于是用(5 imes 5)的矩阵加速递推,初始状态为

[egin{bmatrix} 0&0&1&2&2 end{bmatrix} quad ]

递推矩阵为

[egin{bmatrix} 0&1&0&0&0\ 1&1&0&0&0\ 0&0&0&1&0\ 0&2&1&1&0\ 0&-1&0&0&0 end{bmatrix} quad ]

这些东西都出来了,就已经可以套板子了。

注意本题有模数,这意味着我们得做出一些改变,即把-1改为mod-1.

(mathrm{Code:})

#include <bits/stdc++.h>
#define mod    1000000007
#define int    long long
using namespace std;
int n, T;
int read()
{
    int  s = 0, w = 1;
    char c = getchar();
    while ((c < '0' || c > '9') && c != '-')
        c = getchar();
    if (c == '-')
        w = -1, c = getchar();
    while (c <= '9' && c >= '0')
        s = (s << 3) + (s << 1) + c - '0', c = getchar();
    return s * w;
}
int mul(int a, int b)
{
    return 1LL * a * b % mod;
}
int add(int a, int b)
{
    return a + b >= mod ? a + b - mod : a + b;
}
struct martix
{
    int a[5][5];
    inline martix()
    {
        memset(a, 0, sizeof(a));
    }
    inline martix operator *(martix b)
    {
        martix c;
        for (int i = 0; i < 5; ++i)
            for (int j = 0; j < 5; ++j)
                for (int k = 0; k < 5; ++k)
                    c.a[i][j] = add(c.a[i][j], mul(a[i][k], b.a[k][j]));
        return c;
    }
}   f, d;
int fx[5]    = { 0, 0, 1, 2, 2 };
int dx[5][5] =
{
    { 0,       1, 0, 0, 0 },
    { 1,       1, 0, 0, 0 },
    { 0,       0, 0, 1, 0 },
    { 0,       2, 1, 1, 0 },
    { 0, mod - 1, 0, 0, 1 }
};
martix power(martix a, int b)
{
    martix s;
    for (int i = 0; i < 5; ++i)
        s.a[i][i] = 1;
    for (; b; b >>= 1)
    {
        if (b & 1)
            s = s * a;
        a = a * a;
    }
    return s;
}
signed main()
{
    T = read();
    while (T--)
    {
        memcpy(f.a, fx, sizeof(fx));
        memcpy(d.a, dx, sizeof(dx));
        n = read();
        if (n == 0 || n == 1 || n == 2)
        {
            printf("0
"); continue;
        }
        f = f * power(d, n - 1);
        printf("%lld
", f.a[0][0]);
    }
    return 0;
}

这题在调的过程中给了我很多启发,纠正了很多我在理解上的错误。

hxc神仙给的题就是不一般啊。。。

感觉和黑题的差距还有有一些的,即使标签是黑的。


暂时咕到这里,先去写题了。

哈哈我回来了.


CF1182E Product Oriented Recurrence

开第二题的坑。又一道黑题,然而神仙:CF那题很水的。

这题莫得部分分,所以直接说正解,主要在矩乘的操作上。

100%

有了连乘的恶心操作,我们改咋搞呢?

每项(f_i),显然可以表示成:

[f_i=f_1^{x} imes f_2^{y} imes f_3^{z} imes c^{w} ]

手模前几项:

x y z w
1 0 0 0
0 1 0 0
0 0 1 0
1 1 1 1
1 2 2 6
2 3 4 14
4 6 7 30
...

发现规律了吗?(反正我没有)

总之就可以发现:

[mathrm{x}_i=mathrm{x}_{i-1}+mathrm{x}_{i-2}+mathrm{x}_{i-3}\ mathrm{y}_i=mathrm{y}_{i-1}+mathrm{y}_{i-2}+mathrm{y}_{i-3}\ mathrm{z}_i=mathrm{z}_{i-1}+mathrm{z}_{i-2}+mathrm{z}_{i-3}\ mathrm{w}_i=mathrm{w}_{i-1}+mathrm{w}_{i-2}+mathrm{w}_{i-3}+2 imes mathrm{i} -6 ]

  • 1.对于x、y、z

x、y、z只是初值不同,而递推式相同。

可知初始矩阵为:

[egin{bmatrix} 1&0&0\ 0&1&0\ 0&0&1 end{bmatrix} ]

递推矩阵为

[egin{bmatrix} 0&1&0\ 0&0&1\ 1&1&1 end{bmatrix} ]

这样最后一行即为x、y、z。

可以化简,将初始矩阵往后推一位,得

[egin{bmatrix} 1&1&1 end{bmatrix} ]

直接递推即可。

  • 2.对于w

这个稍微麻烦一点。

(5 imes 5)的矩阵,因为(2i-6)部分要递推。

初始矩阵:

[egin{bmatrix} w_{i-3}&w_{i-2}&w_{i-1}&2i-2&2 end{bmatrix} egin{aligned} &= end{aligned} egin{bmatrix} 2&6&14&8&2 end{bmatrix} ]

递推矩阵

[egin{bmatrix} 0&0&1&0&0\ 1&0&1&0&0\ 0&1&1&0&0\ 0&0&1&1&0\ 0&0&0&1&1 end{bmatrix} ]

然后就可以随便推了。

  • 3.还原答案

我们递推的是指数,想求出最终答案,加个快速幂还原一下就行了。
有了x、y、z、w,可以根据

[f_i=f_1^{x} imes f_2^{y} imes f_3^{z} imes c^{w} ]

求出答案。

  • 4.关于模数

对指数取模显然不能直接用模数。

根据欧拉定理

[a^{m} equiv a^{m\%varphi (p)}(mod p) ]

由于p为质数,(φ(p)=p-1=1e9+6)

至此,本题完美解决。


(mathrm{Code:})

#include <bits/stdc++.h>
#define mod    1000000007
#define md     1000000006
#define N      10
#define int    long long
using namespace std;
int n, m;
int read()
{
    int  s = 0, w = 1;
    char c = getchar();
    while ((c < '0' || c > '9') && c != '-')
        c = getchar();
    if (c == '-')
        w = -1, c = getchar();
    while (c <= '9' && c >= '0')
        s = (s << 3) + (s << 1) + c - '0', c = getchar();
    return s * w;
}
inline int mul(int a, int b)
{
    return 1LL * a * b % mod;
}
inline int add(int a, int b)
{
    return a + b >= mod ? a + b - mod : a + b;
}
int power(int a, int b)
{
    int ans = 1 % mod;
    for (; b; b >>= 1)
    {
        if (b & 1)
            ans = mul(ans, a);
        a = mul(a, a);
    }
    return ans;
}
struct martix
{
    int a[N][N];
    int n, m;
    inline martix()
    {
        n = m = 0;
        memset(a, 0, sizeof(a));
    }
    inline martix operator *(martix b)
    {
        martix c;
        c.n = n; c.m = b.m;
        for (int i = 0; i < c.n; ++i)
            for (int j = 0; j < c.m; ++j)
                for (int k = 0; k < m; ++k)
                    c.a[i][j] = (c.a[i][j] + (a[i][k] * b.a[k][j]) % md) % md;
        return c;
    }
} f1, f2, d1, d2;
martix martix_power(martix a, int b)
{
    martix s; s.n = s.m = a.n;
    for (int i = 0; i < a.n; ++i)
        s.a[i][i] = 1;
    for (; b; b >>= 1)
    {
        if (b & 1)
            s = s * a;
        a = a * a;
    }
    return s;
}
signed main()
{
    n = read();
    int a1 = read(), a2 = read(), a3 = read(), c = read();
    f1.n       = 1; f1.m = 3;
    f1.a[0][0] = f1.a[0][1] = f1.a[0][2] = 1;

    d1.n       = 3; d1.m = 3;
    d1.a[0][1] = d1.a[1][2] = d1.a[2][0] = d1.a[2][1] = d1.a[2][2] = 1;

    f2.n       = 1; f2.m = 5;
    f2.a[0][0] = 2; f2.a[0][1] = 6; f2.a[0][2] = 14;
    f2.a[0][3] = 8; f2.a[0][4] = 2;

    d2.n       = 5; d2.m = 5;
    d2.a[0][2] = d2.a[1][0] = d2.a[1][2] = d2.a[2][1] = d2.a[2][2] = d2.a[3][2] = d2.a[3][3] = d2.a[4][3] = d2.a[4][4] = 1;
    f1         = f1 * martix_power(d1, n - 4);
    f2         = f2 * martix_power(d2, n - 4);
    int ans = 1;
    ans = mul(ans, power(a1, f1.a[0][0]));
    ans = mul(ans, power(a2, f1.a[0][1]));
    ans = mul(ans, power(a3, f1.a[0][2]));
    ans = mul(ans, power(c, f2.a[0][0]));
    printf("%lld
", ans);
    return 0;
}


你以为我会写下一题?

咕了咕了~


G - New Year and Old Subsequence - CodeForces - 750E

ACM的G题,听说cz讲过。

以这个稍微了解一下广义矩乘吧。

线段树优化广义矩乘,其实就是加速dp的递推。

我们将几种状态进行编号:(emptyset,2,20,201,2017)分别编号01234.

这样对于每一个递推,就有一种编号间的递推转移了。

f[i][0/1/2/3/4]表示算到第i位,目前状态为0/1/2/3/4的最小删除次数。

原dp为

[egin{aligned} & f[i][0]=0\ &f[i][1]=min(f[i-1][0],f[i-1][1]+1)[s_i=='2']\ & f[i][2]=min(f[i-1][1],f[i-1][2]+1)[s_i=='0']\ &f[i][3]=min(f[i-1][2],f[i-1][3]+1)[s_i=='1'] \ &f[i][4]=min(f[i-1][3],f[i-1][4]+1)[s_i=='7']\ &f[i][3]=f[i-1][3]+1 , f[i][4]=f[i-1][4]+1[s_i=='6'] end{aligned} ]

重定义乘法为dp的递推,在本题中为加法变取min,乘法变加法。

考虑构造转移矩阵。

对于当前字符为'2'时:若要保持空集,则需删去,代价为1,若不变成空集,则可以不花费代价。其他每个状态保持都无需花费,故转移为:

[mathrm{ egin{bmatrix} 0&1&inf&inf&inf\ inf&0&inf&inf&inf\ inf&inf&0&inf&inf\ inf&inf&inf&0&inf\ inf&inf&inf&inf&0 end{bmatrix} } ]

其他字符类似。

对于字符'6',因为无论如何不能要它,所以从(201、2017)的转移保持代价为1,因为题目为子序列,所以2017这种状态也不能加6

[egin{bmatrix} 0&inf&inf&inf&inf\ inf&0&inf&inf&inf\ inf&inf&0&inf&inf\ inf&inf&inf&1&inf\ inf&inf&inf&inf&1 end{bmatrix} ]

用线段树维护运算,就像正常维护加法一样就行了。

注意左右矩阵相乘顺序,即使这题中没问题,但是其他的广义运算是要注意的。

(mathrm{Code:})

#include <bits/stdc++.h>
#define int    long long
#define N      200010
using namespace std;
int n, m;
struct matrix
{
    int a[5][5];
    int n, m;
    matrix()
    {
        n = m = 0;
        memset(a, 0x3f, sizeof(a));
    }
    inline int *operator [](int x)
    {
        return a[x];
    }
    matrix operator *(matrix b)
    {
        matrix c;
        c.n = n; c.m = b.m;
        for (int i = 0; i < c.n; ++i)
            for (int j = 0; j < c.m; ++j)
                for (int k = 0; k < m; ++k)
                    c.a[i][j] = min(c.a[i][j], a[i][k] + b.a[k][j]);
        return c;
    }
};//矩阵类
struct segment_tree
{
    int    l, r;
    matrix k;
} tr[N << 2];
int read()
{
    int  s = 0, w = 1;
    char c = getchar();
    while ((c < '0' || c > '9') && c != '-')
        c = getchar();
    if (c == '-')
        w = -1, c = getchar();
    while (c <= '9' && c >= '0')
        s = (s << 3) + (s << 1) + c - '0', c = getchar();
    return s * w;
}
char s[N] = {};
inline void push(int p)
{
    tr[p].k = tr[p << 1].k * tr[p << 1 | 1].k;
}
void build(int p, int l, int r)
{
    tr[p].l = l; tr[p].r = r;
    if (l == r)
    {
        tr[p].k.n = tr[p].k.m = 5;
        for (int i = 0; i < 5; ++i)
            tr[p].k[i][i] = 0;
        if (s[l] == '2')
            tr[p].k[0][0] = 1, tr[p].k[0][1] = 0;
        if (s[l] == '0')
            tr[p].k[1][1] = 1, tr[p].k[1][2] = 0;
        if (s[l] == '1')
            tr[p].k[2][2] = 1, tr[p].k[2][3] = 0;
        if (s[l] == '7')
            tr[p].k[3][3] = 1, tr[p].k[3][4] = 0;
        if (s[l] == '6')
            tr[p].k[3][3] = 1, tr[p].k[4][4] = 1;
        return;
    }
    int mid = (l + r) >> 1;
    build(p << 1, l, mid);
    build(p << 1 | 1, mid + 1, r);
    push(p);
}//构造矩阵&线段树
matrix ask(int p, int l, int r)
{
    if (l <= tr[p].l && tr[p].r <= r)
        return tr[p].k;
    int mid = (tr[p].l + tr[p].r) >> 1;
    if (r <= mid)
        return ask(p << 1, l, r);
    if (l > mid)
        return ask(p << 1 | 1, l, r);
    return ask(p << 1, l, r) * ask(p << 1 | 1, l, r);
}//常规询问
signed main()
{
    n = read();
    m = read();
    cin >> (s + 1);
    build(1, 1, n);
    for (int i = 1; i <= m; ++i)
    {
        int    x = read(), y = read();
        matrix ans = ask(1, x, y);
        printf("%lld
", ans[0][4] > n ? -1 : ans[0][4]);
    }//输出答案0转移到4的最小花费
    return 0;
}

大概就这样了。


下一题估计近期不会写,因为要放假了,先鸽着。

原文地址:https://www.cnblogs.com/yywxdgy/p/12196331.html