浅谈矩阵乘法维护线段树标记的技巧

George_Plover · 2022-01-16 23:11:07 · 算法·理论

一、引述

懒标记在一类数据结构问题中的应用是非常普遍的。在数据结构的简单运用中，懒标记的维护通常非常简洁。例如经典的线段树区间加法与区间求和操作，通过懒标记维护区间未向下处理的加法信息，这样在修改的时候不必向下递归，等到查询向下递归的时候再同步下传标记。这个简单的例子中，我们只需要简单维护一个加法标记就可以了。

然而，在近期的竞赛中，这样一类只需维护单个标记的基础题目已经十分少见了。以线段树为例，在提高题目难度的时候，除了考察线段树的变种，还有一种方式就是增加需要支持的操作种类，使得线段树需要维护两个或者多个不同的懒标记。当标记大于一个的时候，我们就需要思考不同标记在下传过程中的相互影响，然而这会极大地增加临场的思考时间。对此，掌握一些线性代数方面的技巧，通过构造向量和矩阵乘法，能够更迅速准确地抓住维护这些标记的方法。

本文将从一些题目入手，以线段树为例，介绍几种维护区间线性变换的问题场景，谈一谈利用矩阵维护懒惰标记的技巧。

二、技巧分析

（〇）区间加乘线段树

在进入主题之前，不妨先来回忆一个最基础的双标记模型：区间加法乘法线段树。

在这个题目中，给出了一个长度为 n 的整数序列 a_1,a_2,\cdots a_n 的序列，要求支持以下操作：

我们都知道，这道题目可以通过维护乘法标记 mul_tag 和加法标记 add_tag 来实现。对于一个暂时没有被下传更新的区间和 a 来说，它的真实值 a^\prime = mul\_tag \cdot a\ +\ add\_tag\cdot len ，（ len 表示区间长度）因此我们在下传的过程中，需要先下传乘法的标记，再下传加法标记，才能够保证正确性。这样两个标记的维护方式是如何设计出来的呢？我们考虑一下线性代数的分析方式。

（1）区间信息的向量表示

我们假设，线段树的每一个叶子结点，维护的不仅仅是 a_i 这个数值，而是一个向量：

\vec{a}= \begin{bmatrix} a_i \\ 1 \end{bmatrix}

（这里向量中的 a_i 即表示题目中的 a_i，向量中的 1 是常数 1）那么对于线段树的非叶子结点 u ，我们定义其维护的向量 \vec {a}_u = \vec {a}_l+\vec {a}_r 。其中 l,r 分别表示结点 u 的左右儿子结点。

这样的定义下，结点向量 \vec a 有着这样的含义：它的第一位维度的值，是其子树中所有叶子的 a_i 的和（sum）；它的第二维度的值，是一个常数，等于其子树中叶子结点的个数，也就是该结点表示的线段区间长度（len）。

接下来，我们考虑对一个结点进行区间的加法操作，不妨设该结点上，\vec a =\begin{bmatrix} sum \\ len \end{bmatrix}。我们对这个结点所表示的区间的每个数字加上常数 c ，可以知到，操作之后有：\vec a^\prime =\begin{bmatrix} sum + len\cdot c \\ len \end{bmatrix} 。对此，我们不妨把这个操作，看作矩阵的乘法：

\vec a^\prime =M_c\vec a=\begin{bmatrix}1& c \\ 0 & 1 \end{bmatrix}\begin{bmatrix} sum \\ len \end{bmatrix}=\begin{bmatrix} sum + len\cdot c \\ len \end{bmatrix}

我们把 M_c=\begin{bmatrix}1& c \\ 0 & 1 \end{bmatrix} 叫做加法操作的操作矩阵。对一个线段树结点维护的区间，执行一次区间加法操作，其本质是对该结点的子树的所有结点上的向量都左乘一个矩阵。矩阵的乘法满足结合律，故而可以用懒标记维护。

接着，我们考虑对一个结点进行区间的乘法操作，不妨设该结点上，\vec a =\begin{bmatrix} sum \\ len \end{bmatrix}。我们对这个结点所表示的区间的每个数字乘上常数 k ，可以知道，操作之后有：\vec a^\prime =\begin{bmatrix} sum \cdot k \\ len \end{bmatrix} 。对此，我们不妨把这个操作，也写作矩阵的乘法：

