[Python] How to find all grid points that correspond to non-reduced fractions in a square?

Given a positive integer N, we can label all grid points in the square N x N, starting at 1, the total number of grid points is N x N, and the grid points are list(itertools.product(range(1, N + 1), repeat=2)).

Now, I want to find all tuples (x, y) that satisfies the condition x/y is a non-reduced fraction, the following is a bruteforce implementation that is guaranteed to be correct, but it is very inefficient:

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
]


Now the next function is slightly smarter, but it generates duplicates and as a result is only noticeably faster but not by much:

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)
}


In [255]: %timeit find_complex_points(1024)
233 ms ± 4.44 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [256]: %timeit find_complex_points_1(1024)
194 ms ± 1.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)


Is there a better way to accomplish this?

(My goal is simple, I want to create a NumPy 2D array of uint8 type with shape (N, N), fill it with 255, and make all pixels (x, y) 0 if (x+1)/(y+1) is a non-reduced fraction)

Continue reading...