This file contains the commented procedures to calculate the list of
CY 3-folds on ~miles/CYlist_output.txt

(The list was calculated with repetitions and in no particular order, and then
sorted and repeated lines removed using a text editor.)

First round: the simple sum of reciprocals

Quasismooth hypersurface CYs.
This was the K3 algorithm that gave the famous 95

Look for d and a[1..5] such that

d = sum_5 a[i]

for each i either a[i] | d or a[i] | d-a[j] for some j

other conditions as appropriate (irreducible, quasismooth)

a[i] | d for all i

Then the question is the well-known question
  1 = sum (reciprocal integers)

We solve this by the tree search

n = 5,
i=1,
a[1],.. a[i-1] already fixed recursively
 remainder = 1 - sum 1/a[j]  (summed over j in [1,.. i-1])
1/a[i] >= remainder/(5-i) so a[i] <= (5-i)/remainder.
fill in a[i] in [a[i-1], (5-i)/remainder]
next i
stop when i = 5 and backtrack.


n := 6;
i := 1;
a := [];
remainder := 1;

// going forward, fill up a[i..5] with largest value
while i lt n do
 a[i] := Floor((n-i)/remainder);
 remainder := remainder - 1/a[i];
 i := i+1;
end while;

MyList := [];
MyOtherList := [];
range := [[2..5]];
for m1 in range[1] do
 range[2] := [Max(-Floor(-1/(1-1/m1))+1,m1)..Floor((5-1)/(1-1/m1))];
 for m2 in range[2] do
  range[3] := [Max(-Floor(-1/(1-1/m1-1/m2))+1,m2)..Floor((5-2)/(1-1/m1-1/m2))];
  for m3 in range[3] do
  range[4] := [Max(-Floor(-1/(1-1/m1-1/m2-1/m3))+1,m3)..Floor((5-3)/(1-1/m1-1/m2-1/m3))];
   for m4 in range[4] do
    x := 1-1/m1-1/m2-1/m3-1/m4;
    N := Numerator(x);
    D := Denominator(x);
    if N eq 1 then
     Append(~MyList,[m1,m2,m3,m4,D]);
     h := LeastCommonMultiple([m1,m2,m3,m4,D]);
     Append(~MyOtherList,[h/D+h/m4+h/m3+h/m2+h/m1,h/D,h/m4,h/m3,h/m2,h/m1]);
    end if;
   end for;
  end for;
 end for;
end for;

for Y in MyOtherList[100..141] do
printf "\tY(%o) in PP(%o,%o,%o,%o,%o)\n",Y[1],Y[2],Y[3],Y[4],Y[5],Y[6];

for Y in MyOtherList do

la := [Q!1 : i in [1..5]]; // The coefficients in sum la[i]*a[i] = d.
range := [2..4];           // The range of values for m1.
for m1 in range[1] do for j1 in [1..5] do
 la1 := la;                // Make a copy of la for editing.
 la1[j1] := la1[j1]-1/m1;
  // New sequence of coeffs la1 for relation sum_>=2 la1[i]*a[i] = d-1/mi.
 top:=; bottom:= Max(m1,),
 range[2] := [bottom, top];
  for m2 in range[2] do
  for j2 in [1..5] do
 

R<x1,x2,x3,x4,x5> := PolynomialRing(Q,5);
TM := []; // "tangent" monomials

 if m1 eq 2 and j1 ne 1 then; else
 TM1 := TM; Append(~TM1,x1^(m1-1)*R.j1); // remember the tangent monomial
 la1 := la;
  if j1 eq 1 then d1 := (1-1/m1);
  else d1 := (1-1/(m1-1)); la1[j1] := la1[j1]-1/(m1-1);
  end if;
 end if;
