分类目录归档:算法

线性代数相关OJ题选讲

这篇文章是关于线性代数的OJ题选讲,这几道题目都不难,主要是考察思维的灵活程度。


Fast Matrix Calculation

2014 Multi-University Training Contest 9, available at HDUOJ 4965

给定 \(N\times k\) 的矩阵 \(A\),\(k \times N\) 的矩阵 \(B\) \((2 \leq k \leq 6, 4 \leq N \leq 1000)\):
Step 1. 计算 \(N \times N\) 的矩阵 \(C = AB\).
Step 2. 计算 \(M = C^{N \times N} \bmod 6\).
Step 3. 输出 \(M\) 中所有元素的和.

题解

根据矩阵乘法的结合律,

\[M = (AB)^{N \times N} = A (BA)^{N \times N – 1} B\]

时间复杂度由裸的快速幂\(O(N^3 \log N)\)变成了\(O(N^2k + Nk^2 + k^3 \log N)\)。


Matrix Multiplication

2014 Multi-University Training Contest 5, available at HDUOJ 4920

给定\(n \times n\)的矩阵\(A, B\),求模3意义下的乘积。 \((1 \leq n \leq 800)\)

题解

令 \(M_x\): \(M_x(i, j) = [M(i, j) = x]\).

\[ AB = (A_1 – A_2)(B_1 – B_2) = A_1B_1 – A_1B_2 – A_2B_1 + A_2B_2 \]

其中 \(A_1, A_2, B_1, B_2\) 是 01矩阵. 使用 std::bitset,可以做做到\(O(4n^3/wordsize)\)。

(当然,使用其他一些优秀的卡常方法也能过这题)


Matrix Revolution

HDUOJ 1759

给定一个 \(n \times n\) 的矩阵 \(A\),计算下述矩阵的非零元素的数量

\[ A + A^2 + \cdots + A^k\]

\(A\) 是稀疏矩阵,用 \(m\) 个三元组 \((a, b, c)\) 表示,每个三元组代表 \(A_{ab} = c\)。此外,保证主对角线上的元素全为 1。其他元素均为0。

\((0 \leq n \leq 1000, 0 \leq m \leq 10n, n < k < 10^{100}, 0 < c < 10^9)\)

题解

利用图的邻接矩阵表示,可知答案就是距离不超过 \(k\) 的点对数。


Matrix Multiplication

POJ 3318

给定 \(n \times n\) 的矩阵 \(A, B, C\)。判断\(A B = C\)是否成立。\((n \leq 500, |A_{ij}|, |B_{ij}| \leq 100, |C_{ij}| \leq 10,000,000)\)

提示:有多组测试数据,\(O(n^3)\) 的算法会 TLE。

题解

随机生成 \(n\) 维列向量 \(\vec{x}\),判断 \(AB\vec{x} = C\vec{x}\) 是否成立。可在 \(O(n^2)\) 时间内完成。

算法出错当且仅当 \((AB-C)\vec{x} = 0\) 而 \(AB – C \neq 0\)。 假设我们在 \(\mathbb{F}_p\) 里做。令 \(M = AB – C\), 则 \(\vec{x}\) 在 \(M\) 的零空间中。仔细想一想 \(M\) 的零空间最多 \(n-1\) 维,因此出错的概率不会超过\(1/p\)。

(当然,使用三次方的算法加一些优秀的卡常,也许也能过)

Remark: 

这一算法被称为Freivalds算法,是一种单边错Monte Carlo随机算法。对于整数矩阵,存在确定性的\(O(n^2)\)的算法判断乘法是否正确。不过由于Freivalds算法非常易于实现,目前实际上仍然大量使用该算法 。


Square

UVA 11542, also available at vjudge

给定 \(n\) 个整数的集合你可以生成 \(2^n-1\) 个非空子集。求有多少个子集,其中的数之积为完全平方数。

\(1 \leq T \leq 30\), \(1 \leq n \leq 100\), 所有整数在 1 和 \(10^{15}\) 之间。整数的素因子最多为 500。

