echo = 2;

LIB("syzextra.so");

noop();

system("--min-time", "0.01");
system("--ticks-per-sec", 100);

LIB "poly.lib"; // for numerator & denominator

option(redSB); option(redTail); // assumed for the results of kStd!
option(prot);
option(mem);
option(notWarnSB);


// option(noloadLib); option(noredefine);

LIB "teachstd.lib";

int reduce_time = 0;

///////////////////////////////////////////////////////////////////////////////
proc MyReduce(f, G)
{
  int t = timer;
  def g = reduce_syz(f, G, 0);
  int tt = timer;

  reduce_time = reduce_time + (tt - t);

  return(g);
};

proc separateSyzGB( module J, int c )
{
  module II, G; vector v; int i;

  J = simplify(J, 2);

  for( i = ncols(J); i > 0; i-- )
  {
    v = J[i];
    if( leadcomp(v) > c )
    {
      II[i] = v;
    } else
    {
      G[i] = v; // leave only gen(i): i <= c
    }
  }

  II = simplify(II, 2);
  G = simplify(G, 2);

  return (list(G, II));
}



proc splitSyzGB( module J, int c )
{
  module JJ; vector v, vv; int i;

  for( i = ncols(J); i > 0; i-- )
  {
    v = J[i];

    vv = 0;
    
    while( leadcomp(v) <= c )
    {
      vv = vv + lead(v);
      v  = v  - lead(v);
    }

    J[i] = vv;
    JJ[i] = v;
  }

  J = simplify(J, 2);
  JJ = simplify(JJ, 2);

  return (list(J, JJ));
}