end for; end for;

 Y(1806) in PP(1,42,258,602,903)
 Y(924) in PP(1,21,132,308,462)
 Y(630) in PP(1,14,90,210,315)
 Y(966) in PP(2,21,138,322,483)
 Y(336) in PP(1,7,48,112,168)
 Y(294) in PP(1,6,42,98,147)
 Y(714) in PP(3,14,102,238,357)
 Y(378) in PP(2,7,54,126,189)
 Y(168) in PP(1,3,24,56,84)
 Y(420) in PP(3,7,60,140,210)
 Y(126) in PP(1,2,18,42,63)
 Y(210) in PP(2,3,30,70,105)
 Y(546) in PP(6,7,78,182,273)
 Y(84) in PP(1,1,12,28,42)
 Y(600) in PP(1,24,75,200,300)
 Y(312) in PP(1,12,39,104,156)
 Y(216) in PP(1,8,27,72,108)
 Y(168) in PP(1,6,21,56,84)
 Y(120) in PP(1,4,15,40,60)
 Y(96) in PP(1,3,12,32,48)
 Y(264) in PP(3,8,33,88,132)
 Y(72) in PP(1,2,9,24,36)
 Y(120) in PP(2,3,15,40,60)
 Y(168) in PP(3,4,21,56,84)
 Y(48) in PP(1,1,6,16,24)
 Y(342) in PP(1,18,38,114,171)
 Y(180) in PP(1,9,20,60,90)
 Y(126) in PP(1,6,14,42,63)
 Y(198) in PP(2,9,22,66,99)
 Y(72) in PP(1,3,8,24,36)
 Y(54) in PP(1,2,6,18,27)
 Y(90) in PP(2,3,10,30,45)
 Y(36) in PP(1,1,4,12,18)
 Y(240) in PP(1,15,24,80,120)
 Y(90) in PP(1,5,9,30,45)
 Y(60) in PP(1,3,6,20,30)
 Y(120) in PP(3,5,12,40,60)
 Y(30) in PP(1,1,3,10,15)
 Y(330) in PP(3,22,30,110,165)
 Y(66) in PP(2,3,6,22,33)
 Y(156) in PP(1,12,13,52,78)
 Y(84) in PP(1,6,7,28,42)
 Y(60) in PP(1,4,5,20,30)
 Y(48) in PP(1,3,4,16,24)
 Y(36) in PP(1,2,3,12,18)
 Y(60) in PP(2,3,5,20,30)
 Y(84) in PP(3,4,7,28,42)
 Y(24) in PP(1,1,2,8,12)
 Y(78) in PP(1,6,6,26,39)
 Y(42) in PP(1,3,3,14,21)
 Y(210) in PP(6,14,15,70,105)
 Y(42) in PP(2,2,3,14,21)
 Y(30) in PP(1,2,2,10,15)
 Y(60) in PP(3,3,4,20,30)
 Y(48) in PP(2,3,3,16,24)
 Y(18) in PP(1,1,1,6,9)
 Y(420) in PP(1,20,84,105,210)
 Y(220) in PP(1,10,44,55,110)
 Y(120) in PP(1,5,24,30,60)
 Y(100) in PP(1,4,20,25,50)
 Y(140) in PP(2,5,28,35,70)
 Y(60) in PP(1,2,12,15,30)
 Y(180) in PP(4,5,36,45,90)
 Y(40) in PP(1,1,8,10,20)
 Y(156) in PP(1,12,26,39,78)
 Y(84) in PP(1,6,14,21,42)
 Y(60) in PP(1,4,10,15,30)
 Y(48) in PP(1,3,8,12,24)
 Y(36) in PP(1,2,6,9,18)
 Y(60) in PP(2,3,10,15,30)
 Y(84) in PP(3,4,14,21,42)
 Y(24) in PP(1,1,4,6,12)
 Y(84) in PP(2,7,12,21,42)
 Y(28) in PP(1,2,4,7,14)
 Y(72) in PP(1,8,9,18,36)
 Y(40) in PP(1,4,5,10,20)
 Y(24) in PP(1,2,3,6,12)
 Y(16) in PP(1,1,2,4,8)
 Y(36) in PP(1,4,4,9,18)
 Y(36) in PP(2,3,4,9,18)
 Y(20) in PP(1,2,2,5,10)
 Y(60) in PP(4,5,6,15,30)
 Y(12) in PP(1,1,1,3,6)
 Y(110) in PP(1,10,22,22,55)
 Y(60) in PP(1,5,12,12,30)
 Y(70) in PP(2,5,14,14,35)
 Y(30) in PP(1,2,6,6,15)
 Y(20) in PP(1,1,4,4,10)
 Y(90) in PP(2,10,15,18,45)
 Y(30) in PP(1,3,5,6,15)
 Y(60) in PP(3,5,10,12,30)
 Y(30) in PP(2,2,5,6,15)
 Y(40) in PP(2,5,5,8,20)
 Y(10) in PP(1,1,1,2,5)
 Y(42) in PP(1,6,7,7,21)
 Y(24) in PP(1,3,4,4,12)
 Y(18) in PP(1,2,3,3,9)
 Y(30) in PP(2,3,5,5,15)
 Y(12) in PP(1,1,2,2,6)
 Y(24) in PP(2,3,3,4,12)
 Y(18) in PP(2,2,2,3,9)
 Y(14) in PP(1,2,2,2,7)
 Y(8) in PP(1,1,1,1,4)
 Y(156) in PP(1,12,39,52,52)
 Y(84) in PP(1,6,21,28,28)
 Y(60) in PP(1,4,15,20,20)
 Y(48) in PP(1,3,12,16,16)
 Y(36) in PP(1,2,9,12,12)
 Y(60) in PP(2,3,15,20,20)
 Y(84) in PP(3,4,21,28,28)
 Y(24) in PP(1,1,6,8,8)
 Y(45) in PP(1,5,9,15,15)
 Y(30) in PP(1,3,6,10,10)
 Y(60) in PP(3,5,12,20,20)
 Y(15) in PP(1,1,3,5,5)
 Y(42) in PP(1,6,7,14,14)
 Y(24) in PP(1,3,4,8,8)
 Y(18) in PP(1,2,3,6,6)
 Y(30) in PP(2,3,5,10,10)
 Y(12) in PP(1,1,2,4,4)
 Y(21) in PP(1,3,3,7,7)
 Y(24) in PP(2,3,3,8,8)
 Y(9) in PP(1,1,1,3,3)
 Y(84) in PP(2,12,21,21,28)
 Y(24) in PP(1,3,6,6,8)
 Y(36) in PP(2,4,9,9,12)
 Y(60) in PP(4,6,15,15,20)
 Y(12) in PP(1,1,3,3,4)
 Y(60) in PP(3,10,12,15,20)
 Y(12) in PP(1,2,2,3,4)
 Y(24) in PP(3,3,4,6,8)
 Y(15) in PP(1,3,3,3,5)
 Y(30) in PP(3,5,6,6,10)
 Y(6) in PP(1,1,1,1,2)
 Y(20) in PP(1,4,5,5,5)
 Y(12) in PP(1,2,3,3,3)
 Y(8) in PP(1,1,2,2,2)
 Y(20) in PP(2,4,4,5,5)
 Y(12) in PP(2,2,2,3,3)
 Y(5) in PP(1,1,1,1,1)


