{bounds=[6,4,3,3];tot=0;

for(r1=1,6,for(r2=0,r1,for(r3=0,r2,for(r4=0,r3,b1=[r1,r2,r3,r4];
   if(b1*b1~<=bounds[1],print(b1);\\bound on |b1|^2<=6
      bs1=b1;   \\beginning of Gram-Schmidt
      if(b1[2]==0,q2=0,q2=1);if(b1[3]==0,q3=0,q3=1);if(b1[4]==0,q4=0,q4=1); \\If a coordinate of the first vector is 0, this guarantees that the corresponding coordinate of the second vector is non-negative
      for(s1=-sqrtint(bounds[2]),sqrtint(bounds[2]),for(s2=-q2*sqrtint(bounds[2]-s1^2),sqrtint(bounds[2]-s1^2),for(s3=-q3*sqrtint(bounds[2]-s1^2-s2^2),sqrtint(bounds[2]-s1^2-s2^2),for(s4=-q4*sqrtint(bounds[2]-s1^2-s2^2-s3^2),sqrtint(bounds[2]-s1^2-s2^2-s3^2),b2=[s1,s2,s3,s4];
          if(b2*b2~<=bounds[2],\\bound on |b2|^2<=4
             mu21=(b2*bs1~)/(bs1*bs1~);bs2=b2-mu21*bs1;\\compute mu21 and second gram-schmidt orthogonalized vector.
             if(bs2*bs2~<=3*2^(4-2-2),\\check bound on bs2 from corollary 4.4
             if(abs(mu21)<=1/2,\\check condition 1 of LLL-reduction
                check=bs2+mu21*bs1;
                if(check*check~>=3/4*(bs1*bs1~),\\check condition 2 of LLL-reduction
                   for(t1=-sqrtint(bounds[3]),sqrtint(bounds[3]),for(t2=-sqrtint(bounds[3]),sqrtint(bounds[3]),for(t3=-sqrtint(bounds[3]),sqrtint(bounds[3]),for(t4=-sqrtint(bounds[3]),sqrtint(bounds[3]),b3=[t1,t2,t3,t4];
                       if(b3*b3~<=bounds[3],\\bound on |b3|^2<=3
                          mu32=(b3*bs2~)/(bs2*bs2~);mu31=(b3*bs1~)/(bs1*bs1~);
                          bs3=b3-mu32*bs2-mu31*bs1;\\compute third vector in gram-schmidt
                          if(abs(mu31)<=1/2,\\check first condition in LLL-reduced
                             if(abs(mu32)<=1/2,\\check first condition in LLL-reduced
                                check=bs3+mu32*bs2;
                                if(check*check~>=3/4*(bs2*bs2~),\\check second condition in LLL-reduced
                                if(bs3*bs3~<=3*2^(4-3-2),\\check bound on bs3 from corollary 4.4
                                   for(u1=-sqrtint(bounds[4]),sqrtint(bounds[4]),for(u2=-sqrtint(bounds[4]),sqrtint(bounds[4]),for(u3=-sqrtint(bounds[4]),sqrtint(bounds[4]),for(u4=-sqrtint(bounds[4]),sqrtint(bounds[4]),b4=[u1,u2,u3,u4];
                                       if(b4*b4~<=bounds[4],
                                          mu43=(b4*bs3~)/(bs3*bs3~);mu42=(b4*bs2~)/(bs2*bs2~);mu41=(b4*bs1~)/(bs1*bs1~);bs4=b4-mu43*bs3-mu42*bs2-mu41*bs1;
                                          if(abs(mu43)<=1/2,
                                             if(abs(mu42)<=1/2,
                                                if(abs(mu41)<=1/2,
                                                   check=bs4+mu43*bs3;
                                                   if(check*check~>=3/4*(bs3*bs3~),
                                                      if(bs4*bs4~<=1,
                                                      M=[b1[1],b1[2],b1[3],b1[4];b2[1],b2[2],b2[3],b2[4];b3[1],b3[2],b3[3],b3[4];b4[1],b4[2],b4[3],b4[4]];
                                                      if(abs(matdet(M))>1,
                                                         tot=tot+1;print(tot,M);
                                                         N=M^(-1);
                                                         P=matrix(4,4,i,j,floor(abs(N[i,j])));
                                                         S=(P[,1]~*P[,1])*(P[,2]~*P[,2])*(P[,3]~*P[,3])*(P[,4]~*P[,4]);
                                                         if(S==0,,print("counterexample: ",M)))))))))))))))))))))))))))))))))));tot
}