Q := Rationals();
RR<x,y,z,t,u> := PolynomialRing(Q,5);

ListY := [ Parent([Integers()|1,2,3]) | ];
ListTM := [ Parent([RR|x^2]) | ];

TM := [RR|];

function Pass(TM) /* Checks that no 3 of TM are divisible by x,y,z,t,u. Proof:
If x divides 3 of the tangent monomials then Y must be singular on x = 0. */
 return &and
  [
  &and[#[ m : m in TM | IsDivisibleBy(m,RR.i)] le 2] : i in [1..5]
  ];
end function;

procedure CYloop(matrix, monrange, var_i, n, ~ListY, ~ListTM, ~TM)
 // var_i represents which variable we are dealing with, out of n

for mon in monrange do
 mi := mon[1]; ji := mon[2];
 TM[var_i] := RR.var_i^(mi-1)*RR.ji;
 adjoin := ZeroMatrix(Integers(),1,n+1);
 adjoin[1,var_i] := mi-1;
 adjoin[1,ji] := adjoin[1,ji]+1;
 adjoin[1,n+1] := -1;
 M := VerticalJoin(matrix, adjoin);
 N := EchelonForm(M);
 if N[var_i+1,n+1] ge 0 then
  max := -1; else
  max := Floor((&+[N[var_i+1,i] : i in [var_i+1..n]])/(-N[var_i+1,n+1])) +1;
 end if;
 range := [ [m,j] : j in [1..n], m in [2 .. max] ];
 
 //printf "Variable: %o, max = %o.\n", var_i, max;

if max eq -1 then
  if var_i eq n-2 then 
    //printf "Bust after variable %o: TMons = %o.\n", n-2, TM[1..n-2]; 
  end if;
else
 if var_i ne (n-1) then
  CYloop(N, range, var_i+1, n, ~ListY, ~ListTM, ~TM);
 else
// We are in the innermost loop (var_i = n-1)
    if IsDivisibleBy(N[n,n],N[n,n+1]) then // pure power works
      TM[n] := u^(N[n,n] div -N[n,n+1]);
      if Pass(TM) then
        Y := Eltseq(Basis(Kernel(Transpose(N)))[1]);
        // this Y is a solution; now test order
          if &and[Y[i] ge Y[i+1] : i in [1..n-1]] and Y[n] ge 0 then
           Append(~ListY,Y); Append(~ListTM,TM); // this is the output
       end if; // test of order
      end if; // refers to Pass(TM)
    else
      Y := Eltseq(Basis(Kernel(Transpose(N)))[1]);
      // debug: if Y[6] eq 1386 then N; Y; end if;
      // we're treating this as a potential solution, but
      // TM[n] not yet decided.
      // if Min(Y) gt 0 then
     if &and[Y[i] ge Y[i+1] : i in [1..n-1]] and Y[n] ge 0 then
      for j in [1..n-1] do //loop over j = jn
       if IsDivisibleBy(Y[n+1]-Y[j], Y[n]) then
         mn := ((Y[n+1]-Y[j]) div Y[n]) + 1;
         TM[n] := u^(mn-1)*RR.j;
         if Pass(TM) then
           Append(~ListY,Y); Append(~ListTM,TM); // this is the output
         end if;
       end if; // ends test of mn
      end for; // ends loop over j = jn
      end if; // ends check of smallest weight being greater than 0
     end if; // ends pure power case division
  end if; // ends test of var_i != n-1
end if; // if 'max', the upper bound, is -1, then we do nothing -- incompatible 

// does this free up memory? probably not, since we just reassign to them anyway
// delete N; delete M; delete adjoin; delete range;

end for;

end procedure;

/* n = 5 variables
n := 5;
range1 := [
  [5,1], // [5,2], [5,3], [5,4], [5,5]  how to know that these are all the same?
  [4,1],[4,2],[4,3],[4,4],[4,5],
  [3,1],[3,2],[3,3],[3,4],[3,5],
  [2,1]];

M := Matrix(1,6,[1,1,1,1,1,-1]);
*/

n := 4;
range1 := [
  [4,1],[4,2],[4,3],[4,4],
  [3,1],[3,2],[3,3],[3,4],
  [2,1]];

M := Matrix(1,n+1,[1,1,1,1,-1]);


CYloop(M, range1, 1, n, ~ListY, ~ListTM, ~TM);

print "Sample: ";
for i in [1000..1020] do
ListY[i];
end for;

#ListY; // 30528, 47265, 52844, 17206, 39337

LListY := Sort(SetToSequence(SequenceToSet([Sort(Y) : Y in ListY])));
#LListY; // 7810
