source: git/Singular/LIB/standard.lib @ d939c1

//////////////////////////////////////////////////////////////////////////////
version="$Id: standard.lib,v 1.83 2006-11-13 10:15:34 Singular Exp $";
category="Miscellaneous";
info="
LIBRARY: standard.lib   Procedures which are always loaded at Start-up

PROCEDURES:
 stdfglm(ideal[,ord])   standard basis of ideal via fglm [and ordering ord]
 stdhilb(ideal[,h])     standard basis of ideal using the Hilbert function
 hilbstd(ideal[,h])     standard basis using (weighted) Hilbert function
 groebner(ideal/module) standard basis using a heuristically chosen method
 res(ideal/module,[i])  free resolution of ideal or module
 sprintf(fmt,...)       returns fomatted string
 fprintf(link,fmt,..)   writes formatted string to link
 printf(fmt,...)        displays formatted string
 weightKB(stc,dd,vl)    degree dd part of a kbase wrt. some weigths
";

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

proc stdfglm (ideal i, list #)
"SYNTAX: @code{stdfglm (} ideal_expression @code{)} @*
         @code{stdfglm (} ideal_expression@code{,} string_expression @code{)}
TYPE:    ideal
PURPOSE: computes the standard basis of the ideal in the basering
         via @code{fglm} (from the ordering given as the second argument
         to the ordering of the basering).@*
         If no second argument is given, \"dp\" is used.
SEE ALSO: fglm, groebner, std, stdhilb
KEYWORDS: fglm
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 #)
"SYNTAX: @code{stdhilb (} ideal_expression @code{)} @*
         @code{stdhilb (} ideal_expression@code{,} intvec_expression @code{)}
TYPE:    ideal
PURPOSE: computes the standard basis of the homogeneous ideal in the basering,
         via a Hilbert driven standard basis computation.@*
         An optional second argument will be used as 1st Hilbert function.
ASSUME:  The optional second argument is the first Hilbert series as computed
         by @code{hilb}.
SEE ALSO: stdfglm, std, groebner
KEYWORDS: 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=ncols(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),dp;
   ideal i=y3+x2,x2y+x2,x3-x2,z4-x2-y;
   ideal i1=stdhilb(i); i1;
   // the latter computation is equivalent to:
   intvec v=hilb(i,1);
   ideal i2=stdhilb(i,v); i2;
}
//////////////////////////////////////////////////////////////////////////

