How to find multiple solutions for 3 variable, 2 degree Diophantine equation?

Question

I have the equation $x^2+xy+y^2=z^2$ to solve it in natural numbers.

The discriminant of it $D=4z^2-3y^2$ and must be perfect square.

I wrote Python program to get solutions for $1<x<100$ by enumeration.

def Solution():
    A=[]
    nMaximum=10**2
    for x in range(1,nMaximum):
        dTemp1a=3*x**2
        for z in range(x+1, nMaximum):
            dDiscriminant=4*z**2-dTemp1a
            dTemp5=int(dDiscriminant**0.5)
            if dTemp5**2!=dDiscriminant:
                continue
            dTemp6=(-1*x+dTemp5)/2
            y=int(dTemp6)
            if not CheckIfExists(A, z):
                A.append([x,y,z])
    return A
def CheckIfExists(arr, z):
    bResult=False
    for s in arr:
        if s[2]==z:
            bResult=True
            break
    return bResult
a = Solution()
print(len(a))
print(a)
[3, 5, 7], [5, 16, 19], [6, 10, 14], [7, 8, 13] ...

Three variable, second degree diophantine equation doesn't explain how to get other solutions when we know the first solution $(3,5,7)$

Could you give me a hint ?

UPDATE asked by Shean: I need to get all solutions based on $(3,5,7)$. See my question as the example of what I am looking for:

First 30 solutions of Pell's equation.

or see these formulas for all solutions of Pell's equation

http://mathworld.wolfram.com/PellEquation.html

Answer 1

Can't fit this into a comment so I'll make it an answer.

Firstly, your question isn't clear.

Three variable, second degree diophantine equation doesn't explain how to get other solutions when we know the first solution $(3,5,7)$

If you are expecting to generate ALL triples from $(3,5,7)$ then I don't think this is possible (actually it may be possible, see Will's comment). If you simply want more solutions then take the triple $(3k,5k,7k)$ for any positive integer $k$.

Now if you're wondering why the solutions to formulas such as \begin{equation} x=m^2-n^2\\ y=2mn+n^2\\\tag{1} z=m^2+mn+n^2 \end{equation} where $n<m$ and $\gcd(m,n)=1$ produce less results to: simply plugging in values for $x$ and $y$ and checking if $x^2+xy+y^2$ is a square, is because $(1)$ generates primitive solutions, that is $\gcd(x,y,z)=1$.

(Note: $(1)$ sometimes produces non-primitive solutions, to remedy this, simply divide $x,y,z$ by $\gcd(x,y,z)$)

Some example code using the system $(1)$ producing only primitive triples:

import math
limit = 10
triples = []
for n in range(1, limit):
    for m in range(n + 1, limit):
        if math.gcd(m, n) == 1:
            x = m2 - n2
            y = 2 * m * n + n2
            z = m2 + m * n + n**2
            if y<x:
                x_temp = x
                x=y
                y=x_temp
            d=math.gcd(x,math.gcd(y,z)) #removing non-primitives
            triples.append((int(x/d), int(y/d), int(z/d)))
#remove duplicates
triples= list(dict.fromkeys(triples))
Sort the triples based on x and then y
sorted_triples = sorted(triples, key=lambda triple: (triple[0], triple[1]))
for triple in sorted_triples:
    print("x =", triple[0])
    print("y =", triple[1])
    print("z =", triple[2])
    print()

Some example code plugging values for $x$ and $y$ and simply verifying the equation (This is what your code above was doing):

limit = 100
x=1
y=1
while x<limit:
    y=x
    while y<limit:
        z=(x**2+x*y+y**2)**0.5
        if z==int(z):
            print("x="+str(x)+"\ny="+str(y)+"\nz="+str(int(z))+"\n")
        y+=1
    x+=1

You can verify for yourself that any solution from the second code is a non-primitive. Or, you can run a variant of the first code which multiplies each primitive by $k=2,3,4,5$ and notice that the outputs are the same. Simply replace:

d=math.gcd(x,math.gcd(y,z)) #removing non-primitives
triples.append((int(x/d), int(y/d), int(z/d)))

with

for k in range(2,5):
    triples.append((x*k,y*k,z*k))