暑假第二周笔记

· · 个人记录

Week2笔记

本周学习了两个知识点：分块思想+莫队与矩阵乘法

Part 1分块思想

1.分块的概念

何为分块？个人对于它的理解就是加了稍微优化的暴力。那这种思想能干啥呢？它是基于根号思想维护"区间"之类的问题。但我们回顾"区间"问题，我们学过了暴力、树状数组、线段树等众多算法。比如给定一个长度为 n 的数列，做 m 次区间修改和区间查询，每次只涉及部分区间。暴力法只是简单的从整体上进行修改和查询，复杂度为 O(mn) ，效率很差，。树状数组和线段树都用到二分的思想，以O(\log_{2}{n}) 的复杂度组织数据结构，每次只处理涉及的区间，从而实现了O(m\log_{2}{n}) 的高效复杂度。

虽然暴力法只能解决小规模的问题，但是它的代码非常简单。

有一种代码比树状数组和线段树简单，效率比暴力法高的算法，称为"分块" ，它能以 O(m\sqrt{n}) 的复杂度解决区间修改+区间查询的问题。简单的说，分块是线段树的分区思想改良的暴力法，它能把数列分成很多"块"，对涉及的块做整体性的维护操作（类似于线段树的Lazy-Tag）,而不是像普通暴力法那样处理整个数列，从而提高了效率

2.分块的操作

分块的基本操作如下

  1. 块的大小用block表示
  2. 块的数量用t表示
  3. 块的左右边界。定义数组st[]和ed[],用st[i]和ed[i]表示块i的第i个和最后一个元素的位置，有st[1]=1,ed[1]=block;st[2]=block+1,ed[2]=2\times block;...;st[i]=(i-1) \times block+1,ed[i]=i\times block...
  4. 每个元素所属的块。定义pos[],pos[i]表示第i个元素所在的块，pos[i]=(i-1)\div block+1

可以结合这张图进行理解

分块的实现

具体看下面的代码，其中block的值等于\sqrt{n}取整，因为此时的时间复杂度是最优的。如果\sqrt{n}的结果不是整数，那么最后是一个较小的碎块，代码中最重要的就是处理这个问题。

/*
分块
*/
int block=sqrt(n);
int t=n/block;
if(n%block)t++;
for(int i=1;i<=t;i++){
    st[i]=(i-1)*block+1;
    ed[i]=i*block;
}
ed[t]=n;
for(int i=1;i<=n;i++)
    pos[i]=(i-1)/block+1;
for(int i=1;i<=t;i++){
    for(int j=st[i];j<=ed[i];j++)
        sum[i]+=a[j];
}
void change(int L,int R,int d){
    int p=pos[L],q=pos[R];
    if(p==q){
        for(int i=L;i<=R;i++)a[i]+=d;
        sum[p]+=d*(R-L+1);
    }
    else{
        for(int i=p+1;i<=q-1;i++)add[i]+=d;
        for(int i=L;i<=ed[p];i++)a[i]+=d;
        sum[p]+=d*(ed[p]-L+1);
        for(int i=st[q];i<=R;i++)a[i]+=d;
        sum[q]+=d*(R-st[q]+1);
    }
}
long long ask(int L,int R){
    int p=pos[L],q=pos[R];
    long long ans=0;
    if(p==q){
        for(int i=L;i<=R;i++)ans+=a[i];
        ans+=add[p]*(R-L+1);
    }
    else{
        for(int i=p+1;i<=q-1;i++)ans+=sum[i]+add[i]*(ed[i]-st[i]+1);
        for(int i=L;i<=ed[p];i++)ans+=a[i];
        ans+=add[p]*(ed[p]-L+1);
        for(int i=st[q];i<=R;i++)ans+=a[i];
        ans+=add[q]*(R-st[q]+1);
    }
    return ans;
}

4.块长与时间效率：