proc hilbstd(ideal i,list #)
"SYNTAX: @code{stdhilb (} ideal_expression @code{)} @*
         @code{stdhilb (} ideal_expression@code{,} intvec_expression @code{)}
TYPE:    ideal
PURPOSE: Computes a Groebner basis of the homogeneous ideal in the basering.
THEORY:  Compute first a standard basis of the (weighted) homogenization
         of the ideal, then the (weighted) Hilbert function and finally
         a Groebner basis in the original ring by using the computed Hilbert
         function@*
         An optional second argument will be used as 1st Hilbert function.
ASSUME:  The optional second argument is the first Hilbert series as computed
         by @code{hilb}.
NOTE:    This procedure makes stdhilb obsolet since it is a generalization
         to rings with some variables having weights >1. Parameters are kept.
SEE ALSO: stdfglm, std, groebner
KEYWORDS: Hilbert function
EXAMPLE: example hilbstd;  shows an example"
{
   def R=basering;
   list lR = ringlist(R);
   intvec v = size(lR[1..size(lR)]);   //size of lists lR[i] in lR
   int n = nvars(R);
   intvec w;                 //ringweights
   int ii, neg;
   for(ii=1; ii<=n; ii++)
   {
     w[ii]=deg(var(ii));
     if(w[ii] <= 0)
     {neg = 1;}
   }
//---------- exclude cases to which hilbstd should no be applied  ----------
   if( homog(i) || find(ordstr(R),"s") || find(ordstr(R),"M")
       || find(ordstr(R),"a") || (neg > 0) )
   {
      if((size(#)!=0) && homog(i))
      {
         return(std(i,#[1]));
      }
      return(std(i));
   }

//----------- create ring for fast computation of hilbert series --------
   list lS=lR;
   lS[2]=insert(lR[2],"@t",v[2]);
   lS[3]=lR[3][1],lR[3][size(lR[3])];
   intvec ww=w,1;
   if(w==1)
   {
     lS[3][1]=list("dp",ww);
   }
   else
   {
     lS[3][1]=list("wp",ww);
   }

   def S = ring(lS);              //ring with one weighted block of variables
   setring S;
   ideal i = homog(imap(R,i),@t); //weighted homog of i
   string s_opt = option();
   int p_opt=(find(s_opt, "prot"));
   if (p_opt) {"std in " + string(S);}
   intvec h = hilb(std(i),1);     //compute weighted hilbert series of i

//------------- use hilbert driven std with original ordering  ------------
   setring R;                     //can access to lR only in R
   lR[2]=lS[2];
   lR[3]=insert(lR[3],list("dp",1),v[3]-1);
   //insert a last block for homogenizing variabble

   def T = ring(lR);             //T = R with 1 homogenizing variable @t
   setring T;                    //added to last block with weight 1
   ideal i=fetch(S,i);           //homogenized i in T

   if (p_opt) {"std with hilb in " + string(T);}
   ideal a=std(i,h,ww);           //use h from S and Hilbert driven std in T

//-------------------- dehomogenize and simplify -------------------------
   a=subst(a,@t,1);              //dehomogenize in T (do not use map!)
   a=simplify(a,34);             //keep only a[j] with different leading terms

//-------------------- back to original ring -------------------------
   setring R;
   ideal a = fetch(T,a);
   attrib(a,"isSB",1);
   return(a);
}
example
{ "EXAMPLE:"; echo = 2;
   ring  r=0,(x,y,z),(wp(43),wp(49,56));
   ideal i=y3+x2,x2y+x2,x3-x2,z4-x2-y;
   ideal i1=hilbstd(i); i1;
   // the latter computation is equivalent to:
   ring r1=0,(x,y,z),wp(43,49,56);
   ideal i = imap(r,i);
   intvec v=hilb(std(i),1);
   setring r;
   ideal i2 = hilbstd(i,v);
}
//////////////////////////////////////////////////////////////////////////

proc groebner(def i, list #)
"SYNTAX: @code{groebner (} ideal_expression @code{)} @*
         @code{groebner (} module_expression @code{)} @*
         @code{groebner (} ideal_expression@code{,} int_expression @code{)} @*
         @code{groebner (} module_expression@code{,} int_expression @code{)}
TYPE:    type of the first argument
PURPOSE: computes the standard basis of the first argument @code{I}
         (ideal or module), by a heuristically chosen method:
         possiblities are @code{std}, @code{slimgb} and/or conversions
         based on @code{fglm}, @code{stdhilb} etc.
         @code{option(prot)} tells about the chosen way.
NOTE: If a 2nd argument @code{wait} is given, then the computation proceeds
      at most @code{wait} seconds. That is, if no result could be computed in
      @code{wait} seconds, then the computation is interrupted, 0 is returned,
      a warning message is displayed, and the global variable
      @code{Standard::groebner_error} is defined.
SEE ALSO: stdhilb, stdfglm, std
KEYWORDS: time limit on computations; MP, groebner basis computations
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);
  int nvarP = nvars(P);
  intvec w;                     //ringweights
  int neg,k;
  for(k=1;  k<=nvarP; k++)
  {
     w[k]=deg(var(k));
     if(w[k] <= 0) {neg=1;}
  }
  if ( ( find(ordstr_P,"s") > 0)
  ||(find(ordstr_P,"M") > 0)
  ||(find(ordstr_P,"a") > 0)
  ||(neg>0) )
  {
    //spaeter den lokalen fall ueber lp oder aehnlich behandeln
    return(std(i));
  }

  if (typeof(basering)=="ring") // slimgb does not know qrings
  {
     //if ordering is global, there are parameters and minpoly is 0
     if (((npars(basering)>0) &&(minpoly==0)))
     { return(slimgb(i)); }
     // ordering dp, char 0
     if((char(P)==0) && (system("nblocks") <= 2) && (find(ordstr_P,"dp")>0))
     { return(slimgb(i)); }
  }
  // for module case, not already hadled by slimgb:
  if (nrows(i)>1)
  { 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); }
  if (find(s_opt, "intStrategy"))
  { option(intStrategy); }

  // construct ring in which first std computation is done
  // CL: 21/09/05 for Singular 3-0 with ringlists....

  list BRlist = ringlist(P);
  int add_vars = 0;
  ideal Qideal = ideal(P);

  if (npars_P > 0)
  {
    for(k=ncols(i); k>0; k--) { i[k]=cleardenom(i[k]); }
  }
  // more than one parameters are converted to ring variables
  if (npars_P > 1)
  {
    for (k=1; k<=npars_P; k++)
    {
      BRlist[2][nvarP+k] = BRlist[1][2][k];
    }
    BRlist[1]=BRlist[1][1];
    add_vars = npars_P;
  }


  // for Hilbert driven approach, Qring structure is removed (defining ideal
  // will be added to the ideal under consideration in the process).
  if (size(BRlist[4])>0)
  {
    int was_qring = 1;
    BRlist[4] = ideal(0);
  }

  // a homogenizing variable is added, if necessary
  int is_homog = (homog(i) && (npars_P <= 1));
  if (! is_homog)
  {
    add_vars = add_vars + 1;
    BRlist[2][nvarP+add_vars] = "@t@";
  }
  // save ordering for later use
  list ord_P = BRlist[3];   // should be ring independent

  //ordering is set to (wp(w,1..1), C) where w are the ringweights
    intvec weight_wp = w;
  for(k=nvarP+1;  k<=nvarP+add_vars; k++)
  {
     weight_wp[k]=1;
  }

  BRlist[3] = list();
  if(weight_wp==1)
  {
    BRlist[3][1]=list("dp",weight_wp);
  }
  else
  {
    BRlist[3][1]=list("wp",weight_wp);
  }
  BRlist[3][2]=list("C",intvec(0));

  //------------ change the ring
  def Phelp = ring(BRlist);
  kill BRlist;
  setring Phelp;

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

  // compute std and hilbert series
  if (p_opt)
  {
    "std in " + string(Phelp);
  }
  intvec hi=hilb(std(qh),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";
    }
    if ( w==1 ) { i = std(i,hi); }
    else { i = std(i,hi,w); }

  }
  else
  {
    // additional variables were introduced
    // need another intermediate ring
    list BRlist = ringlist(Phelp);
    BRlist[3] = list();
    for (k=1; k<=size(ord_P)-1; k++)
    {
      BRlist[3][k] = ord_P[k];
    }

    if( IsSimple_P && (add_vars==1) && (size(ord_P)==2) && (ord_P[1][1]=="lp"))
    {
      // for lp with at most one parameter, we do not need a block ordering
      intvec OW = BRlist[3][1][2];
      OW = OW,1;
      BRlist[3][1][2] = OW; // extend block1 by 1
      BRlist[3][2]=ord_P[2]; // copy block 2
    }
    else
    {
      if( IsSimple_P && (add_vars==1) && (size(ord_P)==2)&&(ord_P[2][1]=="lp"))
      {
        // for lp with at most one parameter, we do not need a block ordering
        intvec OW = ord_P[2][2];
        OW = OW,1;
        BRlist[3][2]=ord_P[2];
        BRlist[3][2][2] = OW; // extend block 2 by 1
      }
      else
      {
        intvec OW = 1;
        for (k=2; k<=add_vars; k++) { OW = OW,1; }
        BRlist[3][size(ord_P)] = list("dp",OW);
        BRlist[3][size(BRlist[3])+1]=ord_P[size(ord_P)];
      }
    }
    // change to intermediate ring
    def Phelp1 = ring(BRlist);
    setring Phelp1;

    ideal qh = imap(Phelp, qh);
    kill Phelp;
    if (p_opt)
    {
      "std with hilb in " + string(Phelp1);;
    }
    // compute std with Hilbert series
    if (weight_wp==1) { qh = std(qh, hi);}
    else { qh = std(qh, hi, weight_wp);}

    // subst 1 for homogenizing var
    if (!is_homog)
    {
      if (p_opt)
      {
        "dehomogenization";
      }
      qh = subst(qh, @t@, 1);
    }

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

  // clean-up time
  option(set, opt);
  if (find(s_opt, "redSB") > 0)
  {
    if (p_opt)
    {
      "interreduction";
    }
    i=interred(i);
  }
  attrib(i, "isSB", 1);
  return (i);
}
example
{ "EXAMPLE: "; echo=2; // LIB "./standard.lib";
  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);
  groebner(i);
  option(noprot);
  if (system("with","MP")) {groebner(i,0);}
  defined(Standard::groebner_error);
}

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

proc res(list #)
"@c we do texinfo here:
@cindex resolution, computation of
@table @code
@item @strong{Syntax:}
@code{res (} ideal_expression@code{,} int_expression @code{[,} any_expression @code{])}
@*@code{res (} module_expression@code{,} int_expression @code{[,} any_expression @code{])}
@item @strong{Type:}
resolution
@item @strong{Purpose:}
computes a (possibly minimal) free resolution of an ideal or module using
a heuristically chosen method.
@* The second (int) argument (say, @code{k}) specifies the length of
the resolution. If it is not positive then @code{k} is assumed to be the
number of variables of the basering.
@* If a third argument is given, the returned resolution is minimized.

Depending on the input, the returned resolution is computed using the
following methods:
@table @asis
@item @strong{quotient rings:}
@code{nres} (classical method using syzygies) , see @ref{nres}.

@item @strong{homogeneous ideals and k=0:}
@code{lres} (La'Scala's method), see @ref{lres}.

@item @strong{not minimized resolution and (homogeneous input with k not 0, or local rings):}
@code{sres} (Schreyer's method), see @ref{sres}.

@item @strong{all other inputs:}
@code{mres} (classical method), see @ref{mres}.
@end table
@item @strong{Note:}
Accessing single elements of a resolution may require that some partial
computations have to be finished and may therefore take some time.
@end table
@c ref
See also
@ref{betti};
@ref{ideal};
@ref{minres};
@ref{module};
@ref{mres};
@ref{nres};
@ref{lres};
@ref{hres};
@ref{sres}.
@ref{resolution}
@c ref
"
{
   def P=basering;
   if (size(#) < 2)
   {
     ERROR("res: need at least two arguments: ideal/module, int");
   }

   def m=#[1]; //the ideal or module
   int i=#[2]; //the length of the resolution
   if (i< 0) { i=0;}

   string varstr_P = varstr(P);

   int p_opt;
   string s_opt = option();
   // set p_opt, if option(prot) is set
   if (find(s_opt, "prot"))
   {
     p_opt = 1;
   }

   if(size(ideal(basering)) > 0)
   {
     // the quick hack for qrings - seems to fit most needs
     // (lres is not implemented for qrings, sres is not so efficient)
     if (p_opt) { "using nres";}
     return(nres(m,i));
   }

   if(homog(m)==1)
   {
      resolution re;
      if (((i==0) or (i>=nvars(basering))) && typeof(m) != "module")
      {
        //LaScala for the homogeneous case and i == 0
        if (p_opt) { "using lres";}
        re=lres(m,i);
        if(size(#)>2)
        {
           re=minres(re);
        }
      }
      else
      {
        if(size(#)>2)
        {
          if (p_opt) { "using mres";}
          re=mres(m,i);
        }
        else
        {
          if (p_opt) { "using sres";}
          re=sres(std(m),i);
        }
      }
      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);
      if (p_opt) { "using mres in another ring";}
      list re=mres(m,i);
      setring P;
      resolution result=imap(Phelp,re);
      if (size(#) > 2) {result = minres(result);}
      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);
      if (p_opt) { "using sres in another ring";}
      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);
    if (p_opt) { "using mres in another ring";}
   list re=mres(m,i);
   setring P;
   resolution result=imap(Phelp,re);
   result = minres(result);
   return(result);
}
example
{"EXAMPLE:"; echo = 2;
  ring r=0,(x,y,z),dp;
  ideal i=xz,yz,x3-y3;
  def l=res(i,0); // homogeneous ideal: uses lres
  l;
  print(betti(l), "betti"); // input to betti may be of type resolution
  l[2];         // element access may take some time
  i=i,x+1;
  l=res(i,0);   // inhomogeneous ideal: uses mres
  l;
  ring rs=0,(x,y,z),ds;
  ideal i=imap(r,i);
  def l=res(i,0); // local ring not minimized: uses sres
  l;
  res(i,0,0);     // local ring and minimized: uses mres
}
/////////////////////////////////////////////////////////////////////////

proc quot (m1,m2,list #)
"SYNTAX: @code{quot (} module_expression@code{,} module_expression @code{)}
         @*@code{quot (} module_expression@code{,} module_expression@code{,}
            int_expression @code{)}
         @*@code{quot (} ideal_expression@code{,} ideal_expression @code{)}
         @*@code{quot (} ideal_expression@code{,} ideal_expression@code{,}
            int_expression @code{)}
TYPE:    ideal
SYNTAX:  @code{quot (} module_expression@code{,} ideal_expression @code{)}
TYPE:    module
PURPOSE: computes the quotient of the 1st and the 2nd argument.
         If a 3rd argument @code{n} is given the @code{n}-th method is used
         (@code{n}=1...5).
SEE ALSO: quotient
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(quotient(m1,m2));
  }
  if (size(#)>0)
  {
    if (typeof(#[1])=="int" )
    {
      return(quot1(m1,m2,#[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 id3=quotient(id1,id2);
  id3;
  ideal id4=quot(id1,id2,1);
  id4;
  ideal id5=quot(id1,id2,2);
  id5;
}

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 quot1; 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=quotient(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=lres(mm,0);
  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=quotient(a,module(B[1]));
  for(i=2;i<=ncols(B);i++)
  {
     re=intersect1(re,quotient(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=ncols(bb);
  ideal re=quotient(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,quotient(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=ncols(b);
  vector @b;
  for(i=1;i<=j;i++)
  {
     @b=@b+@t^(i-1)*@w^(j-i+1)*b[i];
  }
  ideal re=quotient(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=ncols(b);
  vector @b;
  for(i=1;i<=j;i++)
  {
     @b=@b+@t^(i-1)*b[i];
  }
  @b=homog(@b,@w);
  ideal re=quotient(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<=ncols(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));
}

//////////////////////////////////////////////////////////////////
///
/// sprintf, fprintf printf
///
proc sprintf(string fmt, list #)
"SYNTAX:  @code{sprintf (} string_expression @code{[,} any_expressions
               @code{] )}
RETURN:   string
PURPOSE:  @code{sprintf(fmt,...);} performs output formatting. The first
          argument is a format control string. Additional arguments may be
          required, depending on the content of the control string. A series
          of output characters is generated as directed by the control string;
          these characters are returned as a string. @*
          The control string @code{fmt} is simply text to be copied,
          except that the string may contain conversion specifications.@*
          Do @code{help print;} for a listing of valid conversion
          specifications. As an addition to the conversions of @code{print},
          the @code{%n} and @code{%2} conversion specification does not
          consume an additional argument, but simply generates a newline
          character.
NOTE:     If one of the additional arguments is a list, then it should be
          enclosed once more into a @code{list()} command, since passing a list
          as an argument flattens the list by one level.
SEE ALSO: fprintf, printf, print, string
EXAMPLE : example sprintf; shows an example
"
{
  int sfmt = size(fmt);
  if (sfmt  <= 1)
  {
    return (fmt);
  }
  int next, l, nnext;
  string ret;
  list formats = "%l", "%s", "%2l", "%2s", "%t", "%;", "%p", "%b", "%n", "%2";
  while (1)
  {
    if (size(#) <= 0)
    {
      return (ret + fmt);
    }
    nnext = 0;
    while (nnext < sfmt)
    {
      nnext = find(fmt, "%", nnext + 1);
      if (nnext == 0)
      {
        next = 0;
        break;
      }
      l = 1;
      while (l <= size(formats))
      {
        next = find(fmt, formats[l], nnext);
        if (next == nnext) break;
        l++;
      }
      if (next == nnext) break;
    }
    if (next == 0)
    {
      return (ret + fmt);
    }
    if (formats[l] != "%2" && formats[l] != "%n")
    {
      ret = ret + fmt[1, next - 1] + print(#[1], formats[l]);
      # = delete(#, 1);
    }
    else
    {
      ret = ret + fmt[1, next - 1] + print("", "%2s");
    }
    if (size(fmt) <= (next + size(formats[l]) - 1))
    {
      return (ret);
    }
    fmt = fmt[next + size(formats[l]), size(fmt)-next-size(formats[l]) + 1];
  }
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),dp;
  module m=[1,y],[0,x+z];
  intmat M=betti(mres(m,0));
  list l = r, m, M;
  string s = sprintf("s:%s,%n l:%l", 1, 2); s;
  s = sprintf("s:%n%s", l); s;
  s = sprintf("s:%2%s", list(l)); s;
  s = sprintf("2l:%n%2l", list(l)); s;
  s = sprintf("%p", list(l)); s;
  s = sprintf("%;", list(l)); s;
  s = sprintf("%b", M); s;
}

proc printf(string fmt, list #)
"SYNTAX:  @code{printf (} string_expression @code{[,} any_expressions@code{] )}
RETURN:   none
PURPOSE:  @code{printf(fmt,...);} performs output formatting. The first
          argument is a format control string. Additional arguments may be
          required, depending on the content of the control string. A series
          of output characters is generated as directed by the control string;
          these characters are displayed (i.e., printed to standard out). @*
          The control string @code{fmt} is simply text to be copied, except
          that the string may contain conversion specifications. @*
          Do @code{help print;} for a listing of valid conversion
          specifications. As an addition to the conversions of @code{print},
          the @code{%n} and @code{%2} conversion specification does not
          consume an additional argument, but simply generates a newline
          character.
NOTE:     If one of the additional arguments is a list, then it should be
          enclosed once more into a @code{list()} command, since passing a
          list as an argument flattens the list by one level.
SEE ALSO: sprintf, fprintf, print, string
EXAMPLE : example printf; shows an example
"
{
  write("", sprintf(fmt, #));
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),dp;
  module m=[1,y],[0,x+z];
  intmat M=betti(mres(m,0));
  list l=r,m,matrix(M);
  printf("s:%s,l:%l",1,2);
  printf("s:%s",l);
  printf("s:%s",list(l));
  printf("2l:%2l",list(l));
  printf("%p",matrix(M));
  printf("%;",matrix(M));
  printf("%b",M);
}


proc fprintf(link l, string fmt, list #)
"SYNTAX:  @code{fprintf (} link_expression@code{,} string_expression @code{[,}
                any_expressions@code{] )}
RETURN:   none
PURPOSE:  @code{fprintf(l,fmt,...);} performs output formatting.
          The second argument is a format control string. Additional
          arguments may be required, depending on the content of the
          control string. A series of output characters is generated as
          directed by the control string; these characters are
          written to the link l.
          The control string @code{fmt} is simply text to be copied, except
          that the string may contain conversion specifications.@*
          Do @code{help print;} for a listing of valid conversion
          specifications. As an addition to the conversions of @code{print},
          the @code{%n} and @code{%2} conversion specification does not
          consume an additional argument, but simply generates a newline
          character.
NOTE:     If one of the additional arguments is a list, then it should be
          enclosed once more into a @code{list()} command, since passing
          a list as an argument flattens the list by one level.
SEE ALSO: sprintf, printf, print, string
EXAMPLE : example fprintf; shows an example
"
{
  write(l, sprintf(fmt, #));
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),dp;
  module m=[1,y],[0,x+z];
  intmat M=betti(mres(m,0));
  list l=r,m,M;
  link li="";   // link to stdout
  fprintf(li,"s:%s,l:%l",1,2);
  fprintf(li,"s:%s",l);
  fprintf(li,"s:%s",list(l));
  fprintf(li,"2l:%2l",list(l));
  fprintf(li,"%p",list(l));
  fprintf(li,"%;",list(l));
  fprintf(li,"%b",M);
}

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

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

}
*/
///////////////////////////////////////////////////////////////////////////////

proc weightKB(def stc, int dd, list wim)
"SYNTAX: @code{weightKB (} module_expression@code{,} int_expression @code{,}
            list_expression @code{)}@*
         @code{weightKB (} ideal_expression@code{,} int_expression@code{,}
            list_expression @code{)}
RETURN:  the same as the input type of the first argument
PURPOSE:  Computes the part of a vector space basis of the respective quotient with
          weighted degree of the monomials equal to the second argument. The last
          argument contains the information about the weights as a list of intvec:
            @code{wim[1]} for all variables (positive),
            @code{wim[2]} only for module for the generators.
NOTE:     This is a generalisation for the command @code{kbase} with the same first
          two arguments.
SEE ALSO: kbase
EXAMPLE: example weightKB; shows an example
"
{
  if(checkww(wim)){ERROR("wrong weights";);}
  kbclass();
  wwtop=wim[1];
  stc=interred(lead(stc));
  if(typeof(stc)=="ideal")
  {
    stdtop=stc;
    ideal out=widkbase(dd);
    delkbclass();
    return(out);
  }
  list mbase=kbprepare(stc);
  module mout;
  int im,ii;
  if(size(wim)>1){mmtop=wim[2];}
  else{mmtop=0;}
  for(im=size(mbase);im>0;im--)
  {
    stdtop=mbase[im];
    if(im>size(mmtop)){ii=dd;}
    else{ii=dd-mmtop[im];}
    mout=mout+widkbase(ii)*gen(im);
  }
  delkbclass();
  return(mout);
}
///////////////////////////////////////////////////////////////////////////////
// construct global values
static proc kbclass()
{
  intvec wwtop,mmtop;
  export (wwtop,mmtop);
  ideal stdtop,kbtop;
  export (stdtop,kbtop);
  return();
}
// delete global values
static proc delkbclass()
{
  kill wwtop,mmtop;
  kill stdtop,kbtop;
  return();
}
//  select parts of the modul
static proc kbprepare(module mstc)
{
  list rr;
  ideal kk;
  int i1,i2;
  mstc=transpose(mstc);
  for(i1=ncols(mstc);i1>0;i1--)
  {
    kk=0;
    for(i2=nrows(mstc[i1]);i2>0;i2--)
    {
      kk=kk+mstc[i1][i2];
    }
    rr[i1]=kk;
  }
  return(rr);
}
//  check for weights
static proc checkww(list vv)
{
  if(typeof(vv[1])!="intvec"){return(1);}
  intvec ww=vv[1];
  int mv=nvars(basering);
  if(size(ww)<mv){return(1);}
  while(mv>0)
  {
    if(ww[mv]<=0){return(1);}
    mv--;
  }
  if(size(vv)>1)
  {
    if(typeof(vv[2])!="intvec"){return(1);}
  }
  return(0);
}
// The "Caller" for ideals
//    dd   - the degree of the result
static proc widkbase(int dd)
{
  if((size(stdtop)==1)&&(deg(stdtop[1])==0)){return(0);}
  if(dd<=0)
  {
    if(dd<0){return(0);}
    else{return(1);}
  }
  int m1,m2;
  m1=nvars(basering);
  while(wwtop[m1]>dd)
  {
    m1--;
    if(m1==0){return(0);}
  }
  attrib(stdtop,"isSB",1);
  poly mo=1;
  if(m1==1)
  {
    m2=dd/wwtop[1];
    if((m2*wwtop[1])==dd)
    {
      mo=var(1)^m2;
      if(reduce(mo,stdtop)==mo){return(mo);}
      else{return(0);}
    }
  }
  kbtop=0;
  m2=dd;
  weightmon(m1-1,m2,mo);
  while(m2>=wwtop[m1])
  {
    m2=m2-wwtop[m1];
    mo=mo*var(m1);
    if(m2==0)
    {
      if(reduce(mo,stdtop)==mo)
      {
        kbtop=kbtop+mo;
        return(kbtop);
      }
    }
    weightmon(m1-1,m2,mo);
  }
  return(kbtop);
}
// the recursive prozedur
//    va     - number of the variable
//    drest  - rest of the degree
//    mm     - the candidate
static proc weightmon(int va, int drest, poly mm)
{
  while(wwtop[va]>drest)
  {
    va--;
    if(va==0){return();}
  }
  int m2;
  if(va==1)
  {
    m2=drest/wwtop[1];
    if((m2*wwtop[1])==drest)
    {
      mm=mm*var(1)^m2;
      if(reduce(mm,stdtop)==mm){kbtop=kbtop+mm;}
    }
    return();
  }
  m2=drest;
  weightmon(va-1,m2,mm);
  while(m2>=wwtop[va])
  {
    m2=m2-wwtop[va];
    mm=mm*var(va);
    if(m2==0)
    {
      if(reduce(mm,stdtop)==mm)
      {
        kbtop=kbtop+mm;
        return();
      }
    }
    weightmon(va-1,m2,mm);
  }
  return();
}
example
{ "EXAMPLE:"; echo=2;
  ring r=0,(x,y,z),dp;
  ideal i = x6,y4,xyz;
  intvec w = 2,3,6;
  weightKB(i, 12, list(w));
}
///////////////////////////////////////////////////////////////////////////////