题解-CF575A

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

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

题面大意

给你一个数列 $S'$，你需要求出数列 $F$ 的第 $k$ 项。

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

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

解题思路

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

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

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

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

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

我们定义 $P_{i,j}$ 为从 $S_i$ 开始连乘，总长度为 $2^j$ 的矩阵。

那么我们的转移公式就应该是 $P_{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;
}