题解

每个整数的唯一分解序列都看成 \(\mathbb{F}_2\) 上的一个的向量,就是要求这些向量构成矩阵的零空间的大小。

如何编写缓存友好的代码——从一道卡常数的题目说起

WC2017的第二题是一道卡常三合一题。其中,第一个子任务要求对最多2e8个32位无符号整数排序,时限3s,内存2G。输入采用伪随机数,输出采用hash,从而避免输入输出的影响。

看到题目第一反应是暴力。采用标准库中的sort函数,可以很容易地写出下面的代码

#include <cstdio>
#include <algorithm>
using namespace std;

typedef unsigned uint32_t;

inline uint32_t next_integer(uint32_t &x) {
  x = x ^ (x << 13);
  x = x ^ (x >> 17);
  x = x ^ (x << 5);
  return x;
}

const int n = 200000000;
uint32_t arr[n];

void genarray() {
  uint32_t x = 0x98765432;
  for (int i=0; i<n; i++) arr[i] = next_integer(x);
}

uint32_t output(uint32_t* ptr) {
  uint32_t ret = n * 4;
  uint32_t x = 23333333;
  for (int i=0; i<n; i++) {
    ret ^= ptr[i] + x;
    next_integer(x);
  }
  return ret;
}

int main() {
  genarray();
  sort(arr, arr+n);
  printf("%x\n", output(arr));
  return 0;
}

显然,这段代码是会TLE的。在我的电脑配置下(Intel Core i5-7200U, 8G, Ubuntu 16.04 LTS, 编译参数 g++ -O2 -m32 <source>,下同)一共跑了22.4s。

众所周知,标准库中的sort的复杂度是\(O(n \log n)\)。看来,必须使用线性时间排序才行。计数排序在这边用不了,因为32位无符号整数的取值范围太大了。那么,使用基数排序如何?将每个整数拆分为前16位和后16位,然后先以后16位为key使用计数排序,再以前16位为key使用计数排序,就可以在2倍线性时间内完成整个排序。

// ...
typedef unsigned uint32_t;
typedef unsigned short uint16_t;
typedef unsigned char uint8_t;

union {
  uint32_t u32;
  uint16_t word[2];
} arr[n], arr2[n];

uint32_t cnt[2][65536];

void counting() {
  for (int i=0; i<n; i++) {
    cnt[0][arr[i].word[0]]++;
    cnt[1][arr[i].word[1]]++;
  }
  for (int j=0; j<2; j++) partial_sum(cnt[j], cnt[j]+65536, cnt[j]);
}

void sort0() {
  for (int i=n-1; i>=0; i--)
    arr2[--cnt[0][arr[i].word[0]]].u32 = arr[i].u32;
}

void sort1() {
  for (int i=n-1; i>=0; i--)
    arr[--cnt[1][arr2[i].word[1]]].u32 = arr2[i].u32;
}

int main() {
  genarray();
  counting();
  sort0(); sort1();
  printf("%x\n", output((uint32_t*)arr));
  return 0;
}

这回运行时间为5.2s,比之前快了许多,但仍然不能在时间限制内通过。使用perf工具监测cpu用时,发现sort1和sort0函数是瓶颈,时间占比分别达到了34.37%和33.91%。再监测cache miss数量,发现程序总共发生约6.38亿次cache miss,而sort1和sort0中的占比分别为44.88%和43.74%。考虑基数排序的原理:每一趟基数排序是一次计数排序,对于相同key的元素,它们是在内存中连续放置的,而key不同的元素,则需要随机访问。通常,一趟基数排序的cache hit率可以如下估算:

