寻找N×N整数格点中最优裁剪圆的高效求解方法

给定一个10×10的整数格点矩阵V，每个格点(i,j)对应值V(i,j)：

V = [
 [ 1,  1, -3,  0,  0,  3, -1,  3, -3,  2],
 [-2, -1,  0,  1,  0, -2,  0,  0,  1, -3],
 [ 2,  2, -3,  2, -2, -1,  2,  2, -2,  0],
 [-2,  0, -3,  3,  0,  2, -1,  1,  3,  3],
 [-1, -2, -1,  2,  3,  3, -3, -3,  2,  0],
 [-3,  3,  2,  0, -3, -2, -1, -3,  0, -3],
 [ 3,  2,  2, -1,  0, -3,  1,  1, -2,  2],
 [-3,  1,  3,  3,  0, -3, -3,  2, -2,  1],
 [ 0, -3,  0,  3,  2, -2,  3, -2,  3,  3],
 [-1,  3, -3, -2,  0, -1, -2, -1, -1,  2]
]

目标与约束 #

需找到最大化内部格点值之和的裁剪圆，约束条件：

• 圆心固定在(0,0)
• 半径R可为任意正实数（无需整数或有理数）
• 圆可被斜率为-45°的对角线裁剪，该对角线由两个格点确定

额外规则 #

• 裁剪圆的分数：所有同时满足「在圆内（含边界）」且「位于对角线包含(0,0)一侧」的格点值之和。
• 仅需考虑恰好经过格点的圆（共n²个相关半径）；裁剪圆可通过圆与对角线的一个交点确定（交点无需为整数坐标）。
• 若最优解无需裁剪，仅返回对应半径即可。

现有实现 #

已实现无裁剪的最优圆求解算法，时间复杂度与输入规模成正比：

import numpy as np
from math import sqrt
np.random.seed(40)

def find_max(A):
    n = A.shape[0]
    sum_dist = np.zeros(2 * n * n, dtype=np.int32)
    for i in range(n):
        for j in range(n):
            dist = i**2 + j**2
            sum_dist[dist] += A[i, j]
    cusum = np.cumsum(sum_dist)
    # 返回最优半径及其分数
    return sqrt(np.argmax(cusum)), np.max(cusum)

A = np.random.randint(-3, 4, (10, 10))
print(find_max(A))

咨询问题 #

如何高效求解最优裁剪圆？其求解效率如何？

核心求解思路 #

要高效找到最优裁剪圆，需同时结合半径筛选和对角线侧格点过滤两个维度，步骤如下：

1. 格点预处理：对每个格点(i,j)，计算其到原点的距离平方d=i²+j²，以及i+j的值（用于判断是否在对角线包含原点的一侧），将所有格点按d从小到大排序。
2. 对角线参数枚举：斜率为-45°的对角线可由x+y=k（k为整数）表示，k的取值范围为0到2n-2（n为矩阵边长）。对每个k，计算所有满足i+j≤k且在圆内的格点值之和。
3. 动态累加计算：遍历排序后的格点，维护一个按i+j分组的累加数组，实时更新不同k值下的格点和前缀和，记录每个半径对应的最大分数。
4. 全局最优筛选：对比所有k值对应的最优分数与无裁剪的最优分数，取最大值对应的半径（及k值，若需裁剪）作为最终结果。

优化技巧 #

• 预计算格点属性：提前存储每个格点的d和i+j值，避免重复计算。
• 排序复用：按d排序格点后，可动态累加不同k值下的前缀和，无需重复遍历所有格点。
• 内存复用：使用同一个累加数组处理不同k值，减少内存开销。

代码实现示例 #

import numpy as np
from math import sqrt