To solve for qsmooth CYs in wPP^4 I want a1,..a5 in QQ with ai > 0 satisfying

(1) sum ai = 1                        [1,1,1,1,1,-1]

(2_i) for each i there exists mi in NN and j in [1..5] such that

  ai(mi-1) + aj = 1               e.g. [m1-1,.. 1,0,.. -1]

Condition (2_i) means that there is a monomial xi^mi or xi^{mi-1}xj of
degree 1. (We treat the weights and degree as homogeneous quantities.)

Plus whatever other conditions are necessary to ensure irreducibility,
well-formed, nonsingularity in codim 1.

Assume
 (3) a1 >= a2 >= .. >= a5 > 0.
We find all solutions according to the following tree-searching scheme:
 (a) make a choice of the monomial (m1,j1) for (2_1) in some range;
 (b) once the first j-1 monomials (m1,j1),.. (m_{j-1},j_{j-1}) are chosen,
I take the j linear relations (1) and (2_j) and reduce them to row echelon
form to determine a range for the next monomial (ai,ji).

First step
==========
The range of (m1,j1): by (1) and (3) it follows that a1 >= 1/5.
This implies m1 in [2..5]. In detail, for the upper bound, if j1 = 1 (so the
monomial is x1^m1), then (2_1) gives
  m1 = 1/a1 <= 5;

if j1 <> 1 (so the monomial is x1^{m1-1}xj), then (2_1) gives
  m1 = (1-aj)/a1 + 1 < 5 + 1, so m1 <= 5.

For the lower bound, m1 = 1 would give that one of the aj = 1, so all the
others have to be zero; therefore m1 >= 2.