修改查询操作的时间复杂度都是O(N/S+S)级别。根据均值不等式，块长S的值取 \sqrt{n}时总时间复杂度最低，所以一般默认情况下，块大小都取\sqrt{n}。

每次操作的时间复杂度为O(\sqrt{n})级别。m次操作总的时间复杂度为O(m\sqrt{n})

5.分块的感想

分块算法代码简单粗暴，没有复杂数据结构和复杂逻辑，很容易编码。它的思想可以概括为 "整体打包维护，碎片逐个枚举"

Part 2莫队

1.普通莫队

形式

假设 n=m，那么对于序列上的区间询问问题，如果从 [l,r] 的答案能够 O(1) 扩展到 [l-1,r],[l+1,r],[l,r+1],[l,r-1]（即与 [l,r] 相邻的区间）的答案，那么可以在 O(n\sqrt{n}) 的复杂度内求出所有询问的答案。

解释

离线后排序，顺序处理每个询问，暴力从上一个区间的答案转移到下一个区间答案（一步一步移动即可）。

排序方法

对于区间 [l,r], 以 l 所在块的编号为第一关键字，r 为第二关键字从小到大排序。

实现

void move(int pos, int sign) {
  // update nowAns
}

void solve() {
  BLOCK_SIZE = int(ceil(pow(n, 0.5)));
  sort(querys, querys + m);
  for (int i = 0; i < m; ++i) {
    const query &q = querys[i];
    while (l > q.l) move(--l, 1);
    while (r < q.r) move(r++, 1);
    while (l < q.l) move(l++, -1);
    while (r > q.r) move(--r, -1);
    ans[q.id] = nowAns;
  }
}

2.带修莫队

特点

普通莫队是不能带修改的。

我们可以强行让它可以修改，就像 DP 一样，可以强行加上一维 时间维, 表示这次操作的时间。

时间维表示经历的修改次数。

即把询问 [l,r] 变成 [l,r,\text{time}]。

那么我们的坐标也可以在时间维上移动，即 [l,r,\text{time}] 多了一维可以移动的方向，可以变成：

这样的转移也是 O(1) 的，但是我们排序又多了一个关键字，再搞搞就行了。

可以用和普通莫队类似的方法排序转移，做到 O(n^{5/3})。

这一次我们排序的方式是以 n^{2/3} 为一块，分成了 n^{1/3} 块，第一关键字是左端点所在块，第二关键字是右端点所在块，第三关键字是时间。

过程

先考虑普通莫队的做法：

现在再来考虑修改：

怎么强行加上一个修改呢？假设一个修改是修改第 pos 个位置上的颜色，原本 pos 上的颜色为 a，修改后颜色为 b，还假设当前莫队的区间扩展到了 [l,r]。

实现

#include <bits/stdc++.h>
#define SZ (10005)
using namespace std;

template <typename _Tp>
void IN(_Tp& dig) {
  char c;
  dig = 0;
  while (c = getchar(), !isdigit(c))
    ;
  while (isdigit(c)) dig = dig * 10 + c - '0', c = getchar();
}

int n, m, sqn, c[SZ], ct[SZ], c1, c2, mem[SZ][3], ans, tot[1000005], nal[SZ];

struct query {
  int l, r, i, c;

  bool operator<(const query another) const {
    if (l / sqn == another.l / sqn) {
      if (r / sqn == another.r / sqn) return i < another.i;
      return r < another.r;
    }
    return l < another.l;
  }
} Q[SZ];

void add(int a) {
  if (!tot[a]) ans++;
  tot[a]++;
}

void del(int a) {
  tot[a]--;
  if (!tot[a]) ans--;
}

char opt[10];