proc prepareSyz( module I, list # )
{
  int i;
  int k = 0;
  int r = nrows(I);
  int c = ncols(I);


  if( size(#) > 0 )
  {
    if( typeof(#[1]) == "int" || typeof(#[1]) == "bigint" )
    {
      k = #[1];
    }
  }

  if( k < r )
  {
    "// *** Wrong k: ", k, " < nrows: ", r, " => setting k = r = ", r;
    k = r;
  }

//   "k: ", k;  "c: ", c;   "I: ", I;

  for( i = c; i > 0; i-- )
  {
    I[i] = I[i] + gen(k + i);
  }

//  DetailedPrint(I);

  return(I);
}


/// is p - permissible?
proc myPrimeTest(def I, int p)
{
  int i; def v; number c;

  for(i = ncols(I); i > 0; i--)
  {
    v = I[i];
    while(v != 0)
    {
      c = leadcoef(v);
      
      if( (numerator(c) mod p) == 0) { return(0); }
      if( (denominator(c) mod p) == 0) { return(0); }

      v = v - lead(v);
    }
  }
  return(1);
}

proc myPrimeList(module I, int n, list #)
{
  intvec L;
  int i,p;

  if(size(#) == 0)
  {
    p = 2147483647;

    while(myPrimeTest(I,p)==0)
    {
      p = prime(p-1);
      if(p == 2) { ERROR("no more primes"); }
    }
    L[1] = p;
  }
  else
  {
    L = #[1];
    p = prime(L[size(L)]-1);
    while(!myPrimeTest(I,p))
    {
      p = prime(p-1);
      if(p == 2) { ERROR("no more primes"); }
    }
    L[size(L)+1] = p;
  }

  if(p == 2) { ERROR("no more primes"); }
  for(i = 2; i <= n; i++)
  {
    p = prime(p-1);
    while(!myPrimeTest(I,p))
    {
      p = prime(p-1);
      if(p == 2) { ERROR("no more primes"); }
    }
    L[size(L)+1] = p;
  }

  return(L);
}



////////////////////////////////////////////////////////////////////////////////

// int try = 0;

proc GBCandidate0( module M, def P)
{

  def Q = basering;

  int iComp = attrib(basering, "SyzComp");
  
  setring P;
  module M = imap(Q, M);

  M = idPrepare(M, iComp);

/*  
  try  ++;
  
  if( try != 5 )
  {
    return (0);
  }
*/
  
  M = simplify(M, 2);
  
  setring Q;

  M = imap(P, M);
  
  return (M);
}



proc Noro( module M )
{
  def save = basering;


  if( char(save) != 0 )
  {
    ERROR("Noro: char != 0!");
  }

  int iComp = nrows(M);

  
 def Q = MakeSyzCompOrdering(); setring Q;
//  def Q = MakeInducedSchreyerOrdering(); setring Q;


  SetSyzComp(iComp);
  attrib(Q, "SyzComp", iComp);
  DetailedPrint(basering);

  list L = ringlist(basering);

  def M = imap(save, M);

  M = prepareSyz(M, iComp);

  DetailedPrint(basering);

  intvec pp = myPrimeList(M, 10);
  
//  intvec pp = 99981793+1;
  
//  intvec pp = 32003 + 1; pp = myPrimeList(M, 10, pp);
  
  int p;
  int i = 1;

  int try = 0;
  
  while(1)
  {
    try ++;

    // choose p:
    p = pp[i];

    "Current prime: ", p;
      
    if( typeof(L[1]) == "int" || typeof(L[1]) == "bigint" )
    {
      L[1] = p;
    }
    else
    {
      if( typeof(L[1][1]) == "int" || typeof(L[1][1]) == "bigint" )
      {
        L[1][1] = p;
      }
      else
      {
        ERROR("NORO: cannot create p-ring list: wrong input ring?");
      }
    }

    // new ring over F_p
    def P = ring(L); setring P;

   SetSyzComp(iComp); 
//    DetailedPrint(basering);

    setring Q;

//    M;    print(M);

    // Compute GB Candidate for M using the ring P / F_p
    "Computing GB Candidate: ";
    def result = GBCandidate(M, P);

    if( defined(result) && (size(result) > 0) && (typeof(result) == "list") && (result[1] == "ok") )
    {
      "===================================================================";
      "It is DONE!";
      "try: ", try;

      result[2];
      print(result[2]);
//      DetailedPrint(result[2]);

      list GB_SYZ = separateSyzGB(result[2], iComp);

      module GB  = GB_SYZ[1]; //splitSyzGB(, iComp)[1]; 
      module SYZ = GB_SYZ[2];


      GB;
      SYZ;

      

      setring save;

      
      module GB  = imap(Q, GB);

      if( size(GB) > 0 )
      {
        // need the top part!
        GB = transpose(GB);
        GB = GB[1..iComp];
        GB = transpose(GB);
      }

      
      module SYZ = imap(Q, SYZ);

      if( size(SYZ) > 0 )
      {
        int r = nrows(SYZ);

        // need the bottom part:
        SYZ = transpose(SYZ);
        SYZ = SYZ[iComp+1 .. r];
        SYZ = transpose(SYZ);
      }
      
      return (list((GB), (SYZ)));
    }

    kill P;
    setring Q;

    i ++;

    if( i > size(pp) )
    {
      pp = myPrimeList(M, 10, pp);
      pp;
    }

        
    if( !defined(pp) )
    {
      ERROR("NORO: Sorry no more primes!");
      return(0);
    }
  } // while(1)
}
    

// basering is assumed to be Z[x_1, ..., x_n]

// F is a list of vectors of the same length
proc GBCandidate (module F, def Fp)
{
  int ss = 1; // symmetric S-polynomial
  
  // TODO: check input for errors
  def br = basering;

  int i, j;

  int sizeF = size(F);

  intvec C;
  list D; // list of intvecs of length 2, as indices of F
  for(i = 1; i <= sizeF; i++)
  {
    for(j = i+1; j <= sizeF; j++)
    {
      D = insert(D, intvec(i, j));
    }
  }

  vector h, s;
  number c;
  
  setring Fp;
  vector h, s;
  module F;
  
  int p = char(Fp);
  // TODO: check that characteristic of Fp is permissible for F

  

  while(size(D) > 0)
  {
    
    if(char(basering) != p )
    {
      "ERROR: wrong current ring!";
      $$$
    }
    
    C = D[size(D)];
    D = delete(D, size(D));

//    "Pair: (", C[1], ",", C[2], ")";
    
    F = fetch(br, F);
    s = spoly(F[C[1]], F[C[2]], ss);
    h = MyReduce(s, F);
    
    if(h != 0)
    {
      setring br;

      s = spoly(F[C[1]], F[C[2]], ss);
      h = MyReduce(s, F); // h \in Z_<p>[x] ???

      c = leadcoef(h);

      if(h != 0)
      {
        if( (denominator(c) mod p) == 0 )
        {
          "h: ", h;
          "c: ", c;
          "ERROR: GBCandidate: wrong coeff!";
          $$
        }
      }
      
      if(h != 0 && (numerator(c) mod p) != 0)
      {
        for(i = 1; i <= sizeF; i++)
        {
          D = insert(D, intvec(i, sizeF+1));
        }
        sizeF++;
        F[sizeF] = h;

//        "new F[",sizeF,"] element: ", h;
        
        setring Fp;
        
        h = fetch(br, h);
        F[sizeF] = h;
      }
      else
      {
        /*
        "&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&";
        matrix T; //  = lift(F, s);
        basering; // /Q
        
        test(23+32);
        option(teach);
        
        "i, j: ", C[1], C[2];
        

        F;
        "f: ", F[C[1]];
        "g: ", F[C[2]];
        "s: ", s;
        "ZERO::: h: ", h; // 0!!!

        vector m = NFMora(s, F);
        "MORA: ", m;

//        matrix(s)-matrix(F)*T; print(T);

        MyReduce(s, F);
        "h-test: ", _ == h;
        
        "????";

        setring Fp;

        "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
        basering; // /p
        

        F;
        "F-proj-test: ", size(module(matrix(F) - matrix(fetch(br, F)))) == 0;
        
        "f: ", F[C[1]];
        "g: ", F[C[2]];
        "s: ", s;

        "s-proj-test: ", s == fetch(br, s);
        
        "NON-ZERO::: h: ", h;

//        matrix T = fetch(br, T);        def S = (matrix(F)*T);        print(S);        S;
        
        vector m = NFMora(s, F);
        "MORA: ", m;

        "m-proj-test: ", m == fetch(br, m);
        
        
        MyReduce(s, F);
        "h-test: ", _ == h;

        "????";
        $$$
        */    
        return(list("ng", module()));
      }
    }
  }
  setring br;
  F = interred(F); // ?
  return(list("ok", F));
}




echo = 1;


ring R = 0, (w, x, y, z), dp;
module I = w^2 - x*z,  w*x - y*z,  x^2 - w*y, x*y - z^2, y^2 - w*z;



/*
ring r = 0,(x,y,z),(c,dp);
module I = [x+1, y, 1], [xy, z, z2]; // NO SYZ!
*/


/*
// Test: kotsireas 
ring R = (0),(B, b, D, d, F, f),dp;
module I = B*b-b*D-B*d+D*d-2*F+2,B*b+b*D-B*d-D*d-2*b*F+2*d*F-2*B+2*D,b^2-2*b*d+d^2-2*b-2*d+f+1,B^2*b^3-1,D^2*d^3-1,F^2*f^3-1; 
*/



/*
// Cohn3
ring R = (0),(x1, y, z, t),dp;
module I = -x1^3*y^2+2*x1^2*y^2*z-x1^2*y*z^2-144*x1^2*y^2-207*x1^2*y*z+288*x1*y^2*z+78*x1*y*z^2+x1*z^3-3456*x1^2*y-5184*x1*y^2-9504*x1*y*z-432*x1*z^2-248832*x1*y+62208*x1*z-2985984*x1,-x1^3*z*t^2+x1^2*z^2*t^2-6*x1^3*z*t+4*x1^2*z^2*t+32*x1^3*t^2-72*x1^2*z*t^2-87*x1*z^2*t^2-z^3*t^2-8*x1^3*z-432*x1^2*z*t-414*x1*z^2*t+2592*x1*z*t^2+864*z^2*t^2-1728*x1^2*z-20736*x1*z*t+3456*z^2*t-186624*z*t^2-124416*x1*z-1492992*z*t-2985984*z,x1^2*y*t^3-2*x1*y^2*t^3+y^3*t^3+8*x1^2*y*t^2-12*x1*y^2*t^2+4*y^3*t^2-24*x1*y*t^3+24*y^2*t^3+20*x1^2*y*t-20*x1*y^2*t-160*x1*y*t^2+96*y^2*t^2+128*x1*t^3+16*x1^2*y+96*x1*y*t+2304*x1*t^2+1152*x1*y+13824*x1*t+27648*x1,y^3*t^3-y^2*z*t^3+4*y^3*t^2-2*y^2*z*t^2+72*y^2*t^3+71*y*z*t^3+288*y^2*t^2+360*y*z*t^2+6*z^2*t^2+1728*y*t^3-464*z*t^3+432*y*z*t+8*z^2*t+6912*y*t^2-4320*z*t^2+13824*t^3+z^2-13824*z*t+55296*t^2-13824*z; 
*/



/*
// Cyclic 7: (done)
ring R = (0),(a, b, c, d, e),dp;
module I = a+b+c+d+e,a*b+a*e+b*c+c*d+d*e,a*b*c+a*b*e+a*d*e+b*c*d+c*d*e,a*b*c*d+a*b*c*e+a*b*d*e+a*c*d*e+b*c*d*e,a*b*c*d*e-1; 
*/

/*
// HCyclic 7: (running: 80% of 16GB mem while computing syz!... too slow... )
ring R = (0),(x1, x2, x3, x4, x5, x6, x7, w),dp;
module I = x1+x2+x3+x4+x5+x6+x7,x1*x2+x2*x3+x3*x4+x4*x5+x5*x6+x1*x7+x6*x7,x1*x2*x3+x2*x3*x4+x3*x4*x5+x4*x5*x6+x1*x2*x7+x1*x6*x7+x5*x6*x7,x1*x2*x3*x4+x2*x3*x4*x5+x3*x4*x5*x6+x1*x2*x3*x7+x1*x2*x6*x7+x1*x5*x6*x7+x4*x5*x6*x7,x1*x2*x3*x4*x5+x2*x3*x4*x5*x6+x1*x2*x3*x4*x7+x1*x2*x3*x6*x7+x1*x2*x5*x6*x7+x1*x4*x5*x6*x7+x3*x4*x5*x6*x7,x1*x2*x3*x4*x5*x6+x1*x2*x3*x4*x5*x7+x1*x2*x3*x4*x6*x7+x1*x2*x3*x5*x6*x7+x1*x2*x4*x5*x6*x7+x1*x3*x4*x5*x6*x7+x2*x3*x4*x5*x6*x7,x1*x2*x3*x4*x5*x6*x7-w^7; 
*/

// syz(I); print(_); // 

/*
module M;
M[5]=abcde*gen(1)-gen(1)+gen(6);
M[1]=a*gen(1)+b*gen(1)+c*gen(1)+d*gen(1)+e*gen(1)+gen(2);
M[2]=ab*gen(1)+bc*gen(1)+cd*gen(1)+ae*gen(1)+de*gen(1)+gen(3);
M[3]=abc*gen(1)+bcd*gen(1)+abe*gen(1)+ade*gen(1)+cde*gen(1)+gen(4);
M[4]=abcd*gen(1)+abce*gen(1)+abde*gen(1)+acde*gen(1)+bcde*gen(1)+gen(5);

ring P = (99981794),(a, b, c, d, e),dp;
setring R;

GBCandidate(M, P);$$
*/

int tG = timer;
def G = groebner(I);
int ttG = timer;

"Time GB: ", ttG - tG;


int tN = timer;
def L = Noro( I );
int ttN = timer;

"Time Noro: ", ttN - tN;

// All syzygies?
int tS = timer;
module S = syz(I);
int ttS = timer;

"Time Syz: ", ttS - tS;


module GB = L[1];
module SYZ = L[2];

":::::::::::::::::::::::::::::::::::::::::::: GB :::::::::::::::::::::::::::::::::::: ";
print(GB);

// test GB:
if( size(NF(I, GB)) > 0 )
{
  ERROR("NORO was wrong: I is bigger than GB!");
}

if( size(NF(GB, G)) > 0 )
{
  ERROR("NORO was wrong: GB is bigger than I!");
}

":::::::::::::::::::::::::::::::::::::::::::: SYZ :::::::::::::::::::::::::::::::::::: ";
print(SYZ);


if( size(SYZ) > 0 )
{
  // test syzygy
  if( size( module(transpose(SYZ)*transpose(I)) ) > 0 )
  {
    ERROR("NORO was wrong: SYZ are NOT syzygies of I!");
  }
}



// test SYZ:
if( size(NF(SYZ, groebner(S))) > 0 )
{
  ERROR("NORO was wrong: too much syzygies found!!!");
}

if( size(NF(S, groebner(SYZ))) > 0 )
{
  ERROR("NORO was wrong: too few syzygies found!!!");
}


":::::::::::::::::::::::::::::::::::::::::::: GOOD :::::::::::::::::::::::::::::::::::: ";


"Time GB: ", ttG - tG;
"Time Syz: ", ttS - tS;
"Time Noro: ", ttN - tN, " vs GB+SYZ: ", (ttS - tS) + (ttG - tG), "";

"Time NF: ", reduce_time;



number t = (ttN - tN);
number tt = ((ttS - tS) + (ttG - tG));
"Factor: ", t div tt;
"Factor: ", (ttN - tN) div ((ttS - tS) + (ttG - tG));


$$