如何高效查找N×N方格中对应非最简分数的格点？

问题背景 #

给定正整数N，N×N方格中的格点可表示为list(itertools.product(range(1, N + 1), repeat=2))，总计N²个。我们需要找出所有满足x/y为非最简分数的元组(x,y)。

现有实现方案 #

1. 暴力法find_complex_points #

该方法正确性有保证，但效率极低。遍历所有格点，通过计算最大公约数判断是否为非最简分数：

import math
from itertools import product

def find_complex_points(lim: int) -> list[tuple[int, int]]:
    return [
        (x, y)
        for x, y in product(range(1, lim + 1), repeat=2)
        if math.gcd(x, y) > 1
    ]

测试数据：当lim=1024时，耗时约233ms。

2. 优化版find_complex_points_1 #

通过枚举倍数生成非最简分数对应格点，但会产生重复项，仅实现小幅提速：

def find_complex_points_1(lim: int) -> set[tuple[int, int]]:
    lim += 1
    return {
        (x, y)
        for mult in range(2, lim)
        for x, y in product(range(mult, lim, mult), repeat=2)
    }

测试数据：当lim=1024时，耗时约194ms。

进阶优化尝试 #

自行设计的find_complex_points_2实现了大幅提速（lim=1024时约58.9ms），但仍存在重复计算问题（例如元组(10,10)会被多次生成），需要进一步消除冗余。

最终目标实现：生成对应图像 #

将坐标替换为自然数求和（即N替换为(N²+N)/2），生成指定尺寸的uint8型NumPy二维数组：初始填充255，若(x+1)/(y+1)为非最简分数，则将对应区域像素设为0。以下是完整实现代码：

图像缩放函数 #

使用Numba加速的图像缩放方法，将小尺寸图像按指定比例放大：

import numpy as np
import numba as nb

@nb.njit(cache=True)
def resize_img(img: np.ndarray, h_scale: int, w_scale: int) -> np.ndarray:
    height, width = img.shape
    result = np.empty((height, h_scale, width, w_scale), np.uint8)
    result[...] = img[:, None, :, None]
    return result.reshape((height * h_scale, width * w_scale))

无冗余的非最简分数格点查找 #

通过栈结构避免重复计算，高效生成所有符合条件的格点：

def find_composite_points(lim: int) -> set[tuple[int, int]]:
    stack = dict.fromkeys(range(lim, 1, -1))
    lim += 1
    points = set()
    while stack:
        x, _ = stack.popitem()
        points.add((x, x))
        mults = []
        for y in range(x * 2, lim, x):
            stack.pop(y, None)
            mults.append(y)
            points.update([(x, y), (y, x)])

        for i, x_val in enumerate(mults):
            points.add((x_val, x_val))
            for y_val in mults[i + 1 :]:
                points.update([(x_val, y_val), (y_val, x_val)])

    return points

自然数求和工具函数 #

计算1到n的自然数和，用于坐标转换：

def natural_sum(n: int) -> int:
    return (n + 1) * n // 2

图像生成主函数 #

生成初始图像，根据格点信息修改对应区域像素，最后返回缩放后的图像：

def composite_image(lim: int, scale: int) -> np.ndarray:
    length = natural_sum(lim)
    img = np.full((length, length), 255, dtype=np.uint8)
    for x, y in find_composite_points(lim):
        x1, y1 = natural_sum(x - 1), natural_sum(y - 1)
        img[x1 : x1 + x, y1 : y1 + y] = 0

    return resize_img(img, scale, scale)

内容的提问来源于stack exchange，提问作者Ξένη Γήινος