N.B. Also if m1 = 2 then j1 = 1. For if m1 = 2 and j1 <> 1 then a1 + aj = 1,
so again all the other ai have to be zero. However, this mixes up two steps;
it is more systematic to take a more stupid lower bound, and allow the
EchelonForm argument below to give a contradition.

This gives the outer loop

range := [2,5];
for m1 in range[1] do for j1 in [1..5]
 // in principle I could also have: if m1 eq 2 and j1 gt 1 then; else end if;
 //
 // the internal loops for m2, etc.
 end for; end for;

Second step
===========
I subtract a multiple of (2_1) from (1) to get rid of a1; this gives one
linear relation involving a2,.. a5, and I need to determine the range for
m2, j2. This is the first step of a systematic row echelon procedure. (Later
I will have several relations (2_i).)

// m1 := 2; j1 := 2; leads to a contradiction
M := Matrix(2,6,[1,1,1,1,1,-1,m1-1,0,0,0,0,-1]);
M[2,j1] := M[2,j1] + 1;
EchelonForm(M);

e.g., if the first monomial was x^2z then the second row of M is
 [0,2,1,2,2,-1]; therefore a2 >= 1/7, giving m2 <= 7

upperbound := (&+[EchelonForm(M)[2,i] : i in [2..5]])/(-M[2,6]);

As lowerbound for m2,j2?
Necessary conditions: m2 >= m1 and in case of equality j2 > j1.
There is no need to be any cleverer than that at this stage.

range[2] := [m1,upperbound];

So I write the outer nested 2 loops

 M := Matrix(2,6,[1,1,1,1,1,-1,m1-1,0,0,0,0,-1]);
 M[2,j1] := M[2,j1] + 1;
 N := EchelonForm(M);
 if N[2,6] ge 0 then; else
 upperbound := (&+[N[2,i] : i in [2..5]])/(-N[2,6]);
 range[2] := [m1,upperbound];
 for m2 in range[2] do for j2 in [1..5] do

// the internal loops for m3, etc.

  end for; end for; end if; // loop over m2, j2
 end for; end for; // loop over m1, j1

The framework is now clear: I write the 3 outer loops

 // prepare the range for the loop over m2, j2

 // start of the loop over m2, j2
   // should require m2 > m1 or j2 >= j1 to avoid repetition in answer

  // prepare the range for the loop over m3, j3
  M := Matrix(3,6,[1,1,1,1,1,-1, m1-1,0,0,0,0,-1, 0,m2-1,0,0,0,-1]);
  M[2,j1] := M[2,j1] + 1;
  M[3,j2] := M[2,j2] + 1;
  N := EchelonForm(M);
  if N[3,6] ge 0 then; else
  upperbound := (&+[N[3,i] : i in [3..5]])/(-N[3,6]);
  range[3] := [m2,upperbound];
  for m3 in range[3] do for j3 in [1..5] do
   // should require m3 > m2 or j3 >= j2 to avoid repetition

// the internal loops for m4, etc.

   end for; end for; end if; // ends loop over m3, j3
  end for; end for; end if; // ends loop over m2, j2
 end for; end for; // ends loop over m1, j1

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

Now write out the inner loop for m4, m5

  // prepare the range for the loop over m4, j4
   M := Matrix(4,6,[1,1,1,1,1,-1, m1-1,0,0,0,0,-1,
     0,m2-1,0,0,0,-1, 0,0,m3-1,0,0,-1]);
   M[2,j1] := M[2,j1] + 1;
   M[3,j2] := M[3,j2] + 1;
   M[4,j3] := M[4,j3] + 1;
   N := EchelonForm(M);
   if N[4,6] ge 0 then; else
   upperbound := Floor((&+[N[4,i] : i in [4..5]])/(-N[4,6]));
   range[4] := [m3,upperbound];

   // start of the loop over m4, j4
   for m4 in range[4] do for j4 in [1..5] do
    // should require m3 > m2 or j3 >= j2 to avoid repetition
  
  // whatever we do for m5, etc.

    end for; end for; // ends loop over m4, j4
   end for; end for; // ends loop over m3, j3
  end for; end for; // ends loop over m2, j2

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

ListY := [[1,1,1,1,1,5]];
ListTM := [[x^5,y^5,z^5,t^5,u^5]];