\vec a^\prime =M_k\vec a=\begin{bmatrix}k& 0 \\ 0 & 1 \end{bmatrix}\begin{bmatrix} sum \\ len \end{bmatrix}=\begin{bmatrix} sum \cdot k \\ len \end{bmatrix}

我们把 M_c=\begin{bmatrix}k& 0 \\ 0 & 1 \end{bmatrix} 叫做乘法操作的操作矩阵。对一个线段树结点维护的区间，执行一次区间乘法操作，其本质是对该结点的子树的所有结点上的向量都左乘一个矩阵。与加法操作同样，可以用矩阵乘法的懒标记维护。

（2）向量表示的好处

上面的表示方法看起来似乎有些繁琐，用看起来笨重的矩阵乘法去替代了原来的操作。但是这样的维护方法，把两种原本看起来不同类型的操作，都化作了同一种类型：矩阵乘法。

于是，我们可以把每种操作，都看成区间矩阵乘法操作，我们只需要在每个结点维护一个标记即可（即，一个矩阵）。在标记下传的时候，只需要对子结点的标记和向量都左乘下传矩阵即可。

这样，我们可以用线段树在 O(\log n) 的复杂度维护这两种操作，只不过，二阶矩阵的一次乘法有 8 次运算的常数。

（3）从向量法回推双标记法

注意到上面我们提到的加法矩阵 M_c=\begin{bmatrix}1& c \\ 0 & 1 \end{bmatrix} 和乘法矩阵 M_k=\begin{bmatrix}k& 0 \\ 0 & 1 \end{bmatrix} 都是上三角矩阵，并且无论怎样相乘，左下角始终是 0 ，右下角始终是 1 。这意味着，在维护矩阵标记的时候，无论怎样相乘，标记矩阵一定可以写作这样的形式：

M_{tag}=\begin{bmatrix}x& y \\ 0 & 1 \end{bmatrix}

也就是说，我们只需要维护两个变量 x,y ，即可表示一个矩阵标记。其实，这个 x,y 分别就是我们的乘法标记和加法标记！我们来看 x,y 在矩阵乘法下传时的表现就可以发现这一点。

在标记下传的时候，用 M_{tag} 去左乘结点维护的向量 \vec a=\begin{bmatrix} sum\\ len \end{bmatrix} ，就有：

\begin{bmatrix}x& y \\ 0 & 1 \end{bmatrix}\begin{bmatrix}sum \\ len \end{bmatrix}=\begin{bmatrix}x\cdot sum + y\cdot len\\ len \end{bmatrix}

可以看到，第一维度更新后的形式，恰好就是上文提到的： a^\prime = mul\_tag \cdot a\ +\ add\_tag\cdot len 。正因为

同理，在标记下传的时候，会发生标记相乘的情况，我们考虑 M_{tag}M_{tag_0} ：

\begin{bmatrix}x& y \\ 0 & 1 \end{bmatrix}\begin{bmatrix}x_0& y_0 \\ 0 & 1 \end{bmatrix}=\begin{bmatrix}x\cdot x_0& x\cdot y_0+y\\0& 1 \end{bmatrix}

可以看到更新后的形式，和我们双标记法下传时采用的 “乘法标记相乘，加法标记先乘再加” 的维护方式是完全一致的！

综合来看，我们可以不必维护 2\times 2 的矩阵作为标记，因为我们已经推导出矩阵有效的变量仅有两个，即 x,y ，并且已经推导除了 x,y 这两个标记的维护方法。于是我们可以从矩阵回到标记法，这样能省掉矩阵计算时无意义的常数开销。直接使用的 x,y 标记，事实上就是我们传统做法的乘法、加法标记。

在下文中会给出该分析方法在更多问题上的运用。

（一）序列累加线段树

考虑这样一类模型：给出一个长度为 n 的整数序列 a_1,a_2,\cdots a_n，你需要维护以下类型的区间操作：

可以看到，这个模型在传统的区间修改查询的基础上，新增了“把当前序列的值加入到另一个数列”这样的操作。自然而然，可能会想到通过维护两种不同的懒标记来实现。比如这样的思路：用一个标记维护 a 序列的区间加法，再用第二个标记来维护 b 序列对 a_i 的累加。但实践起来，会发现因为 a 序列的变化，两种标记没有办法很好的同时维护。

