source: git/Singular/LIB/standard.lib @ 4893c14

//////////////////////////////////////////////////////////////////////////////
version="$Id: standard.lib,v 1.67 2004-08-26 09:42:54 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
 timeStd(i,d)           std(i) if the standard basis computation finished after
                        d-1 seconds and i otherwhise
 timeFactorize(p,d)     factorize(p) if the factorization finished after d-1 
                        seconds otherwhise f is considered to be irreducible
 factorH(p)             changes variables to become the last variable the 
                        principal one in the multivariate factorization and 
                        factorizes then the polynomial

";

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

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

  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
  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 =";

  if (npars_P > 0)
  {
    for(int k=ncols(i); k>0; k--) { i[k]=cleardenom(i[k]); }
  }
  // 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];
  }
  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
    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)
    {
      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 timeFactorize(poly i,list #)
"USAGE:  timeFactorize(p,d)  poly p , integer d
RETURN:  factorize(p) if the factorization finished after d-1 
         seconds otherwhise f is considered to be irreducible
EXAMPLE: example timeFactorize; shows an example
"
{
  def P=basering;
  if (size(#) > 0)
  {
    if (system("with", "MP"))
    {
      if ((typeof(#[1]) == "int")&&(#[1]))
      {
        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(timeFactorize(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;
          }
          kill (l_fork);
        }
        else
        {
          list result;
          intvec v=1,1;
          result[1]=list(1,i);
          result[2]=v;
          j = system("sh", "kill " + string(pid));
        }
        return (result);
      }
    }
  }
  return(factorH(i));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=0,(x,y),dp;
   poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
   p=p^2;
   timeFactorize(p,2);
   "ERROR: will not run currently:";
   //timeFactorize(p,20);
}

proc timeStd(ideal i,list #)
"USAGE:  timeStd(i,d), i ideal, d integer
RETURN:  std(i) if the standard basis computation finished after
         d-1 seconds and i otherwhise
EXAMPLE: example timeStd; shows an example
"
{
  def P=basering;
  if (size(#) > 0)
  {
    if (system("with", "MP"))
    {
      if ((typeof(#[1]) == "int")&&(#[1]))
      {
        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(timeStd(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;
          }
        }
        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;
          }
          kill (l_fork);
        }
        else
        {
          ideal result=i;
          j = system("sh", "kill " + string(pid));
        }
        return (result);
      }
    }
  }
  return(std(i));
}
example
{ "EXAMPLE:"; echo = 2;
   ring r=32003,(a,b,c,d,e),dp;
   int n=6;
   ideal i=
   a^n-b^n,
   b^n-c^n,
   c^n-d^n,
   d^n-e^n,
   a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;
   timeStd(i,2);
   timeStd(i,20);
}

proc factorH(poly p)
"USAGE:  factorH(p)  p poly
RETURN:  factorize(p)
NOTE:    changes variables to become the last variable the principal
         one in the multivariate factorization and factorizes then
         the polynomial
EXAMPLE: example factorH; shows an example
"
{
   def R=basering;
   int i,j;
   int n=1;
   int d=nrows(coeffs(p,var(1)));
   for(i=1;i<=nvars(R);i++)
   {
      j=nrows(coeffs(p,var(i)));
      if(d>j)
      {
         n=i;
         d=j;
      }
   }
   ideal ma=maxideal(1);    //die letzte Variable ist die Hauptvariable
   ma[nvars(R)]=var(n);
   ma[n]=var(nvars(R));
   map phi=R,ma;
   list fac=factorize(phi(p));
   list re=phi(fac);
   return(re);
}
example
{ "EXAMPLE:"; echo = 2;
  system("random",992851144);
  ring r=32003,(x,y,z,w,t),lp;
  poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3;
  factorize(p);  //fast
  system("random",992851262);
  //factorize(p);  //slow
  system("random",992851262);
  factorH(p);
}
//////////////////////////////////////////////////////////////////////////

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

}
*/