m1 := 3; j1 := 4;
m2 := 4; j2 := 1;
m3 := 5; j3 := 5;
TM := [RR.1^(m1-1)*RR.j1,RR.2^(m2-1)*RR.j2,RR.3^(m3-1)*RR.j3];

range:=[5..20];
for m4 in range /*[4]*/ do for j4 in [1..5] do
 TM[4] := RR.4^(m4-1)*RR.j4;
 M := Matrix(5,6,[1,1,1,1,1,-1, m1-1,0,0,0,0,-1,
 0,m2-1,0,0,0,-1, 0,0,m3-1,0,0,-1, 0,0,0,m4-1,0,-1]);
 M[2,j1] := M[2,j1] + 1;
 M[3,j2] := M[3,j2] + 1;
 M[4,j3] := M[4,j3] + 1;
 M[5,j4] := M[5,j4] + 1;
  N := EchelonForm(M);
   if N[5,6] ge 0 then // tests monomials so far incompatible
   printf "Bust: %o, TMons = %o,%o,%o,%o.\n", N[5],TM[1],TM[2],TM[3],TM[4];
   else if IsDivisibleBy(N[5,5],N[5,6])
    then // tries pure power
     TM[5] := u^(N[5,5] div -N[5,6]);
     Y := Eltseq(Basis(Kernel(Transpose(N)))[1]);
     Append(~ListY,Y); Append(~ListTM,TM); // this is the output
    else
     for j in [1..4] do //loop over j = j5
      if N[j,j] eq - N[j,6] then;
        Num := -N[j,5]*N[5,6]+N[j,6]*N[5,5]+N[j,j]*N[5,5];
        Denom := -N[j,j]*N[5,6];
        if IsDivisibleBy(Num,Denom) then  // ask: can we define m5 as integer?
          m5 := (Num div Denom) + 1;
          TM[5] := u^(m5-1)*RR.j;
     Append(~ListY,Y); Append(~ListTM,TM);
       end if; // ends test of m5
      end if; //  ends treatment of njj <> -nj6
     end for; // ends loop over j = j5
     end if; // ends pure power case division
     end if; // ends treatment of monomials if they are compatible
end for; end for; // ends loop over m4, j4

Handling output
===============

#ListY;

for i in [100..150] do
Y := ListY[i]; f := &+ListTM[i];
if &or[Y eq ListY[j] : j in [2..i-1]] then; else
printf "\tY(%o) in PP(%o,%o,%o,%o,%o), %o\n",
  Y[6],Y[1],Y[2],Y[3],Y[4],Y[5],f;
end if; end for;

i:=5;
PY<xx,yy,zz,tt,uu> := ProjectiveSpace(Q, [ListY[i][j] : j in [1..5]]);
Y := Scheme(PY, &+ ListTM[i]) ;

i:=5;
PY<xx,yy,zz,tt,uu> := ProjectiveSpace(Q, [ListY[i][j] : j in [1..5]]);
Y := Scheme(PY, &+ ListTM[i]) ;
IsNonSingular(Y);
true

Useful hints
============

Eltseq(Basis(Kernel(Transpose(N))));
Eltseq($1[1]);  

Nt := Transpose(N);
XX,YY:=Solution(Nt, Codomain(Nt)!0);
YY;
RSpace of degree 6, dimension 1 over Integer Ring
Echelonized basis:
( 8  3  4  1  1 17)

List of successes:

Y(81) in PP(36,15,20,9,1), x^2*t + x*y^3 + x*t^5 + z^4*u + u^81
Y33) in PP(15,6,8,3,1), x^2*t + x*y^3 + x*t^6 + z^4*u + u^33
Y129) in PP(57,24,32,15,1), x^2*t + x*y^3 + y*t^7 + z^4*u + u^129
Y65) in PP(29,12,16,7,1), x^2*t + x*y^3 + z^4*u + z*t^7 + u^65
Y49) in PP(22,9,12,5,1), x^2*t + x*y^3 + y*t^8 + z^4*u + u^49
Y148) in PP(67,27,36,14,4), x^2*t + x*y^3 + z^4*u + z*t^8 + u^37
Y17) in PP(8,3,4,1,1), x^2*t + x*y^3 + x*t^9 + z^4*u + u^17
Y180) in PP(81,33,44,18,4), x^2*t + x*y^3 + z^4*u + t^10 + u^45
Y146) in PP(65,27,36,16,2), x^2*t + x*y^3 + z^4*u + t^9*u + u^73
Y184) in PP(85,33,44,14,8), x^2*t + x*y^3 + z^4*u + z*t^10 + u^23
Y82) in PP(37,15,20,8,2), x^2*t + x*y^3 + z^4*u + t^10*u + u^41
Y216) in PP(99,39,52,18,8), x^2*t + x*y^3 + z^4*u + t^12 + u^27
Y50) in PP(23,9,12,4,2), x^2*t + x*y^3 + z^4*u + t^12*u + u^25
Y84) in PP(39,15,20,6,4), x^2*t + x*y^3 + z^4*u + t^14 + u^21
Y256) in PP(121,45,60,14,16), x^2*t + x*y^3 + z^4*u + z*t^14 + u^16
Y288) in PP(135,51,68,18,16), x^2*t + x*y^3 + z^4*u + t^16 + u^18
Y120) in PP(57,21,28,6,8), x^2*t + x*y^3 + z^4*u + t^20 + u^15

The first of these is singular as written because 3 of its monomials are
divisible by x. However, t^9 works instead of x*t^5

Y(81) in PP(36,15,20,9,1), x^2*t + x*y^3 + t^9 + z^4*u + u^81

[i : i in [1..#ListY] | ListY[i] eq [36,15,20,9,1,81]];

=================================

Round 3
-------

To avoid all that messing with conditionals both before after the
choice of monomials. and to avoid repetitions in the output.

Think of [m1,j1] in range1 as one object, where range1 is a list of lists.

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

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;

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

TM := [RR|];

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]];

for mon1 in range1 do
 m1 := mon1[1]; j1 := mon1[2];
 TM[1] := RR.1^(m1-1)*RR.j1;
 M := Matrix(2,6,[1,1,1,1,1,-1, m1-1,0,0,0,0,-1]);
 M[2,j1] := M[2,j1]+1;
 N := EchelonForm(M);
 if N[2,6] ge 0 then
  max2 := -1; else
  max2 := Floor((&+[N[2,i] : i in [2..5]])/(-N[2,6])) +1;
 end if;
 range2 := [[m,j] : j in [1..5], m in [2 .. max2] // | such that more conds
     ];

  for mon2 in range2 do
  m2 := mon2[1]; j2 := mon2[2];
  TM[2] := RR.2^(m2-1)*RR.j2;
  M := Matrix(3,6,[1,1,1,1,1,-1, m1-1,0,0,0,0,-1, 0,m2-1,0,0,0,-1]);
  M[2,j1] := M[2,j1]+1;
  M[3,j2] := M[3,j2]+1;
  N := EchelonForm(M);
  if N[3,6] ge 0 then
    max3 := -1; else
    max3 := Floor((&+[N[3,i] : i in [3..5]])/(-N[3,6])) + 1;
   end if;
  range3 := [[m,j] : j in [1..5], m in [2 .. max3] // | such that more conds
     ];

   for mon3 in range3 do
   m3 := mon3[1]; j3 := mon3[2];
   TM[3] := RR.3^(m3-1)*RR.j3;
   M := Matrix(4,6,[1,1,1,1,1,-1, m1-1,0,0,0,0,-1, 0,m2-1,0,0,0,-1,
    0,0,m3-1,0,0,-1]);
   M[2,j1] := M[2,j1]+1;
   M[3,j2] := M[3,j2]+1;
   M[4,j3] := M[3,j3]+1;
   N := EchelonForm(M); // should check that N[4,6] < 0
   if N[4,6] ge 0 then
    printf "Bust after Mon3: TMons = %o,%o,%o.\n", TM[1],TM[2],TM[3];
    max4 := -1; else
    max4 := Floor((&+[N[4,i] : i in [4..5]])/(-N[4,6])) + 1;
   end if;
   range4 := [[m,j] : j in [1..5], m in [2 .. max4] // | more conds
      ];

    for mon4 in range4 do
    m4 := mon4[1]; j4 := mon4[2];
    TM[4] := RR.4^(m4-1)*RR.j4;
    M := Matrix(5,6,[1,1,1,1,1,-1, m1-1,0,0,0,0,-1,
    0,m2-1,0,0,0,-1, 0,0,m3-1,0,0,-1, 0,0,0,m4-1,0,-1]);
    M[2,j1] := M[2,j1] + 1;
    M[3,j2] := M[3,j2] + 1;
    M[4,j3] := M[4,j3] + 1;
    M[5,j4] := M[5,j4] + 1;
    N := EchelonForm(M);
    if N[5,6] ge 0 then; // tests monomials so far incompatible
 //  printf "Bust: %o, TMons = %o,%o,%o,%o.\n", N[5],TM[1],TM[2],TM[3],TM[4];
    else
     if IsDivisibleBy(N[5,5],N[5,6]) then // pure power works
      TM[5] := u^(N[5,5] div -N[5,6]);
       if Pass(TM) then
        Y := Eltseq(Basis(Kernel(Transpose(N)))[1]);
        Append(~ListY,Y); Append(~ListTM,TM); // this is the output
       end if;
     else
      for j in [1..4] do //loop over j = j5
       if N[j,j] eq - N[j,6] then; // this avoids a trivial false solution
        Num := -N[j,5]*N[5,6]+N[j,6]*N[5,5]+N[j,j]*N[5,5];
        Denom := -N[j,j]*N[5,6];
        if IsDivisibleBy(Num,Denom) then // ask: can we define m5 as integer?
         m5 := (Num div Denom) + 1;
         TM[5] := u^(m5-1)*RR.j;
          if Pass(TM) then
           Y := Eltseq(Basis(Kernel(Transpose(N)))[1]);
           Append(~ListY,Y); Append(~ListTM,TM); // this is the output
          end if;
        end if; // ends test of m5
       end if; // test of njj <> -nj6
      end for; // ends loop over j = j5
     end if; // ends pure power case division
    end if; // ends treatment of monomials if they are compatible

   end for; // ends mon4 in range4
  end for; // ends mon3 in range3
 end for; // ends mon2 in range2