不妨回想上面提到的方法：考虑用向量维护。

（1）向量表示

在线段树每个叶子结点维护一个三维向量 \vec a_i=\begin{bmatrix} b_i \\ a_i \\ 1 \end{bmatrix} 。

同样，对每个非叶子结点，我们在上面维护一个向量，等于左右儿子向量的和。于是线段树每个结点维护的三维向量可以写作如下形式：

\vec a=\begin{bmatrix} B \\ A \\ L \end{bmatrix}

其中， B 表示的是该结点子树中所有叶子的 b_i 值的和，A 表示的该结点子树中所有叶子的 a_i 值的和，L 表示的该结点子树中叶子的个数。

（2）操作矩阵化

我们来分析两种操作的对应的矩阵乘法。

首先，考虑区间加法操作，给区间所有 a_i 值加上一个常数 c ，那么只有 a 序列的区间总和会变化，而且增加了 c\cdot L 。因此这一操作等价于如下矩阵乘法：

M_c \vec a =\begin{bmatrix} 1&0&0 \\ 0&1&c \\ 0&0&1 \end{bmatrix}\begin{bmatrix} B \\ A \\ L \end{bmatrix}=\begin{bmatrix} B \\ A+c\cdot L \\ L \end{bmatrix}

然后，考虑把 a 序列的区间对应加到 b 序列的区间的操作，这给每个 b_i 加上了一个 a_i 。那么只有 b 序列的区间总和会变化，这一操作等价于如下矩阵乘法：

M_b \vec a =\begin{bmatrix} 1&1&0 \\ 0&1&0 \\ 0&0&1 \end{bmatrix}\begin{bmatrix} B \\ A \\ L \end{bmatrix}=\begin{bmatrix} B+A \\ A \\ L \end{bmatrix}

于是，我们把两种操作都转化为了三阶矩阵的左乘，可以用线段树维护矩阵标记来实现了。在查询答案的时候，我们只需要想经典线段树区间查询那样，找到 O(\log n) 个结点，将它们的向量中的对应值相加即可。

（3）从矩阵标记回推多标记

如果我们不能容忍矩阵乘法的 O(27) 的常数，我们可以对矩阵作出分析，两种矩阵在互相乘的时候，显然都是上三角矩阵，并且主对角线永远是 1 ，进一步分析得到矩阵标记的统一形式：

\begin{bmatrix} 1&x&y \\ 0&1&z \\ 0&0&1 \end{bmatrix}

也就是说，可以用三个标记 x,y,z 来表示一个矩阵标记。这些标记在下传的时候该如何维护，可以通过展开矩阵乘法的方式推导出来，这和上文的内容是一样的（推导过程留给读者）。

从这里我们可以看出，通过线性操作的方式来维护标记，至少需要三个标记，这也解释了为什么我们在尝试用仅两个标记解决该问题时会遇到困难。

从上面两个例子我们可以逐渐感受到：一些繁琐的序列操作，只要我们能够构造出一个向量，然后能够把所有的操作化为向量的线性变换，我们就可以把操作转变为矩阵乘法。无论题目给出了多少种不同的线性操作，我们只需要构造出对应的矩阵，就能简单地使用线段树维护矩阵乘法的方式进行维护。

（二）线段树维护二次项求和

考虑这样一类模型：给出一个长度为 n 的整数序列 a_1,a_2,\cdots a_n，你需要维护以下类型的区间操作：

该问题的难点在于如何维护 a_i^2 的区间和，这也是一类经典问题，我们用和上面同理的方法来分析并维护。

（1）向量表示

在线段树每个叶子结点维护一个三维向量 \vec a_i=\begin{bmatrix} a_i^2 \\ a_i \\ 1 \end{bmatrix} 。

同样，对每个非叶子结点，我们在上面维护一个向量，等于左右儿子向量的和。于是线段树每个结点维护的三维向量可以写作如下形式：

\vec a=\begin{bmatrix} A_2 \\ A \\ L \end{bmatrix}

其中， A 表示的是该结点子树中所有叶子的 a_i 值的和，A_2 表示的该结点子树中所有叶子的 a_i^2 值的和，L 表示的该结点子树中叶子的个数。

（2）操作矩阵化

