Issues generating all primitive integer triangles with 60 degree angle (Eisenstein Triples)

Question

I am trying to write a program that generates all primitive Eisenstein Triples, i.e. triangles with integer sides that have a 60 degree angle. These triplets satisfy the equation $a^2-ab+b^2=c^2$. Paraphrasing from this link:

Integer triangles with a 60° angle can also be generated by:
$a = m^2 - mn + n^2$
$b = 2mn - n^2$
$c = m^2 - n^2$
with coprime integers $m,n$ , with $0<n<m$. From here, all primitive solutions can be obtained by dividing a, b, and c by their greatest common divisor.

My problem is that I know that the numbers $3, 8, 7$ are an Eisenstein Triple, but it seems this generator algorithm skips it. Furthermore, setting the formula for $c$ above equal to either $3, 8,$ or $7$ leads to different triples, none of which are $3, 8, 7$:

If $c=3$, then $m=2$ and $n=1$, giving $a=3$ and $b=3$
If $c=8$, then $m=3$ and $n=1$, giving $a=7$ and $b=5$
If $c=7$, then $m=4$ and $n=3$, giving $a=13$ and $b=15$

Am I missing something? Thanks in advance for any insight.

Edit: Well, I answered my own question! I forgot to consider the case where $gcd(a, b, c) = 3$. Setting $c = 9$ requires $m=5, n=4$, which gives $a=21, b=24$, and dividing $a, b,$ and $c$ by their gcd does give $3, 8, 7$.

Edit 2: I realized the solution I found above would be skipped by my algorithm anyway, since I am optimizing by limiting $n \leq \frac{m}{2}$ (see the article above). As it turns out, $m=5$ and $n=1$ satisfies that condition and gives the same triple as well. My program was not properly reducing the triple by the gcd, which is why I did not notice this during debugging.

Answer 1

Below is a Python script using numpy library to generate Eisenstein triples. First thing, so the non-primitive versions are excluded, the list of primes "p" can be generated; if you want all the Eisenstein triples, feel free to skip the prime check. This is a brute force way to collect these. Parameterized solutions are dangerous because it's easy to skip solutions that the equation doesn't cover (e.g. Pythagorean quadruples). There is a way to produce these by matrices/linear transformations, but that's a subject for another time.

>import numpy as np
>level = 50
>sq = 0
># antipythagorean triples, minus & plus
>apall = np.array([0,0,0])
>apminus = np.array([0,0,0])
>applus = np.array([0,0,0])
>for x in range(1,level):
>    for y in range(1,level):
>        c = 0
>        for pr in list(p):  #prime check to guarantee primitive triples
>            #print(x,y,pr,end=" ")
>            if x%pr == 0 and y%pr ==0:
>                c = c +1
>                #print(x,y,"is not primitive")
>                #continue
>        #print(" ")
>        if c > 0:
>            continue
>        apm = x**2 + y**2 - x*y
>        app = x**2 + y**2 + x*y
>        apmsqr = apm**0.5
>        appsqr = app**0.5
>        if apmsqr == int(apmsqr) and x != y:
>            plac = np.array([x,y,int(apmsqr)])
>            print(plac,"-")
>            apminus = np.vstack((apminus,plac))
>            apall = np.vstack((apall,plac))
>        if appsqr == int(appsqr):
>            plac = np.array([x,y,int(appsqr)])
>            print(plac,"+")
>            applus = np.vstack((applus,plac))
>            apall = np.vstack((apall,plac))