Question: Given a NxN integer lattice, I want to find the clipped circle which maximizes the sum of its interior lattice point values. Each lattice point (i,j) has a value V(i,j) and are stored in the following matrix 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]]  The goal is to maximize the sum of values V(i,j) of the lattice points lying on the boundary and within interior of a (clipped) circle with radius R, with the assumptions and conditions:  the circle has center at (0,0) the circle can have any positive radius (not necessarily an integer radius, i.e., rational). the circle may be clipped at two lattice points, resulting in a diagonal line as shown in the picture. This diagonal line has a slope of -45 degrees.  Some additional details: The score for a clipped circle is the sum of all the integers that are both within the circle (or on the border) and on the side of the diagonal line including (0,0). The values on (or near) the border are -3, 1, 3, -1, -3, 3, -1, 2, 0, 3. Even though the circle can have any radius, we need only consider circles that intersect a grid point precisely so there are n^2 different relevant radiuses. Further, we need only record one position where the circle intersects with the diagonal line to fully specify the clipped circle. Note that this intersection with the diagonal does not need to be at an integer coordinate. If the optimal solution doesn't have the diagonal clipping the circle at all then we need only return the radius of the circle. What I have found so far: If we only wanted to find the optimal circle we could do that quickly in time proportional to the input size with: 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) # returns optimal radius with its score return sqrt(np.argmax(cusum)), np.max(cusum) A = np.random.randint(-3, 4, (10, 10)) print(find_max(A))  How quickly can the optimal clipped circle be found?

Answer: Start by creating a cumulative frequency table, or a fenwick tree. You'll have a record for each radius of circle, with value corresponding to explored weights at that distance from the origin. Then, begin a BFS from the origin. For each diagonal "frontier", you'll need to update your table/tree with the radius:weight key-value pair (add weight to existing value). You'll also need to then query the table/tree for the current cumulative sum at each radius just added, noting the maximum and updating a global running maximum accordingly. Once your search terminates, you'll have the maximum sum for your clipped-circle. If you want to reconstruct the circle, just store the max radius and BFS depth along with the global max sum itself. This will give you your solution in O(N^2 log N) time, as there will be N^2 updates and queries, which are O(log N) each. The intuition behind this solution is that by exploring along this diagonal "frontier" outward, you implicitly clip all your circles you query since the weights above/right of it haven't been added yet. By calculating the max (at each search depth) for just the radii that were just updated, you also enforce the constraint that the circles intersect the clipping line at an integer coordinate. Update Here is python code showing this in action. It needs cleaned up, but at least it shows the process. I opted to use cumulative frequency / max arrays, instead of trees, since that'll probably lend itself to vectorization with numpy for OP. def solve(matrix): n = len(matrix) max_radius_sqr = 2 * (n - 1) ** 2 num_bins = max_radius_sqr.bit_length() + 1 frontier = [(0, 0)] csum_arr = [[0] * 2 ** i for i in range(num_bins)[::-1]] cmax_arr = [[0] * 2 ** i for i in range(num_bins)[::-1]] max_csum = -float("inf") max_csum_depth = None max_csum_radius_sqr = None depth = 0 while frontier: next_frontier = [] if depth + 1 < n: # BFS up next_frontier.append((0, depth + 1)) # explore frontier, updating csums and maximums per each for x, y in frontier: if x + 1 < n: # BFS right next_frontier.append((x + 1, y)) index = x ** 2 + y ** 2 # index is initially the radius squared for i in range(num_bins): csum_arr[i][index] += matrix[y][x] # update csums if i != 0: # skip first, since no children to take max of sum_left = csum_arr[i-1][index << 1] # left/right is tree notation of the array max_left = cmax_arr[i-1][index << 1] max_right = cmax_arr[i-1][index << 1 | 1] cmax_arr[i][index] = max(max_left, sum_left + max_right) # update csum maximums index >>= 1 # shift off last bit, update sums/maxs again, log2 times # after entire frontier is explored, query for overall max csum over all radii # update running global max and associated values if cmax_arr[-1][0] > max_csum: max_csum = cmax_arr[-1][0] max_csum_depth = depth index = 0 for i in range(num_bins-1)[::-1]: # reconstruct max radius (this could just as well be stored) sum_left = csum_arr[i][index << 1] max_left = cmax_arr[i][index << 1] max_right = cmax_arr[i][index << 1 | 1] index <<= 1 if sum_left + max_right > max_left: index |= 1 max_csum_radius_sqr = index depth += 1 frontier = next_frontier # total max sum, dx + dy of diagonal cut, radius ** 2 return max_csum, max_csum_depth, max_csum_radius_sqr  Calling this with the given test case produces the expected output:  matrix = [ [-1, 3, -3, -2, 0, -1, -2, -1, -1, 2], [ 0, -3, 0, 3, 2, -2, 3, -2, 3, 3], [-3, 1, 3, 3, 0, -3, -3, 2, -2, 1], [ 3, 2, 2, -1, 0, -3, 1, 1, -2, 2], [-3, 3, 2, 0, -3, -2, -1, -3, 0, -3], [-1, -2, -1, 2, 3, 3, -3, -3, 2, 0], [-2, 0, -3, 3, 0, 2, -1, 1, 3, 3], [ 2, 2, -3, 2, -2, -1, 2, 2, -2, 0], [-2, -1, 0, 1, 0, -2, 0, 0, 1, -3], [ 1, 1, -3, 0, 0, 3, -1, 3, -3, 2], ][::-1] print(solve(matrix)) # output: 13 9 54  In other words, it says the maximum total sum is 13, with a diagonal cut stagger (dx + dy) of 9, and radius squared of 54. If I have some time tonight or this weekend, I'll clean up the code a bit.