我们来分析两种操作的对应的矩阵乘法。

首先，考虑区间加法操作，给区间所有 a_i 值加上一个常数 c ，容易发现 (a_i+c)^2=a_i^2+2c\cdot a_i+c^2 ，对应的就有 A_2^\prime=A_2 + 2c\cdot A +c^2\cdot L 。因此这一操作等价于如下矩阵乘法：

M_c \vec a =\begin{bmatrix} 1&2c&c^2 \\ 0&1&c \\ 0&0&1 \end{bmatrix}\begin{bmatrix} A_2 \\ A \\ L \end{bmatrix}=\begin{bmatrix} A_2 + 2c\cdot A +c^2\cdot L \\ A+c\cdot L \\ L \end{bmatrix}

然后，考虑区间乘法操作，不难求出这一操作等价于如下矩阵乘法：

M_k \vec a =\begin{bmatrix} c^2&0&0 \\ 0&c&0 \\ 0&0&1 \end{bmatrix}\begin{bmatrix} A_2 \\ A \\ L \end{bmatrix}=\begin{bmatrix} c^2A_2 \\ cA \\ L \end{bmatrix}

于是，我们把两种操作都转化为了三阶矩阵的左乘。注意到矩阵标记都可以写作 \begin{bmatrix} x_1&x_2&x_3 \\ 0&x_4&x_5 \\ 0&0&1 \end{bmatrix} 的形式，因此使用五个标记代替矩阵，可以优化常数。值得注意的是，这种方法也可以去维护序列的 k 次方求和，不过需要构造 k+1 阶矩阵，单次操作复杂度 O(k^3\log n)。

从上面的三个部分，我们可以看到，设计良好的向量，可以把多样的操作矩乘化，这是一种实用的技巧，下面我们将从几个题目来分析并巩固这种技巧。

三、题目分析

（1）【ICPC-2021 南京】Paimon Segment Tree

（gym/103470/E）

【题意】

给出一个长度为 n 的序列 A=a_1,a_2\cdots a_n 。我们把 A_0=a_{0,1},a_{0,2}\cdots a_{0,n} 称为 A 序列的初始版本。

接下来会进行 m 次操作，其中第 i 次操作形如 (l_i,r_i,c_i) ，是把当前最新版本的 A_{i-1} 进行修改，把区间 [l_i,r_i] 范围的变量分别加上 c_i，得到新的序列版本 A_{i} 。形式化的，有：\begin{cases} a_{i,j}=a_{i-1,j}+c_i,&l_i\le j \le r_i \\ a_{i,j}=a_{i-1,j},&others \end{cases} 。

然后，会有 q 次查询操作 (x,y,l,r) ，询问版本 A_x,A_{x+1}\cdots A_y 中下标在 [l,r] 范围的变量的平方和。形式化的，查询：

\sum_{i=x}^y\sum_{j=l}^r a_{i,j}^2

n,m,q\le 5\times 10^4

【分析】

这道题目在这场比赛中只有 29/680 支队伍通过，属于金牌题的范畴。然而，它只是序列累加线段树和二次求和线段树的综合运用。

首先我们把询问离线，几个连续版本的变量求和，我们可以转化为前缀版本和的差分。前缀版本和，其实本质就是，在每次修改 A 序列之后，都把当前的 A 序列对应加在 S 序列上，这样 S 就是版本前缀和序列了。

具体而言，对于每个询问 (x,y,l,r) ，我们把它拆分成 (x-1,l,r) 和 (y,l,r) 两个询问部分，单个询问部分 (y,l,r) 询问的是 y 及其之前的版本序列的 [l,r] 区间上的值的平方和。两个部分最后作差即是询问的答案。将所有的询问部分按照第一关键字排序，当我们求出第 u 个版本的序列 A_u 的时候，就把序列 A_u 累加在 S 上，这个时候处理所有形如 (u,l,r) 的询问即可。

构造叶子结点维护的四维向量：\vec a_i=\begin{bmatrix} s_i \\ a_i^2 \\ a_i \\ 1 \end{bmatrix}，定义非叶子结点向量为子节点向量和，于是有线段树结点的向量表达式 \vec a_i=\begin{bmatrix} S\\A_2 \\ A \\ L \end{bmatrix} 。

