\\ Search range for a
aFirst = 0;
aLast  = 57;

\\ Polynomial variable for F_a(t,1)
t = 't;

for(a = aFirst, aLast,

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

  \\ F_a(t,1)
  pol = t^4 - a*t^3 - 6*t^2 + a*t + 1;

  \\ Thue initialization.
  \\ The second argument 1 means unconditional computation.
  th = thueinit(pol, 1);

  rawSolutions = List();

  \\ 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 = -bound, bound,
    if(rhs != 0,

      solutions = thue(th, rhs);

      for(i = 1, #solutions,
        x = solutions[i][1];
        y = solutions[i][2];

        if(gcd(x, y) == 1,
          val = x^4 - a*x^3*y - 6*x^2*y^2 + a*x*y^3 + y^4;
          listput(rawSolutions, [x, y, val])
        )
      )
    )
  );

  \\ Remove duplicates and sort.
  \\ Each element is [x, y, F_a(x,y)].
  allSolutions = Set(Vec(rawSolutions));

  \\ Verification for all primitive solutions.
  for(i = 1, #allSolutions,
    if(abs(allSolutions[i][3]) > bound,
      error("bad bound at a = ", a)
    );

    if(gcd(allSolutions[i][1], allSolutions[i][2]) != 1,
      error("not primitive at a = ", a)
    )
  );

  \\ Extract representatives with x >= 0 and y > 0.
  repList = List();

  for(i = 1, #allSolutions,
    x = allSolutions[i][1];
    y = allSolutions[i][2];

    if(x >= 0 && y > 0,
      listput(repList, allSolutions[i])
    )
  );

  representativeSolutions = Vec(repList);

  \\ Verification for representatives.
  for(i = 1, #representativeSolutions,
    if(representativeSolutions[i][1] < 0,
      error("bad representative x at a = ", a)
    );

    if(representativeSolutions[i][2] <= 0,
      error("bad representative y at a = ", a)
    )
  );

  \\ Compact output for this a.
  print("a = ", a);
  print("bound = ", bound);
  print("[x,y,F_a(x,y)] with x >= 0 and y > 0:");
  print(representativeSolutions);
  print("");
);