乘积为完美立方体的子阵列计数

本文介绍一种针对二维矩阵的算法，用于计算其所有乘积为完美立方体（即能表示为 $a^3$ 的正整数 $a$）的子矩阵个数。

算法概述

该算法基于对完全立方体的因数分解。设矩阵 $M$ 的大小为 $n \times m$，则对于 $M$ 中的每个位置 $(i,j)$，我们可以把其所在行（或列）的元素看做一个完全立方体的分解因子。

具体来说，对于矩阵 $M$ 中的每个位置 $(i,j)$，设其所在行的元素为 $r_1, r_2, \cdots, r_m$，则我们可以将这 $m$ 个数字分解为质因数形式，即 $r_k = p_1^{a_{1,k}} \cdot p_2^{a_{2,k}} \cdots p_s^{a_{s,k}}$，其中 $p_1,p_2,\cdots,p_s$ 是 $r_k$ 的所有质因数，$a_{1,k},a_{2,k},\cdots,a_{s,k}$ 是其对应的质因数指数。

由于我们要求的是所有行乘积为完美立方体的子矩阵个数，不妨先将每行元素的质因数指数按照相同的质因数分组，并将分组后的指数求和，即可得到每一行的 “分组后的” 矩阵。

接下来，我们只要枚举所求子矩阵的左右两列（或上下两行），并逐列（或逐行）求出所求子矩阵所有行的 “分组后的” 矩阵的积，判断该积是否为完美立方体即可。

假设一个 “分组后的” 矩阵的大小为 $1 \times k$，其分组情况为 $(a_1,a_2,\cdots,a_s)$，则它对应的完全立方体为 $p_1^{3a_1} \cdot p_2^{3a_2} \cdots p_s^{3a_s}$。因此，对于每一列的所有行的 “分组后的” 矩阵，我们只需要算它们的完全立方体的乘积，然后判断是否为完美立方体即可。为了优化程序，我们可以预处理所有可能的完美立方体，利用哈希表 $h$ 存储每个完美立方体及其个数，然后统计所有乘积在哈希表中的完美立方体的个数之和。

具体过程如下：

1. 将每行元素的质因数指数按照相同的质因数进行分组，并将分组后的指数求和，得到每一行的 “分组后的” 矩阵
2. 枚举所求子矩阵的左右两列（或上下两行）
3. 逐列（或逐行）处理所有行的 “分组后的” 矩阵，得到该列上（或行上）所有 “分组后的” 矩阵的完全立方体的乘积 $x$
4. 判断 $x$ 是否为哈希表 $h$ 中的完美立方体，如果是，则累加该立方体的个数
5. 统计所有乘积在哈希表中的完美立方体的个数之和即为所求子矩阵的个数

具体实现时，我们需要对于所有可能的质因数指数组合都预处理出其乘积是否为完美立方体及其个数，并存储到哈希表 $h$ 中，以加速判断过程。

代码实现

以下是基于 Python 的实现代码：

import math
from collections import defaultdict

# 预处理所有可能的完美立方体
perfect_cubes = {x ** 3 for x in range(1, int(math.pow(2, 21)))}

def perfect_cube_submatrices(matrix):
    n, m = len(matrix), len(matrix[0])
    cnt = 0 # 计数器
    for i in range(m): # 枚举所有子矩阵的左右两列
        # 对所有行，计算其 “分组后的” 矩阵
        groups = [[0] * len(primes) for _ in range(n)]
        for j in range(n):
            num = matrix[j][i]
            for k, p in enumerate(primes):
                while num % p == 0:
                    num //= p
                    groups[j][k] += 1
        # 计算每列上所有 “分组后的” 矩阵的完全立方体的乘积
        prod = [1] * n
        for j in range(i+1, m):
            for k in range(n):
                prod[k] *= groups[k][prime_idx[k][j]]
            # 判断该乘积是否为完美立方体
            if prod.count(0) > 0: # 直接忽略 “分组后的” 矩阵全为 0 的情况
                continue
            cube_product = 1
            for c in prod:
                cube_product *= pow(primes[c > 0], 3 * c)
            cnt += perfect_cubes.get(cube_product, 0)
    return cnt

# 基于欧拉筛法预处理出所有小于等于 10^7 的质数，并调用 `perfect_cube_submatrices()`
def main():
    sieve = [True] * (int(1e7) + 1)
    primes = []
    for i in range(2, len(sieve)):
        if sieve[i]:
            primes.append(i)
        j = 0
        while j < len(primes) and primes[j] * i <= len(sieve):
            sieve[primes[j] * i] = False
            if i % primes[j] == 0:
                break
            j += 1
    global primes
    primes = [p for p in primes if p > 2] # 只保留大于 2 的质数
    global prime_idx
    prime_idx = [[bisect.bisect(primes, d) - 1 for d in r] for r in matrix] # 预处理每个数的质数指数
    return perfect_cube_submatrices(matrix)

性能分析

假设矩阵大小为 $n \times m$，预处理所有可能的完美立方体需要的时间复杂度为 $O(P)$，其中 $P$ 为小于等于 $2^{21}$ 的完全立方体数，约为 $3.6 \times 10^5$。预处理出每个数的质数指数需要的时间复杂度为 $O(nm \log \log \max a_{i,j})$，其中 $\max a_{i,j}$ 是 $M$ 中的最大值。因此，总时间复杂度为 $O(P + nm \log \log \max a_{i,j} + S)$，其中 $S$ 是所有的乘积为完美立方体的子矩阵的个数。

该算法主要的时间瓶颈在于哈希表查询，可以采用更高效的哈希表实现以提高性能。此外，对于较小的矩阵可以直接暴力穷举子矩阵，而不需要使用本文介绍的算法。