// $Id: standard.lib,v 1.20 1998-06-18 15:45:12 siebert Exp $
//////////////////////////////////////////////////////////////////////////////

version="$Id: standard.lib,v 1.20 1998-06-18 15:45:12 siebert Exp $";
info="
LIBRARY: standard.lib   PROCEDURES WHICH ARE ALWAYS LOADED AT START-UP

 stdfglm(ideal[,ord])   standard basis of the ideal via fglm [and ordering ord]
 stdhilb(ideal)         standard basis of the ideal using the Hilbert function
 groebner(ideal/module) standard basis of ideal or module using a
                        heuristically choosen method
 quotient(any,any[,n])  a general quotient procedure calling several algorithms
                  allows module/module, ideal/ideal, module/ideal and a
                  pre-definition of the algorithm by the parameter n
 quotient1(m1,m2) computes quotients by every vector of m2 and intersects them
 quotient2(m1,m2) a heuristic variant: the quotient is just defined by a
                  (not really) general element of m2 which has to be proved
 quotient3(m1,m2) the homogeneous variant of quotient5(m1,m2)
 quotient4(m1,m2) the same as quotient5(m1,m2) using the modulo-command
                  instead of the quotient-command from the kernel
 quotient5(m1,m2) computes with a real general element of m2 by adjoining
                  a new variable
";

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