end for; // ends mon1 in range1

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

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

Good := [i : i in Pass | &and [ListY[i] ne ListY[j] : j in Pass | j lt i ]];
  // impossibly slow: it is comparing 10045 elements for equality
#Good; // 5244

for i in [100*j + 7 : j in [1..20]] do
 Y := ListY[i]; f := &+ListTM[i];
 printf "\tY(%o) in PP(%o,%o,%o,%o,%o), %o\n",
  Y[6],Y[1],Y[2],Y[3],Y[4],Y[5],f;
end for;

 Y(52) in PP(13,16,11,8,4), x^4 + x*t^7 + y^3*u + z^4*t + z*u^6
 Y(60) in PP(14,19,18,6,3), x^3*z + y^3*u + z^3*t + t^10 + u^20
 Y(40) in PP(13,10,8,8,1), x^3*u + y^4 + z^4*t + z*t^4 + u^40
 Y(27) in PP(9,9,6,2,1), x^3 + x*z^3 + y^3 + y*t^9 + u^27
 Y(357) in PP(119,134,89,14,1), x^3 + x*t^17 + y^2*z + z^4*u + u^357
 Y(336) in PP(112,152,32,37,3), x^3 + x*z^7 + y^2*z + t^9*u + u^112
 Y(36) in PP(12,17,4,2,1), x^3 + x*z^6 + y^2*t + z*t^16 + u^36
 Y(90) in PP(30,43,7,4,6), x^3 + x*t^15 + y^2*t + z^12*u + u^15
 Y(156) in PP(52,72,13,7,12), x^3 + x*z^8 + y^2*u + y*t^12 + u^13
 Y(117) in PP(39,26,36,7,9), x^3 + x*y^3 + y*t^13 + z^3*u + u^13
 Y(90) in PP(30,25,15,14,6), x^3 + x*z^4 + y^3*z + t^6*u + u^15
 Y(36) in PP(12,11,6,4,3), x^3 + x*z^4 + y^3*u + t^9 + u^12
 Y(68) in PP(21,26,16,4,1), x^2*y + y^2*z + z^4*t + t^17 + u^68
 Y(260) in PP(68,124,51,12,5), x^2*y + x*t^16 + y^2*t + z^5*u + u^52
 Y(99) in PP(25,49,15,9,1), x^2*y + y^2*u + z^6*t + t^11 + u^99
 Y(144) in PP(54,36,47,3,4), x^2*y + x*t^30 + y^4 + z^3*t + u^36
 Y(150) in PP(54,42,24,25,5), x^2*y + x*z^4 + y^3*z + t^6 + u^30
 Y(35) in PP(14,7,6,3,5), x^2*y + x*t^7 + y^5 + z^5*u + u^7
 Y(127) in PP(53,37,21,15,1), x^2*z + x*y^2 + y*t^6 + z^6*u + u^127
 Y(672) in PP(253,224,166,21,8), x^2*z + y^3 + z^4*u + t^32 + u^84
 Y(135) in PP(63,45,9,13,5), x^2*z + x*z^8 + y^3 + t^10*u + u^27
 Y(80) in PP(27,20,26,5,2), x^2*z + y^4 + y*t^12 + z^3*u + u^40
 Y(292) in PP(110,93,72,13,4), x^2*z + x*t^14 + y^3*t + z^4*u + u^73
 Y(171) in PP(76,54,19,13,9), x^2*z + x*z^5 + y^3*u + y*t^9 + u^19
 Y(72) in PP(29,15,14,12,2), x^2*z + y^4*t + z^5*u + t^6 + u^36
 Y(54) in PP(26,18,7,2,1), x^2*t + x*z^4 + y^3 + y*t^18 + u^54
 Y(84) in PP(39,28,5,6,6), x^2*t + x*z^9 + y^3 + t^14 + u^14
 Y(52) in PP(22,13,5,8,4), x^2*t + x*z^6 + y^4 + t^6*u + u^13
 Y(256) in PP(120,85,34,16,1), x^2*t + x*z^4 + y^3*u + t^16 + u^256
 Y(37) in PP(17,7,9,3,1), x^2*t + y^4*z + y*t^10 + z^4*u + u^37
 Y(390) in PP(194,130,49,15,2), x^2*u + x*z^4 + y^3 + t^26 + u^195
 Y(159) in PP(79,53,19,7,1), x^2*u + y^3 + z^8*t + z*t^20 + u^159
 Y(120) in PP(55,30,13,12,10), x^2*u + x*z^5 + y^4 + t^10 + u^12
 Y(68) in PP(32,19,11,2,4), x^2*u + x*t^18 + y^3*z + z^6*t + u^17
 Y(105) in PP(51,21,25,5,3), x^2*u + y^5 + z^4*t + z*t^16 + u^35
 Y(336) in PP(168,84,67,16,1), x^2 + x*y^2 + z^5*u + t^21 + u^336
 Y(600) in PP(300,200,75,24,1), x^2 + x*z^4 + y^3 + t^25 + u^600
 Y(1602) in PP(801,534,223,41,3), x^2 + y^3 + z^7*t + t^39*u + u^534
 Y(336) in PP(168,112,41,7,8), x^2 + x*t^24 + y^3 + z^8*u + u^42
 Y(330) in PP(165,110,28,22,5), x^2 + y^3 + z^11*t + t^15 + u^66
 Y(186) in PP(93,31,26,30,6), x^2 + x*y^3 + z^6*t + t^6*u + u^31
 Y(414) in PP(207,115,69,17,6), x^2 + x*z^3 + y^3*z + t^24*u + u^69
 Y(108) in PP(54,33,9,8,4), x^2 + x*z^6 + y^3*z + t^13*u + u^27
 Y(182) in PP(91,54,20,15,2), x^2 + x*t^13 + y^3*z + y*u^12 + z^9*u
 Y(324) in PP(162,105,46,9,2), x^2 + x*t^18 + y^3*t + z^7*u + u^162
 Y(184) in PP(92,61,23,7,1), x^2 + x*z^4 + y^3*u + z*t^23 + u^184
 Y(96) in PP(48,31,6,8,3), x^2 + x*z^8 + y^3*u + t^12 + u^32
 Y(16) in PP(8,2,3,2,1), x^2 + x*y^4 + y*t^7 + z^5*u + u^16
 Y(200) in PP(100,40,33,25,2), x^2 + y^5 + z^6*u + t^8 + u^100
 Y(152) in PP(76,32,19,24,1), x^2 + x*z^4 + y^4*t + y*t^5 + u^152
