Skip to content

题解-CF575A

Posted on:2023年8月12日 at 22:13

本文章遵守知识共享协议 CC-BY-NC-SA,同步发表于洛谷题解区,转载时须在文章的任一位置附上原文链接和作者署名(rickyxrc)。推荐在我的个人博客阅读。

是一道矩阵加速的题目,但是矩阵加速的方法不止快速幂,还有倍增。

题面大意

给你一个数列 SS',你需要求出数列 FF 的第 kk 项。

其中 Fi={0(i=0)Fi1Si1+Fi2+Si2(i0)F_i = \begin{cases}0 & (i=0) \\ F_{i-1}S_{i-1} + F_{i-2}+S_{i-2} & (i \neq 0) \end{cases}(类似于斐波那契数列)。

无限的数列 SS 可以看作一个有限的数列 SS' 的循环,但是 SS 中某些位置的值是特殊的,会专门给入。

解题思路

这是一个类似于斐波那契数列的转移,那么我们不难想到使用矩阵维护转移状态。

具体的,定义状态矩阵 Ti=[Fi+1Fi]T_i = \begin{bmatrix}F_{i+1} & F_i\end{bmatrix},转移矩阵 Si=[Si+11Si0]S_i = \begin{bmatrix}S_{i+1} & 1 \\ S_i & 0\end{bmatrix},那么我们可以得到 TiSi=Ti+1T_i S_i = T_{i+1},这是下文讨论内容的前提。

然后我们就考虑如何加速,不难发现特殊点的数量和值域相比极少,也就是说大多数情况下,我们的计算采取的转移策略是相同的,是许多转移矩阵的乘积。

形式化地说,我们是在求 T0i=0kSiT_0\prod_{i=0}^kS_i,而中间的 SiS_i 连乘中,有许多部分是相同的,我们可以将其打包,加速运算。

这里解法出现了分歧,可以使用线段树加矩阵快速幂(类似于分块的思想),当然也可以使用倍增。

我们定义 Pi,jP_{i,j} 为从 SiS_i 开始连乘,总长度为 2j2^j 的矩阵。

那么我们的转移公式就应该是 Pi,j=Pi,j1Pi+2j1,j1P_{i,j} = P_{i,j-1} P_{i+2^{j-1},j-1}(实际情况需按需取模,详见代码)。

然后逻辑就较为简单了(尽管难调试):对没有改变的部分连乘,对有改变的部分单独乘即可。

#include <stdio.h>
#include <vector>
#include <algorithm>
#include <map>

typedef long long i64;
const i64 maxn = 50007;
i64 k, mod = 998244353, pos;
int n, m;

struct modint
{
    i64 x;
    modint(i64 x_ = 0);
    modint operator+=(modint b);
    modint operator*=(modint b);
};

modint operator+(modint a, modint b);
modint operator*(modint a, modint b);

struct matrix
{
    modint data[2][2];
    matrix() { data[0][0] = data[0][1] = data[1][0] = data[1][1] = 0; }
    modint *operator[](int index) { return data[index]; }
};
matrix operator*(matrix a, matrix b);

matrix intl, f[maxn][71];

i64 s[maxn], u, v, loopid, lastloop, pow2[maxn];
std::pair<i64, i64> spec[maxn];
std::map<i64, i64> mp;

i64 getsval(i64 x)
{
    auto it = mp.find(x);
    if (it == mp.end())
        return s[x % n];
    else
        return it->second;
}

matrix fib_trans(i64 si, i64 si1)
{
    matrix res;
    res[0][0] = si1;
    res[0][1] = 1;
    res[1][0] = si;
    res[1][1] = 0;
    return res;
}

matrix gfib(i64 index) { return fib_trans(getsval(index), getsval(index + 1)); }

void build()
{
    for (int i = 0; i < n; i++)
    {
        f[i][0][0][0] = s[(i + 1) % n];
        f[i][0][0][1] = 1;
        f[i][0][1][0] = s[i];
        f[i][0][1][1] = 0;
    }

    for (int j = 1; j < 63; j++)
        for (int i = 0; i < n; i++)
            f[i][j] = f[i][j - 1] * f[(i + (1ll << (j - 1))) % n][j - 1];
}

matrix to_binary(i64 &from, i64 to)
{
    matrix res;
    res[0][0] = res[1][1] = 1;

    if (from >= to)
        return res;

    for (i64 bit = 62; bit >= 0; bit--)
    {
        if (from + (1ll << bit) <= to)
            res = res * f[from % n][bit],
            from += 1ll << bit;
    }

    return res;
}

#define pos_next intl = intl * gfib(pos++)
#define break_if  \
    if (pos == k) \
        break;

int main()
{
    intl[0][0] = intl[1][1] = 1;

    scanf("%lld%lld", &k, &mod);
    if (k == 0)
        return printf("0") & 0;

    scanf("%d", &n);

    for (int i = 0; i < n; i++)
        scanf("%lld", s + i);
    scanf("%d", &m);
    for (int i = 0; i < m; i++)
        scanf("%lld%lld", &u, &v),
            mp[u] = v,
            spec[i] = std::make_pair(u, v);
    std::sort(spec, spec + m);
    spec[m].first = 1000000000000000007ll;
    m++;

    build();

    for (int i = 0; i < m; i++)
    {
        if (k <= spec[i].first)
        {
            intl = intl * to_binary(pos, k);
            break;
        }
        else
        {
            intl = intl * to_binary(pos, spec[i].first - 1);
            if (spec[i].first + 1 == spec[i + 1].first)
            {
                break_if;
                pos_next;
            }
            else
            {
                break_if;
                pos_next;
                break_if;
                pos_next;
                break_if;
            }
        }
    }

    printf("%lld", intl[0][1].x);
    return 0;
}


在 Rickyxrc's blog 出现的文章,若无特殊注明,均采用 CC BY-NC-SA 4.0 协议共享,也就是转载时需要注明本文章的地址,并且引用本文章的文章也要使用相同的方式共享。