proc stdfglm (ideal i, list #)
"USAGE:   stdfglm(i[,s]); i ideal, s string (any allowed ordstr of a ring)
RETURN:  stdfglm(i): standard basis of i in the basering, calculated via fglm
                     from ordering \"dp\" to the ordering of the basering.
         stdfglm(i,s): standard basis of i in the basering, calculated via
                     fglm from ordering s to the ordering of the basering.
EXAMPLE: example stdfglm; shows an example"
{
   string os;
   def dr= basering;
   if( (size(#)==0) or (typeof(#[1]) != "string") )
   {
     os = "dp(" + string( nvars(dr) ) + ")";
     if ( (find( ordstr(dr), os ) != 0) and (find( ordstr(dr), "a") == 0) )
     {
       os= "Dp";
     }
     else
     {
       os= "dp";
     }
   }
   else { os = #[1]; }
   execute "ring sr=("+charstr(dr)+"),("+varstr(dr)+"),"+os+";";
   ideal i= fetch(dr,i);
   intvec opt= option(get);
   option(redSB);
   i=std(i);
   option(set,opt);
   setring dr;
   return (fglm(sr,i));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r  = 0,(x,y,z),lp;
   ideal i = y3+x2, x2y+x2, x3-x2, z4-x2-y;
   ideal i1= stdfglm(i);         //uses fglm from "dp" to "lp"
   i1;
   ideal i2= stdfglm(i,"Dp");    //uses fglm from "Dp" to "lp"
   i2;
}
/////////////////////////////////////////////////////////////////////////////

proc stdhilb(ideal i,list #)
"USAGE:   stdhilb(i);  i ideal
         stdhilb(i,v); i homogeneous ideal, v intvec (the Hilbert function)
RETURN:  stdhilb(i): a standard basis of i (computing v internally)
         stdhilb(i,v): standard basis of i, using the given Hilbert function
EXAMPLE: example stdhilb; shows an example"
{
   def R=basering;

   if((homog(i)==1)||(ordstr(basering)[1]=="d"))
   {
      if ((size(#)!=0)&&(homog(i)==1))
      {
         return(std(i,#[1]));
      }
      return(std(i));
   }

   execute "ring S = ("+charstr(R)+"),("+varstr(R)+",@t),dp;";
   ideal i=homog(imap(R,i),@t);
   intvec v=hilb(std(i),1);
   execute "ring T = ("+charstr(R)+"),("+varstr(R)+",@t),("+ordstr(R)+");";
   ideal i=fetch(S,i);
   ideal a=std(i,v);
   setring R;
   map phi=T,maxideal(1),1;
   ideal a=phi(a);

   int k,j;
   poly m;
   int c=size(i);

   for(j=1;j<c;j++)
   {
     if(deg(a[j])==0)
     {
       a=ideal(1);
       attrib(a,"isSB",1);
       return(a);
     }
     if(deg(a[j])>0)
     {
       m=lead(a[j]);
       for(k=j+1;k<=c;k++)
       {
          if(size(lead(a[k])/m)>0)
          {
            a[k]=0;
          }
       }
     }
   }
   a=simplify(a,2);
   attrib(a,"isSB",1);
   return(a);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r = 0,(x,y,z),lp;
   ideal i = y3+x2, x2y+x2, x3-x2, z4-x2-y;
   ideal i1= stdhilb(i); i1;
   // is in this case equivalent to:
   intvec v=1,0,0,-3,0,1,0,3,-1,-1;
   ideal i2=stdhilb(i,v);
}
//////////////////////////////////////////////////////////////////////////

proc groebner(def i, list #)
"USAGE: groebner(i[, wait]) i -- ideal/module; wait -- int
RETURNS: Standard basis of ideal or module which is computed using a
         heuristically choosen method:
         If the ordering of the current ring is a local ordering, or
         if it is a non-block ordering and the current ring has no
         parameters, then std(i) is returned.
         Otherwise, i is mapped into a ring with no parameters and
         ordering dp, where its Hilbert series is computed. This is
         followed by a Hilbert-series based std computation in the
         original ring.
NOTE: If a 2nd argument 'wait' is given, then the computation proceeds
      at most 'wait' seconds. That is, if no result could be computed in
      'wait' seconds, then the computation is interrupted, 0 is returned,
      a warning message is displayed, and the global variable
      'groebner_error' is defined.
EXAMPLE: example groebner; shows an example"
{
  def P=basering;

  // we have two arguments -- try to use MPfork links
  if (size(#) > 0)
  {
    if (system("with", "MP"))
    {
      if (typeof(#[1]) == "int")
      {
        int wait = #[1];
        int j = 10;

        string bs = nameof(basering);
        link l_fork = "MPtcp:fork";
        open(l_fork);
        write(l_fork, quote(system("pid")));
        int pid = read(l_fork);
        write(l_fork, quote(groebner(eval(i))));

        // sleep in small intervalls for appr. one second
        if (wait > 0)
        {
          while(j < 1000000)
          {
            if (status(l_fork, "read", "ready", j)) {break;}
            j = j + j;
          }
        }

        // sleep in intervalls of one second from now on
        j = 1;
        while (j < wait)
        {
          if (status(l_fork, "read", "ready", 1000000)) {break;}
          j = j + 1;
        }

        if (status(l_fork, "read", "ready"))
        {
          def result = read(l_fork);
          if (bs != nameof(basering))
          {
            def PP = basering;
            setring P;
            def result = imap(PP, result);
            kill PP;
          }
          if (defined(groebner_error))
          {
            kill(groebner_error);
          }
          kill (l_fork);
        }
        else
        {
          ideal result;
          if (! defined(groebner_error))
          {
            int groebner_error = 1;
            export groebner_error;
          }
          "// ** groebner did not finish";
          j = system("sh", "kill " + string(pid));
        }
        return (result);
      }
      else
      {
        "// ** groebner needs int as 2nd arg";
      }
    }
    else
    {
      "// ** groebner with two args is not supported in this configuration";
    }
  }

  // we are still here -- do the actual computation
  string ordstr_P = ordstr(P);
  if (find(ordstr_P,"s") > 0)
  {
    //spaeter den lokalen fall ueber lp oder aehnlich behandeln
    return(std(i));
  }

  int IsSimple_P;
  if (system("nblocks") <= 2)
  {
    if (find(ordstr_P, "M") <= 0)
    {
      IsSimple_P = 1;
    }
  }
  int npars_P = npars(P);

  // return std if no parameters and (dp or wp)
  if ((npars_P <= 1) && IsSimple_P)
  {
    if (find(ordstr_P, "d") > 0)
    {
      return (std(i));
    }
    if (find(ordstr_P,"w") > 0)
    {
      return (std(i));
    }
  }

  // reset options
  intvec opt=option(get);
  int p_opt;
  string s_opt = option();
  option(none);
  // turn on option(prot) and/or option(mem), if previously set
  if (find(s_opt, "prot"))
  {
    option(prot);
    p_opt = 1;
  }
  if (find(s_opt, "mem"))
  {
    option(mem);
  }

  // construct ring in which first std computation is done
  string varstr_P = varstr(P);
  string parstr_P = parstr(P);
  int is_homog = (homog(i) && (npars_P <= 1));
  int add_vars = 0;
  string ri = "ring Phelp =";
  
  // more than one parameters are converted to ring variables
  if (npars_P > 1)
  {
    ri = ri + string(char(P)) + ",(" + varstr_P + "," + parstr_P;
    add_vars = npars_P;
  }
  else
  {
    ri = ri + "(" + charstr(P) + "),(" + varstr_P;
  }

  // a homogenizing variable is added, if necessary
  if (! is_homog)
  {
    ri = ri + ",@t";
    add_vars = add_vars + 1;
  }
  // ordering is set to (dp, C)
  ri = ri + "),(dp,C);";

  // change the ring
  execute(ri);

  // get ideal from previous ring
  if (is_homog)
  {
    ideal qh = imap(P, i);
  }
  else
  {
    // and homogenize
    ideal qh=homog(imap(P,i),@t);
  }

  // compute std and hilbert series
  if (p_opt)
  {
    "std in " + ri[13, size(ri) - 13];
  }
  ideal qh1=std(qh);
  intvec hi=hilb(qh1,1);

  if (add_vars == 0)
  {
    // no additional variables were introduced
    setring P; // can immediately change to original ring
    // simply compute std with hilbert series in original ring
    if (p_opt)
    {
      "std with hilb in basering";
      i = std(i, hi);
    }
  }
  else
  {
    // additional variables were introduced
    // need another intermediate ring
    ri = "ring Phelp1 = (" + charstr(Phelp) 
      + "),(" + varstr(Phelp) + "),(" + ordstr_P;

    // for lp wit at most one parameter, we do not need a block ordering
    if ( ! (IsSimple_P && (add_vars <2) && find(ordstr_P, "l")))
    {
      // need block ordering
      ri = ri + ", dp(" + string(add_vars) + ")";
    }
    ri = ri + ");";

    // change to intermediate ring
    execute(ri);
    ideal qh = imap(Phelp, qh);
    kill Phelp;
    if (p_opt)
    {
      "std with hilb in " + ri[14,size(ri)-14];
    }
    // compute std with Hilbert series
    qh = std(qh, hi);
    // subst 1 for homogenizing var
    if (!is_homog)
    {
      qh = subst(qh, @t, 1);
    }

    // go back to original ring
    setring P;
    // get ideal, delete zeros and clean SB
    i = imap(Phelp1,qh);
    i = simplify(i, 34);
    kill Phelp1;
  }

  // clean-up time
  option(set, opt);
  if (find(s_opt, "redSB") > 0)
  {
    i=interred(i);
  }
  attrib(i, "isSB", 1);
  return (i);
}
example
{
  "EXAMPLE: "; echo = 2;
  ring r = 0, (a,b,c,d), lp;
  option(prot);
  ideal i = a+b+c+d, ab+ad+bc+cd, abc+abd+acd+bcd, abcd-1; // cyclic 4
  groebner(i);
  ring rp = (0, a, b), (c,d), lp;
  ideal i = imap(r, i);
  ideal j = groebner(i);
  option(noprot);
  j; simplify(j, 1); std(i);
  if (system("with", "MP")) {groebner(i, 0);}
  defined(groebner_error);
}


//////////////////////////////////////////////////////////////////////////
proc res(list #)
{
   def P=basering;
   def m=#[1]; //the ideal or module

   int i=#[2]; //the length of the resolution
               //if size(#)>2 a minimal resolution is computed

   //LaScala for the homogeneous case
   if(homog(m)==1)
   {
      resolution re=lres(m,i);
      if(size(#)>2)
      {
         re=minres(re);
      }
      return(re);
   }

   //mres for the global non homogeneous case
   if(find(ordstr(P),"s")==0)
   {
      string ri= "ring Phelp ="
                  +string(char(P))+",("+varstr_P+"),(dp,C);";
      execute(ri);
      def m=imap(P,m);
      list re=mres(m,i);
      setring P;
      resolution result=imap(Phelp,re);
      return(result);
   }

   //sres for the local case and not minimal resolution
   if(size(#)<=2)
   {
      string ri= "ring Phelp ="
                  +string(char(P))+",("+varstr_P+"),(ls,c);";
      execute(ri);
      def m=imap(P,m);
      m=std(m);
      list re=sres(m,i);
      setring P;
      resolution result=imap(Phelp,re);
      return(result);
   }

   //mres for the local case and minimal resolution
   string ri= "ring Phelp ="
                  +string(char(P))+",("+varstr_P+"),(ls,C);";
   execute(ri);
   def m=imap(P,m);
   list re=mres(m,i);
   setring P;
   resolution result=imap(Phelp,re);
   return(result);
}

proc quotient (any m1,any m2,list L)
USAGE:   quotient(m1, m2[, n]); m1, m2 two submodules of k^s,
         n (optional) integer (1<= n <=5)
RETURN:  the quotient of m1 and m2
EXAMPLE: example quot; shows an example
{
  if (((typeof(m1)!="ideal") and (typeof(m1)!="module"))
     or ((typeof(m2)!="ideal") and (typeof(m2)!="module")))
  {
    "USAGE:   quot(m1, m2[, n]); m1, m2 two submodules of k^s,";
    "         n (optional) integer (1<= n <=5)";
    "RETURN:  the quotient of m1 and m2";
    "EXAMPLE: example quot; shows an example";
    return();
  }
  if (typeof(m1)!=typeof(m2))
  {
    return quot(m1,m2);
  }
  if (size(L)>0)
  {
    if (typeof(L[1])=="int" )
    {
      return quot1(m1,m2,L[1]);
    }
  }
  else
  {
    return quot1(m1,m2,2);
  }
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=181,(x,y,z),(c,ls);
  ideal id1=maxideal(4);
  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
  option(prot);
  ideal id6=quot(id1,id2);
  id6;
  ideal id7=quotient(id1,id2,1);
  id7;
  ideal id8=quotient(id1,id2,2);
  id8;
}

static proc quot1 (module m1, module m2,int n)
USAGE:   quot1(m1, m2, n); m1, m2 two submodules of k^s, 
         n integer (1<= n <=5)
RETURN:  the quotient of m1 and m2
EXAMPLE: example quot; shows an example
{
  if (n==1)
  {
    return(quotient1(m1,m2));
  }
  else 
  {   
    if (n==2)
    {
      return(quotient2(m1,m2));
    }
    else 
    { 
      if (n==3)
      {
        return(quotient3(m1,m2));
      }
      else 
      { 
        if (n==4)
        {
          return(quotient4(m1,m2));
        }
        else 
        { 
          if (n==5)
          {
            return(quotient5(m1,m2));
          }
          else
          {
            return(quotient(m1,m2));
          }
        }
      }
    }
  } 
}
example
{ "EXAMPLE:"; echo = 2;
  ring r=181,(x,y,z),(c,ls);
  ideal id1=maxideal(4);
  ideal id2=x2+xyz,y2-z3y,z3+y5xz;
  option(prot);
  ideal id6=quot(id1,id2);
  id6;
  ideal id7=quot1(id1,id2,1);
  id7;
  ideal id8=quot1(id1,id2,2);
  id8;
}

static proc quotient0(module a,module b) 
{
  module mm=b+a;
  resolution rs=system("LaScala",mm);
  list I=list(rs);
  matrix M=I[2];
  matrix A[1][nrows(M)]=M[1..nrows(M),1];
  ideal i=A;
  return (i);
}
proc quotient1(module a,module b)  //17sec
USAGE:   quotient1(m1, m2); m1, m2 two submodules of k^s, 
RETURN:  the quotient of m1 and m2
{
  int i;
  a=std(a);
  module dummy;
  module B=NF(b,a)+dummy;
  ideal re=quot(a,module(B[1]));
  for(i=2;i<=size(B);i++)
  {
     re=intersect1(re,quot(a,module(B[i])));
  }
  return(re);   
}
proc quotient2(module a,module b)    //13sec
USAGE:   quotient2(m1, m2); m1, m2 two submodules of k^s, 
RETURN:  the quotient of m1 and m2
{
  a=std(a);
  module dummy;
  module bb=NF(b,a)+dummy;
  int i=size(bb);
  ideal re=(quot(a,module(bb[i])));
  bb[i]=0;
  module temp;
  module temp1;
  module bbb;
  int mx;
  i=i-1;
  while (1)
  {
    if (i==0) break;
    temp = a+bb*re;
    temp1 = lead(interred(temp));
    mx=ncols(a);
    if (ncols(temp1)>ncols(a))
    {
      mx=ncols(temp1);
    }
    temp1 = matrix(temp1,1,mx)-matrix(lead(a),1,mx);
    temp1 = dummy+temp1;
    if (deg(temp1[1])<0) break;
    re=(intersect1(re,quot(a,module(bb[i]))));
    bb[i]=0;
    i = i-1;
  }
  return(re);   
}
proc quotient3(module a,module b)   //89sec
USAGE:   quotient3(m1, m2); m1, m2 two submodules of k^s, 
         only for global rings
RETURN:  the quotient of m1 and m2
{
  string s="ring @newr=("+charstr(basering)+
           "),("+varstr(basering)+",@t,@w),dp;";
  def @newP=basering;
  execute s;
  module b=imap(@newP,b);
  module a=imap(@newP,a);
  int i;
  int j=size(b);
  vector @b;
  for(i=1;i<=j;i++)
  {
     @b=@b+@t^(i-1)*@w^(j-i+1)*b[i];
  }
  ideal re=quot(a,module(@b));
  setring @newP;
  ideal re=imap(@newr,re);
  return(re);   
}
proc quotient5(module a,module b)   //89sec
USAGE:   quotient5(m1, m2); m1, m2 two submodules of k^s, 
         only for global rings
RETURN:  the quotient of m1 and m2
{
  string s="ring @newr=("+charstr(basering)+
           "),("+varstr(basering)+",@t),dp;";
  def @newP=basering;
  execute s;
  module b=imap(@newP,b);
  module a=imap(@newP,a);
  int i;
  int j=size(b);
  vector @b;
  for(i=1;i<=j;i++)
  {
     @b=@b+@t^(i-1)*b[i];
  }
  @b=homog(@b,@w);
  ideal re=quot(a,module(@b));
  setring @newP;
  ideal re=imap(@newr,re);
  return(re);   
}
proc quotient4(module a,module b)   //95sec
USAGE:   quotient4(m1, m2); m1, m2 two submodules of k^s, 
         only for global rings
RETURN:  the quotient of m1 and m2
{
  string s="ring @newr=("+charstr(basering)+
           "),("+varstr(basering)+",@t),dp;";
  def @newP=basering;
  execute s;
  module b=imap(@newP,b);
  module a=imap(@newP,a);
  int i;
  vector @b=b[1];
  for(i=2;i<=size(b);i++)
  {
     @b=@b+@t^(i-1)*b[i];
  }
  matrix sy=modulo(@b,a);
  ideal re=sy;
  setring @newP;
  ideal re=imap(@newr,re);
  return(re);   
}
static proc intersect1(ideal i,ideal j)
{
  def R=basering;
  execute "ring gnir = ("+charstr(basering)+"),
                       ("+varstr(basering)+",@t),(C,dp);";
  ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j);
  ideal j=eliminate(i,var(nvars(basering)));
  setring R;
  map phi=gnir,maxideal(1);
  return(phi(j));
}
 
/*
proc minres(list #)
{
  if (size(#) == 2)
  {
    if (typeof(#[1]) == "ideal" || typeof(#[1]) == "module")
    {
      if (typeof(#[2] == "int"))
      {
        return (res(#[1],#[2],1));
      }
    }
  }
  
  if (typeof(#[1]) == "resolution")
  {
    return minimizeres(#[1]);
  }
  else
  {
    return minimizeres(#);
  }
  
}

*/