def find_max_clipped(A):
    n = A.shape[0]
    # 预处理所有格点：(距离平方d, s=i+j, val)
    points = []
    max_d = 0
    min_s, max_s = float('inf'), -float('inf')
    for i in range(n):
        for j in range(n):
            d = i**2 + j**2
            s = i + j
            val = A[i][j]
            points.append((d, s, val))
            max_d = max(max_d, d)
            min_s = min(min_s, s)
            max_s = max(max_s, s)
    
    # 按距离平方d从小到大排序格点
    points.sort(key=lambda x: x[0])
    
    # 初始化s相关的累加数组，映射s到0-based索引
    s_offset = -min_s
    s_count = max_s - min_s + 1
    s_sum = np.zeros(s_count, dtype=np.int32)
    d_to_k_sum = []
    current_total = np.zeros(s_count, dtype=np.int32)
    
    prev_d = -1
    for d, s, val in points:
        if d != prev_d:
            # 记录上一个d对应的所有k的累加和
            if prev_d != -1:
                d_to_k_sum.append((prev_d, current_total.copy()))
            prev_d = d
        # 更新当前s的累加值
        s_idx = s + s_offset
        s_sum[s_idx] += val
        # 计算前缀和，current_total[k_idx]对应s≤(min_s +k_idx)的总和
        current_total = np.cumsum(s_sum)
    
    # 添加最后一个d的记录
    d_to_k_sum.append((prev_d, current_total.copy()))
    
    # 初始化全局最优值
    max_score = -float('inf')
    best_r = 0
    best_k = None
    
    # 先计算无裁剪的最优解
    sum_dist = np.zeros(max_d + 1, dtype=np.int32)
    for d, s, val in points:
        sum_dist[d] += val
    cusum = np.cumsum(sum_dist)
    no_clip_r = sqrt(np.argmax(cusum))
    no_clip_score = np.max(cusum)
    if no_clip_score > max_score:
        max_score = no_clip_score
        best_r = no_clip_r
        best_k = None
    
    # 遍历所有可能的k值，计算带裁剪的最优解
    for k in range(min_s, max_s + 1):
        k_idx = k - min_s
        current_max = -float('inf')
        current_best_d = 0
        # 找到当前k下，所有d对应的最大累加和
        for d, k_sum in d_to_k_sum:
            if k_sum[k_idx] > current_max:
                current_max = k_sum[k_idx]
                current_best_d = d
        # 更新全局最优
        if current_max > max_score:
            max_score = current_max
            best_r = sqrt(current_best_d)
            best_k = k
    
    # 返回结果：无裁剪则返回半径和分数，否则返回半径、k值和分数
    if best_k is None:
        return (best_r, max_score)
    else:
        # k对应对角线方程x+y=k，包含原点的一侧为x+y ≤k
        return (best_r, best_k, max_score)

# 传入给定的矩阵V
V = np.array([
 [ 1,  1, -3,  0,  0,  3, -1,  3, -3,  2],
 [-2, -1,  0,  1,  0, -2,  0,  0,  1, -3],
 [ 2,  2, -3,  2, -2, -1,  2,  2, -2,  0],
 [-2,  0, -3,  3,  0,  2, -1,  1,  3,  3],
 [-1, -2, -1,  2,  3,  3, -3, -3,  2,  0],
 [-3,  3,  2,  0, -3, -2, -1, -3,  0, -3],
 [ 3,  2,  2, -1,  0, -3,  1,  1, -2,  2],
 [-3,  1,  3,  3,  0, -3, -3,  2, -2,  1],
 [ 0, -3,  0,  3,  2, -2,  3, -2,  3,  3],
 [-1,  3, -3, -2,  0, -1, -2, -1, -1,  2]
])

print(find_max_clipped(V))

效率分析 #

• 时间复杂度：O(n² log n)。主要开销来自对n²个格点按距离平方排序的操作（O(n² log n)），其余预处理、动态累加和全局筛选步骤均为线性或线性对数级，整体效率与输入规模的平方对数成正比。
• 空间复杂度：O(n²)。用于存储格点信息和中间累加数组，对于n=10的场景内存占用极小；即使n=100，仅需存储10000个格点及对应数组，完全可控。
• 适用性：该算法可高效处理任意规模的n×n格点矩阵，既能准确找到带裁剪的最优圆，也能兼容无裁剪的最优解场景。

内容的提问来源于stack exchange，提问作者Simd