This problem pops up quite often in the search for primes. The number of integral solutions to an elliptical curve ax^2+by^2=n influences the chances that n is prime or not. This fact is exploited in the Sieve of Atkin – see my post about that. If both coefficients are positive, the curve is an ellipse. If one of the coefficients are negative, the curve is a hyperbola. Let’s look at both cases. We will consider only Quadrant I solutions.

## Ellipses

My first idea was to find the maximum x and y values, then generate the Cartesian product of coordinate pairs. In effect, this overlays a grid of points atop the ellipse. But the vast majority of points don’t lie anywhere near the curve! These points would be wasteful to check. Rather than drawing a box around an ellipse, why not walk down the curve like stairs?

We start at the positive (x=1)-intercept of the ellipse, which is (1,\sqrt{\frac{n-a}{b}}). We then iterate x-values up to the floor of the x-intercept, which is (\sqrt{\frac{n}{a}},0). y-values are generated using the x-values, then rounded to the nearest integer. This set of coordinate pairs “steps” down the ellipse.

We then just test these coordinates by plugging into ax^2+by^2 to see if it spits out n. The ones that do are the integral solutions. Here’s my implementation in Python:

def find_quadrant1_lattice_points_of_ellipse(n, a, b): """Steps down an ellipse, generating rounded-off coordinates, then checks them. Returns a list of lattice points if there are any.""" lattice_points = [] x_intercept = math.sqrt(n/a) for x in range(1, math.floor(x_intercept) + 1): y = round(math.sqrt((n - a*x**2)/b)) if y != 0 and n == a*x**2 + b*y**2: lattice_points.append((x, y)) return lattice_points

## Hyperbolas

Hyperbolas are a bit different. For starts, they look a bit like a butterfly with wings stretching out into infinity. Therefore, finding lattice points need some kind of stopping condition. For the Sieve of Atkin, only integral solutions where x>y are considered. So, let’s do that. We will start at the (y=1)-intercept and climb the steps up by iterating y.

def find_quadrant1_lattice_points_of_hyperbola(n, a, b): """Steps up a hyperbola, generating rounded-off coordinates, then checks them. Returns a list of lattice points satisfying x > y, if there are any""" lattice_points = [] stopping_y = math.sqrt(n/(a+b)) for y in range(1, math.floor(stopping_y) + 1): x = round(math.sqrt((n - b*y**2)/a)) if x > y and n == a*x**2 + b*y**2: lattice_points.append((x, y)) return lattice_points

If you want to use them, fetch them from my GitHub.