int main() {
  IN(n), IN(m), sqn = pow(n, (double)2 / (double)3);
  for (int i = 1; i <= n; i++) IN(c[i]), ct[i] = c[i];
  for (int i = 1, a, b; i <= m; i++)
    if (scanf("%s", opt), IN(a), IN(b), opt[0] == 'Q')
      Q[c1].l = a, Q[c1].r = b, Q[c1].i = c1, Q[c1].c = c2, c1++;
    else
      mem[c2][0] = a, mem[c2][1] = ct[a], mem[c2][2] = ct[a] = b, c2++;
  sort(Q, Q + c1), add(c[1]);
  int l = 1, r = 1, lst = 0;
  for (int i = 0; i < c1; i++) {
    for (; lst < Q[i].c; lst++) {
      if (l <= mem[lst][0] && mem[lst][0] <= r)
        del(mem[lst][1]), add(mem[lst][2]);
      c[mem[lst][0]] = mem[lst][2];
    }
    for (; lst > Q[i].c; lst--) {
      if (l <= mem[lst - 1][0] && mem[lst - 1][0] <= r)
        del(mem[lst - 1][2]), add(mem[lst - 1][1]);
      c[mem[lst - 1][0]] = mem[lst - 1][1];
    }
    for (++r; r <= Q[i].r; r++) add(c[r]);
    for (--r; r > Q[i].r; r--) del(c[r]);
    for (--l; l >= Q[i].l; l--) add(c[l]);
    for (++l; l < Q[i].l; l++) del(c[l]);
    nal[Q[i].i] = ans;
  }
  for (int i = 0; i < c1; i++) printf("%d\n", nal[i]);
  return 0;
}

Part3 矩阵乘法

矩阵乘法

矩阵的乘法是向量内积的推广。

矩阵相乘只有在第一个矩阵的列数和第二个矩阵的行数相同时才有意义。

设 A 为 P \times M 的矩阵，B 为 M \times Q 的矩阵，设矩阵 C 为矩阵 A 与 B 的乘积，

其中矩阵 C 中的第 i 行第 j 列元素可以表示为：

C_{i,j} = \sum_{k=1}^MA_{i,k}B_{k,j}

在矩阵乘法中，结果 C 矩阵的第 i 行第 j 列的数，就是由矩阵 A 第 i 行 M 个数与矩阵 B 第 j 列 M 个数分别 相乘再相加 得到的。这里的 相乘再相加，就是向量的内积。乘积矩阵中第 i 行第 j 列的数恰好是乘数矩阵 A 第 i 个行向量与乘数矩阵 B 第 j 个列向量的内积，口诀为 左行右列。

线性代数研究的向量多为列向量，根据这样的对矩阵乘法的定义方法，经常研究对列向量左乘一个矩阵的左乘运算，同时也可以在这里看出「打包处理」的思想，同时处理很多个向量内积。

矩阵乘法满足结合律，不满足一般的交换律。

利用结合律，矩阵乘法可以利用 快速幂 的思想来优化。

在比赛中，由于线性递推式可以表示成矩阵乘法的形式，也通常用矩阵快速幂来求线性递推数列的某一项。

优化

首先对于比较小的矩阵，可以考虑直接手动展开循环以减小常数。

可以重新排列循环以提高空间局部性，这样的优化不会改变矩阵乘法的时间复杂度，但是会在得到常数级别的提升。矩阵加速递推

应用

斐波那契数列大家应该都非常的熟悉了。在斐波那契数列当中，F_1 = F_2 = 1，F_i = F_{i - 1} + F_{i - 2}(i \geq 3)。

如果有一道题目让你求斐波那契数列第 n 项的值，最简单的方法莫过于直接递推了。但是如果 n 的范围达到了 10^{18} 级别，递推就不行了，稳TLE。考虑矩阵加速递推。

设 Fib(n) 表示一个 1 \times 2 的矩阵 \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]。我们希望根据 Fib(n-1)=\left[ \begin{array}{ccc}F_{n-1} & F_{n-2} \end{array}\right] 推出 Fib(n)。

试推导一个矩阵 \text{base}，使 Fib(n-1) \times \text{base} = Fib(n)，即 \left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \text{base} = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]。

