source: git/Singular/LIB/standard.lib @ 0266ac

//////////////////////////////////////////////////////////////////////////////
version="$Id: standard.lib,v 1.76 2005-11-04 14:23:10 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
 groebner(ideal/module) standard basis using a heuristically chosen method
 quot(any,any[,n])      quotient using 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 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: 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
         @code{std(I)} is returned. Otherwise, @code{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 @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{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);
  if ((find(ordstr_P,"s") > 0)
  ||(find(ordstr_P,"M") > 0)
  ||(find(ordstr_P,"w") > 0)
  ||(find(ordstr_P,"W") > 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)); }
     // oerdering 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;
  int nvarP = nvars(P);
  int k;
  ideal Qideal = ideal(P);
  string ri;

  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 (dp, C)
  intvec weight_dp = 1;
  for (k=2; k<=nvarP+add_vars; k++) { weight_dp = weight_dp,1; }
  BRlist[3] = list();
  BRlist[3][1]=list("dp",weight_dp);
  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)
  {
    ri = string(Phelp);
    "std in " + ri;
  }
  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";
    }
    i = std(i, hi);
  }
  else
  {
    // additional variables were introduced
    // need another intermediate ring
    list BRlist = ringlist(Phelp);
    BRlist[3] = list();
    for (k=1; k<=size(ord_P); 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;
    }
    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 = BRlist[3][2][2];
        OW = OW,1;
        BRlist[3][2][2] = OW;
      }
      else
      {
        intvec OW = 1;
        for (k=2; k<=add_vars; k++) { OW = OW,1; }
        BRlist[3][size(ord_P)+1] = list("dp",OW);
      }
    }

    // change to intermediate ring
    def Phelp1 = ring(BRlist);
    setring Phelp1;

    ideal qh = imap(Phelp, qh);
    kill Phelp;
    if (p_opt)
    {
      ri = string(Phelp1);
      "std with hilb in " + ri;
    }
    // compute std with Hilbert series
    qh = std(qh, hi);
    // 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;
  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 #)
"@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 'n' is given the n-th method is used
         (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,M;
  printf("s:%s,l:%l",1,2);
  printf("s:%s",l);
  printf("s:%s",list(l));
  printf("2l:%2l",list(l));
  printf("%p",list(l));
  printf("%;",list(l));
  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)
"USAGE:   weightKB(stc,dd,wim);
          stc ideal/module (standard basis),
          dd int (degree of the result),
          wim list of intvec (weights of the degree):
            wim[1] for all variables (>0),
            wim[2] only for module for the generators.
RETURN:  same type as stc
         (the terms of the kbase w.r.t. stc with
          weighted degree equals dd)
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();
}
///////////////////////////////////////////////////////////////////////////////