\[ hit\% = \min\{1,  \frac{\#\text{cache line}}{k}\} * \left( 1 – \frac{\text{size of element}}{\text{size of cache line}}\right) \]

其中,k表示基数排序中基的大小,这里为65536。Intel i5-7200U的cache信息如下:

L1 data cache: 2 x 32 KBytes, 8-way set associative, 64-byte line size
L2 cache: 2 x 256 KBytes, 4-way set associative, 64-byte line size
L3 cache: 3 MBytes, 12-way set associative, 64-byte line size

经过简单的计算发现L3 cache一共只有约5万行,也就是说,这样写连3级缓存都不能充分利用,CPU必须经常花费数倍的时间访问主存,因此效率十分低下。

L1 cache一共只有512行,那么我们将32位整数分为4段,每段8位,然后进行基数排序行不行呢?让我们重新改写一下程序。

// ...

uint32_t cnt[4][256];

void counting() {
  for (int i=0; i<n; i++) {
    cnt[0][arr[i].byte[0]]++;
    cnt[1][arr[i].byte[1]]++;
    cnt[2][arr[i].byte[2]]++;
    cnt[3][arr[i].byte[3]]++;
  }
  for (int j=0; j<4; j++)
    partial_sum(cnt[j], cnt[j]+256, cnt[j]);
}

void sort0() {
  for (int i=n-1; i>=0; i--)
    arr2[--cnt[0][arr[i].byte[0]]].u32 = arr[i].u32;
}

void sort1() {
  for (int i=n-1; i>=0; i--)
    arr[--cnt[1][arr2[i].byte[1]]].u32 = arr2[i].u32;
}

void sort2() {
  for (int i=n-1; i>=0; i--)
    arr2[--cnt[2][arr[i].byte[2]]].u32 = arr[i].u32;
}

void sort3() {
  for (int i=n-1; i>=0; i--)
    arr[--cnt[3][arr2[i].byte[3]]].u32 = arr2[i].u32;
}

int main() {
  genarray();
  counting();
  sort0(); sort1(); sort2(); sort3();
  printf("%x\n", output((uint32_t*)arr));
  return 0;
}

跑下来运行时间为2.7s,能够在时限内通过。使用性能分析工具perf,发现总共只发生了2.5亿次cache miss。虽然经过修改后,理论运行时间变为2倍,但实际运行时间反而减少了将近一半。

CPU在访问主存时,往往需要等待数百个时钟周期。由于burst机制的存在,CPU可以一次性与主存交换大量连续数据,使得cache能够被十分有效的实现。要增加cache的利用率,主要是要提高程序访问数据的空间局部性和时间局部性。将 \(O(n \log n)\) 的排序改为线性时间排序,是一种典型的space-time tradeoff, 如果排序的基数太大,那么count数组会变得很大,排序后的数据在内存中的分布会很散,降低了缓存利用率,这时就可以通过适当降低基数的方式,提高cache利用率,从而达到理论运行时间和cache利用率之间的平衡。

k次幂前缀和的线性做法

给定正整数 \(x, k\) 和质数 \(p\),求
\[ \sum_{i=1}^{x} i^k \bmod p \]

这是一道很经典的题目,在杜瑜皓的《多项式及求和》中有详细讨论。比较暴力的做法是对于每个 \(i^k\) 都用快速幂算出其值,然后相加,时间复杂度为 \(O(x \log k) \),当 \(x\) 很大时,就超时了。然而利用拉格朗日插值法,这题最快可以做到 \(O(k)\) 的复杂度,其做法主要由以下几部分构成。

线性预处理逆元

求 1..n 之间每个数的逆元,如果都用费马小定理或者扩展欧几里得算,那么复杂度将会达到 \(O(n \log p)\) 或 \(O(n \log n)\)。利用一些递推式,可以线性地求出 1..n 中每个数的逆元,从而复杂度可以减少一个 log。常用的一个递推公式是

\[ i^{-1} = -\left\lfloor \frac{p}{i} \right\rfloor \cdot (p \bmod i)^{-1} \bmod p \]

这一公式的正确性也很好证明,等式右端可以写成

\[ -\frac{\left\lfloor \frac{p}{i} \right\rfloor}{p-\left\lfloor \frac{p}{i} \right\rfloor i} = \frac{1}{i} \]

利用此公式,可以在 \(O(n)\) 时间内求出 1..n 中每一个数的逆元。

线性筛求积性函数的值

注意到 \(i^k\) 是一个积性函数,这样,利用线性筛,我们只需要用快速幂计算素数 \(i\) 的值,其他点的函数值就可以直接得到。根据素数定理,前 \(n\) 个数中素数的个数约为 \(O(\frac{n}{\log n})\),从而利用快速幂计算出素数处的值,可以做到 \(O(\frac{n}{\log n}\log{k}\) ,再加上筛法本身 \(O(n)\) 的复杂度,我们可以在 \(O(n + \frac{n}{\log n}\log{k})\) 的时间复杂度内计算出 \(i^k\) 的前 \(n\) 项值。

拉格朗日插值法

注意到答案一定是关于 \(x\) 的 \(k+1\) 次多项式。而根据线性代数的知识,我们只要知道这个多项式的任意 \(k+2\) 个点处的值,我们就能确定这个多项式。拉格朗日插值公式就给出了这个多项式

\[ f(x) = \sum_{i=0}^{k+1} y_i l_i(x) \]

其中

\[ l_i(x) = \prod_{0 \leq m \leq k+1, m \neq i} \frac{x – x_m}{x_j – x_m} \]

具体到这道题上,插值公式就是

\[ f(n) = \sum_{i = 0}^{k+1} (-1)^{k+i+1} f(i) \frac{n(n-1) \cdots (n-i+1)(n-i-1) \cdots (n-d)}{(d-i)!i!} \]

其中 \(f(i)\) 就是 \(i^k\) 的前缀和,上式右端求和的每一项,分子都是 n 到 n-d 的前缀积乘后缀积,分母可以用递推线性求出,这样就可以在 \(O(k)\) 时间内解决整个问题。

#include <bits/stdc++.h>
using namespace std;

#define rep(i, n) for (int i=0; i<(n); i++)
#define Rep(i, n) for (int i=1; i<=(n); i++)
#define range(x) (x).begin(), (x).end()
typedef long long LL;
typedef unsigned long long ULL;

#define pow owahgrauhgso

ULL n, k, m;

ULL powmod(ULL b, ULL e) {
  ULL r = 1;
  while (e) {
    if (e & 1) r = r * b % m;
    b = b * b % m;
    e >>= 1;
  }
  return r;
}

const int MAXN = 1000006;
ULL inv[MAXN];

void init_inv() {
  inv[1] = 1;
  for (int i = 2; i <= k+1; i++) { 
    inv[i] = (m - m / i) * inv[m % i] % m; 
    assert(inv[i] * i % m == 1);
  }
}

ULL pow[MAXN];
ULL prime[MAXN], cnt;

void sieve() {
  pow[1] = 1;
  for (int i = 2; i <= k+1; i++) {
    if (!pow[i]) {
      pow[i] = powmod(i, k);
      prime[cnt++] = i;
      for (int j = 0; j < cnt && i*prime[j] <= k+1; j++) {
        pow[i * prime[j]] = pow[i] * pow[prime[j]] % m;
        if (i % prime[j] == 0) break;
      }
    }
  }
}

ULL sum[MAXN];
ULL ans[MAXN];

auto addmod = [](ULL a, ULL b) -> ULL {return (a+b)%m;};

ULL lagrange() {
  ULL p;
  p = 1;
  for (int i=0; i<=k+1; i++) {
    if (i) p = p * inv[i] % m;
    ans[i] = (k+1-i)&1 ? m-sum[i] : sum[i];
    ans[i] = ans[i] * p % m;
    p = p * (m + n - i) % m;
  }
  p = 1;
  for (int i=k+1; i>=0; i--) {
    if (k+1-i) p = p * inv[k+1-i] % m;
    ans[i] = ans[i] * p % m;
    p = p * (m + n - i) % m;
  }
  return accumulate(ans, ans+k+2, 0, addmod); 
}

 

int main() { cin >> n >> k >> m; init_inv(); sieve(); partial_sum(pow, pow+k+2, sum, addmod); cout << lagrange() << endl; return 0; }