于是区间加 c 操作等价于左乘 \begin{bmatrix} 1&0&0&0\\0&1&2c&c^2 \\ 0&0&1&c \\ 0&0&0&1 \end{bmatrix} ，版本累加操作等价于左乘\begin{bmatrix} 1&1&0&0\\0&1&0&0 \\ 0&0&1&0 \\ 0&0&0&1 \end{bmatrix} 。

直接维护矩阵乘法常数较大，总时间复杂度为 O(64n\log n)。可以考虑展开写有效部分的矩阵乘法（毕竟是上三角矩阵而且主对角线全是 1 ），常数更为优秀。

【参考代码】

见附录 code1.cpp

（2）【HDU多校-2021 第二场 1007】I love data structure

（记不清题号了）

【题意】

给出两个长度为 n 的整数序列 A=a_1,a_2\cdots a_n,\ B=b_1,b_2,\cdots b_n 。

接下来会进行 m 次操作，每种操作是以下类型之一：

n,m\le 2\times 10^5

【分析】

这道题目更好地呈现了线性变换操作。无论是交换两个数，还是加上常数，亦或是直接的线性变换。这些都可以看作矩阵的乘法。对两个数的向量而言，交换两个数其实就是左乘 \begin{bmatrix} 0&1 \\1&0 \end{bmatrix} ，题目中第3种操作其实就是 \begin{bmatrix} 3&2 \\3&-2 \end{bmatrix}\begin{bmatrix} a_i \\b_i \end{bmatrix} 。

不过查询的值是 a_ib_i 的求和，因此我们构造的向量还需要一些额外的部分，比如 a_ib_i 。注意到，在执行一次第3种操作后，a_ib_i 的值变成了 (3a_i+2b_i)(3a_i-2b_i)=9a_i^2-4b_i^2 ，也就是说一次操作后，a_ib_i 可能变成了 a_i^2 与 b_i^2 的线性组合。所以我们的向量还必须维护 a_i^2 与 b_i^2 。

具体而言，我们可以构造叶子结点的六维向量：

\text{叶子结点：}\vec a_i=\begin{bmatrix} a_ib_i \\ b_i^2 \\ a_i^2\\b_i\\a_i \\ 1 \end{bmatrix}\ \ \ \text{结点统一表达:} \vec a_i=\begin{bmatrix} M \\ B_2 \\ A_2\\B\\A \\ L \end{bmatrix}

下面直接对每种操作给出矩阵左乘：

\begin{matrix} \begin{bmatrix} 1&0&0&x&0&0 \\ 0&1&0&0&0&0 \\ 0&0&1&0&2x&x^2\\ 0&0&0&1&0&0\\ 0&0&0&0&1&x \\ 0&0&0&0&0&1 \end{bmatrix} \begin{bmatrix} M \\ B_2 \\ A_2\\B\\A \\ L \end{bmatrix} & \begin{bmatrix} 1&0&0&0&0&0 \\ 0&0&1&0&0&0 \\ 0&1&0&0&0&0\\ 0&0&0&0&1&0\\ 0&0&0&1&0&0 \\ 0&0&0&0&0&1 \end{bmatrix} \begin{bmatrix} M \\ B_2 \\ A_2\\B\\A \\ L \end{bmatrix} & \begin{bmatrix} 0&-4&9&0&0&0 \\ -12&4&9&0&0&0 \\ 12&4&9&0&0&0\\ 0&0&0&-2&3&0\\ 0&0&0&2&3&0 \\ 0&0&0&0&0&1 \end{bmatrix} \begin{bmatrix} M \\ B_2 \\ A_2\\B\\A \\ L \end{bmatrix} \\ & & \\ a_i \ \text{add}\ x & \text{swap } a_i,b_i & (a_i,b_i)\to(3a_i+2b_i,3a_i-2b_i) \end{matrix}

直接维护六阶矩阵，在乘法运算的时候具有 O(216) 的大常数，因此要通过此题还需要对这些稀疏矩阵做处理，例如，可以注意到矩阵左下角的 3\times 3 区域是恒为 0 的。

