a := 1;

// Bound for |F_a(x,y)|
bound := 6*a + 7;

// Binary quartic form
FormValue := function(a, x, y)
    return x^4 - a*x^3*y - 6*x^2*y^2 + a*x*y^3 + y^4;
end function;

// F_a(t,1)
P<t> := PolynomialRing(Integers());
f := t^4 - a*t^3 - 6*t^2 + a*t + 1;

// Create Thue object.
T := Thue(f);

// Store all primitive solutions as pairs <x,y>
solutionPairs := { <0, 0> };
Exclude(~solutionPairs, <0, 0>);

// Solve F_a(x,y) = rhs for all nonzero rhs with |rhs| <= bound.
// rhs = 0 is skipped because F_a(x,y)=0 has no primitive integer solution.
for rhs in [-bound..bound] do
    if rhs ne 0 then
        for sol in Solutions(T, rhs) do
            x := Integers()!sol[1];
            y := Integers()!sol[2];

            if GCD(x, y) eq 1 then
                Include(~solutionPairs, <x, y>);
            end if;
        end for;
    end if;
end for;

// Sort all primitive solutions
solutionPairsSeq := SetToSequence(solutionPairs);
Sort(~solutionPairsSeq);

// Convert to triples <x,y,F_a(x,y)>
allSolutions := [
    <pair[1], pair[2], FormValue(a, pair[1], pair[2])>
    : pair in solutionPairsSeq
];

// Extract representatives with x >= 0 and y > 0
representativeSolutions := [
    sol : sol in allSolutions
    | sol[1] ge 0 and sol[2] gt 0
];

// Verification
for sol in allSolutions do
    assert Abs(sol[3]) le bound;
    assert GCD(sol[1], sol[2]) eq 1;
end for;

for sol in representativeSolutions do
    assert sol[1] ge 0;
    assert sol[2] gt 0;
end for;

// Output
printf "a = %o\n", a;
printf "bound = %o\n", bound;
printf "number of all primitive solutions = %o\n\n", #allSolutions;

printf "All primitive solutions (x,y,F_a(x,y)):\n";
for sol in allSolutions do
    printf "(%o, %o, %o)\n", sol[1], sol[2], sol[3];
end for;

printf "\nnumber of representative solutions = %o\n\n", #representativeSolutions;

printf "Representative solutions (x,y,F_a(x,y)) with x >= 0 and y > 0:\n";
for sol in representativeSolutions do
    printf "(%o, %o, %o)\n", sol[1], sol[2], sol[3];
end for;