怎么推呢？因为 F_n=F_{n-1}+F_{n-2}，所以 \text{base} 矩阵第一列应该是 \left[\begin{array}{ccc} 1 \\ 1 \end{array}\right]，这样在进行矩阵乘法运算的时候才能令 F_{n-1} 与 F_{n-2} 相加，从而得出 F_n。同理，为了得出 F_{n-1}，矩阵 \text{base} 的第二列应该为 \left[\begin{array}{ccc} 1 \\ 0 \end{array}\right]。

综上所述：\text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right] 原式化为 \left[\begin{array}{ccc}F_{n-1} & F_{n-2}\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right] = \left[ \begin{array}{ccc}F_n & F_{n-1} \end{array}\right]

转化为代码，应该怎么求呢？

定义初始矩阵 \text{ans} = \left[\begin{array}{ccc}F_2 & F_1\end{array}\right] = \left[\begin{array}{ccc}1 & 1\end{array}\right], \text{base} = \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]。那么，F_n 就等于 \text{ans} \times \text{base}^{n-2} 这个矩阵的第一行第一列元素，也就是 \left[\begin{array}{ccc}1 & 1\end{array}\right] \times \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-2} 的第一行第一列元素。

注意，矩阵乘法不满足交换律，所以一定不能写成 \left[\begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array}\right]^{n-2} \times \left[\begin{array}{ccc}1 & 1\end{array}\right] 的第一行第一列元素。另外，对于 n \leq 2 的情况，直接输出 1 即可，不需要执行矩阵快速幂。

为什么要乘上 \text{base} 矩阵的 n-2 次方而不是 n 次方呢？因为 F_1, F_2 是不需要进行矩阵乘法就能求的。也就是说，如果只进行一次乘法，就已经求出 F_3 了。如果还不是很理解为什么幂是 n-2，建议手算一下。

下面是求斐波那契数列第 n 项对 10^9+7 取模的示例代码（核心部分）。

const int mod = 1000000007;

struct Matrix {
  int a[3][3];

  Matrix() { memset(a, 0, sizeof a); }

  Matrix operator*(const Matrix &b) const {
    Matrix res;
    for (int i = 1; i <= 2; ++i)
      for (int j = 1; j <= 2; ++j)
        for (int k = 1; k <= 2; ++k)
          res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j]) % mod;
    return res;
  }
} ans, base;

void init() {
  base.a[1][1] = base.a[1][2] = base.a[2][1] = 1;
  ans.a[1][1] = ans.a[1][2] = 1;
}

void qpow(int b) {
  while (b) {
    if (b & 1) ans = ans * base;
    base = base * base;
    b >>= 1;
  }
}

int main() {
  int n = read();
  if (n <= 2) return puts("1"), 0;
  init();
  qpow(n - 2);
  println(ans.a[1][1] % mod);
}

这是一个稍微复杂一些的例子。

\begin{gathered} f_{1} = f_{2} = 0\\ f_{n} = 7f_{n-1}+6f_{n-2}+5n+4\times 3^n \end{gathered}

我们发现，f_n 和 f_{n-1}, f_{n-2}, n 有关，于是考虑构造一个矩阵描述状态。

但是发现如果矩阵仅有这三个元素 \begin{bmatrix}f_n& f_{n-1}& n\end{bmatrix} 是难以构造出转移方程的，因为乘方运算和 +1 无法用矩阵描述。

于是考虑构造一个更大的矩阵。

\begin{bmatrix}f_n& f_{n-1}& n& 3^n & 1\end{bmatrix}

我们希望构造一个递推矩阵可以转移到

\begin{bmatrix} f_{n+1}& f_{n}& n+1& 3^{n+1} & 1 \end{bmatrix}

转移矩阵即为

\begin{bmatrix} 7 & 1 & 0 & 0 & 0\\ 6 & 0 & 0 & 0 & 0\\ 5 & 0 & 1 & 0 & 0\\ 12 & 0 & 0 & 3 & 0\\ 5 & 0 & 1 & 0 & 1 \end{bmatrix}