limit = radius ** 2 + epsilon
for x, y, z in it.product(list(range(-radius, radius + 1)), repeat=3):
    if x ** 2 + y ** 2 + z ** 2 <= limit:
        yield x + x0, y + y0, z + z0
понимание python с несколькими предложениями «для» и одним «если»