source: git/Singular/LIB/modwalk.lib @ 4083fa

///////////////////////////////////////////////////////////////////////////////
version="version modwalk.lib 4.0.0.0 Jun_2013 ";
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:
modpWalk(ideal,int,int,int[,int,int,int,int])    standard basis conversion of I in prime characteristic
modWalk(ideal,int,int[,int,int,int,int]);        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, int reduction, 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";
    }
  }
<<<<<<< HEAD
=======

//-------------------------  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);
      }
    }
  }

>>>>>>> f533f6f7667328bccb271b19b2f603aaebe41596
//-------------------------  compute a standard basis mod p  -----------------------------
  if(variant == 1)
  {
    if(size(#)>0)
    {
      i = rwalk(i,radius,pert_deg,reduction,#);
    }
    else
    {
      i = rwalk(i,radius,pert_deg,reduction);
    }
  }
  if(variant == 2)
  {
    if(size(#) == 2)
    {
      i = gwalk(i,reduction,#);
    }
    else
    {
      i = gwalk(i,reduction);
    }
  }
  if(variant == 3)
  {
    if(size(#) == 2)
    {
      i = frandwalk(i,radius,reduction,#);
    }
    else
    {
      i = frandwalk(i,radius,reduction);
    }
  }
  if(variant == 4)
  {
    if(size(#) == 2)
    {
      i=fwalk(i,#);
    }
    else
    {
      i=fwalk(i);
    }
  }
  if(variant == 5)
  {
    if(size(#) == 2)
    {
     i=prwalk(i,radius,pert_deg,pert_deg,reduction,#);
    }
    else
    {
      i=prwalk(i,radius,pert_deg,pert_deg,reduction);
    }
  }
  if(variant == 6)
  {
    if(size(#) == 2)
    {
      i=pwalk(i,pert_deg,pert_deg,reduction,#);
    }
    else
    {
<<<<<<< HEAD
      i=pwalk(i,pert_deg,pert_deg,reduction);
=======
      i=pwalk(i,pert_deg,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;
      }
>>>>>>> f533f6f7667328bccb271b19b2f603aaebe41596
    }
  }

  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,x(1..4),(a(a),lp);
  ideal I = std(cyclic(4));
  int reduction = 1;
  ring rb = 0,x(1..4),(a(b),lp);
  ideal I = imap(ra,I);
  modpWalk(I,p,1,reduction,a,b);
  std(I);
}

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

proc modWalk(def II, int variant, int reduction, 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.
@*       - 2, the Groebner Walk algorithm is applied to II[1] or to II, respectively.
@*       - 3, the Fractal Walk algorithm with random element is applied to II[1] or II,
           respectively.
@*       - 4, the Fractal Walk algorithm is applied to II[1] or II, respectively.
@*       - 5, the Perturbation Walk algorithm with random element is applied to II[1]
           or II, respectively, with radius II[3] and perturbation degree II[2].
@*       - 6, the Perturbation Walk algorithm is applied to II[1] or II, respectively,
           with perturbation degree II[3].
         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("--cpus");
  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,reduction,list(curr_weight,target_weight));
    }
  }
  else
  {
    for(i=1; i<=size(L); i++)
    {
      Arguments[i] = list(II,L[i],variant,reduction);
    }
  }
  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 = primeTest(J, prime(random(1000000000,2134567879)));
    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);
    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[size(Arguments)+1] = list(II,L[i],variant,reduction,list(curr_weight,target_weight));
      }
    }
    else
    {
      for(i=j; i<=size(L); i++)
      {
        Arguments[size(Arguments)+1] = list(II,L[i],variant,reduction);
      }
    }
    PP = parallelWaitAll("modpWalk",Arguments);
    if(printlevel >= 10)
    {
      "parallel modpWalk";
    }
    for(i=1; i<=size(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= y3+xyz+y2z+xz3, 3+xy+x2y+y2z;
  int reduction = 0;
  ideal J = modWalk(I,1,1);
  J;
}

////////////////////////////////////////////////////////////////////////////////
static 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);
    }
  }
}

////////////////////////////////////////////////////////////////////////////////
static 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);
}

////////////////////////////////////////////////////////////////////////////////
static 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);
}

////////////////////////////////////////////////////////////////////////////////
static proc mixedTest()
"USAGE:  mixedTest();
RETURN:  1 if ordering of basering is mixed, 0 else
EXAMPLE: example mixedTest(); shows an example
"
{
   int i,p,m;
   for(i = 1; i <= nvars(basering); i++)
   {
      if(var(i) > 1)
      {
         p++;
      }
      else
      {
         m++;
      }
   }
   if((p > 0) && (m > 0)) { return(1); }
   return(0);
}
example
{ "EXAMPLE:"; echo = 2;
   ring R1 = 0,(x,y,z),dp;
   mixedTest();
   ring R2 = 31,(x(1..4),y(1..3)),(ds(4),lp(3));
   mixedTest();
   ring R3 = 181,x(1..9),(dp(5),lp(4));
   mixedTest();
}