由于矩阵过于庞大，分析其中的有效部分还是非常繁琐的。有一种更好的分析方法：ab,a^2,b^2 的线性变换都是因 a,b 线性变换而起的，因此实际上我们可以只维护 3 阶矩阵来表示 \begin{bmatrix}B\\A\\L\end{bmatrix} 的变换标记，当需要维护 6 维向量向量的时候，根据该 3 阶矩阵就可以临时推导出对应的 6 阶矩阵。这样做，在维护矩阵乘法的时候，只有 O(27) 的常数，在维护 6 维向量的时候，只有 O(36) 的常数。因此复杂度可以看作 O(36n\log n) 。

【参考代码】

（因HDU无法访问）略

四、小结

当遇到操作数量繁多且复杂的数据结构操作时，需要观察其性质。在合适的向量构造下，如果所有的操作都是线性的，我们能够把繁杂的操作统一成矩阵的乘法操作，进而解决。

当然，在初步分析的时候，我们可能会得到一个较大的矩阵，为了优化常数，我们还需要对矩阵做一些处理，从中提取有效的、必要的成分，从而减少维护矩阵乘法的代价，降低常数。

五、附录

code1.cpp

#include <cstdio>//出于个人习惯，本代码用的是： 行向量表示，维护矩阵右乘
#include <cstring>
#include <iostream>
#include <algorithm>
#define MAXN 51000
#define MOD 1000000007
using namespace std;
int n,m,q,b[MAXN];
struct Mat{
    int a[5][5];
    void gen(int k)//初始化加法操作矩阵，参数为k
    {
        memset(a,0,sizeof(a));
        a[1][2]=k;
        a[2][3]=a[2][4]=(k+k)%MOD;
        a[1][1]=a[2][2]=a[3][3]=a[4][4]=a[3][4]=1;
        a[1][3]=a[1][4]=1ll*k*k%MOD;
    }
    void I(){//初始化单位矩阵
        clear();a[1][1]=a[2][2]=a[3][3]=a[4][4]=1;
    }
    bool is_I()//判断是否为单位矩阵
    {
        return a[1][2]==0 && a[1][3] == 0 && a[1][4]==0 &&a[2][3]==0 && a[2][4]==0 && a[3][4]==0;
    }
    void clear(){
        memset(a,0,sizeof(a));
    }
};
Mat operator * (Mat A,Mat B)//展开实现的矩阵乘法
{
    Mat C;
    C.clear();
    C.a[1][1] = C.a[2][2]=C.a[3][3]=C.a[4][4]=1;
    C.a[1][2] = (B.a[1][2] + A.a[1][2])%MOD;
    C.a[1][3] = (B.a[1][3] + 1ll*A.a[1][2]*B.a[2][3] + A.a[1][3])%MOD;
    C.a[1][4] = (B.a[1][4] + 1ll*A.a[1][2]*B.a[2][4] + 1ll*A.a[1][3]*B.a[3][4] + A.a[1][4])%MOD;
    C.a[2][3] = (0 + B.a[2][3] + A.a[2][3]) %MOD;
    C.a[2][4] = (0 + B.a[2][4] + 1ll*A.a[2][3]*B.a[3][4] + A.a[2][4])%MOD;
    C.a[3][4] = (B.a[3][4] + A.a[3][4])%MOD;
    return C;
}
struct Vec{//维护的向量
    int a[5];
    Vec(){memset(a,0,sizeof(a));}
};
Vec operator +(Vec A, Vec B)//向量加法
{
    for(int i=1;i<=4;i++)
        A.a[i]=(A.a[i]+B.a[i])%MOD;
    return A;
}
struct Segment_Tree//线段树
{
    struct node{
        int l,r;
        Mat tag;
        Vec A;
        node(){l=r=0;tag.clear();}
    }tr[MAXN*4];
    int num;
    Segment_Tree(){
        num=1;
    }
    void modify_node(int x,Mat &M)//展开实现的向量左乘
    {
        tr[x].tag = tr[x].tag * M;
        Vec tmp;
        tmp.a[1] = tr[x].A.a[1];
        tmp.a[2] = (1ll * tr[x].A.a[1] * M.a[1][2] + tr[x].A.a[2])%MOD;
        tmp.a[3] = (1ll * tr[x].A.a[1] * M.a[1][3] + 1ll * tr[x].A.a[2] * M.a[2][3] + tr[x].A.a[3])%MOD;
        tmp.a[4] = (1ll * tr[x].A.a[1] * M.a[1][4] + 1ll * tr[x].A.a[2] * M.a[2][4] + 1ll * tr[x].A.a[3] * M.a[3][4] + tr[x].A.a[4])%MOD;
        tr[x].A = tmp;
    }
    void pushdown(int x)//下传标记
    {
        if(!tr[x].tag.is_I())
        {
            modify_node(tr[x].l,tr[x].tag);
            modify_node(tr[x].r,tr[x].tag);
            tr[x].tag.I();
        }
    }
    void pushup(int x)
    {
        tr[x].A = tr[tr[x].l].A+tr[tr[x].r].A;
    }
    void Build(int x,int L,int R)//建树操作，初始化叶子的向量 
    {
        tr[x].tag.I();
        if(L==R)
        {
            tr[x].A.a[1]=1;
            tr[x].A.a[2]=b[L];
            tr[x].A.a[4]=tr[x].A.a[3]=1ll*b[L]*b[L]%MOD;
            return;
        }
        int mid=(L+R)/2;
        tr[x].l=++num;
        Build(num,L,mid);
        tr[x].r=++num;
        Build(num,mid+1,R);
        pushup(x);//pushup初始化非叶子的向量
    }

