Changeset e182c8 in git


Ignore:
Timestamp:
Apr 19, 2005, 5:23:42 PM (18 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
62a35c339c4e8ffbdfed49433b37436dd23ac10e
Parents:
c00e9bd260d0a0003533514dedab696c789b8a08
Message:
*lossen/greuel/hannes: new libs


git-svn-id: file:///usr/local/Singular/svn/trunk@7847 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular/LIB
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/brnoeth.lib

    rc00e9b re182c8  
    1 version="$Id: brnoeth.lib,v 1.13 2004-02-17 16:03:20 Singular Exp $";
     1version="$Id: brnoeth.lib,v 1.14 2005-04-19 15:23:38 Singular Exp $";
    22category="Coding theory";
    33info="
    … …  
    135135    poly A=g1;
    136136    poly B=g2;
    137     ring raux1=char(basering),(x,y,a),lp;
     137    ring raux1=char(basering),(x,y,@a),lp;
    138138    poly G;
    139     ring raux2=(char(basering),a),(x,y),lp;
    140     map psi=BR,x,a;
     139    ring raux2=(char(basering),@a),(x,y),lp;
     140    map psi=BR,x,@a;
    141141    minpoly=number(psi(A));
    142142    poly f=psi(B);
    … …  
    151151      setring raux1;
    152152      execute("G="+sg+";");
    153       G=subst(G,a,y);
     153      G=subst(G,@a,y);
    154154      setring BR;
    155155      map ppssii=raux1,var(1),var(2),0;
    … …  
    374374        // the point is non-rational and a field extension with minpoly=aux
    375375        // is needed
    376         ring r_ext=(char(basering),a),(x,y,z),lp;
     376        ring r_ext=(char(basering),@a),(x,y,z),lp;
    377377        poly F=imap(r_auxz,F);
    378378        poly f_xz=subst(F,y,1);
    379379        poly aux=imap(base_r,aux);
    380         minpoly=number(subst(aux,x,a));
    381         map phi=r_ext,x+a,0,z;
     380        minpoly=number(subst(aux,x,@a));
     381        map phi=r_ext,x+@a,0,z;
    382382        poly f_origin=phi(f_xz);
    383383        if ( subst(subst(diff(f_origin,x),x,0),z,0)==0 &&
    … …  
    769769{
    770770  // writes the generator "a" of a subfield of the coefficients field of
    771   //   basering in terms of the the current generator (also called "a") as a
     771  //   basering in terms of the current generator (also called "a") as a
    772772  //   string sf is an existing ring whose coefficient field is such a subfield
    773773  // warning : in basering there must be a variable called "x" and subfield
    … …  
    11481148  int i,j,k;
    11491149  int m,n;
    1150 //  list L@HNE=essdevelop(CHI);
    1151   list L@HNE=ratdevelop(CHI);
     1150//  list L@HNE=essdevelop(CHI); 
     1151  list LLL=ratdevelop(CHI);
     1152  if (typeof(LLL[1])=="ring") {
     1153    def altring=basering;
     1154    def HNEring = LLL[1]; setring HNEring;
     1155    def L@HNE = hne;
     1156    kill hne;
     1157  }
     1158  else {
     1159    def L@HNE=LLL;
     1160    def HNEring=basering;
     1161    setring HNEring;
     1162  }
     1163  kill LLL;
    11521164  export(L@HNE);
    11531165  int n_branches=size(L@HNE);
    … …  
    20002012          If @code{printlevel>=0} additional comments are displayed (default:
    20012013          @code{printlevel=0}).
     2014WARNING:  The parameter of the needed field extensions is called 'a'. Thus,
     2015          there should be no global object named 'a' when executing Adj_div.
    20022016KEYWORDS: Hamburger-Noether expansions; adjunction divisor
    20032017SEE ALSO: closed_points, NSplaces
    … …  
    20182032  {
    20192033    ERROR("Base field not implemented");
     2034  }
     2035  if (defined(a)==1)
     2036  {
     2037    if ((typeof(a)=="int") or (typeof(a)=="intvec") or (typeof(a)=="intmat") or
     2038        (typeof(a)=="string") or (typeof(a)=="proc") or (typeof(a)=="list"))
     2039    {
     2040      print("WARNING: There is an object named 'a' of the global type "
     2041             + typeof(a)+".");
     2042      print("This will cause problems when "
     2043             + "executing Adj_div, as 'a' is used for "
     2044             + "the name of the parameter of the field extensions needed.");
     2045      ERROR("Please rename or kill the object named 'a'");
     2046    }
    20202047  }
    20212048  intvec opgt=option(get);
    … …  
    21702197          See @ref{Adj_div} for a description of the entries in L.
    21712198NOTE:     The list_expression should be the output of the procedure Adj_div.@*
    2172           If @code{printlevel>=0} additional comments are displayed (default:
    2173           @code{printlevel=0}).
     2199          Raising @code{printlevel}, additional comments are displayed
     2200          (default: @code{printlevel=0}).
     2201WARNING:  The parameter of the needed field extensions is called 'a'. Thus,
     2202          there should be no global object named 'a' when executing NSplaces.
    21742203SEE ALSO: closed_points, Adj_div
    21752204EXAMPLE:  example NSplaces; shows an example
    … …  
    21812210  // notice : if h==0 then nothing is done
    21822211  // warning : non-positive entries are forgotten
     2212  if (defined(a)==1)
     2213  {
     2214    if ((typeof(a)=="int") or (typeof(a)=="intvec") or (typeof(a)=="intmat") or
     2215        (typeof(a)=="string") or (typeof(a)=="proc") or (typeof(a)=="list"))
     2216    {
     2217      print("WARNING: There is an object named 'a' of the global type "
     2218             + typeof(a)+".");
     2219      print("This will cause problems when "
     2220             + "executing Adj_div, as 'a' is used for "
     2221             + "the name of the parameter of the field extensions needed.");
     2222      ERROR("Please rename or kill the object named 'a'");
     2223    }
     2224  }
    21832225  if (h==0)
    21842226  {
    … …  
    22112253          for (j=1;j<=s;j=j+1)
    22122254          {
     2255            dbprint(printlevel-voice-2,"Place No.\ "+string(j));
    22132256            pP[3]=j;
    22142257            update_CURVE=place(pP,0,update_CURVE);
    … …  
    28722915  poly LOC_EQ=phi(CHI);
    28732916  kill(A,B,phi);
    2874   list L@HNE=essdevelop(LOC_EQ)[1];
    2875   export(L@HNE);
    2876   int m=nrows(L@HNE[1]);
    2877   int n=ncols(L@HNE[1]);
    2878   intvec Li2=L@HNE[2];
    2879   int Li3=L@HNE[3];
    2880   setring ES;
    2881   string newa=subfield(HNEring);
    2882   poly paux=importdatum(HNEring,"L@HNE[4]",newa);
    2883   matrix Maux[m][n];
    2884   int i,j;
    2885   for (i=1;i<=m;i=i+1)
    2886   {
    2887     for (j=1;j<=n;j=j+1)
    2888     {
    2889       Maux[i,j]=importdatum(HNEring,"L@HNE[1]["+string(i)+","+
    2890                                                 string(j)+"]",newa);
    2891     }
    2892   }
     2917  list LLL=essdevelop(LOC_EQ);
     2918  if (typeof(LLL[1])=="ring") {
     2919    def altring=basering;
     2920    def HNEring = LLL[1];
     2921    setring HNEring;
     2922    def L@HNE = hne[1];
     2923    kill hne;
     2924    export(L@HNE);
     2925    int m=nrows(L@HNE[1]);
     2926    int n=ncols(L@HNE[1]);
     2927    intvec Li2=L@HNE[2];
     2928    int Li3=L@HNE[3];
     2929    setring ES;
     2930    string newa=subfield(HNEring);
     2931    poly paux=importdatum(HNEring,"L@HNE[4]",newa);
     2932    matrix Maux[m][n];
     2933    int i,j;
     2934    for (i=1;i<=m;i=i+1)
     2935    {
     2936      for (j=1;j<=n;j=j+1)
     2937      {
     2938        Maux[i,j]=importdatum(HNEring,"L@HNE[1]["+string(i)+","+
     2939                                                  string(j)+"]",newa);
     2940      }
     2941    }
     2942  }
     2943  else {
     2944    def L@HNE=LLL[1][1];
     2945    matrix Maux=L@HNE[1];
     2946    intvec Li2=L@HNE[2];
     2947    int Li3=L@HNE[3];
     2948    poly paux=L@HNE[4];
     2949    def HNEring=basering;
     2950    setring HNEring;
     2951  }
     2952  kill LLL;
     2953 
    28932954  list BRANCH=list();
    28942955  BRANCH[1]=Maux;
    … …  
    33213382  matrix H0[1][N];
    33223383  int i,j;
    3323   for (i=1;i<=N;i=i+1)
    3324   {
    3325     H0[1,i]=VmW[k,i];
    3326   }
     3384  for (i=1;i<=N;i=i+1) { H0[1,i]=VmW[k,i]; }
    33273385  poly Ho;
    3328   for (i=1;i<=N;i=i+1)
    3329   {
    3330     Ho=Ho+(H0[1,i]*nFORMS(n)[i]);
    3331   }
     3386  for (i=1;i<=N;i=i+1) { Ho=Ho+(H0[1,i]*nFORMS(n)[i]); }
    33323387  list INTERD=intersection_div(Ho,update_CURVE);
    33333388  intvec NHo=INTERD[1];
    … …  
    47064761static proc ratdevelop (poly chi)
    47074762{
    4708   def r=basering;
     4763  int ring_is_changed;
    47094764  int k1=res_deg();
    4710   list HND=essdevelop(chi);
     4765  list L=essdevelop(chi);
     4766  if (typeof(L[1])=="ring") {
     4767    def altring=basering;
     4768    def HNring = L[1]; setring HNring;
     4769    def HND = hne;
     4770    kill hne;
     4771    ring_is_changed=1;
     4772  }
     4773  else {
     4774    def HND=L[1];
     4775  }
    47114776  int k2=res_deg();
    47124777  if (k1<k2)
    … …  
    47494814    }
    47504815  }
    4751   keepring(HNEring);
    4752   return(HND);
    4753 }
    4754 ;
    4755 
     4816  if (ring_is_changed==0) { return(HND); }
     4817  else { export HND; setring altring; return(HNring); }
     4818}
    47564819
    47574820

  • Singular/LIB/equising.lib

    rc00e9b re182c8  
    1 version="$Id: equising.lib,v 1.9 2005-04-15 15:20:12 Singular Exp $";
     1version="$Id: equising.lib,v 1.10 2005-04-19 15:23:39 Singular Exp $";
    22category="Singularities";
    33info="
    … …  
    619619static proc esComputation (int typ, poly p_F, list #)
    620620{
     621  intvec ov=option(get);  // store options set at beginning
     622  option(redSB);
    621623  // Initialize variables
    622624  int branch=1;
    … …  
    640642    if (typ==1) //isEquising
    641643    {
     644      option(set,ov);
    642645      return(1); 
    643646    }
    644647    else
    645648    {
     649      option(set,ov);
    646650      return(list(ideal(0),0));
    647651    }
    … …  
    660664        if (typ==1) //isEquising
    661665        {
     666          option(set,ov);
    662667          return(1); 
    663668        }
    664669        else
    665670        {
     671          option(set,ov);
    666672          return(list(ideal(0),int(1)));
    667673        }
    … …  
    675681      }     
    676682    }
    677     if (typeof(#[1])=="list")
    678     {
    679       def @L=#[1];  // is assumed to be the Hamburger-Noether matrix
     683    else
     684    {
     685      if (typeof(#)=="list")
     686      {
     687        def @L=#;  // is assumed to be the Hamburger-Noether matrix
     688      }
    680689    }
    681690  }
    … …  
    716725      {
    717726        print("Return value (=0) has no meaning!");
     727        option(set,ov);
    718728        return(0); 
    719729      }
    720730      else
    721731      {
     732        option(set,ov);
    722733        return(list( ideal(0),error));
    723734      }
    … …  
    737748      {
    738749        print("Return value (=0) has no meaning!");
     750        option(set,ov);
    739751        return(0); 
    740752      }
    741753      else
    742754      {
     755        option(set,ov);
    743756        return(list(ideal(0),int(1)));
    744757      }
     758      option(set,ov);
    745759      return(0);
    746760    }
    … …  
    750764        def HNering = LLL[1];
    751765        setring HNering;
    752         def @L=hne;
     766        def @L=stripHNE(hne);
    753767      }
    754768      else {
    755         def @L=LLL;
     769        def @L=stripHNE(LLL);
    756770      }
    757771    }
    … …  
    832846
    833847  J=interred(J);
     848  if (defined(artin_bd)) { J=jet(J,artin_bd-1); }
    834849
    835850// J=std(J);
    … …  
    840855    {
    841856      setring old_ring;
     857      option(set,ov);
    842858      return(0);
    843859    }   
    … …  
    908924  Jnew=coef_Mat[2,1..ncols(coef_Mat)];
    909925
     926  // simplification of Jnew
     927
    910928  if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
    911929                         // deformation parameters can be cutted off
    … …  
    913931    Jnew=jet(Jnew,artin_bd-1);
    914932  }
    915 
    916   // simplification of Jnew
    917933  Jnew=interred(Jnew);
    918 
     934  if (defined(artin_bd)) { Jnew=jet(Jnew,artin_bd-1); }
    919935  J=J,Jnew;
     936
    920937  if (typ==1) // isEquising
    921938  {
    … …  
    923940    {
    924941      setring old_ring;
     942      option(set,ov);
    925943      return(0);
    926944    }   
    927945  }
    928 
    929946
    930947  while (fertig<>soll and blowup<nrows(M[3]))
    … …  
    10381055        if (typ==1) // isEquising
    10391056        {
     1057          if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
    10401058          if(ideal(nselect(J,1,no_b))<>0)
    10411059          {
    10421060            setring old_ring;
     1061            option(set,ov);
    10431062            return(0);
    10441063          }   
    10451064        }
    1046 
    10471065      }
    10481066    }
    … …  
    10521070      if (typ==1) // isEquising
    10531071      {
     1072        if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
    10541073        if(ideal(nselect(J,1,no_b))<>0)
    10551074        {
    10561075          setring old_ring;
     1076          option(set,ov);
    10571077          return(0);
    10581078        }   
    … …  
    10901110  if (typ==1) // isEquising
    10911111  {
     1112    if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
    10921113    if(ideal(nselect(J,1,no_b))<>0)
    10931114    {
    10941115      setring old_ring;
     1116      option(set,ov);
    10951117      return(0);
    10961118    }   
    … …  
    11381160         qring QQ=MM;     
    11391161      }
    1140         else
    1141         {   
    1142            qId=qId,MM;
    1143            qring QQ = std(qId);
    1144         }
     1162      else
     1163      {   
     1164         qId=qId,MM;
     1165         qring QQ = std(qId);
     1166      }
    11451167
    11461168      ideal Shortmap=imap(myRing,Shortmap);
    … …  
    11481170
    11491171      ideal J=phiphi(J);
     1172      option(redSB);
    11501173      J=std(J);
    11511174      J=nselect(J,1,no_b);
    … …  
    11591182      J=J,Jnew;
    11601183      J=interred(J);
     1184      if (defined(artin_bd)){ J=jet(J,artin_bd-1); }
    11611185    }
    11621186    else
    … …  
    11641188      J=std(J);
    11651189      J=nselect(J,1,no_b);
     1190      if (defined(artin_bd)){ J=jet(J,artin_bd-1); }
    11661191    }
    11671192  }
    … …  
    11701195  dbprint(i_print,"// ");
    11711196 
    1172   minPolyStr = "";
     1197  minPolyStr = "";option(set,ov);
    11731198  if (minpoly !=0)
    11741199  {
    … …  
    11801205  if (typ==1) // isEquising
    11811206  {
     1207    if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
    11821208    if(J<>0)
    11831209    {
    11841210      setring old_ring;
     1211      option(set,ov);
    11851212      return(0);
    11861213    }   
    … …  
    11881215    {
    11891216      setring old_ring;
     1217      option(set,ov);
    11901218      return(1);
    11911219    }   
    … …  
    12331261    export(ES_all_triv);
    12341262    setring old_ring;
    1235     dbprint(i_print+1,"
     1263    dbprint(i_print+2,"
    12361264// 'esStratum' created a list M of a ring and an integer.
    12371265// To access the ideal defining the equisingularity stratum, type:
    12381266        def ESSring = M[1]; setring ESSring;  ES; ");
    12391267
     1268    option(set,ov);
    12401269    return(list(ESSring,0));
    12411270  }
    … …  
    12601289//              where all equimultiple sections are trivial
    12611290//    _[2] = 0");
     1291
     1292    option(set,ov);
    12621293    return(list(list(ES,ES_all_triv),0));
    12631294  }
    … …  
    19992030poly F=f+A*y*diff(f,x)+B*x*diff(f,x)+C*diff(f,y);
    20002031list M=esStratum(F,6);
    2001 std(M[1]);  // standard basis of equisingularity ideal
     2032std(M[1][1]);  // standard basis of equisingularity ideal
    20022033
    20032034LIB "equising.lib";
    … …  
    20072038poly f=x20+y7+(x-y)^2*x2y2+x2y4; // Newton non-degenerate
    20082039list K=esIdeal(f);
     2040K;
    20092041
    20102042ring rr=0,(x,y),ls;
    … …  
    20152047poly F=Fs[1,1];
    20162048list M=esStratum(F,2);
     2049
     2050LIB "equising.lib";
     2051ring R=0,(A,B,C,D,x,y),ds;
     2052poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13;
     2053poly F=f+Ax9+Bx7y2+Cx9y+Dx8y2;
     2054list M=esStratum(F,2);
     2055
    20172056
    20182057LIB "equising.lib";
    … …  
    20462085setring Px;
    20472086poly F=Fs[1,1];
    2048 int t=timer;
    20492087list M=esStratum(F);
    20502088timer-t;          //-> 0

  • Singular/LIB/hnoether.lib

    rc00e9b re182c8  
    1 -       version="$Id: hnoether.lib,v 1.45 2005-04-15 13:42:53 Singular Exp $";
     1-       version="$Id: hnoether.lib,v 1.46 2005-04-19 15:23:40 Singular Exp $";
    22category="Singularities";
    33info="
    … …  
    573573 if (defined(HNDebugOn))
    574574 {"received polynomial: ",f,", where x =",namex,", y =",namey;}
    575 
     575 kill m;
    576576 int M,N,Q,R,l,e,hilf,eps,getauscht,Abbruch,zeile,exponent,Ausgabe;
    577577
    … …  
    27692769      "// result: "+string(size(hne))+" branch(es) successfully computed.");
    27702770 }
     2771
     2772 // ----- Lossen 10/02 : the branches have to be resorted to be able to
     2773 // -----                display the multsequence in a nice way
     2774 if (size(hne)>2)
     2775 { 
     2776   int i,j,k,m;
     2777   list dummy;
     2778   int nbsave;
     2779   int no_br = size(hne);
     2780   intmat nbhd[no_br][no_br];
     2781   for (i=1;i<no_br;i++)
     2782   {
     2783     for (j=i+1;j<=no_br;j++) 
     2784     {
     2785       nbhd[i,j]=separateHNE(hne[i],hne[j]);
     2786       k=i+1;
     2787       while ( (nbhd[i,k] >= nbhd[i,j]) and (k<j) )
     2788       {
     2789         k++;
     2790       }
     2791       if (k<j)  // branches have to be resorted
     2792       {
     2793         dummy=hne[j];
     2794         nbsave=nbhd[i,j];
     2795         for (m=k; m<j; m++)
     2796         {
     2797           hne[m+1]=hne[m];
     2798           nbhd[i,m+1]=nbhd[i,m];
     2799         }
     2800         hne[k]=dummy;
     2801         nbhd[i,k]=nbsave;
     2802       }
     2803     }
     2804   }
     2805 }
     2806 // -----
     2807 
    27712808 if (field_ext==1) {
    27722809   dbprint(printlevel-voice+3,"
    … …  
    27912828     list hne=imap(HNEring,hne);
    27922829     setring altring;
    2793      map m=HNhelpring,par(1),var(1),var(2);
    2794      list hne=m(hne);
    2795      kill m,HNhelpring;
     2830     map mmm=HNhelpring,par(1),var(1),var(2);
     2831     list hne=mmm(hne);
     2832     kill mmm,HNhelpring;
    27962833   }
    27972834   else {   

  • Singular/LIB/mregular.lib

    rc00e9b re182c8  
    1 // IB/PG/GMG, last modified:  21.7.2000
     1// IB/PG/GMG, last modified:  15.10.2004
    22//////////////////////////////////////////////////////////////////////////////
    3 version = "$Id: mregular.lib,v 1.6 2001-01-16 13:48:33 Singular Exp $";
     3version = "$Id: mregular.lib,v 1.7 2005-04-19 15:23:41 Singular Exp $";
    44category="Commutative Algebra";
    55info="
    6 LIBRARY: mregular.lib   Castelnuovo-Mumford Regularity of CM-Schemes and Curves
     6LIBRARY: mregular.lib   Castelnuovo-Mumford regularity of homogeneous ideals
    77AUTHORS: I.Bermejo,     ibermejo@ull.es
    88@*       Ph.Gimenez,    pgimenez@agt.uva.es
    … …  
    1010
    1111OVERVIEW:
    12  A library for computing the Castelnuovo-Mumford regularity of a subscheme of
    13  the projective n-space that DOES NOT require the computation of a minimal
    14  graded free resolution of the saturated ideal defining the subscheme.
    15  The procedures are based on two papers by Isabel Bermejo and Philippe Gimenez:
     12 A library for computing the Castelnuovo-Mumford regularity of a homogeneous
     13 ideal that DOES NOT require the computation of a minimal graded free
     14 resolution of the ideal.
     15 It also determines depth(basering/ideal) and satiety(ideal).
     16 The procedures are based on 3 papers by Isabel Bermejo and Philippe Gimenez:
    1617 'On Castelnuovo-Mumford regularity of projective curves' Proc.Amer.Math.Soc.
    17  128(5) (2000), and 'Computing the Castelnuovo-Mumford regularity of some
     18 128(5) (2000), 'Computing the Castelnuovo-Mumford regularity of some
    1819 subschemes of Pn using quotients of monomial ideals', Proceedings of
    19  MEGA-2000, J. Pure Appl. Algebra (to appear).
    20  The algorithm assumes the variables to be in Noether position.
     20 MEGA-2000, J. Pure Appl. Algebra 164 (2001), and 'Saturation and
     21 Castelnuovo-Mumford regularity', Preprint (2004).
    2122
    2223PROCEDURES:
    23  reg_CM(id);         regularity of arith. C-M subscheme V(id_sat) of Pn
    24  reg_curve(id,[,e]); regularity of projective curve V(id_sat) in Pn
    25  reg_moncurve(li);   regularity of projective monomial curve defined by li
     24 regIdeal(id,[,e]);    regularity of homogeneous ideal id
     25 depthIdeal(id,[,e]);  depth of S/id with S=basering, id homogeneous ideal
     26 satiety(id,[,e]);     saturation index of homogeneous ideal id
     27 regMonCurve(li);      regularity of projective monomial curve defined by li
     28 NoetherPosition(id);  Noether normalization of ideal id
     29 is_NP(id);            checks whether variables are in Noether position
     30 is_nested(id);        checks whether monomial ideal id is of nested type
    2631";
    2732
    … …  
    3136LIB "poly.lib";
    3237//////////////////////////////////////////////////////////////////////////////
    33 
    34 proc reg_CM (ideal i)
     38//
     39proc regIdeal (ideal i, list #)
    3540"
    36 USAGE:   reg_CM (i); i ideal
    37 RETURN:  an integer, the Castelnuovo-Mumford regularity of i-sat.
    38 ASSUME:  i is a homogeneous ideal of the basering S=K[x(0)..x(n)] where
    39          the field K is infinite, and S/i-sat is Cohen-Macaulay.
    40          Assume that K[x(n-d),...,x(n)] is a Noether normalization of S/i-sat
    41          where d=dim S/i -1. If this is not the case, compute a Noether
    42          normalization e.g. by using the proc noetherNormal from algebra.lib.
    43 NOTE:    The output is reg(X)=reg(i-sat) where X is the arithmetically
    44          Cohen-Macaulay subscheme of the projective n-space defined by i.
    45          If printlevel > 0 (default = 0) additional information is displayed.
    46          In particular, the value of the regularity of the Hilbert function of
    47          S/i-sat is given.
    48 EXAMPLE: example reg_CM; shows some examples
     41USAGE:   regIdeal (i[,e]); i ideal, e integer
     42RETURN:  an integer, the Castelnuovo-Mumford regularity of i.
     43         (returns -1 if i is not homogeneous)
     44ASSUME:  i is a homogeneous ideal of the basering S=K[x(0)..x(n)].
     45         e=0:  (default)
     46               If K is an infinite field, makes random changes of coordinates.
     47               If K is a finite field, works over a transcendental extension.
     48         e=1:  Makes random changes of coordinates even when K is finite.
     49               It works if it terminates, but may result in an infinite
     50               loop. After 30 loops, a warning message is displayed and
     51               -1 is returned.
     52NOTE:    If printlevel > 0 (default = 0), additional info is displayed:
     53         dim(S/i), depth(S/i) and end(H^(depth(S/i))(S/i)) are computed,
     54         and an upper bound for the a-invariant of S/i is given.
     55         The algorithm also determines whether the regularity is attained
     56         or not at the last step of a minimal graded free resolution of i,
     57         and if the answer is positive, the regularity of the Hilbert
     58         function of S/i is given.
     59EXAMPLE: example regIdeal; shows some examples
    4960"
    5061{
    5162//--------------------------- initialisation ---------------------------------
    52    int ii,H,h,d,time,si;
     63   int e,ii,jj,H,h,d,time,lastv,sat,firstind;
     64   int lastind,ch,nesttest,NPtest,nl,N,acc;
     65   intmat ran;
    5366   def r0 = basering;
    5467   int n = nvars(r0)-1;
     68   if ( size(#) > 0 )
     69   {
     70      e = #[1];
     71   }
    5572   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
    5673   execute(s);
    57    ideal i,j,I,K;
    58    j = fetch(r0,i);
    59 //----- Checks saturated ideal, computes saturation if necessary,
    60 //      and computes the initial ideal of the i-sat w.r.t. dp-ordering
     74   ideal i,sbi,I,J,K,chcoord,m;
     75   poly P;
     76   map phi;
     77   i = fetch(r0,i);
    6178   time=rtimer;
    62    list l=sat(j,maxideal(1));
    63    i=l[1]; si=l[2];
    64    I=lead(i);
     79   sbi=std(i);
     80   ch=char(r1);
     81//----- Check ideal homogeneous
     82   if ( homog(sbi) == 0 )
     83   {
     84       "// WARNING from proc regIdeal from lib mregular.lib:
     85// The ideal is not homogeneous!";
     86       return (-1);
     87   }
     88   I=simplify(lead(sbi),1);
    6589   attrib(I,"isSB",1);
    6690   d=dim(I);
     91//----- If the ideal i is not proper:
    6792   if ( d == -1 )
    68    {
    69    H=maxdeg1(quotient(lead(std(j)),maxideal(1)))+1;
    70        "// WARNING from proc reg_CM from lib mregular.lib:
    71 // your ideal i of S is zero-dimensional!
    72 // The Castelnuovo-Mumford regularity of i coincides with the
    73 // regularity of the Hilbert function of S/i and its value is:";
     93     {
     94       dbprint(printlevel-voice+2,
     95               "// The ideal i is (1)!
     96// Its Castelnuovo-Mumford regularity is:");
     97       return (0);
     98     }
     99//----- If the ideal i is 0:
     100   if ( size(I) == 0 )
     101     {
     102       dbprint(printlevel-voice+2,
     103               "// The ideal i is (0)!
     104// Its Castelnuovo-Mumford regularity is:");
     105       return (0);
     106     }
     107//----- When the ideal i is 0-dimensional:
     108   if ( d == 0 )
     109     {
     110       H=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
     111       time=rtimer-time;
     112       // Additional information:
     113       dbprint(printlevel-voice+2,
     114               "// Dimension of S/i : 0");
     115       dbprint(printlevel-voice+2,
     116               "// Time for computing regularity: " + string(time) + " sec.");
     117       dbprint(printlevel-voice+2,
     118"// The Castelnuovo-Mumford regularity of i coincides with its satiety, and
     119 // with the regularity of the Hilbert function of S/i. Its value is:");
    74120       return (H);
    75    }
    76 //----- Check Noether position
    77    ideal J=I;
    78    for ( ii = n-d+1; ii <= n; ii++ )
    79    {
    80    J=subst(J,x(ii),0);
    81    }
    82    attrib(J,"isSB",1);
    83    int dz=dim(J);
    84    if ( dz != d )
    85    {
    86        "// WARNING from proc reg_CM from lib mregular.lib:
    87 // the variables are not in Noether position!";
     121     }
     122//----- Determine the situation: NT, or NP, or nothing.
     123//----- Choose the method depending on the situation, on the
     124//----- characteristic of the ground field, and on the option argument
     125//----- in order to get the mon. ideal of nested type associated to i
     126   if ( e == 1 )
     127     {
     128       ch=0;
     129     }
     130   NPtest=is_NP(I);
     131   if ( NPtest == 1 )
     132     {
     133       nesttest=is_nested(I);
     134     }
     135   if ( ch != 0 )
     136     {
     137       if ( NPtest == 0 )
     138         {
     139           N=d*n-d*(d-1)/2;
     140           s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
     141           execute(s);
     142           ideal chcoord,m,i,I;
     143           poly P;
     144           map phi;
     145           i=imap(r1,i);
     146           chcoord=select1(maxideal(1),1,(n-d+1));
     147           acc=0;
     148           for ( ii = 1; ii<=d; ii++ )
     149             {
     150               matrix trex[1][n-d+ii+1]=t((1+acc)..(n-d+ii+acc)),1;
     151               m=select1(maxideal(1),1,(n-d+1+ii));
     152               for ( jj = 1; jj<=n-d+ii+1; jj++ )
     153                 {
     154                   P=P+trex[1,jj]*m[jj];
     155                 }
     156               chcoord[n-d+1+ii]=P;
     157               P=0;
     158               acc=acc+n-d+ii;
     159               kill trex;
     160             }
     161               phi=rtr,chcoord;
     162               I=simplify(lead(std(phi(i))),1);
     163               setring r1;
     164               I=imap(rtr,I);
     165               attrib(I,"isSB",1);
     166         }
     167       else
     168         {
     169           if ( nesttest == 0 )
     170             {
     171               N=d*(d-1)/2;
     172               s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
     173               execute(s);
     174               ideal chcoord,m,i,I;
     175               poly P;
     176               map phi;
     177               i=imap(r1,i);
     178               chcoord=select1(maxideal(1),1,(n-d+2));
     179               acc=0;
     180               for ( ii = 1; ii<=d-1; ii++ )
     181                 {
     182                   matrix trex[1][ii+1]=t((1+acc)..(ii+acc)),1;
     183                   m=select1(maxideal(1),(n-d+2),(n-d+2+ii));
     184                   for ( jj = 1; jj<=ii+1; jj++ )
     185                     {
     186                       P=P+trex[1,jj]*m[jj];
     187                     }
     188                   chcoord[n-d+2+ii]=P;
     189                   P=0;
     190                   acc=acc+ii;
     191                   kill trex;
     192                 }
     193               phi=rtr,chcoord;
     194               I=simplify(lead(std(phi(i))),1);
     195               setring r1;
     196               I=imap(rtr,I);
     197               attrib(I,"isSB",1);
     198             }
     199         }
     200     }
     201   else
     202     {
     203       if ( NPtest == 0 )
     204         {
     205           while ( nl < 30 )
     206             {
     207               chcoord=select1(maxideal(1),1,(n-d+1));
     208               nl=nl+1;
     209               for ( ii = 1; ii<=d; ii++ )
     210                 {
     211                   ran=random(100,1,n-d+ii);
     212                   ran=intmat(ran,1,n-d+ii+1);
     213                   ran[1,n-d+ii+1]=1;
     214                   m=select1(maxideal(1),1,(n-d+1+ii));
     215                   for ( jj = 1; jj<=n-d+ii+1; jj++ )
     216                     {
     217                       P=P+ran[1,jj]*m[jj];
     218                     }
     219                   chcoord[n-d+1+ii]=P;
     220                   P=0;
     221                 }
     222               phi=r1,chcoord;
     223               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
     224               I=simplify(lead(std(phi(i))),1);
     225               attrib(I,"isSB",1);
     226               NPtest=is_NP(I);
     227               if ( NPtest == 1 )
     228                 {
     229                   break;
     230                 }
     231             }
     232           if ( NPtest == 0 )
     233             {
     234       "// WARNING from proc regIdeal from lib mregular.lib:
     235// The procedure has entered in 30 loops and could not put the variables
     236// in Noether position: in your example the method using random changes
     237// of coordinates may enter an infinite loop when the field is finite.
     238// Try removing this optional argument.";
    88239       return (-1);
    89    }
    90 //----- Check Cohen-Macaulay property of S/i-sat
     240             }
     241           i=phi(i);
     242           nesttest=is_nested(I);
     243         }
     244       if ( nesttest == 0 )
     245         {
     246           while ( nl < 30 )
     247             {
     248               chcoord=select1(maxideal(1),1,(n-d+2));
     249               nl=nl+1;
     250               for ( ii = 1; ii<=d-1; ii++ )
     251                 {
     252                   ran=random(100,1,ii);
     253                   ran=intmat(ran,1,ii+1);
     254                   ran[1,ii+1]=1;
     255                   m=select1(maxideal(1),(n-d+2),(n-d+2+ii));
     256                   for ( jj = 1; jj<=ii+1; jj++ )
     257                     {
     258                       P=P+ran[1,jj]*m[jj];
     259                     }
     260                   chcoord[n-d+2+ii]=P;
     261                   P=0;
     262                 }
     263               phi=r1,chcoord;
     264               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
     265               I=simplify(lead(std(phi(i))),1);
     266               attrib(I,"isSB",1);
     267               nesttest=is_nested(I);
     268               if ( nesttest == 1 )
     269                 {
     270                   break;
     271                 }
     272             }
     273           if ( nesttest == 0 )
     274             {
     275       "// WARNING from proc regIdeal from lib mregular.lib:
     276// The procedure has entered in 30 loops and could not find a monomial
     277// ideal of nested type with the same regularity as your ideal: in your
     278// example the method using random changes of coordinates may enter an
     279// infinite loop when the field is finite.
     280// Try removing this optional argument.";
     281       return (-1);
     282             }
     283         }
     284     }
     285//
     286// At this stage, we have obtained a monomial ideal I of nested type
     287// such that reg(i)=reg(I). We now compute reg(I).
     288//
     289//----- When S/i is Cohen-Macaulay:
    91290   for ( ii = n-d+2; ii <= n+1; ii++ )
    92    {
    93    K=K+select(I,ii);
    94    }
    95    if ( size(K) != 0 )
    96    {
    97          "// WARNING from proc reg_CM from lib mregular.lib:
    98 // the ring basering/i-sat is NOT Cohen-Macaulay !!";
    99          return (-1);
    100    }
    101 // Now, compute the regularity!
    102    s="ring r2 = ",charstr(r0),",x(0..n-d),dp;";
    103    execute(s);
    104    ideal I,qq;
    105    I = imap(r1,I);
    106    qq=quotient(I,maxideal(1));
    107    H=maxdeg1(qq)+1; // The value of the regularity.
    108    time=rtimer-time;
    109 // Additional information:
    110    dbprint(printlevel-voice+2,
    111 "// Ideal i of S defining an arithm. Cohen-Macaulay subscheme X of P"+
    112            string(n) + ":");
    113    dbprint(printlevel-voice+2,
    114 "//   - dimension of X: "+string(d-1));
    115    if ( si == 0 )
    116    {
    117    dbprint(printlevel-voice+2,"//   - i is saturated: YES");
    118    }
    119    else
    120    {
    121    dbprint(printlevel-voice+2,"//   - i is saturated: NO");
    122    }
    123    dbprint(printlevel-voice+2,
    124 "//   - regularity of the Hilbert function of S/i-sat: " + string(H-d));
    125    dbprint(printlevel-voice+2,
    126 "//   - time for computing reg(X): " + string(time) + " sec.");
    127    dbprint(printlevel-voice+2,
    128 "// Castelnuovo-Mumford regularity of X:");
    129    return(H);
     291     {
     292       K=K+select(I,ii);
     293     }
     294   if ( size(K) == 0 )
     295     {
     296       s="ring nr = ",charstr(r0),",x(0..n-d),dp;";
     297       execute(s);
     298       ideal I;
     299       I = imap(r1,I);
     300       H=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
     301       time=rtimer-time;
     302       // Additional information:
     303       dbprint(printlevel-voice+2,
     304               "// S/i is Cohen-Macaulay");
     305       dbprint(printlevel-voice+2,
     306               "// Dimension of S/i ( = depth(S/i) ): "+string(d));
     307       dbprint(printlevel-voice+2,
     308               "// Regularity attained at the last step of m.g.f.r. of i: YES");
     309       dbprint(printlevel-voice+2,
     310               "// Regularity of the Hilbert function of S/i: " + string(H-d));
     311       dbprint(printlevel-voice+2,
     312               "// Time for computing regularity: " + string(time) + " sec.");
     313       dbprint(printlevel-voice+2,
     314               "// The Castelnuovo-Mumford regularity of i is:");
     315       return(H);
     316     }
     317//----- When d=1:
     318   if ( d == 1 )
     319     {
     320       H=maxdeg1(simplify(reduce(quotient(I,maxideal(1)),I),2))+1;
     321       sat=H;
     322       J=subst(I,x(n),1);
     323       s = "ring nr = ",charstr(r0),",x(0..n-1),dp;";
     324       execute(s);
     325       ideal J=imap(r1,J);
     326       attrib(J,"isSB",1);
     327       h=maxdeg1(minbase(quotient(J,maxideal(1))))+1;
     328       time=rtimer-time;
     329       if ( h > H )
     330         {
     331           H=h;
     332         }
     333       // Additional information:
     334       dbprint(printlevel-voice+2,
     335               "// Dimension of S/i: 1");
     336       dbprint(printlevel-voice+2,
     337               "// Depth of S/i: 0");
     338       dbprint(printlevel-voice+2,
     339               "// Satiety of i: "+string(sat));
     340       dbprint(printlevel-voice+2,
     341               "// Upper bound for the a-invariant of S/i: end(H^1(S/i)) <= "+
     342               string(h-2));
     343       if ( H == sat )
     344         {
     345           dbprint(printlevel-voice+2,
     346                   "// Regularity attained at the last step of m.g.f.r. of i: YES");
     347           dbprint(printlevel-voice+2,
     348                   "// Regularity of the Hilbert function of S/i: "+string(H));
     349         }
     350       else
     351         {
     352           dbprint(printlevel-voice+2,
     353                   "// Regularity attained at the last step of m.g.f.r. of i: NO");
     354         }
     355           dbprint(printlevel-voice+2,
     356                   "// Time for computing regularity: "+ string(time) + " sec.");
     357           dbprint(printlevel-voice+2,
     358                   "// The Castelnuovo-Mumford regularity of i is:");
     359           return(H);
     360     }
     361//----- Now d>1 and S/i is not Cohen-Macaulay:
     362//
     363//----- First, determine the last variable really occuring
     364       lastv=n-d;
     365       h=n;
     366       while ( lastv == n-d and h > n-d )
     367         {
     368           K=select(I,h+1);
     369           if ( size(K) == 0 )
     370             {
     371               h=h-1;
     372             }
     373           else
     374             {
     375               lastv=h;
     376             }
     377         }
     378//----- and compute Castelnuovo-Mumford regularity:
     379       s = "ring nr = ",charstr(r0),",x(0..lastv),dp;";
     380       execute(s);
     381       ideal I,K,KK,LL;
     382       I=imap(r1,I);
     383       attrib(I,"isSB",1);
     384       K=simplify(reduce(quotient(I,maxideal(1)),I),2);
     385       H=maxdeg1(K)+1;
     386       firstind=H;
     387       KK=minbase(subst(I,x(lastv),1));
     388       for ( ii = n-lastv; ii<=d-2; ii++ )
     389         {
     390           LL=minbase(subst(I,x(n-ii-1),1));
     391           attrib(LL,"isSB",1);
     392           s = "ring mr = ",charstr(r0),",x(0..n-ii-1),dp;";
     393           execute(s);
     394           ideal K,KK;
     395           KK=imap(nr,KK);
     396           attrib(KK,"isSB",1);
     397           K=simplify(reduce(quotient(KK,maxideal(1)),KK),2);
     398           h=maxdeg1(K)+1;
     399           if ( h > H )
     400             {
     401               H=h;
     402             }
     403           setring nr;
     404           kill mr;
     405           KK=LL;
     406         }
     407       // We must determine one more sat. index:
     408       s = "ring mr = ",charstr(r0),",x(0..n-d),dp;";
     409       execute(s);
     410       ideal KK,K;
     411       KK=imap(nr,KK);
     412       attrib(KK,"isSB",1);
     413       K=simplify(reduce(quotient(KK,maxideal(1)),KK),2);
     414       h=maxdeg1(K)+1;
     415       lastind=h;
     416       if ( h > H )
     417         {
     418           H=h;
     419         }
     420       setring nr;
     421       kill mr;
     422       time=rtimer-time;
     423       // Additional information:
     424       dbprint(printlevel-voice+2,
     425               "// Dimension of S/i: "+string(d));
     426       dbprint(printlevel-voice+2,
     427               "// Depth of S/i: "+string(n-lastv));
     428       dbprint(printlevel-voice+2,
     429               "// end(H^"+string(n-lastv)+"(S/i)) = "
     430               +string(firstind-n+lastv-1));
     431       dbprint(printlevel-voice+2,
     432               "// Upper bound for the a-invariant of S/i: end(H^"
     433               +string(d)+"(S/i)) <= "+string(lastind-d-1));
     434       if ( H == firstind )
     435         {
     436           dbprint(printlevel-voice+2,
     437                   "// Regularity attained at the last step of m.g.f.r. of i: YES");
     438           dbprint(printlevel-voice+2,
     439                   "// Regularity of the Hilbert function of S/i: "
     440                   +string(H-n+lastv));
     441         }
     442       else
     443         {
     444           dbprint(printlevel-voice+2,
     445                   "// Regularity attained at the last step of m.g.f.r. of i: NO");
     446         }
     447       dbprint(printlevel-voice+2,
     448               "// Time for computing regularity: "+ string(time) + " sec.");
     449       dbprint(printlevel-voice+2,
     450               "// The Castelnuovo-Mumford regularity of i is:");
     451       return(H);
    130452}
    131453example
    132454{ "EXAMPLE:"; echo = 2;
     455   ring r=0,(x,y,z,t,w),dp;
     456   ideal i=y2t,x2y-x2z+yt2,x2y2,xyztw,x3z2,y5+xz3w-x2zw2,x7-yt2w4;
     457   regIdeal(i);
     458   regIdeal(lead(std(i)));
     459// Additional information is displayed if you change printlevel (=1);
     460}
     461////////////////////////////////////////////////////////////////////////////////
     462/*
     463Out-commented examples:
     464//
    133465   ring s=0,x(0..5),dp;
    134466   ideal i=x(2)^2-x(4)*x(5),x(1)*x(2)-x(0)*x(5),x(0)*x(2)-x(1)*x(4),
    135467           x(1)^2-x(3)*x(5),x(0)*x(1)-x(2)*x(3),x(0)^2-x(3)*x(4);
    136    reg_CM(i);
    137 // Additional information can be obtained as follows:
    138    printlevel = 1;
    139    reg_CM(i);
    140 }
    141 
    142 ///////////////////////////////////////////////////////////////////////////////
    143 /*
    144 Out-commented examples:
    145    ring r=0,(x,y,z,t),dp;
    146    ideal j=x17y14-y31, x20y13, x60-y36z24-x20z20t20;
    147    reg_CM(j);
    148 //
    149 // The polynomial ring in the last variables MUST be a Noether
    150 // Normalization of basering/ideal:
    151    ring rr=0,(t,x,y,z),dp;
    152    ideal j=imap(r,i);
    153 // The same ideal as before but with a different order of the var. in ring
    154    reg_CM(j);
    155 //
    156 // When S/i-sat is not Cohen-Macaulay, one gets an error message:
     468   regIdeal(i);
     469   // Our procedure works when a min. graded free resol. can
     470   // not be computed. In this easy example, regularity can also
     471   // be obtained using a m.g.f.r.:
     472   nrows(betti(mres(i,0)));
    157473   ring r1=0,(x,y,z,t),dp;
    158    ideal i=y4-t3z, x3t-y2z2, x3y2-t2z3, x6-tz5;
    159    reg_CM(i);
    160 //
    161 // When i is zero-dimensional, one gets an error message but the regularity
    162 // of the ideal is computed:
    163    reg_CM(maxideal(4));
     474// Ex.2.5 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
     475   ideal i  = x17y14-y31, x20y13, x60-y36z24-x20z20t20;
     476   regIdeal(i);
     477// Ex.2.9 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
     478   int k=43;
     479   ideal j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
     480   regIdeal(j);
     481   k=14;
     482   j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
     483   regIdeal(j);
     484   k=22;
     485   j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
     486   regIdeal(j);
     487   k=315;
     488   j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
     489   regIdeal(j);
     490// Example in Rk.2.10 in [Bermejo-Gimenez], ProcAMS 128(5):
     491   ideal h=x2-3xy+5xt,xy-3y2+5yt,xz-3yz,2xt-yt,y2-yz-2yt;
     492   regIdeal(h);
     493// The initial ideal is not saturated
     494   regIdeal(lead(std(h)));
     495// More examples:
     496   i=y4-t3z, x3t-y2z2, x3y2-t2z3, x6-tz5;
     497   regIdeal(i);
     498//
     499   regIdeal(maxideal(4));
    164500//
    165501   ring r2=0,(x,y,z,t,w),dp;
    166    ideal i=xy-zw,x3-yw2,x2z-y2w,y3-xz2,-y2z3+xw4+tw4+w5,-yz4+x2w3+xtw3+xw4,
    167      -z5+x2tw2+x2w3+yw4;
    168    reg_CM(i);
     502   ideal i = xy-zw,x3-yw2,x2z-y2w,y3-xz2,-y2z3+xw4+tw4+w5,-yz4+x2w3+xtw3+xw4,
     503            -z5+x2tw2+x2w3+yw4;
     504   regIdeal(i);
     505//
    169506   ring r3=0,(x,y,z,t,w,u),dp;
    170507   ideal i=imap(r2,i);
    171    reg_CM(i);
    172 //
     508   regIdeal(i);
    173509// Next example is the defining ideal of the 2nd. Veronesean of P3, a variety
    174510// in P8 which is arithmetically Cohen-Macaulay:
    … …  
    179515   ring r5=0,x(0..9),dp;
    180516   ideal i=imap(r4,ei);
    181    reg_CM(i);
    182 //
    183 // Now let's build a non saturated ideal defining the same subscheme:
    184    ideal mi=intersect(i,maxideal(3));
    185    size(mi);
    186    reg_CM(mi);
    187 // The result is the same since both ideals define the same subscheme.
    188 //
     517   regIdeal(i);
     518// Here is an example where the computation of a m.g.f.r. of I costs:
     519   ring r8=0,(x,y,z,t,u,a,b),dp;
     520   ideal i=u-b40,t-a40,x-a23b17,y-a22b18+ab39,z-a25b15;
     521   ideal ei=eliminate(i,ab); // It takes a few seconds to compute the ideal
     522   ring r9=0,(x,y,z,t,u),dp;
     523   ideal i=imap(r8,ei);
     524   regIdeal(i);   // This is very fast.
     525// Now you can use mres(i,0) to compute a m.g.f.r. of the ideal!
     526//
     527// The computation of the m.g.f.r. of the following example did not succeed
     528// using the command mres:
     529   ring r10=0,(x(0..8),s,t),dp;
     530   ideal i=x(0)-st24,x(1)-s2t23,x(2)-s3t22,x(3)-s9t16,x(4)-s11t14,x(5)-s18t7,
     531           x(6)-s24t,x(7)-t25,x(8)-s25;
     532   ideal ei=eliminate(i,st);
     533   ring r11=0,x(0..8),dp;
     534   ideal i=imap(r10,ei);
     535   regIdeal(i);
     536// More examples where not even sres works:
     537// Be careful: elimination takes some time here, but it succeeds!
     538   ring r12=0,(s,t,u,x(0..14)),dp;
     539   ideal i=x(0)-st6u8,x(1)-s5t3u7,x(2)-t11u4,x(3)-s9t4u2,x(4)-s2t7u6,x(5)-s7t7u,
     540           x(6)-s10t5,x(7)-s4t6u5,x(8)-s13tu,x(9)-s14u,x(10)-st2u12,x(11)-s3t9u3,
     541           x(12)-s15,x(13)-t15,x(14)-u15;
     542   ideal ei=eliminate(i,stu);
     543   size(ei);
     544   ring r13=0,x(0..14),dp;
     545   ideal i=imap(r12,ei);
     546   size(i);
     547   regIdeal(i);
    189548*/
    190 //////////////////////////////////////////////////////////////////////////////
     549///////////////////////////////////////////////////////////////////////////////
     550///////////////////////////////////////////////////////////////////////////////
     551///////////////////////////////////////////////////////////////////////////////
    191552
    192 proc reg_curve (ideal i, list #)
     553proc depthIdeal (ideal i, list #)
    193554"
    194 USAGE:   reg_curve (i[,e]); i ideal, e integer
    195 RETURN:  an integer, the Castelnuovo-Mumford regularity of i-sat.
    196 ASSUME:  i is a homogeneous ideal of the basering S=K[x(0)..x(n)] where
    197          the field K is infinite, and it defines a projective curve C in
    198          the projective n-space (dim(i)=2). We assume that K[x(n-1),x(n)]
    199          is a Noether normalization of S/i-sat.
    200          e=0: (default)
    201               Uses a random choice of an element of K when it is necessary.
    202               This is absolutly safe (if the element is bad, another random
    203               choice will be done until a good element is found).
    204          e=1: Substitutes the random choice of an element of K by a simple
    205               transcendental field extension of K.
    206 NOTE:    The output is the integer reg(C)=reg(i-sat).
    207          If printlevel > 0 (default = 0) additional information is displayed.
    208          In particular, says if C is arithmetically Cohen-Macaulay or not,
    209          determines in which step of a minimal graded free resolution of i-sat
    210          the regularity of C is attained, and sometimes gives the value of the
    211          regularity of the Hilbert function of S/i-sat (otherwise, an upper
    212          bound is given).
    213 EXAMPLE: example reg_curve; shows some examples
     555USAGE:   depthIdeal (i[,e]); i ideal, e integer
     556RETURN:  an integer, the depth of S/i where S=K[x(0)..x(n)] is the basering.
     557         (returns -1 if i is not homogeneous or if i=(1))
     558ASSUME:  i is a proper homogeneous ideal.
     559         e=0:  (default)
     560               If K is an infinite field, makes random changes of coordinates.
     561               If K is a finite field, works over a transcendental extension.
     562         e=1:  Makes random changes of coordinates even when K is finite.
     563               It works if it terminates, but may result in an infinite
     564               loop. After 30 loops, a warning message is displayed and
     565               -1 is returned.
     566NOTE:    If printlevel > 0 (default = 0), dim(S/i) is also displayed.
     567EXAMPLE: example depthIdeal; shows some examples
    214568"
    215569{
    216570//--------------------------- initialisation ---------------------------------
    217    int d,ii,jj,H,HR,e,dd,h,hh,hm,sK,ts,si,time,ran,rant;
     571   int e,ii,jj,h,d,time,lastv,ch,nesttest,NPtest,nl,N,acc;
     572   intmat ran;
    218573   def r0 = basering;
    219574   int n = nvars(r0)-1;
    … …  
    224579   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
    225580   execute(s);
    226    ideal i,j,I,I0,J,K,II,q,qq,m;
     581   ideal i,sbi,I,J,K,chcoord,m;
     582   poly P;
     583   map phi;
    227584   i = fetch(r0,i);
    228 //----- Checks saturated ideal, computes saturation if necessary,
    229 //      and computes the initial ideal of the i-sat w.r.t. dp-ordering
    230585   time=rtimer;
    231    list l=sat(i,maxideal(1));
    232    i=l[1];si=l[2];
    233    I=lead(i);
     586   sbi=std(i);
     587   ch=char(r1);
     588//----- Check ideal homogeneous
     589   if ( homog(sbi) == 0 )
     590   {
     591       "// WARNING from proc depthIdeal from lib mregular.lib:
     592// The ideal is not homogeneous!";
     593       return (-1);
     594   }
     595   I=simplify(lead(sbi),1);
    234596   attrib(I,"isSB",1);
    235597   d=dim(I);
    236 //----- Check if the ideal defines a curve
    237    if ( d != 2 )
    238    {
    239    "// WARNING from proc reg_curve from lib mregular.lib:
    240 // your ideal does not define a projective curve!";
    241    return (-1);
    242    }
    243 //----- Check Noether position
    244    J=subst(I,x(n),0);
    245    K=subst(J,x(n-1),0);
    246    attrib(K,"isSB",1);
    247    int dz=dim(K);
    248    if ( dz != 2 )
    249    {
    250        "// WARNING from proc reg_curve from lib mregular.lib:
    251 // the variables are not in Noether position!";
     598//----- If the ideal i is not proper:
     599   if ( d == -1 )
     600     {
     601       "// WARNING from proc depthIdeal from lib mregular.lib:
     602// The ideal i is (1)!";
    252603       return (-1);
    253    }
    254 //--------- When S/i-sat is Cohen-Macaulay we can compute regularity:
    255    K=select(I,n+1);
    256    K=K+select(I,n);
    257    sK=size(K);
    258    s="ring r2 = ",charstr(r0),",x(0..n-2),dp;";
    259    if (sK == 0)
    260    {
    261    execute(s);
    262    ideal I,qq;
    263    I = imap(r1,I);
    264    qq=quotient(I,maxideal(1));
    265    H=maxdeg1(qq)+1; // this is the value of the regularity.
    266    time=rtimer-time;
    267 // Additional information:
    268    dbprint(printlevel-voice+2,
    269 "// Ideal i of S defining a projective curve C in P" + string(n) + ":");
    270    if ( si == 0 )
    271    {
    272    dbprint(printlevel-voice+2,"//   - i is saturated: YES");
    273    }
     604     }
     605//----- If the ideal i is 0:
     606   if ( size(I) == 0 )
     607     {
     608       dbprint(printlevel-voice+2,
     609               "// The ideal i is (0)!
     610// The depth of S/i is:");
     611       return (d);
     612     }
     613//----- When the ideal i is 0-dimensional:
     614   if ( d == 0 )
     615     {
     616       time=rtimer-time;
     617       // Additional information:
     618       dbprint(printlevel-voice+2,
     619               "// Dimension of S/i : 0 (S/i is Cohen-Macaulay)");
     620       dbprint(printlevel-voice+2,
     621               "// Time for computing the depth: " + string(time) + " sec.");
     622       dbprint(printlevel-voice+2,
     623               "// The depth of S/i is:");
     624       return (0);
     625     }
     626//----- Determine the situation: NT, or NP, or nothing.
     627//----- Choose the method depending on the situation, on the
     628//----- characteristic of the ground field, and on the option argument
     629//----- in order to get the mon. ideal of nested type associated to i
     630   if ( e == 1 )
     631     {
     632       ch=0;
     633     }
     634   NPtest=is_NP(I);
     635   if ( NPtest == 1 )
     636     {
     637       nesttest=is_nested(I);
     638     }
     639   if ( ch != 0 )
     640     {
     641       if ( NPtest == 0 )
     642         {
     643           N=d*n-d*(d-1)/2;
     644           s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
     645           execute(s);
     646           ideal chcoord,m,i,I;
     647           poly P;
     648           map phi;
     649           i=imap(r1,i);
     650           chcoord=select1(maxideal(1),1,(n-d+1));
     651           acc=0;
     652           for ( ii = 1; ii<=d; ii++ )
     653             {
     654               matrix trex[1][n-d+ii+1]=t((1+acc)..(n-d+ii+acc)),1;
     655               m=select1(maxideal(1),1,(n-d+1+ii));
     656               for ( jj = 1; jj<=n-d+ii+1; jj++ )
     657                 {
     658                   P=P+trex[1,jj]*m[jj];
     659                 }
     660               chcoord[n-d+1+ii]=P;
     661               P=0;
     662               acc=acc+n-d+ii;
     663               kill trex;
     664             }
     665               phi=rtr,chcoord;
     666               I=simplify(lead(std(phi(i))),1);
     667               setring r1;
     668               I=imap(rtr,I);
     669               attrib(I,"isSB",1);
     670         }
     671       else
     672         {
     673           if ( nesttest == 0 )
     674             {
     675               N=d*(d-1)/2;
     676               s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
     677               execute(s);
     678               ideal chcoord,m,i,I;
     679               poly P;
     680               map phi;
     681               i=imap(r1,i);
     682               chcoord=select1(maxideal(1),1,(n-d+2));
     683               acc=0;
     684               for ( ii = 1; ii<=d-1; ii++ )
     685                 {
     686                   matrix trex[1][ii+1]=t((1+acc)..(ii+acc)),1;
     687                   m=select1(maxideal(1),(n-d+2),(n-d+2+ii));
     688                   for ( jj = 1; jj<=ii+1; jj++ )
     689                     {
     690                       P=P+trex[1,jj]*m[jj];
     691                     }
     692                   chcoord[n-d+2+ii]=P;
     693                   P=0;
     694                   acc=acc+ii;
     695                   kill trex;
     696                 }
     697               phi=rtr,chcoord;
     698               I=simplify(lead(std(phi(i))),1);
     699               setring r1;
     700               I=imap(rtr,I);
     701               attrib(I,"isSB",1);
     702             }
     703         }
     704     }
    274705   else
    275    {
    276    dbprint(printlevel-voice+2,"//   - i is saturated: NO");
    277    }
    278    dbprint(printlevel-voice+2,
    279 "//   - C is arithm. Cohen-Macaulay: YES");
    280    dbprint(printlevel-voice+2,
    281 "//   - reg(C) attained at the last step of a m.g.f.r. of i-sat: YES");
    282    dbprint(printlevel-voice+2,
    283 "//   - regularity of the Hilbert function of S/i-sat: " + string(H-2));
    284    dbprint(printlevel-voice+2,
    285 "//   - time for computing reg(C): " + string(time) + " sec.");
    286    dbprint(printlevel-voice+2,
    287 "// Castelnuovo-Mumford regularity of C:");
    288    return(H);
    289    }
    290 //----- Not Cohen-Macaulay case:
    291 //----- First, determine the associated monomial ideal.
    292    l=sat(I,maxideal(1));
    293    ts=l[2];
    294    if (ts != 0)
    295    {
    296      if (e != 0)
    297      {
    298      s= "ring r4 = (char(r0),a),x(0..n),dp;";
    299      execute(s);
    300      ideal i,I,j,m,K;
    301      i=imap(r1,i);
    302      m=nselect(maxideal(1),n+1);
    303      m[n+1]=a*x(n-1);
    304      map phi=r4,m;
    305      j=phi(i);
    306      I=lead(std(j));
    307      K=normalize(I);
    308      setring r1;
    309      I=imap(r4,K);
    310      }
    311    else
    312      {
    313      rant=size(select(I,n+1));
    314      while (rant != 0)
    315        {
    316        ran=random(-1000,1000);
    317        m=nselect(maxideal(1),n+1);
    318        m[n+1]=ran*x(n-1);
    319        map phi=r1,m;
    320        j=phi(i);
    321        I=lead(std(j));
    322        rant=size(select(I,n+1));
    323        dbprint(printlevel-voice+2,"// (random choice of an element of K)");
    324        }
    325      }
    326    }
    327    else
    328    {
    329    I=subst(I,x(n),x(n-1));
    330    }
    331    I0 = subst(I,x(n-1),0);    // Generators of I which are x(n-1)-free;
    332    K= select(I,n);            // The other generators;
    333 //--------------- Computation of H=H(E) --------------------
    334    s="ring r2 = ",charstr(r0),",x(0..n-2),dp;";
    335    execute(s);
    336    ideal I0,qq,ki,mov;
    337    I0 = imap(r1,I0);
    338    qq=quotient(I0,maxideal(1));
    339    H=maxdeg1(qq)+1;
    340 //------------ Computation of HR=H(R)+1 -----------------
    341 //First, order elements in K
    342    s="ring r3 = ",charstr(r0),",(x(n-1),x(n),x(0..n-2)),lp;";
    343    execute(s);
    344    ideal K,KK,ki;
    345    K=imap(r1,K);
    346    KK=sort(K)[1];
    347 //The first step is different to avoid to compute quotient(I0,max) twice
    348    ki=subst(KK[1],x(n-1),1);
    349    dd=leadexp(KK[1])[1];
    350    setring r2;
    351    ki=imap(r3,ki);
    352    attrib(ki,"isSB",1);
    353    for    (jj=1; jj<= size(qq); jj++)
    354        {
    355        if ( reduce(qq[jj],ki)== 0 )
    356           { hm=deg(qq[jj]);
    357             if ( hm > hh)
    358             { hh = hm; }
    359           }
    360        }
    361    HR=hh+dd;
    362 //If K has more than 1 element, recursive steps to compute HR
    363    setring r1;
    364    if (size(K) != 1)
    365    {
    366    setring r2;
    367    mov=I0+ki;
    368    setring r3;
    369    for (ii=2; ii<= size(KK); ii++)
    370      {
    371        ki=subst(KK[ii],x(n-1),1);
    372        dd=leadexp(KK[ii])[1];
    373        setring r2;
    374        qq=quotient(mov,maxideal(1));
    375        ki=imap(r3,ki);
    376        attrib(ki,"isSB",1);
    377        hh=0;
    378        for    (jj=1; jj<= size(qq); jj++)
    379        {
    380            if ( reduce(qq[jj],ki)==0 )
    381               { hm=deg(qq[jj]);
    382                 if ( hm > hh)
    383                 { hh = hm; }
    384               }
    385        }
    386        hh=hh+dd;
    387        if ( hh > HR )
    388        { HR=hh; }
    389        mov=mov+ki;
    390        setring r3;
    391      }
    392    }
    393 //Now one has HR=H(R)+1 and H=H(E) and one can conclude:
    394   time=rtimer-time;
    395   if( HR > H )
    396   {
    397 // Additional information:
    398    dbprint(printlevel-voice+2,
    399 "// Ideal i of S defining a projective curve C in P" + string(n) + ":");
    400    if ( si == 0 )
    401    {
    402    dbprint(printlevel-voice+2,"//   - i is saturated: YES");
    403    }
    404    else
    405    {
    406    dbprint(printlevel-voice+2,"//   - i is saturated: NO");
    407    }
    408    dbprint(printlevel-voice+2,
    409 "//   - C is arithm. Cohen-Macaulay: NO");
    410    dbprint(printlevel-voice+2,
    411 "//   - reg(C) attained at the last step of a m.g.f.r. of i-sat: YES");
    412    dbprint(printlevel-voice+2,
    413 "//   - regularity of the Hilbert function of S/i-sat: " + string(HR-1));
    414    dbprint(printlevel-voice+2,
    415 "//   - time for computing reg(C): "+ string(time) + " sec.");
    416    dbprint(printlevel-voice+2,
    417 "// Castelnuovo-Mumford regularity of C:");
    418    return(HR);
    419   }
    420   if( HR < H )
    421   {
    422 // Additional information:
    423    dbprint(printlevel-voice+2,
    424 "// Ideal i of S defining a projective curve C in P" + string(n) + ":");
    425    if ( si == 0 )
    426    {
    427    dbprint(printlevel-voice+2,"//   - i is saturated: YES");
    428    }
    429    else
    430    {
    431    dbprint(printlevel-voice+2,"//   - i is saturated: NO");
    432    }
    433    dbprint(printlevel-voice+2,
    434 "//   - C is arithm. Cohen-Macaulay: NO");
    435    dbprint(printlevel-voice+2,
    436 "//   - reg(C) attained at the last step of a m.g.f.r. of i-sat: NO");
    437    dbprint(printlevel-voice+2,
    438 "//   - reg(C) attained at the second last step of a m.g.f.r. of i-sat: YES");
    439    dbprint(printlevel-voice+2,
    440 "//   - regularity of the Hilbert function of S/i-sat: strictly smaller than "
    441            + string(H-1));
    442    dbprint(printlevel-voice+2,
    443 "//   - time for computing reg(C): " + string(time) + " sec.");
    444    dbprint(printlevel-voice+2,
    445 "// Castelnuovo-Mumford regularity of C:");
    446    return(H);
    447   }
    448   if( HR == H )
    449   {
    450 // Additional information:
    451    dbprint(printlevel-voice+2,
    452 "// Ideal i of S defining a projective curve C in P" + string(n) + ":");
    453    if ( si == 0 )
    454    {
    455    dbprint(printlevel-voice+2,"//   - i is saturated: YES");
    456    }
    457    else
    458    {
    459    dbprint(printlevel-voice+2,"//   - i is saturated: NO");
    460    }
    461    dbprint(printlevel-voice+2,
    462 "//   - C is arithm. Cohen-Macaulay: NO");
    463    dbprint(printlevel-voice+2,
    464 "//   - reg(C) attained at the last step of a m.g.f.r. of i-sat: YES");
    465    dbprint(printlevel-voice+2,
    466 "//   - reg(C) attained at the second last step of a m.g.f.r. of i-sat: YES");
    467    dbprint(printlevel-voice+2,
    468 "//   - regularity of the Hilbert function of S/i-sat: " + string(HR-1));
    469    dbprint(printlevel-voice+2,
    470 "//   - time for computing reg(C): " + string(time) + " sec.");
    471    dbprint(printlevel-voice+2,
    472 "// Castelnuovo-Mumford regularity of C:");
    473    return(HR);
    474   }
     706     {
     707       if ( NPtest == 0 )
     708         {
     709           while ( nl < 30 )
     710             {
     711               chcoord=select1(maxideal(1),1,(n-d+1));
     712               nl=nl+1;
     713               for ( ii = 1; ii<=d; ii++ )
     714                 {
     715                   ran=random(100,1,n-d+ii);
     716                   ran=intmat(ran,1,n-d+ii+1);
     717                   ran[1,n-d+ii+1]=1;
     718                   m=select1(maxideal(1),1,(n-d+1+ii));
     719                   for ( jj = 1; jj<=n-d+ii+1; jj++ )
     720                     {
     721                       P=P+ran[1,jj]*m[jj];
     722                     }
     723                   chcoord[n-d+1+ii]=P;
     724                   P=0;
     725                 }
     726               phi=r1,chcoord;
     727               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
     728               I=simplify(lead(std(phi(i))),1);
     729               attrib(I,"isSB",1);
     730               NPtest=is_NP(I);
     731               if ( NPtest == 1 )
     732                 {
     733                   break;
     734                 }
     735             }
     736           if ( NPtest == 0 )
     737             {
     738       "// WARNING from proc depthIdeal from lib mregular.lib:
     739// The procedure has entered in 30 loops and could not put the variables
     740// in Noether position: in your example the method using random changes
     741// of coordinates may enter an infinite loop when the field is finite.
     742// Try removing this optional argument.";
     743       return (-1);
     744             }
     745           i=phi(i);
     746           nesttest=is_nested(I);
     747         }
     748       if ( nesttest == 0 )
     749         {
     750           while ( nl < 30 )
     751             {
     752               chcoord=select1(maxideal(1),1,(n-d+2));
     753               nl=nl+1;
     754               for ( ii = 1; ii<=d-1; ii++ )
     755                 {
     756                   ran=random(100,1,ii);
     757                   ran=intmat(ran,1,ii+1);
     758                   ran[1,ii+1]=1;
     759                   m=select1(maxideal(1),(n-d+2),(n-d+2+ii));
     760                   for ( jj = 1; jj<=ii+1; jj++ )
     761                     {
     762                       P=P+ran[1,jj]*m[jj];
     763                     }
     764                   chcoord[n-d+2+ii]=P;
     765                   P=0;
     766                 }
     767               phi=r1,chcoord;
     768               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
     769               I=simplify(lead(std(phi(i))),1);
     770               attrib(I,"isSB",1);
     771               nesttest=is_nested(I);
     772               if ( nesttest == 1 )
     773                 {
     774                   break;
     775                 }
     776             }
     777           if ( nesttest == 0 )
     778             {
     779       "// WARNING from proc depthIdeal from lib mregular.lib:
     780// The procedure has entered in 30 loops and could not find a monomial
     781// ideal of nested type with the same depth as your ideal: in your
     782// example the method using random changes of coordinates may enter an
     783// infinite loop when the field is finite.
     784// Try removing this optional argument.";
     785       return (-1);
     786             }
     787         }
     788     }
     789//
     790// At this stage, we have obtained a monomial ideal I of nested type
     791// such that depth(S/i)=depth(S/I). We now compute depth(I).
     792//
     793//----- When S/i is Cohen-Macaulay:
     794   for ( ii = n-d+2; ii <= n+1; ii++ )
     795     {
     796       K=K+select(I,ii);
     797     }
     798   if ( size(K) == 0 )
     799     {
     800       time=rtimer-time;
     801       // Additional information:
     802       dbprint(printlevel-voice+2,
     803               "// Dimension of S/i: "+string(d)+" (S/i is Cohen-Macaulay)");
     804       dbprint(printlevel-voice+2,
     805               "// Time for computing depth: " + string(time) + " sec.");
     806       dbprint(printlevel-voice+2,
     807               "// The depth of S/i is:");
     808       return(d);
     809     }
     810//----- When d=1 (and S/i is not Cohen-Macaulay) ==> depth =0:
     811   if ( d == 1 )
     812     {
     813       time=rtimer-time;
     814       // Additional information:
     815       dbprint(printlevel-voice+2,
     816               "// Dimension of S/i: 1");
     817       dbprint(printlevel-voice+2,
     818               "// Time for computing depth: "+ string(time) + " sec.");
     819       dbprint(printlevel-voice+2,
     820               "// The depth of S/i is:");
     821       return(0);
     822
     823     }
     824//----- Now d>1 and S/i is not Cohen-Macaulay:
     825//
     826//----- First, determine the last variable really occuring
     827       lastv=n-d;
     828       h=n;
     829       while ( lastv == n-d and h > n-d )
     830         {
     831           K=select(I,h+1);
     832           if ( size(K) == 0 )
     833             {
     834               h=h-1;
     835             }
     836           else
     837             {
     838               lastv=h;
     839             }
     840         }
     841//----- and compute the depth:
     842       time=rtimer-time;
     843       // Additional information:
     844       dbprint(printlevel-voice+2,
     845               "// Dimension of S/i: "+string(d));
     846       dbprint(printlevel-voice+2,
     847               "// Time for computing depth: "+ string(time) + " sec.");
     848       dbprint(printlevel-voice+2,
     849               "// The depth of S/i is:");
     850       return(n-lastv);
    475851}
    476852example
    477853{ "EXAMPLE:"; echo = 2;
    478    ring s = 0,(x,y,z,t),dp;
    479 // 1st example is Ex.2.5 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
     854   ring r=0,(x,y,z,t,w),dp;
     855   ideal i=y2t,x2y-x2z+yt2,x2y2,xyztw,x3z2,y5+xz3w-x2zw2,x7-yt2w4;
     856   depthIdeal(i);
     857   depthIdeal(lead(std(i)));
     858// Additional information is displayed if you change printlevel (=1);
     859}
     860////////////////////////////////////////////////////////////////////////////////
     861/*
     862Out-commented examples:
     863   ring s=0,x(0..5),dp;
     864   ideal i=x(2)^2-x(4)*x(5),x(1)*x(2)-x(0)*x(5),x(0)*x(2)-x(1)*x(4),
     865           x(1)^2-x(3)*x(5),x(0)*x(1)-x(2)*x(3),x(0)^2-x(3)*x(4);
     866   depthIdeal(i);
     867   // Our procedure works when a min. graded free resol. can
     868   // not be computed. In this easy example, depth can also
     869   // be obtained using a m.g.f.r. (Auslander-Buchsbaum formula):
     870   nvars(s)-ncols(betti(mres(i,0)))+1;
     871   ring r1=0,(x,y,z,t),dp;
     872// Ex.2.5 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
    480873   ideal i  = x17y14-y31, x20y13, x60-y36z24-x20z20t20;
    481    reg_curve(i);
    482 // 2nd example is Ex.2.9 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
     874   depthIdeal(i);
     875// Ex.2.9 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
    483876   int k=43;
    484877   ideal j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
    485    reg_curve(j);
    486 // Additional information can be obtained as follows:
    487    printlevel = 1;
    488    reg_curve(j);
    489 }
    490 
    491 ///////////////////////////////////////////////////////////////////////////////
    492 /*
    493 Out-commented examples:
    494 //
    495 // More examples are obtained changing the value of k in the previous example:
    496    k=14;
    497    j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
    498    reg_curve(j);
    499    k=22;
    500    j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
    501    reg_curve(j);
    502    k=315;
    503    j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41*z^k-y40*z^(k+1);
    504    reg_curve(j);
    505 // If the ideal does not define a curve, one gets an error message:
    506    ring s1=0,(a,b,c,d,x(0..9)),dp;
     878   depthIdeal(j);
     879// Example in Rk.2.10 in [Bermejo-Gimenez], ProcAMS 128(5):
     880   ideal h=x2-3xy+5xt,xy-3y2+5yt,xz-3yz,2xt-yt,y2-yz-2yt;
     881   depthIdeal(h);
     882// The initial ideal is not saturated
     883   depthIdeal(lead(std(h)));
     884// More examples:
     885   i=y4-t3z, x3t-y2z2, x3y2-t2z3, x6-tz5;
     886   depthIdeal(i);
     887//
     888   depthIdeal(maxideal(4));
     889//
     890   ring r2=0,(x,y,z,t,w),dp;
     891   ideal i = xy-zw,x3-yw2,x2z-y2w,y3-xz2,-y2z3+xw4+tw4+w5,-yz4+x2w3+xtw3+xw4,
     892            -z5+x2tw2+x2w3+yw4;
     893   depthIdeal(i);
     894//
     895   ring r3=0,(x,y,z,t,w,u),dp;
     896   ideal i=imap(r2,i);
     897   depthIdeal(i);
     898// Next example is the defining ideal of the 2nd. Veronesean of P3, a variety
     899// in P8 which is arithmetically Cohen-Macaulay:
     900   ring r4=0,(a,b,c,d,x(0..9)),dp;
    507901   ideal i= x(0)-ab,x(1)-ac,x(2)-ad,x(3)-bc,x(4)-bd,x(5)-cd,
    508902            x(6)-a2,x(7)-b2,x(8)-c2,x(9)-d2;
    509903   ideal ei=eliminate(i,abcd);
    510    ring s2=0,x(0..9),dp;
    511    ideal i=imap(s1,ei);
    512    reg_curve(i);
    513 // Since this example is Cohen-Macaulay, use reg_CM instead of reg_curve.
    514 //
    515 // Be carefull!: the polynomial ring in the last variables MUST be a Noether
    516 // Normalization of basering/ideal:
    517    ring s3=0,(t,x,y,z),dp;
    518    ideal j=imap(s,j);
    519    reg_curve(j);
    520 //
    521 // The following is example in Rk.2.10 in [Bermejo-Gimenez], ProcAMS 128(5):
    522    setring s;
    523    ideal h=x2-3xy+5xt,xy-3y2+5yt,xz-3yz,2xt-yt,y2-yz-2yt;
    524    reg_curve(h);
    525 // The initial ideal is not saturated but the regularity of the subscheme
    526 // it defines is the regularity of the saturation and can be computed:
    527    reg_curve(lead(std(h)));
    528 //
    529 // Here is an example where the computation of a m.g.f.r. of I costs a lot:
    530    ring s4=0,(x,y,z,t,u,a,b),dp;
     904   ring r5=0,x(0..9),dp;
     905   ideal i=imap(r4,ei);
     906   depthIdeal(i);
     907// Here is an example where the computation of a m.g.f.r. of I costs:
     908   ring r8=0,(x,y,z,t,u,a,b),dp;
    531909   ideal i=u-b40,t-a40,x-a23b17,y-a22b18+ab39,z-a25b15;
    532910   ideal ei=eliminate(i,ab); // It takes a few seconds to compute the ideal
    533    ring s5=0,(x,y,z,t,u),dp;
    534    ideal i=imap(s4,ei);
    535    reg_curve(i);   // This is very fast.
     911   ring r9=0,(x,y,z,t,u),dp;
     912   ideal i=imap(r8,ei);
     913   depthIdeal(i);   // This is very fast.
    536914// Now you can use mres(i,0) to compute a m.g.f.r. of the ideal!
    537915//
    538 // The computation of the m.g.f.r. of the following example did not succeed
    539 // using the command mres:
    540    ring s6=0,(x(0..8),s,t),dp;
     916// Another one:
     917   ring r10=0,(x(0..8),s,t),dp;
    541918   ideal i=x(0)-st24,x(1)-s2t23,x(2)-s3t22,x(3)-s9t16,x(4)-s11t14,x(5)-s18t7,
    542919           x(6)-s24t,x(7)-t25,x(8)-s25;
    543920   ideal ei=eliminate(i,st);
    544    ring s7=0,x(0..8),dp;
    545    ideal i=imap(s6,ei);
    546    reg_curve(i);
     921   ring r11=0,x(0..8),dp;
     922   ideal i=imap(r10,ei);
     923   depthIdeal(i);
     924// More examples where not even sres works:
     925// Be careful: elimination takes some time here, but it succeeds!
     926ring r12=0,(s,t,u,x(0..14)),dp;
     927ideal i=x(0)-st6u8,x(1)-s5t3u7,x(2)-t11u4,x(3)-s9t4u2,x(4)-s2t7u6,x(5)-s7t7u,
     928        x(6)-s10t5,x(7)-s4t6u5,x(8)-s13tu,x(9)-s14u,x(10)-st2u12,x(11)-s3t9u3,
     929        x(12)-s15,x(13)-t15,x(14)-u15;
     930ideal ei=eliminate(i,stu);
     931size(ei);
     932ring r13=0,x(0..14),dp;
     933ideal i=imap(r12,ei);
     934size(i);
     935depthIdeal(i);
    547936//
    548937*/
    549 //////////////////////////////////////////////////////////////////////////////
     938///////////////////////////////////////////////////////////////////////////////
     939///////////////////////////////////////////////////////////////////////////////
     940///////////////////////////////////////////////////////////////////////////////
    550941
    551 proc reg_moncurve (list #)
     942proc satiety (ideal i, list #)
    552943"
    553 USAGE:   reg_moncurve (a0,...,an) ; ai integers with a0=0 < a1 < ... < an=:d
    554 RETURN:  an integer, the Castelnuovo-Mumford regularity of the projective
    555          monomial curve C in Pn parametrically defined by:
    556                x(0)=t^d , x(1)=s^(a1)t^(d-a1), ... , x(n)=s^d.
    557 ASSUME:  a0=0 < a1 < ... < an are integers and the base field is infinite.
    558 NOTE:    The defining ideal I(C) in S is determined using elimination.
    559          The procedure reg_curve is improved in this case since one
    560          knows beforehand that the dimension is 2, that the variables are
    561          in Noether position, that I(C) is prime.
    562          If printlevel > 0 (default = 0) additional information is displayed.
    563          In particular, says if C is arithmetically Cohen-Macaulay or not,
    564          determines in which step of a minimal graded free resolution of I(C)
    565          the regularity is attained, and sometimes gives the value of the
    566          regularity of the Hilbert function of S/I(C) (otherwise, an upper
    567          bound is given).
    568 EXAMPLE: example reg_moncurve; shows some examples
     944USAGE:   satiety (i[,e]); i ideal, e integer
     945RETURN:  an integer, the satiety of i.
     946         (returns -1 if i is not homogeneous)
     947ASSUME:  i is a homogeneous ideal of the basering S=K[x(0)..x(n)].
     948         e=0:  (default)
     949               The satiety is computed determining the fresh elements in the
     950               socle of i. It works over arbitrary fields.
     951         e=1:  Makes random changes of coordinates to find a monomial ideal
     952               with same satiety. It works over infinite fields only. If K
     953               is finite, it works if it terminates, but may result in an
     954               infinite loop. After 30 loops, a warning message is displayed
     955               and -1 is returned.
     956THEORY:  The satiety, or saturation index, of a homogeneous ideal i is the
     957         least integer s such that, for all d>=s, the degree d part of the
     958         ideals i and isat=sat(i,maxideal(1))[1] coincide.
     959NOTE:    If printlevel > 0 (default = 0), dim(S/i) is also displayed.
     960EXAMPLE: example satiety; shows some examples
    569961"
    570962{
    571963//--------------------------- initialisation ---------------------------------
    572    int ii,jj,H,HR,dd,h,hh,hm,sK,time,ttime;
     964   int e,ii,jj,h,d,time,lastv,nesttest,NPtest,nl,sat;
     965   intmat ran;
     966   def r0 = basering;
     967   int n = nvars(r0)-1;
     968   if ( size(#) > 0 )
     969   {
     970      e = #[1];
     971   }
     972   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
     973   execute(s);
     974   ideal i,sbi,I,K,chcoord,m,KK;
     975   poly P;
     976   map phi;
     977   i = fetch(r0,i);
     978   time=rtimer;
     979   sbi=std(i);
     980//----- Check ideal homogeneous
     981   if ( homog(sbi) == 0 )
     982   {
     983       "// WARNING from proc satiety from lib mregular.lib:
     984// The ideal is not homogeneous!";
     985       return (-1);
     986   }
     987   I=simplify(lead(sbi),1);
     988   attrib(I,"isSB",1);
     989   d=dim(I);
     990//----- If the ideal i is not proper:
     991   if ( d == -1 )
     992     {
     993       dbprint(printlevel-voice+2,
     994               "// The ideal i is (1)!
     995// Its satiety is:");
     996       return (0);
     997     }
     998//----- If the ideal i is 0:
     999   if ( size(I) == 0 )
     1000     {
     1001       dbprint(printlevel-voice+2,
     1002               "// The ideal i is (0)!
     1003// Its satiety is:");
     1004       return (0);
     1005     }
     1006//----- When the ideal i is 0-dimensional:
     1007   if ( d == 0 )
     1008     {
     1009       sat=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
     1010       time=rtimer-time;
     1011       // Additional information:
     1012       dbprint(printlevel-voice+2,
     1013               "// Dimension of S/i: 0");
     1014       dbprint(printlevel-voice+2,
     1015               "// Time for computing the satiety: " + string(time) + " sec.");
     1016       dbprint(printlevel-voice+2,
     1017               "// The satiety of i is:");
     1018       return (sat);
     1019     }
     1020//----- When one has option e=1:
     1021//
     1022//----- Determine the situation: NT, or NP, or nothing.
     1023//----- Choose the method depending on the situation in order to
     1024//----- get the mon. ideal of nested type associated to i
     1025   if ( e == 1 )
     1026     {
     1027       NPtest=is_NP(I);
     1028       if ( NPtest == 0 )
     1029         {
     1030           while ( nl < 30 )
     1031             {
     1032               chcoord=select1(maxideal(1),1,(n-d+1));
     1033               nl=nl+1;
     1034               for ( ii = 1; ii<=d; ii++ )
     1035                 {
     1036                   ran=random(100,1,n-d+ii);
     1037                   ran=intmat(ran,1,n-d+ii+1);
     1038                   ran[1,n-d+ii+1]=1;
     1039                   m=select1(maxideal(1),1,(n-d+1+ii));
     1040                   for ( jj = 1; jj<=n-d+ii+1; jj++ )
     1041                     {
     1042                       P=P+ran[1,jj]*m[jj];
     1043                     }
     1044                   chcoord[n-d+1+ii]=P;
     1045                   P=0;
     1046                 }
     1047               phi=r1,chcoord;
     1048               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
     1049               I=simplify(lead(std(phi(i))),1);
     1050               attrib(I,"isSB",1);
     1051               NPtest=is_NP(I);
     1052               if ( NPtest == 1 )
     1053                 {
     1054                   break;
     1055                 }
     1056             }
     1057           if ( NPtest == 0 )
     1058             {
     1059       "// WARNING from proc satiety from lib mregular.lib:
     1060// The procedure has entered in 30 loops and could not put the variables
     1061// in Noether position: in your example the method using random changes
     1062// of coordinates may enter an infinite loop when the field is finite.
     1063// Try removing the optional argument.";
     1064       return (-1);
     1065             }
     1066           i=phi(i);
     1067         }
     1068       nesttest=is_nested(I);
     1069       if ( nesttest == 0 )
     1070         {
     1071           while ( nl < 30 )
     1072             {
     1073               chcoord=select1(maxideal(1),1,(n-d+2));
     1074               nl=nl+1;
     1075               for ( ii = 1; ii<=d-1; ii++ )
     1076                 {
     1077                   ran=random(100,1,ii);
     1078                   ran=intmat(ran,1,ii+1);
     1079                   ran[1,ii+1]=1;
     1080                   m=select1(maxideal(1),(n-d+2),(n-d+2+ii));
     1081                   for ( jj = 1; jj<=ii+1; jj++ )
     1082                     {
     1083                       P=P+ran[1,jj]*m[jj];
     1084                     }
     1085                   chcoord[n-d+2+ii]=P;
     1086                   P=0;
     1087                 }
     1088               phi=r1,chcoord;
     1089               dbprint(printlevel-voice+2,"// (1 random change of coord.)");
     1090               I=simplify(lead(std(phi(i))),1);
     1091               attrib(I,"isSB",1);
     1092               nesttest=is_nested(I);
     1093               if ( nesttest == 1 )
     1094                 {
     1095                   break;
     1096                 }
     1097             }
     1098           if ( nesttest == 0 )
     1099             {
     1100       "// WARNING from proc satiety from lib mregular.lib:
     1101// The procedure has entered in 30 loops and could not find a monomial
     1102// ideal of nested type with the same satiety as your ideal: in your
     1103// example the method using random changes of coordinates may enter an
     1104// infinite loop when the field is finite.
     1105// Try removing the optional argument.";
     1106       return (-1);
     1107             }
     1108         }
     1109//
     1110// At this stage, we have obtained a monomial ideal I of nested type
     1111// such that depth(S/i)=depth(S/I). We now compute depth(I).
     1112//
     1113//----- When S/i is Cohen-Macaulay:
     1114//
     1115       for ( ii = n-d+2; ii <= n+1; ii++ )
     1116         {
     1117           K=K+select(I,ii);
     1118         }
     1119       if ( size(K) == 0 )
     1120         {
     1121           time=rtimer-time;
     1122           // Additional information:
     1123           dbprint(printlevel-voice+2,
     1124                   "// Dimension of S/i: "+string(d));
     1125           dbprint(printlevel-voice+2,
     1126                   "// Time for computing satiety: " + string(time) + " sec.");
     1127           dbprint(printlevel-voice+2,
     1128                   "// The satiety of i is:");
     1129           return(0);
     1130         }
     1131//----- When d=1 (and S/i is not Cohen-Macaulay) ==> depth =0:
     1132       if ( d == 1 )
     1133         {
     1134           KK=simplify(reduce(quotient(I,maxideal(1)),I),2);
     1135           sat=maxdeg1(KK)+1;
     1136           time=rtimer-time;
     1137           // Additional information:
     1138           dbprint(printlevel-voice+2,
     1139                   "// Dimension of S/i: 1");
     1140           dbprint(printlevel-voice+2,
     1141                   "// Time for computing satiety: "+ string(time) + " sec.");
     1142           dbprint(printlevel-voice+2,
     1143                   "// The satiety of i is:");
     1144           return(sat);
     1145         }
     1146//----- Now d>1 and S/i is not Cohen-Macaulay:
     1147//
     1148//----- First, determine the last variable really occuring
     1149       lastv=n-d;
     1150       h=n;
     1151       while ( lastv == n-d and h > n-d )
     1152         {
     1153           K=select(I,h+1);
     1154           if ( size(K) == 0 )
     1155             {
     1156               h=h-1;
     1157             }
     1158           else
     1159             {
     1160               lastv=h;
     1161             }
     1162         }
     1163//----- and compute the satiety:
     1164       sat=0;
     1165       if ( lastv == n )
     1166         {
     1167           KK=simplify(reduce(quotient(I,maxideal(1)),I),2);
     1168           sat=maxdeg1(KK)+1;
     1169         }
     1170       time=rtimer-time;
     1171       // Additional information:
     1172       dbprint(printlevel-voice+2,
     1173               "// Dimension of S/i: "+string(d));
     1174       dbprint(printlevel-voice+2,
     1175               "// Time for computing satiety: "+ string(time) + " sec.");
     1176       dbprint(printlevel-voice+2,
     1177               "// The satiety of i is:");
     1178       return(sat);
     1179     }
     1180//---- If no option: direct computation
     1181   sat=maxdeg1(reduce(quotient(i,maxideal(1)),sbi))+1;
     1182   time=rtimer-time;
     1183       // Additional information:
     1184   dbprint(printlevel-voice+2,
     1185           "// Dimension of S/i: "+string(d)+";");
     1186   dbprint(printlevel-voice+2,
     1187           "// Time for computing satiety: "+ string(time) + " sec.");
     1188   dbprint(printlevel-voice+2,
     1189           "// The satiety of i is:");
     1190   return(sat);
     1191}
     1192example
     1193{ "EXAMPLE:"; echo = 2;
     1194   ring r=0,(x,y,z,t,w),dp;
     1195   ideal i=y2t,x2y-x2z+yt2,x2y2,xyztw,x3z2,y5+xz3w-x2zw2,x7-yt2w4;
     1196   satiety(i);
     1197   ideal I=lead(std(i));
     1198   satiety(I);   // First  method: direct computation
     1199   satiety(I,1); // Second method: doing changes of coordinates
     1200// Additional information is displayed if you change printlevel (=1);
     1201}
     1202////////////////////////////////////////////////////////////////////////////////
     1203/*
     1204Out-commented examples:
     1205   ring s1=0,(x,y,z,t),dp;
     1206   ideal I=zt3,z2t2,yz2t,xz2t,xy2t,x3y;
     1207   satiety(I);
     1208   satiety(I,1);
     1209// Another example:
     1210   ring s2=0,(z,y,x),dp;
     1211   ideal I=z38,z26y2,z14y4,z12x,z10x5,z8x9,z6x16,z4x23,z2y6,y32;
     1212   satiety(I);
     1213   satiety(I,1);
     1214// One more:
     1215   ring s3=0,(s,t,u,x(0..8)),dp;
     1216   ideal i=x(0)-st6u8,x(1)-s5t3u7,x(2)-t11u4,x(3)-s9t4u2,
     1217           x(4)-s2t7u6,x(5)-s7t7u,x(6)-s15,x(7)-t15,x(8)-u15;
     1218   ideal ei=eliminate(i,stu);
     1219   size(ei);
     1220   ring s4=0,x(0..8),dp;
     1221   ideal i=imap(s3,ei);
     1222   ideal m=maxideal(1);
     1223   m[8]=m[8]+m[7];
     1224   map phi=m;
     1225   ideal phii=phi(i);
     1226   ideal nI=lead(std(phii));
     1227   ring s5=0,x(0..7),dp;
     1228   ideal nI=imap(s4,nI);
     1229   satiety(nI);
     1230   satiety(nI,1);
     1231   ideal I1=subst(nI,x(7),1);
     1232   ring s6=0,x(0..6),dp;
     1233   ideal I1=imap(s5,I1);
     1234   satiety(I1);
     1235   satiety(I1,1);
     1236   ideal I2=subst(I1,x(6),1);
     1237   ring s7=0,x(0..5),dp;
     1238   ideal I2=imap(s6,I2);
     1239   satiety(I2);
     1240   satiety(I2,1);
     1241//
     1242*/
     1243///////////////////////////////////////////////////////////////////////////////
     1244///////////////////////////////////////////////////////////////////////////////
     1245///////////////////////////////////////////////////////////////////////////////
     1246
     1247proc regMonCurve (list #)
     1248"
     1249USAGE:   regMonCurve (a0,...,an) ; ai integers with a0=0 < a1 < ... < an=:d
     1250RETURN:  an integer, the Castelnuovo-Mumford regularity of the projective
     1251         monomial curve C in Pn(K) parametrically defined by
     1252              x(0) = t^d , x(1) = s^(a1)t^(d-a1) , ..... , x(n) = s^d
     1253         where K is the field of complex numbers.
     1254         (returns -1 if a0=0 < a1 < ... < an is not satisfied)
     1255ASSUME:  a0=0 < a1 < ... < an are integers.
     1256NOTES:   1. The defining ideal of the curve C, I in S=K[x(0),...,x(n)], is 
     1257            determined by elimination.
     1258         2. The procedure regIdeal has been improved in this case since one
     1259            knows beforehand that the monomial ideal J=lead(std(I)) is of
     1260            nested type if the monomial ordering is dp, and that
     1261            reg(C)=reg(J) (see preprint 'Saturation and Castelnuovo-Mumford
     1262            regularity' by Bermejo-Gimenez, 2004).
     1263         3. If printlevel > 0 (default = 0) additional info is displayed:
     1264            - It says whether C is arithmetically Cohen-Macaulay or not.
     1265            - If C is not arith. Cohen-Macaulay, end(H^1(S/I)) is computed
     1266              and an upper bound for the a-invariant of S/I is given.
     1267            - It also determines one step of the minimal graded free
     1268              resolution (m.g.f.r.) of I where the regularity is attained
     1269              and gives the value of the regularity of the Hilbert function
     1270              of S/I when reg(I) is attained at the last step of a m.g.f.r.
     1271EXAMPLE: example regMonCurve; shows some examples
     1272"
     1273{
     1274//--------------------------- initialisation ---------------------------------
     1275   int ii,H,h,hh,time,ttime,firstind,lastind;
    5731276   int n = size(#)-1;
    5741277//------------------  Check assumptions on integers  -------------------------
    5751278   if ( #[1] != 0 )
    576    {"// WARNING from proc reg_moncurve from lib mregular.lib:
     1279   {"// WARNING from proc regMonCurve from lib mregular.lib:
    5771280// USAGE: your input must be a list of integers a0,a1,...,an such that
    5781281// a0=0 < a1 < a2 < ... < an";
    … …  
    5831286      if ( #[ii] >= #[ii+1] )
    5841287      {
    585       "// WARNING from proc reg_moncurve from lib mregular.lib:
     1288      "// WARNING from proc regMonCurve from lib mregular.lib:
    5861289// USAGE: your input must be a list of integers a0,a1,...,an such that
    5871290// a0=0 < a1 < a2 < ... < an";
    … …  
    5891292      }
    5901293   }
    591    ring r4=0,(x(0..n),s,t),dp;
     1294   ring R=0,(x(0..n),s,t),dp;
    5921295   ideal param,m,i;
    5931296   poly f(0..n);
    … …  
    6021305   ttime=rtimer;
    6031306   i=eliminate(i,st);
    604    ring r1=0,(x(1..n),x(0)),dp;
    605    ideal i,I,I0,K;
    606    i=imap(r4,i);
     1307   ring r=0,(x(1..n),x(0)),dp;
     1308   ideal i,I;
     1309   i=imap(R,i);
     1310   I=minbase(lead(std(i)));
     1311   attrib(I,"isSB",1);
    6071312   ttime=rtimer-ttime;
    6081313   time=rtimer;
    609    I=lead(std(i));
     1314   ring nr=0,x(1..n),dp;
     1315   ideal I,K,KK,J;
     1316   I=imap(r,I);
    6101317   attrib(I,"isSB",1);
    611    I0 = subst(I,x(n),0);
    612    K= select(I,n);
    613    sK=size(K);
    614    ring r2=0,x(1..n-1),dp;
    615    ideal I0,qq,ki,mov;
    616    I0 = imap(r1,I0);
    617    qq=quotient(I0,maxideal(1));
    618    H=maxdeg1(qq)+1;
     1318   K=select(I,n);
    6191319//------------------ Cohen-Macaulay case ------------
    620    if (sK == 0)
    621    {
    622    time=rtimer-time;
    623 // Additional information:
    624    dbprint(printlevel-voice+2,
    625 "// Sequence of integers defining a monomial curve C in P" + string(n) + ":");
    626    dbprint(printlevel-voice+2,
    627 "//   - time for computing ideal I(C) of S (elimination): "
    628              + string(ttime) + " sec.");
    629    dbprint(printlevel-voice+2,
    630 "//   - C is arithm. Cohen-Macaulay: YES");
    631    dbprint(printlevel-voice+2,
    632 "//   - reg(C) attained at the last step of a m.g.f.r. of I(C): YES");
    633    dbprint(printlevel-voice+2,
    634 "//   - regularity of the Hilbert function of S/I(C): " + string(H-2));
    635    dbprint(printlevel-voice+2,
    636 "//   - time for computing reg(C): " + string(time) + " sec.");
    637    dbprint(printlevel-voice+2,
    638 "// Castelnuovo-Mumford regularity of C:");
    639    return(H);
    640    }
    641 //------------ non Cohen-Macaulay case : computation of HR=H(R)+1 -------
    642 //First, order elements in K
    643    ring r3=0,(x(n),x(0),x(1..n-1)),lp;
    644    ideal K,KK,ki;
    645    K=imap(r1,K);
    646    KK=sort(K)[1];
    647 //The first step is different to avoid to compute quotient(I0,max) twice
    648    ki=subst(KK[1],x(n),1);
    649    dd=leadexp(KK[1])[1];
    650    setring r2;
    651    ki=imap(r3,ki);
    652    attrib(ki,"isSB",1);
    653    for    (jj=1; jj<= size(qq); jj++)
    654        {
    655        if ( reduce(qq[jj],ki)== 0 )
    656           { hm=deg(qq[jj]);
    657             if ( hm > hh)
    658             { hh = hm; }
    659           }
    660        }
    661    HR=hh+dd;
    662 //If K has more than 1 element, recursive steps to compute HR
    663    setring r1;
    664    if (size(K) != 1)
    665    {
    666    setring r2;
    667    mov=I0+ki;
    668    setring r3;
    669    for (ii=2; ii<= size(KK); ii++)
    670      {
    671        ki=subst(KK[ii],x(n),1);
    672        dd=leadexp(KK[ii])[1];
    673        setring r2;
    674        qq=quotient(mov,maxideal(1));
    675        ki=imap(r3,ki);
    676        attrib(ki,"isSB",1);
    677        hh=0;
    678        for    (jj=1; jj<= size(qq); jj++)
    679        {
    680            if ( reduce(qq[jj],ki)==0 )
    681               { hm=deg(qq[jj]);
    682                 if ( hm > hh)
    683                 { hh = hm; }
    684               }
    685        }
    686        hh=hh+dd;
    687        if ( hh > HR )
    688        { HR=hh; }
    689        mov=mov+ki;
    690        setring r3;
    691      }
    692    }
    693 //Now one has HR=H(R)+1 and H=H(E) and one can conclude:
    694   time=rtimer-time;
    695   if( HR > H )
    696   {
    697 // Additional information:
    698    dbprint(printlevel-voice+2,
    699 "// Sequence of integers defining a monomial curve C in P"+string(n)+":");
    700    dbprint(printlevel-voice+2,
    701 "//   - time for computing ideal I(C) of S (elimination): "
    702             + string(ttime) + " sec.");
    703    dbprint(printlevel-voice+2,
    704 "//   - C is arithm. Cohen-Macaulay: NO");
    705    dbprint(printlevel-voice+2,
    706 "//   - reg(C) attained at the last step of a m.g.f.r. of I(C): YES");
    707    dbprint(printlevel-voice+2,
    708 "//   - regularity of the Hilbert function of S/I(C): " + string(HR-1));
    709    dbprint(printlevel-voice+2,
    710 "//   - time for computing reg(C): "+ string(time) + " sec.");
    711    dbprint(printlevel-voice+2,
    712 "// Castelnuovo-Mumford regularity of C:");
    713    return(HR);
    714   }
    715   if( HR < H )
    716   {
    717 // Additional information:
    718    dbprint(printlevel-voice+2,
    719 "// Sequence of integers defining a monomial curve C in P"+string(n)+":");
    720    dbprint(printlevel-voice+2,
    721 "//   - time for computing ideal I(C) of S (elimination): "
    722         + string(ttime) + " sec.");
    723    dbprint(printlevel-voice+2,
    724 "//   - C is arithm. Cohen-Macaulay: NO");
    725    dbprint(printlevel-voice+2,
    726 "//   - reg(C) attained at the last step of a m.g.f.r. of I(C): NO");
    727    dbprint(printlevel-voice+2,
    728 "//   - reg(C) attained at the second last step of a m.g.f.r. of I(C): YES");
    729    dbprint(printlevel-voice+2,
    730 "//   - regularity of the Hilbert function of S/I(C): striclty smaller than "
    731           + string(H-1));
    732    dbprint(printlevel-voice+2,
    733 "//   - time for computing reg(C): " + string(time) + " sec.");
    734    dbprint(printlevel-voice+2,
    735 "// Castelnuovo-Mumford regularity of C:");
    736    return(H);
    737   }
    738   if( HR == H )
    739   {
    740 // Additional information:
    741    dbprint(printlevel-voice+2,
    742 "// Sequence of integers defining a monomial curve C in P"+string(n)+":");
    743    dbprint(printlevel-voice+2,
    744 "//   - time for computing ideal I(C) of S (elimination): "
    745           + string(ttime) + " sec.");
    746    dbprint(printlevel-voice+2,
    747 "//   - C is arithm. Cohen-Macaulay: NO");
    748    dbprint(printlevel-voice+2,
    749 "//   - reg(C) attained at the last step of a m.g.f.r. of I(C): YES");
    750    dbprint(printlevel-voice+2,
    751 "//   - reg(C) attained at the second last step of a m.g.f.r. of I(C): YES");
    752    dbprint(printlevel-voice+2,
    753 "//   - regularity of the Hilbert function of S/I(C): " + string(HR-1));
    754    dbprint(printlevel-voice+2,
    755 "//   - time for computing reg(C): " + string(time) + " sec.");
    756    dbprint(printlevel-voice+2,
    757 "// Castelnuovo-Mumford regularity of C:");
    758    return(HR);
    759   }
     1320   if ( size(K) == 0 )
     1321     {
     1322       ring mr=0,x(1..n-1),dp;
     1323       ideal I=imap(nr,I);
     1324       H=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
     1325       time=rtimer-time;
     1326       // Additional information:
     1327       dbprint(printlevel-voice+2,
     1328               "// The sequence of integers defines a monomial curve C in P"
     1329               + string(n));
     1330       dbprint(printlevel-voice+2,
     1331               "//    C is arithmetically Cohen-Macaulay");
     1332       dbprint(printlevel-voice+2,
     1333               "//    Regularity attained at the last step of a m.g.f.r. of I(C)");
     1334       dbprint(printlevel-voice+2,
     1335               "//    Regularity of the Hilbert function of S/I(C): "
     1336               + string(H-2));
     1337       dbprint(printlevel-voice+2,
     1338               "//    Time for computing ideal I(C) (by elimination): "
     1339               + string(ttime) + " sec.");
     1340       dbprint(printlevel-voice+2,
     1341               "//    Time for computing reg(C) once I(C) has been determined: "
     1342               + string(time) + " sec.");
     1343       dbprint(printlevel-voice+2,
     1344               "// The Castelnuovo-Mumford regularity of C is:");
     1345       return(H);
     1346     }
     1347   else
     1348     {
     1349       KK=simplify(reduce(quotient(I,maxideal(1)),I),2);
     1350       firstind=maxdeg1(KK)+1;
     1351       J=subst(I,x(n),1);
     1352       ring mr=0,x(1..n-1),dp;
     1353       ideal J=imap(nr,J);
     1354       lastind=maxdeg1(minbase(quotient(J,maxideal(1))))+1;
     1355       H=firstind;
     1356       if ( lastind > H )
     1357         {
     1358           H=lastind;
     1359         }
     1360       time=rtimer-time;
     1361       // Additional information:
     1362       dbprint(printlevel-voice+2,
     1363               "// The sequence of integers defines a monomial curve C in P"
     1364               + string(n));
     1365       dbprint(printlevel-voice+2,
     1366               "//    C is not arithmetically Cohen-Macaulay");
     1367       dbprint(printlevel-voice+2,
     1368               "//    end(H^1(S/I(C))) = "
     1369               +string(firstind-2));
     1370       dbprint(printlevel-voice+2,
     1371               "//    Upper bound for the a-invariant of S/I(C): end(H^2(S/I(C))) <= "
     1372               +string(lastind-3));
     1373       if ( H == firstind )
     1374         {
     1375           dbprint(printlevel-voice+2,
     1376                   "//    Regularity attained at the last step of a m.g.f.r. of I(C)");
     1377           dbprint(printlevel-voice+2,
     1378                   "//    Regularity of the Hilbert function of S/I(C): "
     1379                   + string(H-1));
     1380         }
     1381       else
     1382         {
     1383           dbprint(printlevel-voice+2,
     1384                   "//    Regularity attained at the second last step of a m.g.f.r. of I(C)");
     1385           dbprint(printlevel-voice+2,
     1386                   "//    (and not attained at the last step)");
     1387         }
     1388       dbprint(printlevel-voice+2,
     1389               "//    Time for computing ideal I(C) (by elimination): "
     1390               + string(ttime) + " sec.");
     1391       dbprint(printlevel-voice+2,
     1392               "//    Time for computing reg(C) once I(C) has been determined: "
     1393               + string(time) + " sec.");
     1394       dbprint(printlevel-voice+2,
     1395               "// The Castelnuovo-Mumford regularity of C is:");
     1396       return(H);
     1397     }
    7601398}
    7611399example
    7621400{ "EXAMPLE:"; echo = 2;
    7631401// The 1st example is the twisted cubic:
    764    reg_moncurve(0,1,2,3);
     1402   regMonCurve(0,1,2,3);
    7651403// The 2nd. example is the non arithm. Cohen-Macaulay monomial curve in P4
    7661404// parametrized by: x(0)-s6,x(1)-s5t,x(2)-s3t3,x(3)-st5,x(4)-t6:
    767    reg_moncurve(0,1,3,5,6);
    768 // Additional information can be obtained as follows:
    769    printlevel = 1;
    770    reg_moncurve(0,1,3,5,6);
     1405   regMonCurve(0,1,3,5,6);
     1406// Additional information is displayed if you change printlevel (=1);
    7711407}
    772 ///////////////////////////////////////////////////////////////////////////////
     1408////////////////////////////////////////////////////////////////////////////////
    7731409/*
    7741410Out-commented examples:
    7751411//
    7761412// The sequence of integers must be strictly increasing
    777 // and the first integer is 0:
    778    reg_moncurve(1,4,6,9);
    779    reg_moncurve(0,3,8,5,23);
    780    reg_moncurve(0,4,7,7,9);
     1413   regMonCurve(1,4,6,9);
     1414   regMonCurve(0,3,8,5,23);
     1415   regMonCurve(0,4,7,7,9);
    7811416//
    7821417// A curve in P3 s.t. the regularity is attained at the last step:
    783    reg_moncurve(0,2,12,15);
     1418   regMonCurve(0,2,12,15);
    7841419//
    7851420// A curve in P4 s.t. the regularity attained at the last but one
    786 // but NOT at the last step:
    787    reg_moncurve(0,5,9,11,20);
    788 //
    789 // A curve in P8 s.t. the m.g.f.r. of the defining ideal could not be obtained
    790 // by the command mres:
    791    reg_moncurve(0,1,2,3,9,11,18,24,25);
     1421// but NOT at the last step (Ex. 3.3 Preprint 2004):
     1422   regMonCurve(0,5,9,11,20);
     1423//
     1424// A curve in P8 s.t. the m.g.f.r. of the defining ideal is not easily
     1425// obtained through m.g.f.r.:
     1426   regMonCurve(0,1,2,3,9,11,18,24,25);
    7921427//
    7931428// A curve in P11 of degree 37:
    794    reg_moncurve(0,1,2,7,16,17,25,27,28,30,36,37);
     1429   regMonCurve(0,1,2,7,16,17,25,27,28,30,36,37);
    7951430// It takes some time to compute the eliminated ideal; the computation of
    796 // the regularity is then rather fast as one can check using proc_curve:
     1431// the regularity is then rather fast as one can check using proc regIdeal:
    7971432   ring q=0,(s,t,x(0..11)),dp;
    798    ideal i=x(0)-st36,x(1)-s2t35,x(2)-s7t30,x(3)-s16t21,x(4)-s17t20,x(5)-s25t12
     1433   ideal i=x(0)-st36,x(1)-s2t35,x(2)-s7t30,x(3)-s16t21,x(4)-s17t20,x(5)-s25t12,
    7991434           x(6)-s27t10,x(7)-s28t9,x(8)-s30t7,x(9)-s36t,x(10)-s37,x(11)-t37;
    8001435   ideal ei=eliminate(i,st);
    8011436   ring qq=0,x(0..11),dp;
    8021437   ideal i=imap(q,ei);
    803    reg_curve(i,1);
     1438   regIdeal(i);
    8041439//
    8051440// A curve in P14 of degree 55:
    806    reg_moncurve(0,1,2,7,16,17,25,27,28,30,36,37,40,53,55);
    807 // In the last three examples, the m.g.f.r. could not be obtained using mres.
     1441   regMonCurve(0,1,2,7,16,17,25,27,28,30,36,37,40,53,55);
    8081442//
    8091443*/
    810 //////////////////////////////////////////////////////////////////////////////
     1444///////////////////////////////////////////////////////////////////////////////
     1445///////////////////////////////////////////////////////////////////////////////
     1446///////////////////////////////////////////////////////////////////////////////
    8111447
     1448proc NoetherPosition (ideal i)
     1449"
     1450USAGE:   NoetherPosition (i); i ideal
     1451RETURN:  ideal such that, for the homogeneous linear transformation
     1452         map phi=S,NoetherPosition(i);
     1453         one has that K[x(n-d+1),...,x(n)] is a Noether normalization of
     1454         S/phi(i) where S=K[x(0),...x(n)] is the basering and d=dim(S/i).
     1455         (returns -1 if i = (0) or (1)).
     1456ASSUME:  The field K is infinite and i is a nonzero proper ideal.
     1457NOTES    1. It works also if K is a finite field if it terminates, but
     1458            may result in an infinite loop. If the procedure enters more
     1459            than 30 loops, -1 is returned and a warning message is displayed.
     1460         2. If printlevel > 0 (default = 0), additional info is displayed:
     1461            dim(S/i) and K[x(n-d+1),...,x(n)] are given.
     1462EXAMPLE: example NoetherPosition; shows some examples
     1463"
     1464{
     1465//--------------------------- initialisation ---------------------------------
     1466   int ii,jj,d,time,nl,NPtest;
     1467   intmat ran;
     1468   def r0 = basering;
     1469   ideal K,chcoord;
     1470   int n = nvars(r0)-1;
     1471   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
     1472   execute(s);
     1473   ideal i,sbi,I,K,chcoord,m;
     1474   poly P;
     1475   map phi;
     1476   i = fetch(r0,i);
     1477   time=rtimer;
     1478   sbi=std(i);
     1479   I=simplify(lead(sbi),1);
     1480   attrib(I,"isSB",1);
     1481   d=dim(I);
     1482//----- If the ideal i is not proper:
     1483   if ( d == -1 )
     1484     {
     1485       "// WARNING from proc NoetherPosition from lib mregular.lib:
     1486// The ideal i is (1)!";
     1487       return (-1);
     1488     }
     1489//----- If the ideal i is 0:
     1490   if ( size(I) == 0 )
     1491     {
     1492       "// WARNING from proc NoetherPosition from lib mregular.lib:
     1493// The ideal i is (0)!";
     1494       return (-1);
     1495     }
     1496//----- When the ideal i is 0-dimensional:
     1497   if ( d == 0 )
     1498     {
     1499       time=rtimer-time;
     1500       // Additional information:
     1501       dbprint(printlevel-voice+2,
     1502               "// Dimension of S/i: 0");
     1503       dbprint(printlevel-voice+2,
     1504               "// Time for computing a Noether normalization: "
     1505               + string(time) + " sec.");
     1506       dbprint(printlevel-voice+2,
     1507               "// K is a Noether normalization of S/phi(i)");
     1508       dbprint(printlevel-voice+2,
     1509               "// where the map phi: S --> S is:");
     1510       setring r0;
     1511       return (maxideal(1));
     1512     }
     1513   NPtest=is_NP(I);
     1514   if ( NPtest == 1 )
     1515     {
     1516       K=x(n-d+1..n);
     1517       setring r0;
     1518       K=fetch(r1,K);
     1519       time=rtimer-time;
     1520       // Additional information:
     1521       dbprint(printlevel-voice+2,
     1522               "// Dimension of S/i: " + string(d) );
     1523       dbprint(printlevel-voice+2,
     1524               "// Time for computing a Noether normalization: " +
     1525               string(time) + " sec.");
     1526       dbprint(printlevel-voice+2,
     1527               "// K[" + string(K) +
     1528               "] is a Noether normalization of S/phi(i)");
     1529       dbprint(printlevel-voice+2,
     1530               "// where the map phi: S --> S is:");
     1531       return (maxideal(1));
     1532     }
     1533//---- Otherwise, random change of coordinates and
     1534//---- test for Noether normalization.
     1535//---- If we were unlucky, another change of coord. will be done:
     1536   while ( nl < 30 )
     1537   {
     1538     chcoord=select1(maxideal(1),1,(n-d+1));
     1539     nl=nl+1;
     1540     for ( ii = 1; ii<=d; ii++ )
     1541     {
     1542       ran=random(100,1,n-d+ii);
     1543       ran=intmat(ran,1,n-d+ii+1);
     1544       ran[1,n-d+ii+1]=1;
     1545       m=select1(maxideal(1),1,(n-d+1+ii));
     1546       for ( jj = 1; jj<=n-d+ii+1; jj++ )
     1547       {
     1548       P=P+ran[1,jj]*m[jj];
     1549       }
     1550       chcoord[n-d+1+ii]=P;
     1551       P=0;
     1552     }
     1553     phi=r1,chcoord;
     1554     dbprint(printlevel-voice+2,"// (1 random change of coord.)");
     1555     I=simplify(lead(std(phi(i))),1);
     1556     attrib(I,"isSB",1);
     1557     NPtest=is_NP(I);
     1558     if ( NPtest == 1 )
     1559       {
     1560       K=x(n-d+1..n);
     1561       setring r0;
     1562       K=fetch(r1,K);
     1563       chcoord=fetch(r1,chcoord);
     1564       time=rtimer-time;
     1565       // Additional information:
     1566       dbprint(printlevel-voice+2,
     1567               "// Dimension of S/i: " + string(d) );
     1568       dbprint(printlevel-voice+2,
     1569               "// Time for computing a Noether normalization: " +
     1570               string(time) + " sec.");
     1571       dbprint(printlevel-voice+2,
     1572               "// K[" + string(K) +
     1573               "] is a Noether normalization of S/phi(i)");
     1574       dbprint(printlevel-voice+2,
     1575               "// where the map phi: S --> S is:");
     1576       return (chcoord);
     1577       }
     1578   }
     1579       "// WARNING from proc NoetherPosition from lib mregular.lib:
     1580// The procedure has entered in more than 30 loops: in your example
     1581// the method may enter an infinite loop over a finite field!";
     1582       return (-1);
     1583}
     1584example
     1585{ "EXAMPLE:"; echo = 2;
     1586   ring r=0,(x,y,z,t,u),dp;
     1587   ideal i1=y,z,t,u; ideal i2=x,z,t,u; ideal i3=x,y,t,u; ideal i4=x,y,z,u;
     1588   ideal i5=x,y,z,t; ideal i=intersect(i1,i2,i3,i4,i5);
     1589   map phi=r,NoetherPosition(i);
     1590   phi;
     1591   ring r5=5,(x,y,z,t,u),dp;
     1592   ideal i=imap(r,i);
     1593   map phi=r5,NoetherPosition(i);
     1594   phi;
     1595// Additional information is displayed if you change printlevel (=1);
     1596}
     1597///////////////////////////////////////////////////////////////////////////////
     1598///////////////////////////////////////////////////////////////////////////////
     1599
     1600proc is_NP (ideal i)
     1601"
     1602USAGE:   is_NP (i); i ideal
     1603RETURN:  1  if K[x(n-d+1),...,x(n)] is a Noether normalization of
     1604            S/i where S=K[x(0),...x(n)] is the basering, and d=dim(S/i),
     1605         0  otherwise.
     1606         (returns -1 if i=(0) or i=(1)).
     1607ASSUME:  i is a nonzero proper homogeneous ideal.
     1608NOTE:    1. If i is not homogeneous and is_NP(i)=1 then K[x(n-d+1),...,x(n)]
     1609            is a Noether normalization of S/i. The converse may be wrong if
     1610            the ideal is not homogeneous.
     1611         2. is_NP is used in the procedures regIdeal, depthIdeal, satiety,
     1612            and NoetherPosition.
     1613EXAMPLE: example is_NP; shows some examples
     1614"
     1615{
     1616//--------------------------- initialisation ---------------------------------
     1617   int ii,d,dz;
     1618   def r0 = basering;
     1619   int n = nvars(r0)-1;
     1620   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
     1621   execute(s);
     1622   ideal i,sbi,I,J;
     1623   i = fetch(r0,i);
     1624   sbi=std(i);
     1625   I=simplify(lead(sbi),1);
     1626   attrib(I,"isSB",1);
     1627   d=dim(I);
     1628//----- If the ideal i is not proper:
     1629   if ( d == -1 )
     1630     {
     1631       "// WARNING from proc is_NP from lib mregular.lib:
     1632// The ideal i is (1)!";
     1633       return (-1);
     1634     }
     1635//----- If the ideal i is 0:
     1636   if ( size(I) == 0 )
     1637     {
     1638       "// WARNING from proc is_NP from lib mregular.lib:
     1639// The ideal i is (0)!";
     1640       return (-1);
     1641     }
     1642//----- When the ideal i is 0-dimensional:
     1643   if ( d == 0 )
     1644     {
     1645       return (1);
     1646     }
     1647//----- Check Noether position
     1648   J=I;
     1649   for ( ii = n-d+1; ii <= n; ii++ )
     1650     {
     1651       J=subst(J,x(ii),0);
     1652     }
     1653   attrib(J,"isSB",1);
     1654   dz=dim(J);
     1655   if ( dz == d )
     1656     {
     1657       return (1);
     1658     }
     1659   else
     1660     {
     1661       return(0);
     1662     }
     1663}
     1664example
     1665{ "EXAMPLE:"; echo = 2;
     1666   ring r=0,(x,y,z,t,u),dp;
     1667   ideal i1=y,z,t,u; ideal i2=x,z,t,u; ideal i3=x,y,t,u; ideal i4=x,y,z,u;
     1668   ideal i5=x,y,z,t; ideal i=intersect(i1,i2,i3,i4,i5);
     1669   is_NP(i);
     1670   ideal ch=x,y,z,t,x+y+z+t+u;
     1671   map phi=ch;
     1672   is_NP(phi(i));
     1673}
     1674///////////////////////////////////////////////////////////////////////////////
     1675///////////////////////////////////////////////////////////////////////////////
     1676
     1677proc is_nested (ideal i)
     1678"
     1679USAGE:   is_nested (i); i monomial ideal
     1680RETURN:  1 if i is of nested type, 0 otherwise.
     1681         (returns -1 if i=(0) or i=(1)).
     1682ASSUME:  i is a nonzero proper monomial ideal.
     1683NOTES:   1. The ideal must be monomial, otherwise the result has no meaning
     1684            (so check this before using this procedure).
     1685         2. is_nested is used in procedures depthIdeal, regIdeal and satiety.
     1686         3. When i is a monomial ideal of nested type of S=K[x(0)..x(n)],
     1687            the a-invariant of S/i coincides with the upper bound obtained
     1688            using the procedure regIdeal with printlevel > 0.
     1689THEORY:  A monomial ideal is of nested type if its associated primes are all
     1690         of the form (x(0),...,x(i)) for some i<=n.
     1691         (see definition and effective criterion to check this property in
     1692         the preprint 'Saturation and Castelnuovo-Mumford regularity' by
     1693         Bermejo-Gimenez, 2004).
     1694EXAMPLE: example is_nested; shows some examples
     1695"
     1696{
     1697//--------------------------- initialisation ---------------------------------
     1698   int ii,d,tev,lastv,h,NPtest;
     1699   def r0 = basering;
     1700   int n = nvars(r0)-1;
     1701   string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
     1702   execute(s);
     1703   ideal I,K,KK,LL;
     1704   I = fetch(r0,i);
     1705   I=minbase(I);
     1706   attrib(I,"isSB",1);
     1707   d=dim(I);
     1708//----- If the ideal i is not proper:
     1709   if ( d == -1 )
     1710   {
     1711       "// WARNING from proc is_nested from lib mregular.lib:
     1712// The ideal i is (1)!";
     1713       return (-1);
     1714   }
     1715//----- If the ideal i is 0:
     1716   if ( size(I) == 0 )
     1717     {
     1718       "// WARNING from proc is_nested from lib mregular.lib:
     1719// The ideal i is (0)!";
     1720       return (-1);
     1721     }
     1722//----- When the ideal i is 0-dimensional:
     1723   if ( d == 0 )
     1724     {
     1725       return (1);
     1726     }
     1727//----- Check Noether position
     1728   NPtest=is_NP(I);
     1729   if ( NPtest != 1 )
     1730   {
     1731       return (0);
     1732   }
     1733//----- When ideal is 1-dim. + var. in Noether position -> Nested Type
     1734   if ( d == 1 )
     1735     {
     1736       return (1);
     1737     }
     1738//----- Determ. of the last variable really occuring
     1739   lastv=n-d;
     1740   h=n;
     1741   while ( lastv == n-d and h > n-d )
     1742     {
     1743       K=select(I,h+1);
     1744       if ( size(K) == 0 )
     1745         {
     1746           h=h-1;
     1747         }
     1748       else
     1749         {
     1750           lastv=h;
     1751         }
     1752     }
     1753//----- Check the second property by evaluation when NP + d>1
     1754   KK=subst(I,x(lastv),1);
     1755   for ( ii = n-lastv; ii<=d-2; ii++ )
     1756   {
     1757     LL=minbase(subst(I,x(n-ii-1),1));
     1758     attrib(LL,"isSB",1);
     1759     tev=size(reduce(KK,LL));
     1760     if ( tev > 0 )
     1761       {
     1762         return(0);
     1763       }
     1764     KK=LL;
     1765   }
     1766   return(1);
     1767}
     1768example
     1769{ "EXAMPLE:"; echo = 2;
     1770   ring s=0,(x,y,z,t),dp;
     1771   ideal i1=x2,y3; ideal i2=x3,y2,z2; ideal i3=x3,y2,t2;
     1772   ideal i=intersect(i1,i2,i3);
     1773   is_nested(i);
     1774   ideal ch=x,y,z,z+t;
     1775   map phi=ch;
     1776   ideal I=lead(std(phi(i)));
     1777   is_nested(I);
     1778}
     1779///////////////////////////////////////////////////////////////////////////////
     1780///////////////////////////////////////////////////////////////////////////////
     1781

  • Singular/LIB/zeroset.lib

    rc00e9b re182c8  
    1 // Last change 12.02.2001 (Eric Westenberger)
    2 ///////////////////////////////////////////////////////////////////////////////
    3 version="$Id: zeroset.lib,v 1.9 2005-02-22 10:36:38 Singular Exp $";
     1// changed 12.02.2001 (Eric Westenberger)
     2// last change 29.10.2004 (Bayer Thomas, Quotient, QuotientMain)
     3///////////////////////////////////////////////////////////////////////////////
     4version="$Id: zeroset.lib,v 1.10 2005-04-19 15:23:42 Singular Exp $";
    45category="Symbolic-numerical solving";
    56info="
    67LIBRARY:  zeroset.lib      Procedures For Roots and Factorization
    7 AUTHOR:   Thomas Bayer,    email: tbayer@mathematik.uni-kl.de
    8           http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/
     8AUTHOR:   Thomas Bayer,    email: bayert@web.de
     9          http://www14.informatik.tu-muenchen.de/personen/bayert/
    910          Current Adress: Institut fuer Informatik, TU Muenchen
    1011
    … …  
    1415 where a is an algebraic number. Written in the frame of the
    1516 diploma thesis (advisor: Prof. Gert-Martin Greuel) 'Computations of moduli
    16  spaces of semiquasihomogeneous singularities and an implementation in Singular'.
     17 spaces of semiquasihomogeneous singularities and an implementation in
     18 Singular'.
    1719 This library is meant as a preliminary extension of the functionality
    1820 of Singular for univariate factorization of polynomials over simple algebraic
    … …  
    2325
    2426PROCEDURES:
    25  Quotient(f, g)    quotient q  of f w.r.t. g (in f = q*g + remainder)
    26  Remainder(f,g)    remainder of the division of f by g
    27  Roots(f)    computes all roots of f in an extension field of Q
    28  SQFRNorm(f)    norm of f (f must be squarefree)
    29  ZeroSet(I)    zero-set of the 0-dim. ideal I
     27 EGCD(f, g)      gcd over an algebraic extension field of Q
     28 Factor(f)       factorization of f over an algebraic extension field
     29 Quotient(f, g)  quotient q  of f w.r.t. g (in f = q*g + remainder)
     30 Remainder(f,g)  remainder of the division of f by g
     31 Roots(f)        computes all roots of f in an extension field of Q
     32 SQFRNorm(f)     norm of f (f must be squarefree)
     33 ZeroSet(I)      zero-set of the 0-dim. ideal I
    3034
    3135AUXILIARY PROCEDURES:
     36 EGCDMain(f, g)       gcd over an algebraic extension field of Q
     37 FactorMain(f)        factorization of f over an algebraic extension field
    3238 InvertNumberMain(c)  inverts an element of an algebraic extension field
    33  QuotientMain(f, g)  quotient of f w.r.t. g
    34  RemainderMain(f,g)  remainder of the division of f by g
    35  RootsMain(f)    computes all roots of f, might extend the ground field
    36  SQFRNormMain(f)  norm of f (f must be squarefree)
     39 QuotientMain(f, g)   quotient of f w.r.t. g
     40 RemainderMain(f,g)   remainder of the division of f by g
     41 RootsMain(f)         computes all roots of f, might extend the ground field
     42 SQFRNormMain(f)      norm of f (f must be squarefree)
    3743 ContainedQ(data, f)  f in data ?
    38  SameQ(a, b)    a == b (list a,b)
     44 SameQ(a, b)          a == b (list a,b)
    3945";
    40 // Factor(f)    factorization of f over an algebraic extension field
    41 // EGCD(f, g)    gcd over an algebraic extension field of Q
    42 // EGCDMain(f, g)    gcd over an algebraic extension field of Q
    43 // FactorMain(f)    factorization of f over an algebraic extension field
    4446
    4547LIB "primitiv.lib";
    … …  
    4749
    4850// note : return a ring : ring need not be exported !!!
    49 
    50 // Artihmetic in Q(a)[x] without built-in procedures
     51// Arithmetic in Q(a)[x] without built-in procedures
    5152// assume basering = Q[x,a] and minpoly is represented by mpoly(a).
    52 // the algorithms are taken from "Polynomial Algorithms in Computer Algebra",
     53// The algorithms are taken from "Polynomial Algorithms in Computer Algebra",
    5354// F. Winkler, Springer Verlag Wien, 1996.
    5455
    … …  
    6869// extension field is represented as Q[x...,a] together with the ideal
    6970// 'mpoly' (attribute "isSB"). The arithmetic in the extension field is
    70 // implemented in the procedures in the procedures 'MultPolys' (multiplication)
     71// implemented in the procedures 'MultPolys' (multiplication)
    7172// and 'InvertNumber' (inversion). After addition and substraction one should
    7273// apply 'SimplifyPoly' to the result to reduce the result w.r.t. 'mpoly'.
    … …  
    105106  def ROR = TransferRing(basering);
    106107  setring ROR;
    107   export(ROR);
    108108
    109109  // get the polynomial f and find the roots
    … …  
    173173"
    174174{
     175  def altring=basering;
    175176  int i, linFactors, nlinFactors, dbPrt;
    176177  intvec wt = 1,0;    // deg(a) = 0
    … …  
    338339  def ZSR = TransferRing(basering);
    339340  setring ZSR;
    340   export(ZSR);
    341341
    342342  // get ideal I and find the zero-set
    … …  
    467467
    468468proc Quotient(poly f, poly g)
    469 "USAGE:   Quotient(f, g); where f,g are polynomials;
    470 PURPOSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < deg(g)
     469"USAGE:   Quotient(f, g); poly f, g;
     470PUROPSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < g
    471471RETURN:  list of polynomials
    472   @format
    473   _[1] = quotient  q
    474   _[2] = remainder r
    475   @end format
     472         _[1] = quotient  q
     473         _[2] = remainder r
    476474ASSUME:  basering = Q[x] or Q(a)[x]
    477475EXAMPLE: example  Quotient; shows an example
    478476"
    479477{
    480   def QUOB = basering;
    481   def QUOR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
    482   setring QUOR;
    483   export(QUOR);
    484   poly f = imap(QUOB, f);
    485   poly g = imap(QUOB, g);
    486   list result = QuotientMain(f, g);
    487 
    488   setring(QUOB);
    489   list result = imap(QUOR, result);
    490   kill(QUOR);
    491   return(result);
    492 }
     478        def QUOB = basering;
     479        def QUOR = TransferRing(basering);      // new ring with parameter 'a' replaced by a variable
     480        setring QUOR;
     481        poly f = imap(QUOB, f);
     482        poly g = imap(QUOB, g);
     483        list result = QuotientMain(f, g);
     484
     485        setring(QUOB);
     486        list result = imap(QUOR, result);
     487        kill(QUOR);
     488        return(result);
     489}
     490
    493491example
    494492{"EXAMPLE:";  echo = 2;
    … …  
    503501
    504502proc QuotientMain(poly f, poly g)
    505 "USAGE:   QuotientMain(f, g); where f,g are polynomials
    506 PURPOSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < deg(g)
     503"USAGE:   QuotientMain(f, g); poly f, g
     504PUROPSE: compute the quotient q and remainder r s.t. f = g*q + r, deg(r) < g
    507505RETURN:  list of polynomials
    508   @format
    509   _[1] = quotient  q
    510   _[2] = remainder r
    511   @end format
    512 ASSUME:  basering = Q[x,a] and ideal mpoly is defined (it might be 0),
    513          this represents the ring Q(a)[x] together with its minimal polynomial.
     506         _[1] = quotient  q
     507         _[2] = remainder r
     508ASSUME:  basering = Q[x,a] and ideal 'mpoly' is defined (it might be 0),
     509         this represents the ring Q(a)[x] together with 'minpoly'.
    514510EXAMPLE: example  Quotient; shows an example
    515511"
    516512{
    517   if(g == 0) { ERROR("Division by zero !");}
    518 
    519   def QMB = basering;
    520   def QMR = NewBaseRing();
    521   setring QMR;
    522   poly f, g, h;
    523   h = imap(QMB, f) / imap(QMB, g);
    524   setring QMB;
    525   return(list(imap(QMR, h), 0));
     513        def altring=basering;
     514        int dg;
     515        intvec wt = 1,0;
     516        list ltList;
     517        poly fn, q, p, lcgInv, lmonog;
     518
     519        if(g == 0) { ERROR("Division by zero !");}
     520
     521        ltList = LeadTerm(g, 1);
     522        dg = deg(ltList[2]);
     523        if(dg == 0) {                   // g is a constant
     524                result[1] = MultPolys(f, InvertNumberMain(g));
     525                result[2] = 0;
     526                return(result);
     527        }
     528        fn = f;
     529        q = 0;
     530        lcgInv = InvertNumberMain(ltList[3]);   // inverse of leading coef of g
     531        lmonog = ltList[2];                     // leading monomial of g
     532
     533        while(deg(fn, wt) >= dg) {              // degree of fn > degree of g
     534                p = MultPolys(LeadTerm(fn, 1)[1]/lmonog, lcgInv);
     535                q = SimplifyPoly(q + p);
     536                fn = SimplifyPoly(fn - MultPolys(p, g));
     537        }
     538        return(list(q, fn));
    526539}
    527540
    … …  
    539552  def REMR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
    540553  setring(REMR);
    541   export(REMR);
    542554  poly f = imap(REMB, f);
    543555  poly g = imap(REMB, g);
    … …  
    566578"
    567579{
     580  def altring=basering;
    568581  int dg;
    569582  intvec wt = 1,0;;
    … …  
    587600ASSUME:  basering = Q(a)[t]
    588601EXAMPLE: example  EGCD; shows an example
    589 NOTE: obsolete: use gcd
    590602"
    591603{
    … …  
    593605  def GCDR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
    594606  setring GCDR;
    595   export(GCDR);
    596607  poly f = imap(GCDB, f);
    597608  poly g = imap(GCDB, g);
    … …  
    663674  execute("poly @f = " + @sf + ";");
    664675  execute("poly @g = " + @sg + ";");
    665   export(@RGCD);
    666676  poly @h = EGCD(@f, @g);
    667677  setring(@RGCDB);
    … …  
    690700  def SNR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
    691701  setring SNR;
    692   export(SNR);
    693702  poly f = imap(SNB, f);
    694703  list result = SQFRNormMain(f);  // squarefree norm of f
    … …  
    722731"
    723732{
     733  def altring=basering;
    724734  int s = 0;
    725735  intvec wt = 1,0;
    … …  
    770780         is Q or a simple extension of Q given by a minpoly.
    771781NOTE:    if basering = Q[t] then this is the built-in @code{factorize}
    772 NOTE:    obsolete: use factorize
    773782EXAMPLE: example  Factor; shows an example
    774783"
    … …  
    777786  def FRR = TransferRing(basering);  // new ring with parameter 'a' replaced by a variable
    778787  setring FRR;
    779   export(FRR);
    780788  poly f = imap(FRB, f);
    781789  list result = FactorMain(f);  // squarefree norm of f
    … …  
    11321140
    11331141  def NLZR = basering;
    1134   export(NLZR);
    11351142
    11361143  n = nvars(basering) - 1;
    … …  
    11611168      list roots;
    11621169      poly f = imap(NLZR, f);
    1163       export(RIS1);
    11641170      export(mpoly);
    11651171      roots = RootsMain(f);
    1166 
    11671172      // get "old" basering with new minpoly
    11681173
Note: See TracChangeset for help on using the changeset viewer.