///////////////////////////////////////////////////////////////////////////////
version="version modwalk.lib 4.0.0.0 Jun_2013 "; // $Id$
category = "Commutative Algebra";
info="
LIBRARY:  modwalk.lib      Groebner basis convertion

AUTHORS:  S. Oberfranz    oberfran@mathematik.uni-kl.de

OVERVIEW:

  A library for converting Groebner bases of an ideal in the polynomial
  ring over the rational numbers using modular methods. The procedures are
  inspired by the following paper:
  Elizabeth A. Arnold: Modular algorithms for computing Groebner bases.
  Journal of Symbolic Computation 35, 403-419 (2003).

PROCEDURES:
 modWalk(I);        standard basis conversion of I using modular methods (chinese remainder)
";

LIB "poly.lib";
LIB "ring.lib";
LIB "parallel.lib";
LIB "rwalk.lib";
LIB "grwalk.lib";
LIB "modstd.lib";


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

proc modpWalk(def II, int p, int variant, list #)
"USAGE:  modpWalk(I,p,#); I ideal, p integer, variant integer
ASSUME:  If size(#) > 0, then 
           #[1] is an intvec describing the current weight vector
           #[2] is an intvec describing the target weight vector
RETURN:  ideal - a standard basis of I mod p, integer - p
NOTE:    The procedure computes a standard basis of the ideal I modulo p and
         fetches the result to the basering. 
EXAMPLE: example modpWalk; shows an example
"
{
  option(redSB);
  int k,nvar@r;
  def R0 = basering;
  string ordstr_R0 = ordstr(R0);
  list rl = ringlist(R0);
  int sizerl = size(rl);
  int neg = 1 - attrib(R0,"global");
  if(typeof(II) == "ideal")
  {
    ideal I = II;
    int radius = 2;
    int pert_deg = 2;
  }
  if(typeof(II) == "list" && typeof(II[1]) == "ideal")
  {
    ideal I = II[1];
    if(size(II) == 2)
    {
      int radius = II[2];
      int pert_deg = 2;
    }
    if(size(II) == 3)
    {
      int radius = II[2];
      int pert_deg = II[3];
    }
  }
  rl[1] = p;
  int h = homog(I);
  def @r = ring(rl);
  setring @r;
  ideal i = fetch(R0,I);
  string order;
  if(system("nblocks") <= 2)
  {
    if(find(ordstr_R0, "M") + find(ordstr_R0, "lp") + find(ordstr_R0, "rp") <= 0)
    {
      order = "simple";
    }
  }

//-------------------------  make i homogeneous  -----------------------------

  if(!mixedTest() && !h)
  {
    if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg))
    {
      if(!((order == "simple") || (sizerl > 4)))
      {
        list rl@r = ringlist(@r);
        nvar@r = nvars(@r);      
        intvec w;
        for(k = 1; k <= nvar@r; k++)
        {
          w[k] = deg(var(k));
        }
        w[nvar@r + 1] = 1;
        rl@r[2][nvar@r + 1] = "homvar";
        rl@r[3][2][2] = w;
        def HomR = ring(rl@r);
        setring HomR;
        ideal i = imap(@r, i);
        i = homog(i, homvar);
      }
    }
  }

//-------------------------  compute a standard basis mod p  -----------------------------

  if(variant == 1)
  {
    if(size(#)>0)
    {
      i = rwalk(i,radius,pert_deg,#);
     // rwalk(i,radius,pert_deg,#); std(i);
    }
    else
    {
      i = rwalk(i,radius,pert_deg);
    }
  }
  if(variant == 2)
  {
    if(size(#) == 2)
    {
      i = gwalk(i,#);
    }
    else
    {
      i = gwalk(i);
    }
  }
  if(variant == 3)
  {
    if(size(#) == 2)
    {
      i = frandwalk(i,radius,#);
    }
    else
    {
      i = frandwalk(i,radius);
    }
  }
  if(variant == 4)
  {
    if(size(#) == 2)
    {
      trwalk(i,radius,pert_deg,#);
    }
    else
    {
      trwalk(i,radius,pert_deg);
    }
  }
  if(!mixedTest() && !h)
  {
    if(!((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg))
    {
      if(!((order == "simple") || (sizerl > 4)))
      {
        i = subst(i, homvar, 1);
        i = simplify(i, 34);
        setring @r;
        i = imap(HomR, i);
        i = interred(i);
        kill HomR;
      }
    }
  }
  setring R0;
  return(list(fetch(@r,i),p));
}
example
{
  "EXAMPLE:"; echo = 2;
  option(redSB);

  int p = 181;
  intvec a = 2,1,3,4;
  intvec b = 1,9,1,1;
  ring ra = 0,(w,x,y,z),(a(a),lp);
  ideal I = std(cyclic(4));
  ring rb = 0,(w,x,y,z),(a(b),lp);
  ideal I = imap(ra,I);
  modpWalk(I,p,1,a,b);
  std(I);
/*  ~
  list P = modpWalk(I,p,2,a,b);
  P;
  list P1 = modpWalk(I,p,1,a,b);
  P1;
  ring R =181,(w,x,y,z),(a(b),lp);
  std(cyclic(4));
  ~
  a = 1,3,5,7,9;
  b = 5,2,7,2,6;
  ring raa = 0,(v,w,x,y,z),(a(a),lp);
  ideal I = std(cyclic(5));
  ring rbb =0,(v,w,x,y,z),(a(b),lp);
  ideal I = imap(raa,I);
  int q = 32003;
  list Q = modpWalk(I,q,1,a,b);
  ideal J = rwalk(I,2,3,a,b);
  J;
  std(I);
  Q;
  ring r = 32003,(v,w,x,y,z),(a(b),lp);
  std(cyclic(5));*/
}

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

proc modWalk(def II, int variant, list #)
"USAGE:  modWalk(II); II ideal or list(ideal,int)
ASSUME:  If variant = 1 the random walk algorithm with radius II[2] is applied
         to II[1] if II = list(ideal, int). It is applied to II with radius 2
         if II is an ideal. If variant = 2, the Groebner walk algorithm is
         applied to II[1] or to II, respectively.
         If size(#) > 0, then # contains either 1, 2 or 4 integers such that
@*       - #[1] is the number of available processors for the computation,
@*       - #[2] is an optional parameter for the exactness of the computation,
                if #[2] = 1, the procedure computes a standard basis for sure,
@*       - #[3] is the number of primes until the first lifting,
@*       - #[4] is the constant number of primes between two liftings until
           the computation stops.
RETURN:  a standard basis of I if no warning appears.
NOTE:    The procedure converts a standard basis of I (over the rational
         numbers) from the ordering \"a(v),lp\", "dp\" or \"Dp\" to the ordering
         \"(a(w),lp\" or \"a(1,0,...,0),lp\" by using modular methods.
         By default the procedure computes a standard basis of I for sure, but
         if the optional parameter #[2] = 0, it computes a standard basis of I
         with high probability.
EXAMPLE: example modWalk; shows an example
"
{
  int TT = timer;
  int RT = rtimer;
  int i,j,pTest,sizeTest,weighted,n1;
  bigint N;

  def R0 = basering;
  list rl = ringlist(R0);
  if((npars(R0) > 0) || (rl[1] > 0))
  {
    ERROR("Characteristic of basering should be zero, basering should have no parameters.");
  }

  if(typeof(II) == "ideal")
  {
    ideal I = II;
    kill II;
    list II;
    II[1] = I;
    II[2] = 2;
    II[3] = 2;
  }
  else
  {
    if(typeof(II) == "list" && typeof(II[1]) == "ideal")
    {
      ideal I = II[1];
      if(size(II) == 1)
      {
        II[2] = 2;
        II[3] = 2;
      }
      if(size(II) == 2)
      {
        II[3] = 2;
      }
    
    }
    else
    {
      ERROR("Unexpected type of input.");
    }
  }

//--------------------  Initialize optional parameters  ------------------------
  n1 = system("cpu");
  if(size(#) == 0)
  {
    int exactness = 1;
    int n2 = 10;
    int n3 = 10;
  }
  else
  {
    if(size(#) == 1)
    {
      if(typeof(#[1]) == "int")
      {
        if(#[1] < n1)
        {
          n1 = #[1];
        }
        int exactness = 1;
        if(n1 >= 10)
        {
          int n2 = n1 + 1;
          int n3 = n1;
        }
        else
        {
          int n2 = 10;
          int n3 = 10;
        }
      }
      else
      {
        ERROR("Unexpected type of input.");
      }
    }
    if(size(#) == 2)
    {
      if(typeof(#[1]) == "int" && typeof(#[2]) == "int")
      {
        if(#[1] < n1)
        {
          n1 = #[1];
        }
        int exactness = #[2];
        if(n1 >= 10)
        {
          int n2 = n1 + 1;
          int n3 = n1;
        }
        else
        {
          int n2 = 10;
          int n3 = 10;
        }
      }
      else
      {
        if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec")
        {
          intvec curr_weight = #[1];
          intvec target_weight = #[2];
          weighted = 1;
          int n2 = 10;
          int n3 = 10;
          int exactness = 1;
        }
        else
        {
          ERROR("Unexpected type of input.");
        }
      }  
    }
    if(size(#) == 3)
    {
      if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int")
      {
        intvec curr_weight = #[1];
        intvec target_weight = #[2];
        weighted = 1;
        n1 = #[3];
        int n2 = 10;
        int n3 = 10;
        int exactness = 1;
      }
      else
      {
        ERROR("Unexpected type of input.");
      }
    }
    if(size(#) == 4)
    {
      if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int" && typeof(#[4]) == "int")
      {
        intvec curr_weight = #[1];
        intvec target_weight = #[2];
        weighted = 1;
        if(#[1] < n1)
        {
          n1 = #[3];
        }
        int exactness = #[4];
        if(n1 >= 10)
        {
          int n2 = n1 + 1;
          int n3 = n1;
        }
        else
        {
          int n2 = 10;
          int n3 = 10;
        }
      }
      else
      {
        if(typeof(#[1]) == "int" && typeof(#[2]) == "int" && typeof(#[3]) == "int" && typeof(#[4]) == "int")
        {
          if(#[1] < n1)
          {
            n1 = #[1];
          }
          int exactness = #[2];
          if(n1 >= #[3])
          {
            int n2 = n1 + 1;
          }
          else
          {
            int n2 = #[3];
          }
          if(n1 >= #[4])
          {
            int n3 = n1;
          }
          else
          {
            int n3 = #[4];
          }
        }
        else
        {
          ERROR("Unexpected type of input.");
        }
      }
    }
    if(size(#) == 6)
    {
      if(typeof(#[1]) == "intvec" && typeof(#[2]) == "intvec" && typeof(#[3]) == "int" && typeof(#[4]) == "int" && typeof(#[5]) == "int" && typeof(#[6]) == "int")
      {
        intvec curr_weight = #[1];
        intvec target_weight = #[2];
        weighted = 1;
        if(#[3] < n1)
        {
          n1 = #[3];
        }
        int exactness = #[4];
        if(n1 >= #[5])
        {
          int n2 = n1 + 1;
        }
        else
        {
          int n2 = #[5];
        }
        if(n1 >= #[6])
        {
          int n3 = n1;
        }
        else
        {
          int n3 = #[6];
        }
      }
      else
      {
        ERROR("Expected list(intvec,intvec,int,int,int,int) as optional parameter list.");
      }
    }
    if(size(#) == 1 || size(#) == 5 || size(#) > 6)
    {
      ERROR("Expected 0,2,3,4 or 5 optional arguments.");
    }
  }
  if(printlevel >= 10)
  {
  "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)+", exactness = "+string(exactness);
  }

//-------------------------  Save current options  -----------------------------
  intvec opt = option(get);
  option(redSB);

//--------------------  Initialize the list of primes  -------------------------
  int tt = timer;
  int rt = rtimer;
  int en = 2134567879;
  int an = 1000000000;
  intvec L = primeList(I,n2);
  if(n2 > 4)
  {
  //  L[5] = prime(random(an,en));
  }
  if(printlevel >= 10)
  {
    "CPU-time for primeList: "+string(timer-tt)+" seconds.";
    "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
  }
  int h = homog(I);
  list P,T1,T2,LL,Arguments,PP;
  ideal J,K,H;

//-------------------  parallelized Groebner Walk in positive characteristic  --------------------

  if(weighted)
  {
    for(i=1; i<=size(L); i++)
    {
      Arguments[i] = list(II,L[i],variant,list(curr_weight,target_weight));
    }
  }
  else
  {
    for(i=1; i<=size(L); i++)
    {
      Arguments[i] = list(II,L[i],variant);
    }
  }
  P = parallelWaitAll("modpWalk",Arguments);
  for(i=1; i<=size(P); i++)
  {
    T1[i] = P[i][1];
    T2[i] = bigint(P[i][2]);
  }

  while(1)
  {
    LL = deleteUnluckyPrimes(T1,T2,h);
    T1 = LL[1];
    T2 = LL[2];
//-------------------  Now all leading ideals are the same  --------------------
//-------------------  Lift results to basering via farey  ---------------------

    tt = timer; rt = rtimer;
    N = T2[1];
    for(i=2; i<=size(T2); i++)
    {
      N = N*T2[i];
    }
    H = chinrem(T1,T2);
    //J = parallelFarey(H,N,n1);
    J=farey(H,N);
    if(printlevel >= 10)
    {
      "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
      "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
    }

//----------------  Test if we already have a standard basis of I --------------

    tt = timer; rt = rtimer;
    //pTest = pTestSB(I,J,L,variant);
    pTest = primeTestSB(I,J,L,variant);
    if(printlevel >= 10)
    {
      "CPU-time for pTest is "+string(timer - tt)+" seconds.";
      "Real-time for pTest is "+string(rtimer - rt)+" seconds.";
    }
    if(pTest)
    {
      if(printlevel >= 10)
      {
        "CPU-time for computation without final tests is "+string(timer - TT)+" seconds.";
        "Real-time for computation without final tests is "+string(rtimer - RT)+" seconds.";
      }
      attrib(J,"isSB",1);
      if(exactness == 0)
      {
        option(set, opt);
        return(J);
      }
      else
      {
        tt = timer;
        rt = rtimer;
        sizeTest = 1 - isIdealIncluded(I,J,n1);
        if(printlevel >= 10)
        {
          "CPU-time for checking if I subset <G> is "+string(timer - tt)+" seconds.";
          "Real-time for checking if I subset <G> is "+string(rtimer - rt)+" seconds.";
        }
        if(sizeTest == 0)
        {
          tt = timer;
          rt = rtimer;
          K = std(J);
          if(printlevel >= 10)
          {
            "CPU-time for last std-computation is "+string(timer - tt)+" seconds.";
            "Real-time for last std-computation is "+string(rtimer - rt)+" seconds.";
          }
          if(size(reduce(K,J)) == 0)
          {
            option(set, opt);
            return(J);
          }
        }
      }
    }
//--------------  We do not already have a standard basis of I, therefore do the main computation for more primes  --------------

    T1 = H;
    T2 = N;
    j = size(L)+1;
    tt = timer; rt = rtimer;
    L = primeList(I,n3,L,n1);
L;
    if(printlevel >= 10)
    {
      "CPU-time for primeList: "+string(timer-tt)+" seconds.";
      "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
    }
    Arguments = list();
    PP = list();
    if(weighted)
    {
      for(i=j; i<=size(L); i++)
      {
        //Arguments[i-j+1] = list(II,L[i],variant,list(curr_weight,target_weight));
        Arguments[size(Arguments)+1] = list(II,L[i],variant,list(curr_weight,target_weight));
      }
    }
    else
    {
      for(i=j; i<=size(L); i++)
      {
        //Arguments[i-j+1] = list(II,L[i],variant);
        Arguments[size(Arguments)+1] = list(II,L[i],variant);
      }
    }
    PP = parallelWaitAll("modpWalk",Arguments);
    if(printlevel >= 10)
    {
      "parallel modpWalk";
     // ~
    }
    for(i=1; i<=size(PP); i++)
    {
      //P[size(P) + 1] = PP[i];
      T1[size(T1) + 1] = PP[i][1];
      T2[size(T2) + 1] = bigint(PP[i][2]);
    }
  }
  if(printlevel >= 10)
  {
    "CPU-time for computation with final tests is "+string(timer - TT)+" seconds.";
    "Real-time for computation with final tests is "+string(rtimer - RT)+" seconds.";
  }
}

example
{
  "EXAMPLE:";
  echo = 2;
  ring R=0,(x,y,z),lp;
  ideal I=-x+y2z-z,xz+1,x2+y2-1;
  // I is a standard basis in dp
  modWalk(I,1);
  modWalk(I,2,2,0);
  modWalk(I,3,system("cpu"),0);
  std(I);

  ring r0=0,x(1..6),dp;
  ideal i0=std(cyclic(6));
  ring r=0,x(1..6),lp;
  ideal i=fetch(r0,i0);
  modWalk(i,1,system("cpu"),0);
  modWalk(i,3);
}

////////////////////////////////////////////////////////////////////////////////
proc isIdealIncluded(ideal I, ideal J, int n1)
"USAGE:  isIdealIncluded(I,J,int n1); I ideal, J ideal, n1 integer
"
{
  if(n1 > 1)
  {
    int k;
    list args,results;
    for(k=1; k<=size(I); k++)
    {
      args[k] = list(ideal(I[k]),J,1);
    }
    results = parallelWaitAll("reduce",args);
    for(k=1; k<=size(results); k++)
    {
      if(results[k] == 0)
      {
        return(1);
      }
    }
    return(0);
  }
  else
  {
    if(reduce(I,J,1) == 0)
    {
      return(1);
    }
    else
    {
      return(0);
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
proc parallelChinrem(list T1, list T2, int n1)
"USAGE:  parallelChinrem(T1,T2); T1 list of ideals, T2 list of primes, n1 integer"
{
  int i,j,k;

  ideal H,J;

  list arguments_chinrem,results_chinrem;
  for(i=1; i<=size(T1); i++)
  {
    J = ideal(T1[i]);
    attrib(J,"isSB",1);
    arguments_chinrem[size(arguments_chinrem)+1] = list(list(J),T2);
  }
  results_chinrem = parallelWaitAll("chinrem",arguments_chinrem);
    for(j=1; j <= size(results_chinrem); j++)  
    {
      J = results_chinrem[j];
      attrib(J,"isSB",1);
      if(isIdealIncluded(J,H,n1) == 0)
      {
        if(H == 0)
        {
          H = J;
        }
        else
        {
          H = H,J;
        }
      }
    }
  return(H);
}

////////////////////////////////////////////////////////////////////////////////
proc parallelFarey(ideal H, bigint N, int n1)
"USAGE:  parallelFarey(H,N,n1); H ideal, N bigint, n1 integer
"
{
  int i,j;
  int ii = 1;
  list arguments_farey,results_farey;
  for(i=1; i<=size(H); i++)
  {
    for(j=1; j<=size(H[i]); j++)
    {
      arguments_farey[size(arguments_farey)+1] = list(H[i][j],N);
    }
  }
  results_farey = parallelWaitAll("farey",arguments_farey);
  ideal J,K;
  poly f_farey;
  while(ii<=size(results_farey))
  {
    for(i=1; i<=size(H); i++)
    {
      f_farey = 0;
      for(j=1; j<=size(H[i]); j++)
      {
        f_farey = f_farey + results_farey[ii][1];
        ii++;
      }
      K = ideal(f_farey);
      attrib(K,"isSB",1);
      attrib(J,"isSB",1);
      if(isIdealIncluded(K,J,n1) == 0)
      {
        if(J == 0)
        {
          J = K;
        }
        else
        {
          J = J,K;
        }
      }
    }
  }
  return(J);
}

proc primeTestSB(def II, ideal J, list L, int variant, list #)
"USAGE:  primeTestSB(I,J,L,variant,#); I,J ideals, L intvec of primes, variant int
RETURN:  1 (resp. 0) if for a randomly chosen prime p that is not in L
         J mod p is (resp. is not) a standard basis of I mod p
EXAMPLE: example primeTestSB; shows an example
"
{
if(typeof(II) == "ideal")
  {
  ideal I = II;
  int radius = 2;
  }
if(typeof(II) == "list")
  {
  ideal I = II[1];
  int radius = II[2];
  }

int i,j,k,p;
def R = basering;
list r = ringlist(R);

while(!j)
  {
  j = 1;
  p = prime(random(1000000000,2134567879));
  for(i = 1; i <= size(L); i++)
    {
    if(p == L[i])
      {
      j = 0;
      break;
      }
    }
  if(j)
    {
    for(i = 1; i <= ncols(I); i++)
      {
      for(k = 2; k <= size(I[i]); k++)
        {
        if((denominator(leadcoef(I[i][k])) mod p) == 0)
          {
          j = 0;
          break;
          }
        }
      if(!j)
        {
        break;
        }
      }
    }
  if(j)
    {
    if(!primeTest(I,p))
      {
      j = 0;
      }
    }
  }
r[1] = p;
def @R = ring(r);
setring @R;
ideal I = imap(R,I);
ideal J = imap(R,J);
attrib(J,"isSB",1);

int t = timer;
j = 1;
if(isIncluded(I,J) == 0)
  {
  j = 0;
  }
if(printlevel >= 11)
  {
  "isIncluded(I,J) takes "+string(timer - t)+" seconds";
  "j = "+string(j);
  }
t = timer;
if(j)
  {
  if(size(#) > 0)
    {
    ideal K = modpWalk(I,p,variant,#)[1];
    }
  else
    {
    ideal K = modpWalk(I,p,variant)[1];
    }
  t = timer;
  if(isIncluded(J,K) == 0)
    {
    j = 0;
    }
  if(printlevel >= 11)
    {
    "isIncluded(K,J) takes "+string(timer - t)+" seconds";
    "j = "+string(j);
    }
  }
setring R;

return(j);
}
example
{ "EXAMPLE:"; echo = 2;
   intvec L = 2,3,5;
   ring r = 0,(x,y,z),lp;
   ideal I = x+1,x+y+1;
   ideal J = x+1,y;
   primeTestSB(I,I,L,1);
   primeTestSB(I,J,L,1);
}