    void Update(int x,int L,int R,int al,int ar,Mat &M)//区间矩阵乘法
    {
        if(R<al||ar<L)return;
        if(al<=L&&R<=ar)
        {
            modify_node(x,M);
            return;
        }
        int mid=(L+R)/2;
        pushdown(x);
        Update(tr[x].l,L,mid,al,ar,M);
        Update(tr[x].r,mid+1,R,al,ar,M);
        pushup(x);
    }

    int Query(int x,int L,int R,int al,int ar)//查询版本前缀和
    {
        if(R<al||ar<L)return 0;
        if(al<=L&&R<=ar)
            return tr[x].A.a[4];
        pushdown(x);
        int mid=(L+R)/2;
        return (Query(tr[x].l,L,mid,al,ar)+Query(tr[x].r,mid+1,R,al,ar))%MOD;
    }
}seg;

int ans[MAXN];//对询问离线的一些辅助变量
struct Op{
    int l,r,k;
}op[MAXN];
struct Q{
    int t,l,r;
    bool is_l;
    int belong;
}p[MAXN*2];
bool cmp(Q x,Q y){
    return x.t<y.t;
}
int cnt;

int main()
{
    scanf("%d%d%d",&n,&m,&q);
    for(int i=1;i<=n;i++)
    {
        scanf("%d",&b[i]);
        b[i]=(b[i]+MOD)%MOD;
    }
    seg.Build(1,1,n);
    for(int i=1;i<=m;i++)
    {
        scanf("%d%d%d",&op[i].l,&op[i].r,&op[i].k);
        if(op[i].k<0)op[i].k=(op[i].k+MOD)%MOD;
    }
    for(int i=1;i<=q;i++)
    {
        int x,y,l,r;
        scanf("%d%d%d%d",&l,&r,&x,&y);
        p[++cnt].belong=i;
        p[cnt].l=l;p[cnt].r=r;
        p[cnt].t=x-1;
        p[cnt].is_l=1;
        if(x==0)cnt--;

        p[++cnt].belong=i;
        p[cnt].l=l;p[cnt].r=r;
        p[cnt].is_l=0;
        p[cnt].t=y;
    }
    sort(p+1,p+cnt+1,cmp);
    for(int i=0,j=1;i<=m;i++)
    {
        if(i>0)
        {
            Mat tmp;
            Mat nsh;nsh.gen(0);
            tmp.gen(op[i].k);
            seg.Update(1,1,n,op[i].l,op[i].r,tmp);
            if(op[i].l>1)
                seg.Update(1,1,n,1,op[i].l-1,nsh);
            if(op[i].r<n)
                seg.Update(1,1,n,op[i].r+1,n,nsh);

        }
        while(j<=cnt && p[j].t==i)
        {
            int tmp = seg.Query(1,1,n,p[j].l,p[j].r);
            int u=p[j].belong;
            if(p[j].is_l)
                ans[u]=(ans[u]-tmp+MOD)%MOD;
            else
                ans[u]=(ans[u]+tmp)%MOD;
            j++;
        }
    }
    for(int i=1;i<=q;i++)
        printf("%d\n",ans[i]);
    return 0;
}