Changeset f6e355 in git for Singular/LIB/equising.lib


Ignore:
Timestamp:
Apr 15, 2005, 5:20:12 PM (19 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
d7952be20b58d1b7b3d37a3c0f03393781afc143
Parents:
20f0951437446030fdebc2cf2b571587c25443e6
Message:
*lossen:  new lib


git-svn-id: file:///usr/local/Singular/svn/trunk@7836 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/equising.lib

    r20f095 rf6e355  
    1 version="$Id: equising.lib,v 1.8 2001-08-27 14:47:49 Singular Exp $";
     1version="$Id: equising.lib,v 1.9 2005-04-15 15:20:12 Singular Exp $";
    22category="Singularities";
    33info="
    44LIBRARY:  equising.lib  Equisingularity Stratum of a Family of Plane Curves
    5 AUTHOR:   Andrea Mindnich, mindnich@mathematik.uni-kl.de
    6 
    7 PROCEDURES:
    8  esStratum(F[,m]);   computes the equisingularity stratum of a deformation
    9  isEquising(F[,m]);  tests if a given deformation is equisingular
     5AUTHOR:   Christoph Lossen, lossen@mathematik.uni-kl.de
     6          Andrea Mindnich, mindnich@mathematik.uni-kl.de
     7
     8MAIN PROCEDURES:
     9 tau_es(f);             codim of mu-const stratum in semi-universal def. base
     10 esIdeal(f);            (Wahl's) equisingularity ideal of f
     11 esStratum(F[,m,L]);    equisingularity stratum of a family F
     12 isEquising(F[,m,L]);   tests if a given deformation is equisingular
     13
     14AUXILIARY PROCEDURE:
     15 control_Matrix(M);     computes list of blowing-up data
    1016";
    1117
     18LIB "hnoether.lib";
    1219LIB "poly.lib";
    1320LIB "elim.lib";
    14 LIB "hnoether.lib";
    15 ///////////////////////////////////////////////////////////////////////////////
     21LIB "deform.lib";
     22LIB "sing.lib";
     23
     24////////////////////////////////////////////////////////////////////////////////
     25//
     26//  The following (static) procedures are used by esComputation
     27//
     28////////////////////////////////////////////////////////////////////////////////
    1629// COMPUTES  a weight vector. x and y get weight 1 and all other
    1730//           variables get weight 0.
    … …  
    2437  return (iv);
    2538}
    26 ///////////////////////////////////////////////////////////////////////////////
    27 // exchanges the variables x and y in the polynomial p_f
     39
     40////////////////////////////////////////////////////////////////////////////////
     41// exchanges the variables x and y in the polynomial f
    2842static proc swapXY(poly f)
    2943{
    … …  
    3650  return (f);
    3751}
    38 ///////////////////////////////////////////////////////////////////////////////
    39 // ASSUME:   p_mon is a monomial
    40 //            and p_g is the product of the variables in p_mon.
    41 // COMPUTES the coefficient of p_mon in p_h.
    42 static proc coefficient(poly p_h, poly p_mon, poly p_g)
     52
     53////////////////////////////////////////////////////////////////////////////////
     54// computes m-jet w.r.t. the variables x,y (other variables weighted 0
     55static proc m_Jet(poly F,int m);
    4356{
    44   matrix coefMatrix = coef(p_h,p_g);
    45   int nc = ncols(coefMatrix);
    46   int ii=1;
    47   poly p_c=0;
    48   while(ii<=nc)
    49   {
    50     if (coefMatrix[1,ii] == p_mon)
    51     {
    52       p_c = coefMatrix[2,ii];
    53       break;
    54     }
    55     ii++;
    56   }
    57   return (p_c);
     57  intvec w=xyVector(); 
     58  poly Fd=jet(F,m,w); 
     59  return(Fd);
    5860}
    59 ///////////////////////////////////////////////////////////////////////////////
    60 // in p_F the variable p_vari is substituted by the polynomial p_g.
    61 static proc eSubst(poly p_F, poly p_vari, poly p_g)
     61
     62
     63////////////////////////////////////////////////////////////////////////////////
     64// computes the 4 control matrices (input is multsequence(L))
     65proc control_Matrix(list M);
     66"USAGE:   control_Matrix(L); L list
     67ASSUME:  L is the output of multsequence(hnexpansion(f)).
     68RETURN:  list M of 4 intmat's:
     69@format
     70  M[1] contains the multiplicities at the respective infinitely near points
     71       p[i,j] (i=step of blowup+1, j=branch) -- if branches j=k,...,k+m pass
     72       through the same p[i,j] then the multiplicity is stored in M[1][k,j],
     73       while M[1][k+1]=...=M[1][k+m]=0.   
     74  M[2] contains the number of branches meeting at p[i,j] (again, the information
     75       is stored according to the above rule)   
     76  M[3] contains the information about the splitting of M[1][i,j] with respect to
     77       different tangents of branches at p[i,j] (information is stored only for
     78       minimal j>=k corresponding to a new tangent direction).
     79       The entries are the sum of multiplicities of all branches with the
     80       respective tangent.
     81  M[4] contains the maximal sum of higher multiplicities for a branch passing
     82       through p[i,j] ( = degree Bound for blowing up) 
     83@end format
     84NOTE:    the branches are ordered in such a way that only consecutive branches
     85         can meet at an infinitely near point. @*
     86         the final rows of the matrices M[1],...,M[3] is (1,1,1,...,1), and
     87         correspond to infinitely near points such that the strict transforms
     88         of the branches are smooth and intersect the exceptional divisor
     89         transversally.
     90SEE ALSO: multsequence
     91EXAMPLE: example control_Matrix; shows an example
     92"
    6293{
    63   def r_base = basering;
    64   ideal MI;
    65   map phi;
    66   int ii = rvar(p_vari);
    67   if (ii != 0)
    68   {
    69     MI = maxideal(1);
    70     MI[ii] = p_g;
    71     phi = r_base, MI;
    72     p_F = phi(p_F);
    73   }
    74   return (p_F);
     94  int i,j,k,dummy;
     95
     96  dummy=0;
     97  for (j=1;j<=ncols(M[2]);j++)
     98  {
     99    dummy=dummy+M[1][nrows(M[1])-1,j]-M[1][nrows(M[1]),j];
     100  }
     101  intmat S[nrows(M[1])+dummy][ncols(M[1])];
     102  intmat T[nrows(M[1])+dummy][ncols(M[1])];
     103  intmat U[nrows(M[1])+dummy][ncols(M[1])];
     104  intmat maxDeg[nrows(M[1])+dummy][ncols(M[1])];
     105 
     106  for (i=1;i<=nrows(M[2]);i++)
     107  {
     108    dummy=1;
     109    for (j=1;j<=ncols(M[2]);j++)
     110    {
     111      for (k=dummy;k<dummy+M[2][i,j];k++)
     112      {
     113        T[i,dummy]=T[i,dummy]+1;
     114        S[i,dummy]=S[i,dummy]+M[1][i,k];
     115        if (i>1)
     116        {
     117          U[i-1,dummy]=U[i-1,dummy]+M[1][i-1,k];
     118        }
     119      }
     120      dummy=k;
     121    }
     122  }
     123
     124  // adding an extra row (in some cases needed to control ES-Stratum
     125  // computation)
     126  for (i=nrows(M[1]);i<=nrows(S);i++)
     127  {
     128    for (j=1;j<=ncols(M[2]);j++)
     129    {
     130      S[i,j]=1; 
     131      T[i,j]=1; 
     132      U[i,j]=1;
     133    }
     134  }
     135 
     136  // Computing the degree Bounds to be stored in M[4]:
     137  for (i=1;i<=nrows(S);i++)
     138  {
     139    dummy=1;
     140    for (j=1;j<=ncols(S);j++)
     141    {
     142      for (k=dummy;k<dummy+T[i,j];k++)
     143      {
     144        maxDeg[i,k]=S[i,dummy];  // multiplicity at i-th blowup
     145      }
     146      dummy=k;
     147    }
     148  }
     149  // adding up multiplicities
     150  for (i=nrows(S);i>=2;i--) 
     151  {
     152    for (j=1;j<=ncols(S);j++)
     153    {
     154      maxDeg[i-1,j]=maxDeg[i-1,j]+maxDeg[i,j];
     155    }
     156  }
     157 
     158  list L=S,T,U,maxDeg;
     159  return(L);
    75160}
    76 ///////////////////////////////////////////////////////////////////////////////
    77 // All ring variables of p_F which occur in (the set of generators of) the
    78 // ideal  Id, are substituted by 0
    79 static proc SimplifyF(poly p_F, ideal Id)
     161
     162
     163////////////////////////////////////////////////////////////////////////////////
     164//  matrix of higher tangent directions:
     165//  returns list: 1) tangent directions
     166//                2) swapping information (x <--> y)
     167static proc inf_Tangents(list L,int s); // L aus hnexpansion,
    80168{
    81   int i=1;
    82   int si = size(Id);
    83   for (i=1; i <= si; i++)
    84   {
    85     if (rvar(Id[i]))
    86     {
    87       p_F = subst(p_F, Id[i], 0);
    88     }
    89   }
    90   return(p_F);
     169  int nv=nvars(basering);
     170  matrix M;
     171  matrix B[s][size(L)];
     172  intvec V;
     173  intmat Mult=multsequence(L)[1];
     174
     175  int i,j,k,counter,e;
     176  for (k=1;k<=size(L);k++)
     177  {
     178    V[k]=L[k][3];  // switch: 0 --> tangent 2nd parameter
     179                   //         1 --> tangent 1st parameter
     180    e=0;   
     181    M=L[k][1];
     182    counter=1;
     183    B[counter,k]=M[1,1];
     184   
     185    for (i=1;i<=nrows(M);i++)
     186    {
     187      for (j=2;j<=ncols(M);j++)
     188      {
     189        counter=counter+1;
     190        if (M[i,j]==var(nv-1))
     191        {
     192          if (i<>nrows(M))
     193          {
     194            B[counter,k]=M[i,j];
     195            j=ncols(M)+1; // goto new row of HNmatrix...
     196            if (counter<>s)
     197            {
     198              if (counter+1<=nrows(Mult))
     199              {
     200                e=Mult[counter-1,k]-Mult[counter,k]-Mult[counter+1,k];
     201              }
     202              else
     203              {
     204                e=Mult[counter-1,k]-Mult[counter,k]-1;
     205              }
     206            }
     207          }
     208          else
     209          {
     210            B[counter,k]=0;
     211            j=ncols(M)+1; // goto new row of HNmatrix...
     212          }
     213        }
     214        else
     215        {
     216          if (e<=0)
     217          {
     218            B[counter,k]=M[i,j];
     219          }
     220          else  // point is still proximate to an earlier point
     221          {
     222            B[counter,k]=y; // marking proximity (without swap....)
     223            if (counter+1<=nrows(Mult))
     224            {
     225              e=e-Mult[counter+1,k];
     226            }
     227            else
     228            {
     229              e=e-1;
     230            }
     231          }
     232        }
     233       
     234        if (counter==s) // given number of points determined
     235        {
     236            j=ncols(M)+1; 
     237            i=nrows(M)+1;
     238            // leave procedure
     239        }
     240      }
     241    }
     242  }
     243  L=B,V;
     244  return(L);
    91245}
    92246
    93 
    94 ///////////////////////////////////////////////////////////////////////////////
    95 
    96 
    97 // Checks, if the basering is admissible.
    98 static proc checkBasering ()
     247////////////////////////////////////////////////////////////////////////////////
     248// compute "good" upper bound for needed number of help variables
     249//
     250static proc Determine_no_b(intmat U,matrix B)
     251// U is assumed to be 3rd output of control_Matrix
     252// B is assumed to be 1st output of inf_Tangents
    99253{
    100   int error = 0;
    101   if(find(charstr(basering),"real"))
    102   {
    103     ERROR ("cannot compute esStratum with 'real' as coefficient field");
    104   }
    105   if (nvars(basering) <= 2)
    106   {
    107    ERROR ("there are to less ringvariables to compute esStratum")
    108   }
    109   error = checkQIdeal(ideal(basering));
    110   return(error);
     254  int nv=nvars(basering);
     255  int i,j,counter;
     256  for (j=1;j<=ncols(U);j++)
     257  {
     258    for (i=1;i<=nrows(U);i++)
     259    {
     260      if (U[i,j]>1)
     261      {
     262        if (B[i,j]<>var(nv-1) and B[i,j]<>var(nv))
     263        {
     264          counter=counter+1;
     265        }
     266      }
     267     
     268    }
     269  }
     270  counter=counter+ncols(U);
     271  return(counter);
    111272}
    112 ///////////////////////////////////////////////////////////////////////////////
    113 static  proc getInput (list #)
     273
     274////////////////////////////////////////////////////////////////////////////////
     275// compute number of infinitely near free points corresponding to non-zero
     276// entries in control_Matrix[1] (except first row)
     277//
     278static proc no_freePoints(intmat Mult,matrix B)
     279// Mult is assumed to be 1st output of control_Matrix
     280// U is assumed to be 3rd output of control_Matrix
     281// B is assumed to be 1st output of inf_Tangents
    114282{
    115   def r_base = basering;
    116 
    117   int maxStep = -1;
    118 
    119   if (size(#) >= 1)
    120   {
    121     if (typeof(#[1]) == "int")
    122     {
    123       maxStep = #[1];
    124     }
    125 
    126     else
    127     {
    128       ERROR("expected esStratum('poly', 'int') ");
    129     }
    130   }
    131 
    132   return(maxStep);
    133 }
    134 //////////////////////////////////////////////////////////////////////////////
    135 // RETURNS: 0, if the ideal cond only depends on the deformation parameters
    136 //          1, otherwise.
    137 static proc checkQIdeal (ideal cond)
    138 {
    139   def r_base = basering;
    140 
    141   int i_print = printlevel-voice + 4;
    142   int i_nvars = nvars(basering);
    143 
    144   ideal id_help = subst(cond,var(i_nvars),0,var(i_nvars-1),0) - cond;
    145 
    146   // cond depends only on the first i_nvars-2 variables <=>
    147   //  id_help == <0>
    148   if ( simplify(id_help, 2) != 0)
    149   {
    150     dbprint(i_print,
    151     "ideal(basering) must only depend on the deformation parameters");
    152     return(1);
    153   }
    154 
    155   return(0);
     283  int i,j,k,counter;
     284  for (j=1;j<=ncols(Mult);j++)
     285  {
     286    for (i=2;i<=nrows(Mult);i++)
     287    {
     288      if (Mult[i,j]>=1)
     289      {
     290        if (B[i-1,j]<>x and B[i-1,j]<>y)
     291        {
     292          counter=counter+1;
     293        }
     294      }
     295    }
     296  }
     297  return(counter);
    156298}
    157299
    … …  
    188330
    189331///////////////////////////////////////////////////////////////////////////////
    190 // Defines a new ring without deformation-parameters.
    191 static proc createHNERing()
     332//
     333// DEFINES: A new basering, "myRing",
     334//          with new names for the parameters and variables.
     335//          The new names for the parameters are a(1..k),
     336//          and t(1..s),x,y for the variables
     337//          The ring ordering is ordStr.
     338// NOTE:    This proc uses 'execute'.
     339static proc createMyRing_new(poly p_F, string ordStr,
     340                                string minPolyStr, int no_b)
    192341{
    193   string str;
    194   string minpolyStr = string(minpoly);
    195 
    196   str = " ring HNERing = (" + charstr(basering) + "), (x,y), ls;";
    197   execute (str);
    198 
    199   str = "minpoly ="+ minpolyStr+";";
    200   execute(str);
    201 
    202   keepring(HNERing);
     342  def r_old = basering;
     343
     344  int chara = char(basering);
     345  string charaStr;
     346  int i;
     347  string helpStr;
     348  int nDefParams = nvars(r_old)-2;
     349
     350  ideal qIdeal = ideal(basering);
     351
     352  if ((npars(basering)==0) and (minPolyStr==""))
     353  {
     354    helpStr = "ring myRing1 ="
     355              + string(chara)+ ", (t(1..nDefParams), x, y),("+ ordStr +");";
     356    execute(helpStr);
     357  }
     358  else
     359  {
     360    charaStr = charstr(basering);
     361    if (charaStr == string(chara) + "," + parstr(basering) or minPolyStr<>"")
     362    {
     363      if (minPolyStr<>"")
     364      {
     365        helpStr = "ring myRing1 =
     366                 (" + string(chara) + ",a),
     367                 (t(1..nDefParams), x, y),(" + ordStr + ");";
     368        execute(helpStr);
     369
     370        execute (minPolyStr);
     371      }
     372      else // no minpoly given
     373      {
     374        helpStr = "ring myRing1 =
     375                  (" + string(chara) + ",a(1..npars(basering)) ),
     376                  (t(1..nDefParams), x, y),(" + ordStr + ");";
     377        execute(helpStr);
     378      }
     379    }
     380    else
     381    {
     382      // ground field is of type (p^k,a)....
     383      i = find (charaStr,",");
     384      helpStr = "ring myRing1 = (" + charaStr[1,i] + "a),
     385              (t(1..nDefParams), x, y),(" + ordStr + ");";
     386      execute (helpStr);
     387    }
     388  }
     389
     390  ideal mIdeal = maxideal(1);
     391  ideal qIdeal = fetch(r_old, qIdeal);
     392  poly p_F = fetch(r_old, p_F);
     393  export p_F,mIdeal;
     394
     395  // Extension by no_b auxiliary variables
     396  if (no_b>0)
     397  {
     398    if (npars(basering) == 0)
     399    {
     400      ordStr = "(dp("+string(no_b)+"),"+ordStr+")";
     401      helpStr = "ring myRing ="
     402                + string(chara)+ ", (b(1..no_b), t(1..nDefParams), x, y),"
     403                + ordStr +";";
     404      execute(helpStr);
     405    }
     406    else
     407    {
     408      charaStr = charstr(basering);
     409      if (charaStr == string(chara) + "," + parstr(basering))
     410      {
     411        if (minpoly !=0)
     412        {
     413          ordStr = "(dp(" + string(no_b) + ")," + ordStr + ")";
     414          minPolyStr = makeMinPolyString("a");
     415          helpStr = "ring myRing =
     416                   (" + string(chara) + ",a),
     417                   (b(1..no_b), t(1..nDefParams), x, y)," + ordStr + ";";
     418          execute(helpStr);
     419
     420          helpStr = "minpoly =" + minPolyStr + ";";
     421          execute (helpStr); 
     422        }
     423        else // no minpoly given
     424        { 
     425          ordStr = "(dp(" + string(no_b) + ")," + ordStr + ")";
     426          helpStr = "ring myRing =
     427                    (" + string(chara) + ",a(1..npars(basering)) ),
     428                    (b(1..no_b), t(1..nDefParams), x, y)," + ordStr + ";";
     429          execute(helpStr);
     430        }
     431      }
     432      else
     433      {
     434        i = find (charaStr,",");
     435        ordStr = "(dp(" + string(no_b) + ")," + ordStr + ")";
     436        helpStr = "ring myRing =
     437                (" + charaStr[1,i] + "a),
     438                (b(1..no_b), t(1..nDefParams), x, y)," + ordStr + ";";
     439        execute (helpStr);
     440      }
     441    }
     442    ideal qIdeal = imap(myRing1, qIdeal);
     443 
     444    if(qIdeal != 0)
     445    {
     446      def r_base = basering;
     447      setring r_base;
     448      kill myRing;
     449      qring myRing = std(qIdeal);
     450    }
     451
     452    poly p_F = imap(myRing1, p_F);
     453    ideal mIdeal = imap(myRing1, mIdeal);
     454    export p_F,mIdeal;
     455    kill myRing1;
     456  }
     457  else
     458  {
     459    if(qIdeal != 0)
     460    {
     461      def r_base = basering;
     462      setring r_base;
     463      kill myRing1;
     464      qring myRing = std(qIdeal);
     465      poly p_F = imap(myRing1, p_F);
     466      ideal mIdeal = imap(myRing1, mIdeal);
     467      export p_F,mIdeal;
     468    }
     469    else
     470    {
     471      def myRing=myRing1;
     472    }
     473    kill myRing1;   
     474  }
     475 
     476  setring r_old;
     477  return(myRing);
    203478}
     479
     480////////////////////////////////////////////////////////////////////////////////
     481// returns list of coef, leadmonomial 
     482//
     483static proc determine_coef (poly Fm)
     484{
     485  def r_base = basering; // is assumed to be the result of createMyRing
     486
     487  int chara = char(basering);
     488  string charaStr;
     489  int i;
     490  string minPolyStr = "";
     491  string helpStr = "";
     492
     493  if (npars(basering) == 0)
     494  {
     495    helpStr = "ring myRing1 ="
     496              + string(chara)+ ", (y,x),ds;";
     497    execute(helpStr);
     498  }
     499  else
     500  {
     501    charaStr = charstr(basering);
     502    if (charaStr == string(chara) + "," + parstr(basering))
     503    {
     504      if (minpoly !=0)
     505      {
     506        minPolyStr = makeMinPolyString("a");
     507        helpStr = "ring myRing1 = (" + string(chara) + ",a), (y,x),ds;";
     508        execute(helpStr);
     509
     510        helpStr = "minpoly =" + minPolyStr + ";";
     511        execute (helpStr);
     512      }
     513      else // no minpoly given
     514      {
     515        helpStr = "ring myRing1 =
     516                  (" + string(chara) + ",a(1..npars(basering)) ), (y,x),ds;";
     517        execute(helpStr);
     518      }
     519    }
     520    else
     521    {
     522      i = find (charaStr,",");
     523
     524      helpStr = " ring myRing1 = (" + charaStr[1,i] + "a), (y,x),ds;";
     525      execute (helpStr);
     526    }
     527  }
     528  poly f=imap(r_base,Fm);
     529  poly g=leadmonom(f);
     530  setring r_base;
     531  poly g=imap(myRing1,g);
     532  kill myRing1;
     533  def M=coef(Fm,xy);
     534 
     535  for (i=1; i<=ncols(M); i++)
     536  {
     537    if (M[1,i]==g)
     538    {
     539      poly h=M[2,i];  // determine coefficient of leading monomial (in K[t])
     540      i=ncols(M)+1;
     541    }
     542  }
     543  return(list(h,g));
     544}
     545
    204546///////////////////////////////////////////////////////////////////////////////
    205547// RETURNS: 1, if p_f = 0 or char(basering) divides the order of p_f
    … …  
    213555  if (p_f == 0)
    214556  {
    215     dbprint(i_print,"The Input is a 'deformation'  of the zero polynomial");
     557    print("Input is a 'deformation'  of the zero polynomial!");
    216558    return(1);
    217559  }
    … …  
    221563  if (number(i_ord) == 0)
    222564  {
    223     dbprint(i_print,
    224             "The characteristic of the coefficient field
    225              divides the order of the equation");
     565    print("Characteristic of coefficient field "
     566            +"divides order of zero-fiber !");
    226567    return(1);
    227568  }
    … …  
    229570  if (squarefree(p_f) != p_f)
    230571  {
    231     dbprint(i_print, "The curve is reducible");
     572    print("Original polynomial (= zero-fiber) is not reduced!");
    232573    return(1);
    233574  }
    … …  
    235576  return(0);
    236577}
     578
     579////////////////////////////////////////////////////////////////////////////////
     580static proc make_ring_small(ideal J)
     581// returns varstr for new ring, the map and the number of vars
     582{
     583  attrib(J,"isSB",1);
     584  int counter=0;
     585  ideal newmap;
     586  string newvar="";
     587  for (int i=1; i<=nvars(basering); i++)
     588  {
     589    if (reduce(var(i),J)<>0)
     590    {
     591      newmap[i]=var(i);
     592     
     593      if (newvar=="")
     594      {
     595        newvar=newvar+string(var(i));
     596        counter=counter+1;       
     597      }
     598      else
     599      {
     600        newvar=newvar+","+string(var(i));
     601        counter=counter+1;       
     602      }
     603    }
     604    else
     605    {
     606      newmap[i]=0;
     607    } 
     608  }
     609  list L=newvar,newmap,counter;
     610  attrib(J,"isSB",0);
     611  return(L);
     612}
     613
    237614///////////////////////////////////////////////////////////////////////////////
    238 // COMPUTES  the multiplicity sequence of p_f
    239 static proc calcMultSequence (poly p_f)
     615//  The following procedure is called by esStratum (typ=0), resp. by
     616//  isEquising (typ=1)
     617///////////////////////////////////////////////////////////////////////////////
     618
     619static proc esComputation (int typ, poly p_F, list #)
    240620{
    241   int i_print = printlevel-voice + 3;
    242   intvec multSeq=0;
    243   list hneList;
    244   int xNotTransversal;
    245   int fIrreducible = 1;
    246 
    247   // if basering =  (p,a) or (p,a(1..s)),
    248   //            p prime, a algebraic, a(1..s) transcendent use reddevelop
    249   // otherwise use develop
    250   if (char(basering) != 0
    251       && npars(basering) !=0
    252       && charstr(basering) == string(char(basering)) + "," + parstr(basering))
    253   {
    254     hneList = reddevelop(p_f, -1);
    255     if ( size(hneList)>=2)
    256     {
    257       fIrreducible = 0;
    258       dbprint(i_print, "The curve is reducible");
    259       return(multSeq, xNotTransversal, fIrreducible);
    260     }
    261 
    262     hneList = hneList[1];
    263 
    264     xNotTransversal= hneList[3];
    265   }
    266   else
    267   {
    268     hneList = develop(p_f, -1);
    269 
    270     xNotTransversal= hneList[3];
    271     fIrreducible = hneList[5];
    272   }
    273 
    274   if ( ! fIrreducible)
    275   {
    276     dbprint(i_print, "The curve is reducible");
    277     return(multSeq, xNotTransversal, fIrreducible);
    278   }
    279 
    280   multSeq = multsequence (hneList);
    281   return(multSeq, xNotTransversal, fIrreducible);
    282 }
    283 ///////////////////////////////////////////////////////////////////////////////
    284 // ASSUME:  The basering is no qring and has at least 3 variables
    285 // DEFINES: A new basering, "myRing",
    286 //          with new names for the parameters and variables.
    287 //          The new names for the parameters are a(1..k),
    288 //          and t(1..s),x,y for the variables
    289 //          The ring ordering is ordStr.
    290 // NOTE:    This proc uses 'execute'.
    291 static proc createMyRing(poly p_F, string ordStr )
    292 {
    293   def r_old = basering;
    294 
    295   int chara = char(basering);
    296   string charaStr;
    297   int i;
     621  // Initialize variables
     622  int branch=1;
     623  int blowup=1;
     624  int auxVar=1;
     625  int nVars;
     626 
     627  intvec upper_bound, upper_bound_old, fertig, soll;
     628  list blowup_string;
     629  int i_print= printlevel-voice+2;
     630
     631  int no_b, number_of_branches, swapped;
     632  int i,j,k,m, counter, dummy;
     633  string helpStr = "";
     634  string ordStr = "";
     635  string MinPolyStr = "";
     636
     637  if (nvars(basering)<=2)
     638  {
     639    print("family is trivial (no deformation parameters)!");
     640    if (typ==1) //isEquising
     641    {
     642      return(1); 
     643    }
     644    else
     645    {
     646      return(list(ideal(0),0));
     647    }
     648  }
     649
     650  if (size(#)>0)
     651  {
     652    if (typeof(#[1])=="int")
     653    {
     654      def artin_bd=#[1];  // compute modulo maxideal(artin_bd)
     655      if (artin_bd <= 1)
     656      {
     657        print("Do you really want to compute over Basering/maxideal("
     658              +string(artin_bd)+") ?");
     659        print("No computation performed !");
     660        if (typ==1) //isEquising
     661        {
     662          return(1); 
     663        }
     664        else
     665        {
     666          return(list(ideal(0),int(1)));
     667        }
     668      }
     669      if (size(#)>1)
     670      {
     671        if (typeof(#[2])=="list")
     672        {
     673          def @L=#[2];  // is assumed to be the Hamburger-Noether matrix
     674        }
     675      }     
     676    }
     677    if (typeof(#[1])=="list")
     678    {
     679      def @L=#[1];  // is assumed to be the Hamburger-Noether matrix
     680    }
     681  }
     682  int ring_is_changed;
     683  def old_ring=basering;
     684  if(defined(@L)<=0)
     685  {
     686    // define a new ring without deformation-parameters and change to it:
     687    string str;
     688    string minpolyStr = string(minpoly);
     689    str = " ring HNERing = (" + charstr(basering) + "), (x,y), ls;";
     690    execute (str);
     691    str = "minpoly ="+ minpolyStr+";";
     692    execute(str);
     693    ring_is_changed=1;
     694    // Basering changed to HNERing (variables x,y, with ls ordering)
     695
     696    k=nvars(old_ring);
     697    matrix Map_Phi[1][k];
     698    Map_Phi[1,k-1]=x;
     699    Map_Phi[1,k]=y;
     700    map phi=old_ring,Map_Phi;
     701    poly f=phi(p_F);
     702
     703    // Heuristics: if x,y are transversal parameters then computation of HNE
     704    // can be much faster when exchanging variables...!
     705    if (2*size(coeffs(f,x))<size(coeffs(f,y)))
     706    {
     707      swapped=1;
     708      f=swapXY(f);
     709    }
     710   
     711    int error=checkPoly(f);
     712    if (error)
     713    {
     714      setring old_ring;
     715      if (typ==1) //isEquising
     716      {
     717        print("Return value (=0) has no meaning!");
     718        return(0); 
     719      }
     720      else
     721      {
     722        return(list( ideal(0),error));
     723      }
     724    }
     725   
     726    dbprint(i_print,"// ");
     727    dbprint(i_print,"// Compute HN expansion");
     728    dbprint(i_print,"// ---------------------");
     729    i=printlevel;
     730    printlevel=printlevel-5;
     731    list LLL=hnexpansion(f);
     732
     733    if (size(LLL)==0) { // empty list returned by hnexpansion
     734      setring old_ring;
     735      print(i_print,"Unable to compute HN expansion !");     
     736      if (typ==1) //isEquising
     737      {
     738        print("Return value (=0) has no meaning!");
     739        return(0); 
     740      }
     741      else
     742      {
     743        return(list(ideal(0),int(1)));
     744      }
     745      return(0);
     746    }
     747    else
     748    {
     749      if (typeof(LLL[1])=="ring") {
     750        def HNering = LLL[1];
     751        setring HNering;
     752        def @L=hne;
     753      }
     754      else {
     755        def @L=LLL;
     756      }
     757    }
     758    printlevel=i;
     759    dbprint(i_print,"// finished");
     760    dbprint(i_print,"// ");
     761  }
     762  def HNEring=basering;
     763  list M=multsequence(@L);
     764  M=control_Matrix(M);     // this returns the 4 control matrices 
     765  def maxDeg=M[4];
     766
     767  list L1=inf_Tangents(@L,nrows(M[1]));
     768  matrix B=L1[1];
     769  intvec V=L1[2];
     770  kill L1;
     771
     772  // if we have computed the HNE for f after swapping x and y, we have
     773  // to reinterprete the (swap) matrix V:
     774  if (swapped==1)
     775  {
     776    for (i=1;i<=size(V);i++) { V[i]=V[i]-1; } // turns 0 into -1, 1 into 0
     777  }
     778
     779  // Determine maximal number of needed auxiliary parameters (free tangents):
     780  no_b=Determine_no_b(M[3],B); 
     781
     782  // test whether HNexpansion needed field extension....
    298783  string minPolyStr = "";
    299   string helpStr;
    300   int nDefParams = nvars(r_old)-2;
    301 
    302 
    303   ideal qIdeal = ideal(basering);
    304 
    305   if (npars(basering) == 0)
    306   {
    307     helpStr = "ring myRing ="
    308               + string(chara)+ ", (t(1..nDefParams), x, y),"+ ordStr +";";
    309     execute(helpStr);
    310   }
    311 
    312   else
    313   {
    314     charaStr = charstr(basering);
    315     if (charaStr == string(chara) + "," + parstr(basering))
    316     {
     784  if (minpoly !=0)
     785  {
     786    minPolyStr = makeMinPolyString("a");
     787    minPolyStr = "minpoly =" + minPolyStr + ";";
     788  }
     789
     790  setring old_ring;
     791
     792  def myRing=createMyRing_new(p_F,"dp",minPolyStr,no_b);
     793  setring myRing;  // comes with mIdeal
     794  map hole=HNEring,mIdeal;
     795  // basering has changed to myRing, in particular, the "old"
     796  // variable names, e.g., A,B,C,z,y are replaced by t(1),t(2),t(3),x,y
     797
     798  // Initialize some variables:
     799  map phi;
     800  poly G, F_save;
     801  poly b_dummy;
     802  ideal J,Jnew,final_Map;
     803  number_of_branches=ncols(M[1]);
     804  for (i=1;i<=number_of_branches;i++)
     805  {
     806    poly F(i);
     807    ideal bl_Map(i); 
     808  }
     809  upper_bound[number_of_branches]=0;
     810  upper_bound[1]=number_of_branches;
     811  upper_bound_old=upper_bound;
     812  fertig[number_of_branches]=0;
     813  for (i=1;i<=number_of_branches;i++){ soll[i]=1; }
     814
     815  // Hole:  B = matrix of blowup points
     816  if (ring_is_changed==0) { matrix B=hole(B); }
     817  else                    { matrix B=imap(HNEring,B); }
     818  m=M[1][blowup,branch];    // multiplicity at 0
     819   
     820  // now, we start by checking equimultiplicity along trivial section
     821  poly Fm=m_Jet(p_F,m-1);
     822
     823  matrix coef_Mat = coef(Fm,xy);
     824  Jnew=coef_Mat[2,1..ncols(coef_Mat)];
     825  J=J,Jnew;
     826
     827  if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
     828                         // deformation parameters can be cutted off
     829  {
     830    J=jet(J,artin_bd-1);
     831  }
     832
     833  J=interred(J);
     834
     835// J=std(J);
     836
     837  if (typ==1) // isEquising
     838  {
     839    if(ideal(nselect(J,1,no_b))<>0)
     840    {
     841      setring old_ring;
     842      return(0);
     843    }   
     844  }
     845 
     846  F(1)=p_F;
     847
     848  // and reduce the remaining terms in F(1):
     849  bl_Map(1)=maxideal(1);
     850
     851  attrib(J,"isSB",1);
     852  bl_Map(1)=reduce(bl_Map(1),J);
     853  attrib(J,"isSB",0);
     854
     855  phi=myRing,bl_Map(1);
     856  F(1)=phi(F(1));
     857
     858  // simplify F(1)
     859  attrib(J,"isSB",1);
     860  F(1)=reduce(F(1),J);
     861  attrib(J,"isSB",0);
     862
     863  // now we compute the m-jet:
     864  Fm=m_Jet(F(1),m);
     865
     866  G=1;
     867  counter=branch;
     868  k=upper_bound[branch];
     869 
     870  F_save=F(1);  // is truncated differently in the following loop
     871 
     872  while(counter<=k)
     873  {
     874    F(counter)=m_Jet(F_save,maxDeg[blowup,counter]);
     875    if (V[counter]==0) // 2nd ring variable is tangent to this branch
     876    {
     877      G=G*(y-(b(auxVar)+B[blowup,counter])*x)^(M[3][blowup,counter]);
     878    }
     879    else // 1st ring variable is tangent to this branch
     880    {
     881      G=G*(x-(b(auxVar)+B[blowup,counter])*y)^(M[3][blowup,counter]);
     882      F(counter)=swapXY(F(counter));
     883    }
     884    bl_Map(counter)=maxideal(1);
     885    bl_Map(counter)[nvars(basering)]=xy+(b(auxVar)+B[blowup,counter])*x;
     886   
     887    auxVar=auxVar+1;
     888    upper_bound[counter]=counter+M[2][blowup+1,counter]-1;
     889    counter=counter+M[2][blowup+1,counter];
     890
     891  }
     892
     893  list LeadDataFm=determine_coef(Fm);
     894  def LeadDataG=coef(G,xy);
     895 
     896  for (i=1; i<=ncols(LeadDataG); i++)
     897  {
     898    if (LeadDataG[1,i]==LeadDataFm[2])
     899    {
     900      poly LeadG = LeadDataG[2,i];  // determine the coefficient of G
     901      i=ncols(LeadDataG)+1;
     902    }
     903  }
     904   
     905  G=LeadDataFm[1]*G-LeadG*Fm;  // leading terms in y should cancel...
     906
     907  coef_Mat = coef(G,xy);
     908  Jnew=coef_Mat[2,1..ncols(coef_Mat)];
     909
     910  if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
     911                         // deformation parameters can be cutted off
     912  {
     913    Jnew=jet(Jnew,artin_bd-1);
     914  }
     915
     916  // simplification of Jnew
     917  Jnew=interred(Jnew);
     918
     919  J=J,Jnew;
     920  if (typ==1) // isEquising
     921  {
     922    if(ideal(nselect(J,1,no_b))<>0)
     923    {
     924      setring old_ring;
     925      return(0);
     926    }   
     927  }
     928
     929
     930  while (fertig<>soll and blowup<nrows(M[3]))
     931  {
     932    upper_bound_old=upper_bound;
     933    dbprint(i_print,"// Blowup Step "+string(blowup)+" completed");
     934    blowup=blowup+1;
     935   
     936    for (branch=1;branch<=number_of_branches;branch=branch+1)
     937    {
     938      Jnew=0;
     939     
     940      // First we check if the branch still has to be considered:
     941      if (branch==upper_bound_old[branch] and fertig[branch]<>1)
     942      {
     943        if (M[3][blowup-1,branch]==1 and
     944               ((B[blowup,branch]<>x and B[blowup,branch]<>y)
     945            or (blowup==nrows(M[3])) ))
     946        {
     947          fertig[branch]=1;
     948          dbprint(i_print,"// 1 branch finished");
     949        }
     950      }
     951     
     952      if (branch<=upper_bound_old[branch] and fertig[branch]<>1)
     953      {
     954        for (i=branch;i>=1;i--)
     955        {
     956          if (M[1][blowup-1,i]<>0)
     957          {
     958            m=M[1][blowup-1,i]; // multiplicity before blowup
     959            i=0;         
     960          }
     961        }
     962         
     963        // we blow up the branch and take the strict transform:
     964        attrib(J,"isSB",1);
     965        bl_Map(branch)=reduce(bl_Map(branch),J);
     966        attrib(J,"isSB",0);
     967
     968        phi=myRing,bl_Map(branch);
     969        F(branch)=phi(F(branch))/x^m;
     970
     971        // simplify F
     972        attrib(Jnew,"isSB",1);
     973
     974        F(branch)=reduce(F(branch),Jnew);
     975        attrib(Jnew,"isSB",0);
     976
     977        m=M[1][blowup,branch]; // multiplicity after blowup
     978        Fm=m_Jet(F(branch),m); // homogeneous part of lowest degree
     979 
     980
     981        // we check for Fm=F[k]*...*F[k+s] where
     982        //
     983        //    F[j]=(y-b'(j)*x)^m(j), respectively F[j]=(-b'(j)*y+x)^m(j)
     984        //
     985        // according to the entries m(j)= M[3][blowup,j] and
     986        //                          b'(j) mod m_A = B[blowup,j]
     987        // computed from the HNE of the special fibre of the family:
     988        G=1;
     989        counter=branch;
     990        k=upper_bound[branch];
     991       
     992        F_save=F(branch);
     993       
     994        while(counter<=k)
     995        {
     996          F(counter)=m_Jet(F_save,maxDeg[blowup,counter]);
     997       
     998          if (B[blowup,counter]<>x and B[blowup,counter]<>y)
     999          {
     1000            G=G*(y-(b(auxVar)+B[blowup,counter])*x)^(M[3][blowup,counter]);
     1001            bl_Map(counter)=maxideal(1);
     1002            bl_Map(counter)[nvars(basering)]=
     1003                                        xy+(b(auxVar)+B[blowup,counter])*x;
     1004            auxVar=auxVar+1;
     1005          }
     1006          else
     1007          {
     1008            if (B[blowup,counter]==x)
     1009            {
     1010              G=G*x^(M[3][blowup,counter]);  // branch has tangent x !!
     1011              F(counter)=swapXY(F(counter)); // will turn x to y for blow up
     1012              bl_Map(counter)=maxideal(1);
     1013              bl_Map(counter)[nvars(basering)]=xy;
     1014            }
     1015            else
     1016            {
     1017              G=G*y^(M[3][blowup,counter]); // tangent has to be y
     1018              bl_Map(counter)=maxideal(1);
     1019              bl_Map(counter)[nvars(basering)]=xy;
     1020            }
     1021          }
     1022          upper_bound[counter]=counter+M[2][blowup+1,counter]-1;
     1023          counter=counter+M[2][blowup+1,counter];
     1024        }
     1025        G=determine_coef(Fm)[1]*G-Fm;  // leading terms in y should cancel
     1026        coef_Mat = coef(G,xy);
     1027        Jnew=coef_Mat[2,1..ncols(coef_Mat)];
     1028        if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
     1029                              // deformation parameters can be cutted off
     1030        {
     1031          Jnew=jet(Jnew,artin_bd-1);
     1032        }
     1033
     1034        // simplification of J
     1035        Jnew=interred(Jnew);
     1036
     1037        J=J,Jnew;
     1038        if (typ==1) // isEquising
     1039        {
     1040          if(ideal(nselect(J,1,no_b))<>0)
     1041          {
     1042            setring old_ring;
     1043            return(0);
     1044          }   
     1045        }
     1046
     1047      }
     1048    }
     1049    if (number_of_branches>=2)
     1050    {
     1051      J=interred(J);
     1052      if (typ==1) // isEquising
     1053      {
     1054        if(ideal(nselect(J,1,no_b))<>0)
     1055        {
     1056          setring old_ring;
     1057          return(0);
     1058        }   
     1059      }
     1060    }
     1061  }
     1062 
     1063  // computation for all equimultiple sections being trivial (I^s(f))
     1064  ideal Jtriv=J;
     1065  for (i=1;i<=no_b; i++)
     1066  {
     1067    Jtriv=subst(Jtriv,b(i),0);
     1068  }
     1069  Jtriv=std(Jtriv);
     1070
     1071  dbprint(i_print,"// ");
     1072  dbprint(i_print,"// Elimination starts:");
     1073  dbprint(i_print,"// -------------------");
     1074
     1075  poly gg;
     1076  int b_left=no_b; 
     1077
     1078  for (i=1;i<=no_b; i++)
     1079  {
     1080    attrib(J,"isSB",1); 
     1081    gg=reduce(b(i),J);
     1082    if (gg==0)
     1083    {
     1084      b_left = b_left-1;  // another b(i) has to be 0
     1085    }                 
     1086    J = subst(J, b(i), gg);
     1087    attrib(J,"isSB",0);
     1088  }
     1089  J=simplify(J,10);
     1090  if (typ==1) // isEquising
     1091  {
     1092    if(ideal(nselect(J,1,no_b))<>0)
     1093    {
     1094      setring old_ring;
     1095      return(0);
     1096    }   
     1097  }
     1098
     1099  ideal J_no_b = nselect(J,1,no_b);
     1100  if (size(J) > size(J_no_b))
     1101  {
     1102    dbprint(i_print,"// std computation started");
     1103    // some b(i) didn't appear in linear conditions and have to be eliminated
     1104    if (defined(artin_bd))
     1105    {
     1106      // first we make the ring smaller (removing variables, which are
     1107      // forced to 0 by J
     1108      list LL=make_ring_small(J);
     1109      ideal Shortmap=LL[2];
     1110      minPolyStr = "";
    3171111      if (minpoly !=0)
    3181112      {
    319         minPolyStr = makeMinPolyString("a");
    320 
    321         helpStr = "ring myRing =
    322                  (" + string(chara) + ",a),
    323                  (t(1..nDefParams), x, y)," + ordStr + ";";
    324         execute(helpStr);
    325 
    326         helpStr = "minpoly =" + minPolyStr + ";";
    327         execute (helpStr);
    328       }
    329       else
     1113        minPolyStr = "minpoly = "+string(minpoly);
     1114      }
     1115      ordStr = "dp(" + string(b_left) + "),dp";
     1116      ideal qId = ideal(basering);
     1117   
     1118      helpStr = "ring Shortring = ("
     1119                + charstr(basering) + "),("+ LL[1] +") , ("+ ordStr  +");";
     1120      execute(helpStr);
     1121      execute(minPolyStr);
     1122      // ring has changed to "Shortring"
     1123
     1124      ideal MM=maxideal(artin_bd);
     1125      MM=subst(MM,x,0);
     1126      MM=subst(MM,y,0);
     1127      MM=simplify(MM,2);
     1128      dbprint(i_print-1,"// maxideal("+string(artin_bd)+") has "
     1129                         +string(size(MM))+" elements");
     1130      dbprint(i_print-1,"//");
     1131
     1132      // we change to the qring mod m^artin_bd
     1133      // first, we have to check if we were in a qring when starting
     1134      ideal qId = imap(myRing, qId); 
     1135      if (qId == 0)
    3301136      {
    331         helpStr = "ring myRing =
    332                   (" + string(chara) + ",a(1..npars(basering)) ),
    333                   (t(1..nDefParams), x, y)," + ordStr + ";";
    334         execute(helpStr);
    335       }
     1137         attrib(MM,"isSB",1);
     1138         qring QQ=MM;     
     1139      }
     1140        else
     1141        {   
     1142           qId=qId,MM;
     1143           qring QQ = std(qId);
     1144        }
     1145
     1146      ideal Shortmap=imap(myRing,Shortmap);
     1147      map phiphi=myRing,Shortmap;
     1148
     1149      ideal J=phiphi(J);
     1150      J=std(J);
     1151      J=nselect(J,1,no_b);
     1152
     1153      setring myRing;
     1154      // back to "myRing"
     1155
     1156      J=nselect(J,1,no_b);
     1157      Jnew=imap(QQ,J);
     1158
     1159      J=J,Jnew;
     1160      J=interred(J);
    3361161    }
    3371162    else
    3381163    {
    339       i = find (charaStr,",");
    340 
    341       helpStr = " ring myRing =
    342               (" + charaStr[1,i] + "a),
    343               (t(1..nDefParams), x, y)," + ordStr + ";";
    344 
    345       execute (helpStr);
    346     }
    347   }
    348 
    349   ideal qIdeal = fetch(r_old, qIdeal);
    350 
    351   if(qIdeal != 0)
    352   {
    353     def r_base = basering;
    354     kill myRing;
    355     qring myRing = std(qIdeal);
    356   }
    357 
    358   poly p_F = fetch(r_old, p_F);
    359   ideal ES;
    360 
    361   keepring(myRing);
     1164      J=std(J);
     1165      J=nselect(J,1,no_b);
     1166    }
     1167  }
     1168 
     1169  dbprint(i_print,"// finished");
     1170  dbprint(i_print,"// ");
     1171 
     1172  minPolyStr = "";
     1173  if (minpoly !=0)
     1174  {
     1175   minPolyStr = "minpoly = "+string(minpoly);
     1176  }
     1177
     1178  kill HNEring;
     1179 
     1180  if (typ==1) // isEquising
     1181  {
     1182    if(J<>0)
     1183    {
     1184      setring old_ring;
     1185      return(0);
     1186    }   
     1187    else
     1188    {
     1189      setring old_ring;
     1190      return(1);
     1191    }   
     1192  }
     1193 
     1194  setring old_ring;
     1195  // we are back in the original ring
     1196
     1197  if (npars(myRing)<>0)
     1198  {
     1199    ideal qIdeal = ideal(basering);
     1200    helpStr = "ring ESSring = ("
     1201                 + string(char(basering))+ "," + parstr(myRing) +
     1202                 ") , ("+ varstr(basering)+") , ("+ ordstr(basering) +");";
     1203    execute(helpStr);
     1204    execute(minPolyStr); 
     1205    // basering has changed to ESSring
     1206
     1207    ideal qIdeal = fetch(old_ring, qIdeal);
     1208    if(qIdeal != 0)
     1209    {
     1210      def r_base = basering;
     1211      kill ESSring;
     1212      qring ESSring = std(qIdeal);
     1213    }
     1214    kill qIdeal;
     1215 
     1216    ideal SSS;
     1217    for (int ii=1;ii<=nvars(basering);ii++)
     1218    {
     1219      SSS[ii+no_b]=var(ii);
     1220    }
     1221    map phi=myRing,SSS;   // b(i) variables are mapped to zero
     1222
     1223    ideal ES=phi(J);
     1224    ideal ES_all_triv=phi(Jtriv);
     1225    kill phi;
     1226 
     1227    if (defined(p_F)<=0)
     1228    {
     1229      poly p_F=fetch(old_ring,p_F);
     1230      export(p_F);
     1231    }
     1232    export(ES);
     1233    export(ES_all_triv);
     1234    setring old_ring;
     1235    dbprint(i_print+1,"
     1236// 'esStratum' created a list M of a ring and an integer.
     1237// To access the ideal defining the equisingularity stratum, type:
     1238        def ESSring = M[1]; setring ESSring;  ES; ");
     1239
     1240    return(list(ESSring,0));
     1241  }
     1242  else
     1243  {
     1244    // no new ring definition necessary
     1245    ideal SSS;
     1246    for (int ii=1;ii<=nvars(basering);ii++)
     1247    {
     1248      SSS[ii+no_b]=var(ii);
     1249    }
     1250    map phi=myRing,SSS;  // b(i) variables are mapped to zero
     1251
     1252    ideal ES=phi(J);
     1253    ideal ES_all_triv=phi(Jtriv);
     1254    kill phi;
     1255
     1256    setring old_ring;
     1257    dbprint(i_print,"// output of 'esStratum' is a list consisting of:
     1258//    _[1][1] = ideal defining the equisingularity stratum
     1259//    _[1][2] = ideal defining the part of the equisingularity stratum
     1260//              where all equimultiple sections are trivial
     1261//    _[2] = 0");
     1262    return(list(list(ES,ES_all_triv),0));
     1263  }
     1264 
    3621265}
    3631266
    364 ///////////////////////////////////////////////////////////////////////////////
    365 /////////// procedures to compute the equisingularity stratum /////////////////
    366 ///////////////////////////////////////////////////////////////////////////////
    367 // DEFINES a new basering, myRing,  which has one variable
    368 //         more than the old ring.
    369 //         The name for the new variable is "H(nhelpV)".
    370 //         p_F and ES are "imaped" into the new ring.
    371 static proc extendRing (poly p_F,  ideal ES, int nHelpV, ideal HCond)
    372 {
    373   def r_old = basering;
    374 
    375   string helpStr;
    376   string minPolyStr = "";
    377 
    378   ideal qIdeal = ideal(basering);
    379 
    380   if (minpoly != 0)
    381   {
    382     if (charstr(basering) == string(char(basering)) + "," + parstr(basering))
    383     {
    384       minPolyStr = string(minpoly);
    385     }
    386   }
    387 
    388   string str =
    389  "ring myRing = (" + charstr(r_old) + "),
    390                 (H(" + string(nHelpV)+ ")," + string(maxideal(1)) + "),
    391                 (dp(" + string(nHelpV) + "),dp);";
    392   execute (str);
    393 
    394   if (minPolyStr != "")
    395   {
    396     helpStr = "minpoly =" + minPolyStr + ";";
    397     execute(helpStr);
    398   }
    399 
    400   ideal qIdeal = imap(r_old, qIdeal);
    401   if(qIdeal != 0)
    402   {
    403     def r_base = basering;
    404 
    405     kill myRing;
    406 
    407     attrib(qIdeal,"isSB",1);
    408     qring myRing = qIdeal;
    409   }
    410 
    411 
    412   poly p_F = imap(r_old, p_F);
    413   ideal ES = imap(r_old, ES);
    414   ideal HCond = imap(r_old, HCond);
    415 
    416   keepring(myRing);
    417 }
    418 ///////////////////////////////////////////////////////////////////////////////
    419 // COMPUTES an ideal equimultCond, such that F_(n-1) mod equimultCond =0,
    420 //          where F_(n-1) is the (nNew-1)-jet of p_F with respect to x,y.
    421 static proc calcEquimultCond(poly p_F, int nNew)
    422 {
    423   ideal equimultCond = 0;
    424   poly p_FnMinus1;
    425   matrix coefMatrix;
    426   int nc;
    427   int ii = 1;
    428 
    429   p_FnMinus1 = jet(p_F, nNew-1, xyVector());
    430 
    431   coefMatrix = coef(p_FnMinus1, xy);
    432 
    433   nc = ncols(coefMatrix);
    434 
    435   for (ii=1; ii<=nc; ii++)
    436   {
    437     equimultCond[ii] = NF(coefMatrix[2,ii],std(0));
    438   }
    439 
    440   p_F = p_F - p_FnMinus1;
    441   p_F = SimplifyF(p_F, equimultCond);
    442 
    443   return(equimultCond, p_F);
    444 }
    445 ///////////////////////////////////////////////////////////////////////////////
    446 // COMPUTES smallest integer >= nNew/nOld -1
    447 static proc calcNZeroSteps (int nOld,int nNew)
    448 {
    449  int nZeroSteps;
    450 
    451  if (nOld mod nNew == 0)
    452  {
    453    nZeroSteps = nOld div nNew -1;
    454  }
    455 
    456  else
    457  {
    458   nZeroSteps  = nOld div nNew;
    459  }
    460 
    461  return(nZeroSteps);
    462 }
    463 ///////////////////////////////////////////////////////////////////////////////
    464 // ASSUME: ord_(X,Y)(F)=nNew
    465 // COMPUTES an ideal I such that (p_F mod I)_nNew  = p_c*y^nNew.
    466 static proc purePowerOfY (poly p_F, int nNew)
    467 {
    468   ideal id_help = 0;
    469   poly  p_Fn;
    470   matrix coefMatrix;
    471   int nc;
    472   poly p_c;
    473   int ii=1;
    474 
    475   p_Fn = jet(p_F, nNew, xyVector());
    476 
    477   coefMatrix = coef(p_Fn, xy);
    478   nc = ncols(coefMatrix);
    479 
    480   p_c = coefMatrix[2,nc];
    481 
    482   for (ii=1; ii <= nc-1; ii++)
    483   {
    484     id_help[ii] = NF(coefMatrix[2,ii],std(0));
    485   }
    486 
    487   p_F = p_F - p_Fn + p_c*y^nNew;
    488 
    489   p_F = SimplifyF(p_F, id_help);
    490 
    491   return(id_help, p_F, p_c);
    492 }
    493 ///////////////////////////////////////////////////////////////////////////////
    494 // ASSUME: ord_(X,Y)(F)=nNew
    495 // COMPUTES an ideal I such that p_Fn mod I = p_c*(y-p_a*x)^nNew,
    496 //          where p_Fn is the homogeneous part of p_F of order nNew.
    497 static proc purePowerOfLin (poly p_F, ideal HCond, int nNew, int nHelpV)
    498 {
    499   ideal id_help = 0;
    500   poly p_Fn;
    501   matrix coefMatrix;
    502   poly p_c;
    503   poly p_ca;
    504   poly p_help;
    505   poly p_a;
    506   int ii;
    507   int jj;
    508   int bico;
    509 
    510   p_Fn = jet(p_F, nNew, xyVector());
    511 
    512   coefMatrix = coeffs(subst(p_Fn,x,1),y);
    513 
    514   p_c = coefMatrix[nNew+1,1];
    515   p_ca = coefMatrix[nNew,1]/(-nNew);
    516 
    517   if (npars(basering)==1 && charstr(basering) != string(char(basering)) + "," + parstr(basering))
    518   {
    519     p_a = H(nHelpV);
    520     HCond = HCond + ideal(p_ca - p_a*p_c);
    521   }
    522   else
    523   {
    524     p_help = p_ca/p_c;
    525     if (p_help * p_c == p_ca)
    526     {
    527       p_a = p_help;
    528     }
    529     else
    530     {
    531       p_a = H(nHelpV);
    532       HCond = HCond + ideal(p_ca - p_a*p_c);
    533     }
    534   }
    535 
    536   bico = (nNew*(nNew-1))/2;
    537 
    538   for (ii = 2; ii <= nNew ; ii++)
    539   {
    540 
    541     if (coefMatrix[nNew+1-ii,1] == 0)
    542     {
    543       if (number(bico) != 0)
    544       // Then a^i=0 since c is a unit
    545       {
    546         id_help = id_help + ideal(NF(p_a^(ii),std(0)));
    547         for (jj = ii+1; jj <= nNew; jj++)
    548         // the remaining coefficients (of y^(nnew-k)*x^k with k>i)
    549         // are also zero.
    550         {
    551           id_help = id_help
    552                     + ideal(NF(coefMatrix[nNew+1-jj,1],std(0)));
    553         }
    554         break;
    555       }
    556     }
    557 
    558     else
    559     {
    560       id_help = id_help +
    561             ideal(NF(coefMatrix[nNew+1-ii,1] - bico*p_c*(-p_a)^ii,std(0)));
    562     }
    563 
    564     bico = (bico*(nNew-ii))/(ii+1);
    565   }
    566 
    567   p_F = SimplifyF(p_F, id_help);
    568 
    569   return(id_help, HCond, p_F, p_c, p_a);
    570 }
    571 ///////////////////////////////////////////////////////////////////////////////
    572 // eliminates the variables H(1),..,H(nHelpV) from the ideal ES + HCond
    573 static proc helpVarElim(ideal ES, ideal HCond, int nHelpV)
    574 {
    575   ES = ES + HCond;
    576   ES = std(ES);
    577   ES = nselect(ES,1,nHelpV);
    578 
    579   return(ES);
    580 }
    581 ///////////////////////////////////////////////////////////////////////////////
    582 // ASSUME:   ord(F)=nNew and p_c(y-p_a*x)^n is the nNew-jet of F with respect
    583 //           to X,Y
    584 // COMPUTES  F(x,yx+a*x)/x^n
    585 static proc formalBlowUp(poly p_F, poly p_c, poly p_a, int nNew)
    586 {
    587 
    588   p_F = p_F - jet(p_F, nNew, xyVector());
    589 
    590   if (p_a != 0)
    591   {
    592     p_F = eSubst(p_F, y , yx + p_a*x);
    593   }
    594   else
    595   {
    596     p_F = subst(p_F, y, xy);
    597   }
    598 
    599   p_F = p_F/(x^nNew);
    600 
    601   p_F = p_F + p_c * y^nNew;
    602 
    603   return (p_F);
    604 }
    605 ///////////////////////////////////////////////////////////////////////////////
    606 // ASSUME:  p_F in K[t(1)..t(s),x,y]
    607 // COMPUTES the minimal ideal ES, such that the deformation p_F mod ES is
    608 //          equisingular.
    609 //          The computation is done up to iteration step 'maxstep'.
    610 // RETURNS: list l, such that
    611 //          l[1]=1 if some error has occured,
    612 //          l[1]=0 otherwise and then l[2] = es_cond.
    613 static proc calcEsCond(poly p_F, intvec multSeq, int maxStep)
    614 {
    615   def r_old = basering;
    616 
    617   ideal ES = 0;
    618 
    619   int ii;
    620   int step = 1;
    621   int nNew = multSeq[step];
    622   int nOld = nNew;
    623   int nZeroSteps;
    624   int nHelpV = 1;
    625   ideal HCond = 0;
    626   int maxDeg = 0;
    627   int z = printlevel - voice + 3;
    628   string str;
    629 
    630   extendRing(p_F, ES, nHelpV, HCond);
    631 
    632   poly p_c;
    633   poly p_a;
    634   ideal I;
    635 
    636   for (ii = 1; ii <= size(multSeq); ii++)
    637   {
    638     maxDeg = maxDeg + multSeq[ii];
    639   }
    640 
    641   while (step <= maxStep)
    642   {
    643 
    644     nOld = nNew;
    645     nNew = multSeq[step];
    646 
    647     p_F = jet(p_F, maxDeg, xyVector());
    648     maxDeg = maxDeg - nNew;
    649 
    650     if (nNew<nOld)
    651     //start a new line in the HNE of F
    652     //               _            _
    653     // for the next | nold/nnew -1 | iteration steps the coefficient 'a'
    654     // in the leading form Fn = c(y-ax)^n should  be zero.
    655     {
    656       p_F = swapXY(p_F);
    657       nZeroSteps = calcNZeroSteps (nOld, nNew);
    658     }
    659 
    660     I, p_F = calcEquimultCond (p_F, nNew);
    661     ES = ES + I;
    662 
    663     if (z>1)
    664     {
    665       "// conditions for equimultiplicity in step:", step;
    666       I;
    667       if (nHelpV >1)
    668       {
    669         str = string(nHelpV);
    670         "// conditions for help variables H(1),..,H("+str+"):";
    671         HCond;
    672       }
    673       pause("press <return> to continue");
    674     }
    675 
    676     if (nZeroSteps > 0)
    677     {
    678       nZeroSteps--;
    679 
    680       // compute conditions, s.th. Fn = c*y^nnew ?
    681       I, p_F, p_c = purePowerOfY (p_F, nNew);
    682       ES = ES + I;
    683 
    684       if (z>1)
    685       {
    686         "// conditions for pure power in step:", step;
    687         I;
    688         if (nHelpV > 1)
    689         {
    690           str = string(nHelpV);
    691           "// conditions for help variables H(1),..,H("+str+"):";
    692           HCond;
    693         }
    694         pause("press <return> to continue");
    695       }
    696       p_a =0;
    697     }
    698 
    699     else
    700     {
    701       I, HCond, p_F, p_c, p_a = purePowerOfLin (p_F, HCond, nNew, nHelpV);
    702 
    703       ES = ES + I;
    704 
    705       if (z>1)
    706       {
    707         "// conditions for pure power in step:", step;
    708         I;
    709         str = string(nHelpV);
    710         "// conditions for help variables H(1),..,H("+str+"):";
    711         HCond;
    712         pause("press <return> to continue");
    713       }
    714     }
    715 
    716     p_F = formalBlowUp (p_F, p_c, p_a, nNew);
    717 
    718     if (p_a == H(nHelpV))
    719     {
    720       nHelpV++;
    721 
    722       def r_base = basering;
    723       kill myRing;
    724 
    725       extendRing(p_F, ES, nHelpV, HCond);
    726 
    727       kill r_base;
    728 
    729       poly p_a;
    730       poly p_c;
    731       ideal I;
    732     }
    733     step++;
    734   }
    735   if (nHelpV > 1)
    736   {
    737     ES = helpVarElim(ES, HCond, nHelpV);
    738   }
    739 
    740   if (nameof(basering)=="myRing")
    741   {
    742     setring r_old;
    743     ES = imap(myRing, ES);
    744   }
    745 
    746   return(ES);
    747 }
    748 
    749 ///////////////////////////////////////////////////////////////////////////////
    750 //                           main procedure
    751 ///////////////////////////////////////////////////////////////////////////////
    752 
    753 proc esStratum (poly p_F, list #)
    754 "USAGE:   esStratum(F[,m]); F poly, m int
    755 ASSUME:  F defines a deformation of an irreducible bivariate polynomial f
    756          and the characteristic of the basering does not divide mult(f). @*
    757          If nv is the number of variables of the basering, then the first nv-2
    758          variables are the deformation parameters. @*
    759          If the basering is a qring, ideal(basering) must only depend
    760          on the deformation parameters.
    761 RETURN:  list l of an ideal and an integer, where
    762 @format
    763   l[1] is the ideal in the deformation parameters, defining the ES-stratum of F,
    764   l[2]=1 if some error has occured,  l[2]=0 otherwise.
    765 @end format
    766 NOTE:    If m is given, the computation stops after m steps of the iteration. @*
    767          printlevel > 0 displays comments and pauses after intermediate
    768          computations (default: printlevel=0) @*
    769          This procedure uses @code{execute} or calls a procedure using
    770          @code{execute}.
    771 EXAMPLE: example esStratum; shows an example
     1267////////////////////////////////////////////////////////////////////////////////
     1268
     1269proc tau_es (poly f,list #)
     1270"USAGE:   tau_es(f); f poly
     1271ASSUME:  f is a reduced bivariate polynomial, the basering has precisely
     1272         two variables, is local and no qring.
     1273RETURN:  int, the codimension of the mu-const stratum in the semi-universal
     1274         deformation base.
     1275NOTE:    printlevel>=1 displays additional information.
     1276         When called with any additional parameter, the computation of the
     1277         Milnor number is avoided (no check for NND).
     1278SEE ALSO: esIdeal, tjurina, invariants
     1279EXAMPLE: example tau_es; shows an example.
    7721280"
    7731281{
    774   def r_user = basering;
    775 
    776   int ii = 1;
    777   int i_nvars = nvars(basering);
    778   int error = 0;
    779   int xNotTransversal;
    780   int fIrreducible;
    781   int maxStep;
    782   int userMaxStep;
    783   ideal cond;
    784   intvec multSeq;
    785   ideal ES = 0;
    786 
    787   error = checkBasering();
    788   if (error)
    789   {
    790     return(list(ES,error));
    791   }
    792 
    793   userMaxStep = getInput(#);
    794 
    795   // define a new basering "myRing" with new names for parameters
    796   // and variables.
    797   // The new names are 'a(1)', ..., 'a(npars)' for the parameters
    798   // and 't(1)', ..., 't(nvars-2)', 'x', 'y' for the variables.
    799   createMyRing(p_F, "dp");
    800 
    801   // define a ring without deformation parameters, to compute the HNE
    802   // of F mod <t_1,..,t_s>
    803   createHNERing();
    804 
    805   poly p_f = imap(myRing,p_F);
    806 
    807   error = checkPoly(p_f);
    808   if (error)
    809   {
    810     setring r_user;
    811     return(list( ideal(0),error));
    812   }
    813 
    814   // compute the multiplicitysequence of p_f.
    815   multSeq, xNotTransversal, fIrreducible = calcMultSequence(p_f);
    816 
    817   if ( ! fIrreducible)
    818   {
    819     setring r_user;
    820     return(list(ideal(0),1));
    821   }
    822 
    823   setring myRing;
    824 
    825   if (xNotTransversal)
    826   {
    827     p_F = swapXY(p_F);
    828   }
    829 
    830   if (userMaxStep != -1 && userMaxStep <  size(multSeq)-1)
    831   {
    832     maxStep = userMaxStep;
    833   }
     1282  int i,j,k,s;
     1283  int slope_x, slope_y, upper;
     1284  int i_print = printlevel - voice + 3;
     1285  string MinPolyStr;
     1286
     1287  // some checks first
     1288  if ( nvars(basering)<>2 )
     1289  {
     1290    print("// basering has not the correct number (two) of variables !");
     1291    print("// computation stopped");
     1292    return(0);
     1293  }
     1294  if ( mult(std(1+var(1)+var(2))) <> 0)
     1295  {
     1296    print("// basering is not local !");
     1297    print("// computation stopped");
     1298    return(0);
     1299  } 
     1300
     1301  if (mult(std(f))<=1)
     1302  {
     1303    // f is rigid
     1304    return(0);
     1305  }
     1306
     1307  if ( deg(squarefree(f))!=deg(f) )
     1308  {
     1309    print("// input polynomial was not reduced");
     1310    print("// try    squarefree(f);   first");
     1311    return(0);
     1312  }
     1313
     1314  def old_ring=basering;
     1315  execute("ring @myRing=("+charstr(basering)+"),("+varstr(basering)+"),ds;");
     1316  poly f=imap(old_ring,f);
     1317
     1318  ideal Jacobi_Id = jacob(f);
     1319
     1320  // check for A_k singularity
     1321  // ----------------------------------------
     1322  if (mult(std(f))==2)
     1323  {
     1324    dbprint(i_print-1,"// ");
     1325    dbprint(i_print-1,"// polynomial defined A_k singularity");
     1326    dbprint(i_print-1,"// ");
     1327    return( vdim(std(Jacobi_Id)) );
     1328  }
     1329
     1330  // check for D_k singularity
     1331  // ----------------------------------------
     1332  if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3)
     1333  {
     1334    dbprint(i_print,"// ");
     1335    dbprint(i_print,"// polynomial defined D_k singularity");
     1336    dbprint(i_print,"// ");   
     1337    ideal ES_Id = f, jacob(f);
     1338    return( vdim(std(Jacobi_Id)));
     1339  }
     1340
     1341
     1342  if (size(#)==0)
     1343  {
     1344    // check if Newton polygon non-degenerate
     1345    // ----------------------------------------
     1346    Jacobi_Id=std(Jacobi_Id);
     1347    int mu = vdim(Jacobi_Id);
     1348    poly f_tilde=f+var(1)^mu+var(2)^mu;  //to obtain convenient Newton-polygon
     1349   
     1350    list NP=newtonpoly(f_tilde);
     1351    dbprint(i_print-1,"// Newton polygon:");
     1352    dbprint(i_print-1,NP);
     1353    dbprint(i_print-1,"");
     1354
     1355    if(is_NND(f,mu,NP))          // f is Newton non-degenerate
     1356    {
     1357      upper=NP[1][2];
     1358      ideal ES_Id= x^k*y^upper;
     1359      dbprint(i_print-1,"polynomial is Newton non-degenerate");
     1360      dbprint(i_print-1,"");
     1361      k=0;
     1362      for (i=1;i<=size(NP)-1;i++)
     1363      {
     1364        slope_x=NP[i+1][1]-NP[i][1];
     1365        slope_y=NP[i][2]-NP[i+1][2];
     1366        for (k=NP[i][1]+1; k<=NP[i+1][1]; k++)
     1367        {
     1368          while ( slope_x*upper + slope_y*k >=
     1369                  slope_x*NP[i][2] + slope_y*NP[i][1])
     1370          {
     1371            upper=upper-1;
     1372          }
     1373          upper=upper+1;
     1374          ES_Id=ES_Id, x^k*y^upper;
     1375        }
     1376      }
     1377      ES_Id=std(ES_Id);
     1378      dbprint(i_print-2,"ideal of monomials above Newton bd. is generated by:");
     1379      dbprint(i_print-2,ES_Id);
     1380      ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
     1381      ES_Id = ES_Id, Jacobi_Id;
     1382      ES_Id = std(ES_Id);
     1383      dbprint(i_print-1,"// ");
     1384      dbprint(i_print-1,"// Equisingularity ideal is computed!");
     1385      dbprint(i_print-1,"");
     1386      return(vdim(ES_Id));
     1387    }
     1388    else
     1389    {
     1390      dbprint(i_print-1,"polynomial is Newton degenerate !");
     1391      dbprint(i_print-1,"");
     1392    }
     1393  }
     1394
     1395  // for Newton degenerate polynomials, we compute the HN expansion, and
     1396  // count the number of free points .....
     1397 
     1398  dbprint(i_print-1,"// ");
     1399  dbprint(i_print-1,"// Compute HN expansion");
     1400  dbprint(i_print-1,"// ---------------------");
     1401  i=printlevel;
     1402  printlevel=printlevel-5;
     1403  if (2*size(coeffs(f,x))<size(coeffs(f,y)))
     1404  {
     1405    f=swapXY(f);
     1406  }   
     1407  list LLL=hnexpansion(f);
     1408  if (size(LLL)==0) { // empty list returned by hnexpansion
     1409    setring old_ring;
     1410    ERROR("Unable to compute HN expansion !");
     1411  }   
    8341412  else
    8351413  {
    836     maxStep = size(multSeq)-1;
    837   }
    838 
    839   ES = calcEsCond(p_F, multSeq, maxStep);
    840 
    841   setring r_user;
    842   ES = fetch(myRing, ES);
    843 
    844   return(list(ES, error));
     1414    if (typeof(LLL[1])=="ring") {
     1415      def HNering = LLL[1];
     1416      setring HNering;
     1417      def @L=hne;
     1418    }
     1419    else {
     1420      def @L=LLL;
     1421    }
     1422  }
     1423  def HNEring=basering;
     1424
     1425  printlevel=i;
     1426  dbprint(i_print-1,"// finished");
     1427  dbprint(i_print-1,"// ");
     1428
     1429  list M=multsequence(@L);
     1430  M=control_Matrix(M);     // this returns the 4 control matrices
     1431  intmat Mult=M[1];
     1432 
     1433  list L1=inf_Tangents(@L,nrows(M[1]));
     1434  matrix B=L1[1];
     1435 
     1436  // determine sum_i m_i(m_i+1)/2 (over inf. near points)
     1437  int conditions=0;
     1438  for (i=1;i<=nrows(Mult);i++)
     1439  {
     1440    for (j=1;j<=ncols(Mult);j++)
     1441    {
     1442      conditions=conditions+(Mult[i,j]*(Mult[i,j]+1)/2);
     1443    }
     1444  }
     1445  int freePts=no_freePoints(M[1],B);
     1446  int taues=conditions-freePts-2;
     1447 
     1448  setring old_ring;
     1449  return(taues);
    8451450}
    846 
    8471451example
    8481452{
    8491453   "EXAMPLE:"; echo=2;
    850    ring r = 11,(a,b,c,d,e,f,g,x,y),ds;
     1454   ring r=32003,(x,y),ds;
     1455   poly f=(x4-y4)^2-x10;
     1456   tau_es(f);
     1457}
     1458
     1459
     1460////////////////////////////////////////////////////////////////////////////////
     1461
     1462proc esIdeal (poly f,list #)
     1463"USAGE:   esIdeal(f[,any]]); f poly
     1464ASSUME:  f is a reduced bivariate polynomial, the basering has precisely
     1465         two variables, is local and no qring, and the characteristic of
     1466         the ground field does not divide mult(f).
     1467RETURN:  if called with only one parameter: list of two ideals,
     1468@format
     1469          _[1]:  equisingularity ideal of f (in sense of Wahl),
     1470          _[2]:  ideal of equisingularity with fixed position of the
     1471                 singularity;
     1472@end format
     1473         if called with more than one parameter: list of three ideals,
     1474@format
     1475          _[1]:  equisingularity ideal of f (in sense of Wahl)
     1476          _[2]:  ideal of equisingularity with fixed position of the
     1477                 singularity;
     1478          _[3]:  ideal of all g such that the deformation defined by f+eg
     1479                 (e^2=0) is isomorphic to an equisingular deformation
     1480                 of V(f) with all equimultiple sections being trivial.
     1481@end format
     1482NOTE:    if some of the above condition is not satisfied then return
     1483         value is list(0,0).
     1484SEE ALSO: tau_es, esStratum
     1485KEYWORDS: equisingularity ideal
     1486EXAMPLE: example esIdeal; shows examples.
     1487"
     1488{
     1489
     1490  int typ;
     1491  if (size(#)>0) { typ=1; }   // I^s is also computed
     1492  int i,k,s;
     1493  int slope_x, slope_y, upper;
     1494  int i_print = printlevel - voice + 3;
     1495  string MinPolyStr;
     1496
     1497  // some checks first
     1498  if ( nvars(basering)<>2 )
     1499  {
     1500    print("// basering has not the correct number (two) of variables !");
     1501    print("// computation stopped");
     1502    return(list(0,0));
     1503  }
     1504  if ( mult(std(1+var(1)+var(2))) <> 0)
     1505  {
     1506    print("// basering is not local !");
     1507    print("// computation stopped");
     1508    return(list(0,0));
     1509  } 
     1510
     1511  if (mult(std(f))<=1)
     1512  {
     1513    // f is rigid
     1514    if (typ==0)
     1515    {
     1516      return(list(ideal(1),ideal(1)));
     1517    }
     1518    else
     1519    {
     1520      return(list(ideal(1),ideal(1),ideal(1)));
     1521    }
     1522  }
     1523
     1524  if ( deg(squarefree(f))!=deg(f) )
     1525  {
     1526    print("// input polynomial was not squarefree");
     1527    print("// try    squarefree(f);   first");
     1528    return(list(0,0));
     1529  }
     1530
     1531  if (char(basering)<>0)
     1532  {
     1533    if (mult(std(f)) mod char(basering)==0)
     1534    {
     1535      print("// characteristic of ground field divides "
     1536            + "multiplicity of polynomial !");
     1537      print("// computation stopped");
     1538      return(list(0,0));   
     1539    }
     1540  }
     1541 
     1542  // check for A_k singularity
     1543  // ----------------------------------------
     1544  if (mult(std(f))==2)
     1545  {
     1546    dbprint(i_print,"// ");
     1547    dbprint(i_print,"// polynomial defined A_k singularity");
     1548    dbprint(i_print,"// ");   
     1549    ideal ES_Id = f, jacob(f);
     1550    ES_Id = interred(ES_Id);
     1551    ideal ESfix_Id = f, maxideal(1)*jacob(f);
     1552    ESfix_Id= interred(ESfix_Id);
     1553    if (typ==0) // only for computation of I^es and I^es_fix
     1554    {
     1555      return( list(ES_Id,ESfix_Id) );
     1556    }
     1557    else
     1558    {
     1559      return( list(ES_Id,ESfix_Id,ES_Id) );
     1560    }
     1561  }
     1562
     1563  // check for D_k singularity
     1564  // ----------------------------------------
     1565  if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3)
     1566  {
     1567    dbprint(i_print,"// ");
     1568    dbprint(i_print,"// polynomial defined D_k singularity");
     1569    dbprint(i_print,"// ");   
     1570    ideal ES_Id = f, jacob(f);
     1571    ES_Id = interred(ES_Id);
     1572    ideal ESfix_Id = f, maxideal(1)*jacob(f);
     1573    ESfix_Id= interred(ESfix_Id);
     1574    if (typ==0) // only for computation of I^es and I^es_fix
     1575    {
     1576      return( list(ES_Id,ESfix_Id) );
     1577    }
     1578    else
     1579    {
     1580      return( list(ES_Id,ESfix_Id,ES_Id) );
     1581    }
     1582  }
     1583
     1584  // check if Newton polygon non-degenerate
     1585  // ----------------------------------------
     1586  int mu = milnor(f);
     1587  poly f_tilde=f+var(1)^mu+var(2)^mu;  //to obtain a convenient Newton-polygon
     1588 
     1589  list NP=newtonpoly(f_tilde);
     1590  dbprint(i_print-1,"// Newton polygon:");
     1591  dbprint(i_print-1,NP);
     1592  dbprint(i_print-1,"");
     1593
     1594  if(is_NND(f,mu,NP))          // f is Newton non-degenerate
     1595  {
     1596    upper=NP[1][2];
     1597    ideal ES_Id= x^k*y^upper;
     1598    dbprint(i_print,"polynomial is Newton non-degenerate");
     1599    dbprint(i_print,"");
     1600    k=0;
     1601    for (i=1;i<=size(NP)-1;i++)
     1602    {
     1603      slope_x=NP[i+1][1]-NP[i][1];
     1604      slope_y=NP[i][2]-NP[i+1][2];
     1605      for (k=NP[i][1]+1; k<=NP[i+1][1]; k++)
     1606      {
     1607        while ( slope_x*upper + slope_y*k >=
     1608                slope_x*NP[i][2] + slope_y*NP[i][1])
     1609        {
     1610          upper=upper-1;
     1611        }
     1612        upper=upper+1;
     1613        ES_Id=ES_Id, x^k*y^upper;
     1614      }
     1615    }
     1616    ES_Id=std(ES_Id);
     1617    dbprint(i_print-1,"ideal of monomials above Newton bd. is generated by:");
     1618    dbprint(i_print-1,ES_Id);
     1619    ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
     1620    ES_Id = ES_Id, f, jacob(f);
     1621    dbprint(i_print,"// ");
     1622    dbprint(i_print,"// equisingularity ideal is computed!");
     1623    if (typ==0)
     1624    {
     1625      return(list(ES_Id,ESfix_Id));
     1626    }
     1627    else
     1628    {
     1629      return(list(ES_Id,ESfix_Id,ES_Id));
     1630    }
     1631  }
     1632  else
     1633  {
     1634    dbprint(i_print,"polynomial is Newton degenerate !");
     1635    dbprint(i_print,"");
     1636  }
     1637 
     1638  def old_ring=basering;
     1639
     1640  dbprint(i_print,"// ");
     1641  dbprint(i_print,"// versal deformation with triv. section");
     1642  dbprint(i_print,"// =====================================");
     1643  dbprint(i_print,"// ");
     1644 
     1645  ideal JJ=maxideal(1)*jacob(f);
     1646  ideal kbase_versal=kbase(std(JJ));
     1647  s=size(kbase_versal);
     1648  string ring_versal="ring @Px = ("+charstr(basering)+"),(t(1.."+string(s)+"),"
     1649                        +varstr(basering)+"),(ds("+string(s)+"),"
     1650                        +ordstr(basering)+");";
     1651  MinPolyStr = string(minpoly);
     1652                           
     1653  execute(ring_versal);
     1654  if (MinPolyStr<>"0")
     1655  {
     1656    MinPolyStr = "minpoly="+MinPolyStr;
     1657    execute(MinPolyStr);
     1658  }
     1659  // basering has changed to @Px
     1660
     1661  poly F=imap(old_ring,f);
     1662  ideal kbase_versal=imap(old_ring,kbase_versal);
     1663  for (i=1; i<=s; i++)
     1664  {
     1665    F=F+var(i)*kbase_versal[i];
     1666  }
     1667  dbprint(i_print-1,F);
     1668  dbprint(i_print-1,"");
     1669
     1670
     1671  ideal ES_Id,ES_Id_all_triv;
     1672  poly Ftriv=F;
     1673 
     1674  dbprint(i_print,"// ");
     1675  dbprint(i_print,"// Compute equisingularity Stratum over Spec(C[t]/t^2)");
     1676  dbprint(i_print,"// ===================================================");
     1677  dbprint(i_print,"// ");
     1678  list M=esStratum(F,2);
     1679  dbprint(i_print,"// finished");
     1680  dbprint(i_print,"// ");
     1681
     1682  if (M[2]==1) // error occured during esStratum computation
     1683  {
     1684    print("Some error has occured during the computation");
     1685    return(list(0,0));
     1686  }
     1687 
     1688  if ( typeof(M[1])=="list" )
     1689  {
     1690    int defpars = nvars(basering)-2;
     1691    poly Fred,Ftrivred;
     1692    poly g;
     1693    F=reduce(F,std(M[1][1]));
     1694    Ftriv=reduce(Ftriv,std(M[1][2]));
     1695 
     1696    for (i=1; i<=defpars; i++)
     1697    {
     1698      Fred=reduce(F,std(var(i)));
     1699      Ftrivred=reduce(Ftriv,std(var(i)));
     1700
     1701      g=subst(F-Fred,var(i),1);
     1702      ES_Id=ES_Id, g;
     1703      F=Fred;
     1704
     1705      g=subst(Ftriv-Ftrivred,var(i),1);
     1706      ES_Id_all_triv=ES_Id_all_triv, g;
     1707      Ftriv=Ftrivred;
     1708    }
     1709
     1710    setring old_ring;
     1711    // back to original ring
     1712
     1713    ideal ES_Id = imap(@Px,ES_Id);   
     1714    ES_Id = interred(ES_Id);
     1715
     1716    ideal ES_Id_all_triv = imap(@Px,ES_Id_all_triv);   
     1717    ES_Id_all_triv = interred(ES_Id_all_triv);
     1718
     1719    ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
     1720    ES_Id = ES_Id, f, jacob(f);
     1721    ES_Id_all_triv = ES_Id_all_triv, f, jacob(f);
     1722
     1723    if (typ==0)
     1724    {
     1725      return(list(ES_Id,ESfix_Id));
     1726    }
     1727    else
     1728    {
     1729      return(list(ES_Id,ESfix_Id,ES_Id_all_triv));
     1730    }
     1731  }
     1732  else
     1733  {
     1734    def AuxRing=M[1];
     1735 
     1736    dbprint(i_print,"// ");
     1737    dbprint(i_print,"// change ring to ESSring");
     1738 
     1739    setring AuxRing;  // contains p_F, ES
     1740
     1741    int defpars = nvars(basering)-2;
     1742    poly Fred,Fredtriv;
     1743    poly g;
     1744    ideal ES_Id,ES_Id_all_triv;
     1745
     1746    poly p_Ftriv=p_F
     1747
     1748    p_F=reduce(p_F,std(ES));
     1749    p_Ftriv=reduce(p_Ftriv,std(ES_all_triv));
     1750    for (i=1; i<=defpars; i++)
     1751    {
     1752      Fred=reduce(p_F,std(var(i)));
     1753      Fredtriv=reduce(p_Ftriv,std(var(i)));
     1754
     1755      g=subst(p_F-Fred,var(i),1);
     1756      ES_Id=ES_Id, g;
     1757      p_F=Fred;
     1758
     1759      g=subst(p_Ftriv-Fredtriv,var(i),1);
     1760      ES_Id_all_triv=ES_Id_all_triv, g;
     1761      p_Ftriv=Fredtriv;
     1762
     1763    }
     1764
     1765    dbprint(i_print,"// ");
     1766    dbprint(i_print,"// back to the original ring");
     1767
     1768    setring old_ring;
     1769    // back to original ring
     1770
     1771    ideal ES_Id = imap(AuxRing,ES_Id);
     1772    ES_Id = interred(ES_Id);
     1773
     1774    ideal ES_Id_all_triv = imap(AuxRing,ES_Id_all_triv);
     1775    ES_Id_all_triv = interred(ES_Id_all_triv);
     1776
     1777    kill @Px;
     1778    kill AuxRing;
     1779
     1780    ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
     1781    ES_Id = ES_Id, f, jacob(f);
     1782    ES_Id_all_triv = ES_Id_all_triv, f, jacob(f);
     1783    dbprint(i_print,"// ");
     1784    dbprint(i_print,"// equisingularity ideal is computed!");
     1785    if (typ==0)
     1786    {
     1787      return(list(ES_Id,ESfix_Id));
     1788    }
     1789    else
     1790    {
     1791      return(list(ES_Id,ESfix_Id,ES_Id_all_triv));
     1792    }
     1793  }
     1794}
     1795example
     1796{
     1797  "EXAMPLE:"; echo=2;
     1798  ring r=0,(x,y),ds;
     1799  poly f=x7+y7+(x-y)^2*x2y2;
     1800  list K=esIdeal(f);
     1801  option(redSB);
     1802  // Wahl's equisingularity ideal:
     1803  std(K[1]);
     1804
     1805  ring rr=0,(x,y),ds;
     1806  poly f=x4+4x3y+6x2y2+4xy3+y4+2x2y15+4xy16+2y17+xy23+y24+y30+y31;
     1807  list K=esIdeal(f);
     1808  vdim(std(K[1]));
     1809  // the latter should be equal to:
     1810  tau_es(f);     
     1811}
     1812
     1813///////////////////////////////////////////////////////////////////////////////
     1814
     1815proc esStratum (poly p_F, list #)
     1816"USAGE:   esStratum(F[,m,L]); F poly, m int, L list
     1817ASSUME:  F defines a deformation of a reduced bivariate polynomial f
     1818         and the characteristic of the basering does not divide mult(f). @*
     1819         If nv is the number of variables of the basering, then the first
     1820         nv-2 variables are the deformation parameters. @*
     1821         If the basering is a qring, ideal(basering) must only depend
     1822         on the deformation parameters.
     1823COMPUTE: equations for the stratum of equisingular deformations with 
     1824         fixed (trivial) section.
     1825RETURN:  list l: either consisting of a list and an integer, where
     1826@format
     1827  l[1][1]=ideal defining the equisingularity stratum
     1828  l[1][2]=ideal defining the part of the equisingularity stratum where all
     1829          equimultiple sections through the non-nodes of the reduced total
     1830          transform are trivial sections
     1831  l[2]=1 if some error has occured,  l[2]=0 otherwise;
     1832@end format
     1833or consisting of a ring and an integer, where
     1834@format
     1835  l[1]=ESSring is a ring extension of basering containing the ideal ES
     1836        (describing the ES-stratum), the ideal ES_all_triv (describing the
     1837        part with trival equimultiple sections) and the poly p_F=F,
     1838  l[2]=1 if some error has occured,  l[2]=0 otherwise.
     1839@end format
     1840NOTE:    L is supposed to be the output of hnexpansion (with the given ordering
     1841         of the variables appearing in f). @*
     1842         If m is given, the ES Stratum over A/maxideal(m) is computed. @*
     1843         This procedure uses @code{execute} or calls a procedure using
     1844         @code{execute}.
     1845         printlevel>=2 displays additional information.
     1846SEE ALSO: esIdeal, isEquising
     1847KEYWORDS: equisingularity stratum
     1848EXAMPLE: example esStratum; shows examples.
     1849"
     1850{
     1851  list l=esComputation (0,p_F,#);
     1852  return(l);
     1853}
     1854example
     1855{
     1856   "EXAMPLE:"; echo=2;
     1857   int p=printlevel;
     1858   printlevel=1;
     1859   ring r = 0,(a,b,c,d,e,f,g,x,y),ds;
    8511860   poly F = (x2+2xy+y2+x5)+ax+by+cx2+dxy+ey2+fx3+gx4;
    852    esStratum(F);
    853    esStratum(F,2);
     1861   list M = esStratum(F);
     1862   M[1][1];
     1863   
     1864   printlevel=3;     // displays additional information
     1865   esStratum(F,2)  ; // ES-stratum over Q[a,b,c,d,e,f,g] / <a,b,c,d,e,f,g>^2
     1866
    8541867   ideal I = f-fa,e+b;
    8551868   qring q = std(I);
    8561869   poly F = imap(r,F);
    8571870   esStratum(F);
     1871   printlevel=p;
    8581872}
    8591873
    8601874///////////////////////////////////////////////////////////////////////////////
    861 //                  procedures for equisingularity test
    862 ///////////////////////////////////////////////////////////////////////////////
    863 
    864 // DEFINES a new basering, myRing,  which has one variable
    865 //         more than the old ring.
    866 //         The name for the new variable is "H(nhelpV)".
    867 static proc T_extendRing(poly p_F, int nHelpV, ideal HCond)
    868 {
    869   def r_old = basering;
    870 
    871   ideal qIdeal = ideal(basering);
    872 
    873   string helpStr;
    874   string minPolyStr = "";
    875 
    876   if(minpoly != 0)
    877   {
    878     if (charstr(basering) == string(char(basering)) + "," + parstr(basering))
    879     {
    880       minPolyStr = string(minpoly);
    881     }
    882   }
    883 
    884   string str = "ring myRing =
    885                (" + charstr(r_old) + "),
    886                (H(" + string( nHelpV)+ ")," + string(maxideal(1)) + "),
    887                (dp(" + string( nHelpV) + "), ds);";
    888   execute (str);
    889 
    890   if (minPolyStr != "")
    891   {
    892     helpStr = "minpoly =" + minPolyStr + ";";
    893     execute(helpStr);
    894   }
    895 
    896   ideal qIdeal = imap(r_old, qIdeal);
    897   if(qIdeal != 0)
    898   {
    899     def r_base = basering;
    900     kill myRing;
    901     qring myRing = std(qIdeal);
    902   }
    903 
    904   poly p_F =imap(r_old, p_F);
    905   ideal HCond = imap(r_old, HCond);
    906 
    907   keepring(myRing);
    908 }
    909 ///////////////////////////////////////////////////////////////////////////////
    910 // tests, if ord p_F = nNew.
    911 static proc equimultTest (poly p_F, int nHelpV, int nNew, ideal HCond)
    912 {
    913   poly p_FnMinus1;
    914   ideal id_help;
    915   matrix coefMatrix;
    916   int i;
    917   int nc;
    918 
    919   p_FnMinus1 = jet(p_F, nNew-1, xyVector());
    920 
    921   coefMatrix = coef(p_FnMinus1, xy);
    922   nc = ncols(coefMatrix);
    923 
    924   for (i=1; i<=nc; i++)
    925   {
    926     id_help[i] = coefMatrix[2,i];
    927   }
    928 
    929   id_help = T_helpVarElim(id_help, HCond, nHelpV);
    930 
    931   if (reduce(id_help, std(0)) !=0 )
    932   {
    933     return(0, p_F);
    934   }
    935 
    936   p_F = p_F - p_FnMinus1;
    937 
    938   return(1, p_F);
    939 }
    940 ///////////////////////////////////////////////////////////////////////////////
    941 // ASSUME: ord(p_F)=nNew
    942 // tests, if p_F =  p_c*y^nNew for some p_c.
    943 static proc pPOfYTest (poly p_F, int nHelpV, int nNew, ideal HCond)
    944 {
    945   poly  p_Fn;
    946   poly p_c;
    947   ideal id_help;
    948   int nc;
    949   int i=1;
    950   matrix coefMatrix;
    951 
    952   p_Fn = jet(p_F, nNew, xyVector());
    953 
    954   coefMatrix = coef(p_Fn, xy);
    955   nc = ncols(coefMatrix);
    956 
    957   p_c = coefMatrix[2,1];
    958 
    959   for (i = 2; i <= nc; i++)
    960   {
    961     id_help[i] = coefMatrix[2,i];
    962   }
    963 
    964   id_help = T_helpVarElim(id_help, HCond, nHelpV);
    965 
    966   if (reduce(id_help, std(0)) !=0 )
    967   {
    968     return(0, p_c);
    969   }
    970 
    971   return(1, p_c);
    972 }
    973 ///////////////////////////////////////////////////////////////////////////////
    974 // ASSUME: ord(p_F)=nNew
    975 // tests, if p_F =  p_c*(y - p_a*x)^nNew for some p_a, p_c.
    976 static proc pPOfLinTest(poly p_F, int nNew, int nHelpV, ideal HCond)
    977 {
    978   poly p_Fn;
    979   poly p_c;
    980   poly p_ca;
    981   poly p_help;
    982   poly p_a;
    983   ideal id_help;
    984 
    985   p_Fn = jet(p_F, nNew, xyVector());
    986 
    987   p_c  = coefficient(p_Fn,y^nNew,y);
    988   p_ca = coefficient(p_Fn,y^(nNew-1)*x,xy)/-nNew;
    989 
    990   if (npars(basering)==1
    991   && charstr(basering) != string(char(basering)) + "," + parstr(basering))
    992   {
    993     p_a = H(nHelpV);
    994     HCond = HCond + ideal(p_ca - p_a*p_c);
    995   }
    996   else
    997   {
    998     p_help = p_ca/p_c;
    999     if (p_help * p_c == p_ca)
    1000     {
    1001       p_a = p_help;
    1002     }
    1003     else
    1004     {
    1005       p_a = H(nHelpV);
    1006       HCond = HCond + ideal(p_ca - p_a*p_c);
    1007     }
    1008   }
    1009 
    1010   id_help = ideal(p_Fn - p_c *(y - p_a *x)^nNew);
    1011   id_help = T_helpVarElim(id_help, HCond, nHelpV);
    1012 
    1013   if (reduce(id_help, std(0)) != 0 )
    1014   {
    1015     return(0, p_F, p_c, p_a, HCond);
    1016   }
    1017 
    1018   return(1, p_F, p_c, p_a, HCond);
    1019 }
    1020 //////////////////////////////////////////////////////////////////////////////
    1021 // eliminates the variables H(1),..,H(nHelpV) from the ideal ES + HCond
    1022 static proc T_helpVarElim(ideal ES, ideal HCond, int nHelpV)
    1023 {
    1024 
    1025   def r_old = basering;
    1026 
    1027   ideal qIdeal = ideal(basering);
    1028 
    1029   string helpStr;
    1030   string minPolyStr = "";
    1031 
    1032   if(minpoly != 0)
    1033   {
    1034     if (charstr(basering) == string(char(basering)) + "," + parstr(basering))
    1035     {
    1036       minPolyStr = string(minpoly);
    1037     }
    1038   }
    1039 
    1040   string str = "ring myRing =
    1041                (" + charstr(r_old) + "),(" + string(maxideal(1)) + "),
    1042                (dp(" + string( nHelpV) + "), dp);";
    1043   execute (str);
    1044 
    1045   if (minPolyStr != "")
    1046   {
    1047     helpStr = "minpoly =" + minPolyStr + ";";
    1048     execute(helpStr);
    1049   }
    1050 
    1051   ideal qIdeal = imap(r_old, qIdeal);
    1052   if(qIdeal != 0)
    1053   {
    1054     def r_base = basering;
    1055     kill myRing;
    1056     qring myRing = std(qIdeal);
    1057   }
    1058 
    1059   ideal ES = imap(r_old, ES);
    1060   ideal HCond = imap(r_old, HCond);
    1061 
    1062   ES = ES + HCond;
    1063   ES = std(ES);
    1064   ES = nselect(ES,1,nHelpV);
    1065 
    1066   setring r_old;
    1067   ES = imap (myRing,ES);
    1068 
    1069   return(ES);
    1070 }
    1071 ///////////////////////////////////////////////////////////////////////////////
    1072 // ASSUME:  F in K[t(1)..t(s),x,y],
    1073 //          the ringordering is ds
    1074 // RETURNS: list l, such that
    1075 //          l[1]=1 if some error has occured,
    1076 //          l[1]=0 otherwise and then
    1077 //          l[2] = 1, if the deformation is equisingular and
    1078 //          l[2] = 0 otherwise.
    1079 static proc equisingTest (poly p_F, intvec multSeq, int maxStep)
    1080 {
    1081   def r_old = basering;
    1082 
    1083   ideal id_Es = 0;
    1084 
    1085   int isES = 1;
    1086   int step = 1;
    1087   int nNew = multSeq[step];
    1088   int nOld = nNew;
    1089   int zeroSteps;
    1090   ideal  HCond = 0;
    1091   int  nHelpV = 1;
    1092 
    1093   T_extendRing (p_F, nHelpV, HCond);
    1094 
    1095   poly p_c;
    1096   poly p_a;
    1097 
    1098   while (step <= maxStep)
    1099   {
    1100     nOld = nNew;
    1101     nNew = multSeq[step];
    1102 
    1103     if (nNew < nOld)
    1104     //start a new line in the HNE of F
    1105     //               _            _
    1106     // for the next | nold/nnew -1 | iteration steps the coefficient 'a'
    1107     // in the leading form Fn = c(y-ax) should  be zero
    1108     {
    1109       p_F = swapXY(p_F);
    1110       zeroSteps = calcNZeroSteps (nOld, nNew);
    1111     }
    1112 
    1113     isES, p_F = equimultTest (p_F, nHelpV, nNew, HCond);
    1114 
    1115     if (! isES)
    1116     {
    1117       return(0);
    1118     }
    1119 
    1120     if (zeroSteps > 0)
    1121     {
    1122       zeroSteps--;
    1123 
    1124       isES, p_c = pPOfYTest (p_F,  nHelpV, nNew, HCond);
    1125       p_a = 0;
    1126     }
    1127     else
    1128     {
    1129       isES, p_F, p_c, p_a, HCond = pPOfLinTest (p_F, nNew, nHelpV, HCond);
    1130     }
    1131 
    1132     if (! isES)
    1133     {
    1134       return(0);
    1135     }
    1136 
    1137     p_F = formalBlowUp (p_F, p_c, p_a, nNew);
    1138 
    1139     if (p_a == H(nHelpV))
    1140     {
    1141       nHelpV++;
    1142 
    1143       def r_base = basering;
    1144       kill myRing;
    1145 
    1146       T_extendRing(p_F, nHelpV, HCond);
    1147 
    1148       kill r_base;
    1149 
    1150       poly p_a;
    1151       poly p_c;
    1152     }
    1153 
    1154     step++;
    1155   }
    1156 
    1157   return(1);
    1158 }
    1159 ///////////////////////////////////////////////////////////////////////////////
    11601875
    11611876proc isEquising (poly p_F, list #)
    1162 "USAGE:   isEquising(F[,m]); F poly, m int
    1163 ASSUME:   F defines a deformation of an irreducible bivariate polynomial f
     1877"USAGE:   isEquising(F[,m,L]); F poly, m int, L list
     1878ASSUME:  F defines a deformation of a reduced bivariate polynomial f
    11641879         and the characteristic of the basering does not divide mult(f). @*
    1165          If nv is the number of variables of the basering, then the first nv-2
    1166          variables are the deformation parameters. @*
     1880         If nv is the number of variables of the basering, then the first
     1881         nv-2 variables are the deformation parameters. @*
    11671882         If the basering is a qring, ideal(basering) must only depend
    11681883         on the deformation parameters.
    1169 RETURN:  list l of two integers, where
    1170 @format
    1171    l[1]=1 if F is an equisingular deformation,  l[1]=0 otherwise.
    1172    l[2]=1 if some error has occured,  l[2]=0 otherwise.
    1173 @end format
    1174 NOTE:    If m is given, the computation stops after m steps of the iteration. @*
    1175          This procedure uses @code{execute} or calls a procedure using
     1884COMPUTE: tests if the given family is equisingular along the trivial
     1885         section.
     1886RETURN:  int: 1 if the family is equisingular, 0 otherwise.
     1887NOTE:    L is supposed to be the output of hnexpansion (with the given ordering
     1888         of the variables appearing in f). @*
     1889         If m is given, the family is considered over A/maxideal(m). @*
     1890         This procedure uses @code{execute} or calls a procedure using 
    11761891         @code{execute}.
    1177 EXAMPLE: example isEquising; shows an example
     1892         printlevel>=2 displays additional information.
     1893EXAMPLE: example isEquising; shows examples.
    11781894"
    11791895{
    1180   def r_user = basering;
    1181 
    1182   int ii = 1;
    1183   int i_nvars = nvars(basering);
    1184   int error = 0;
    1185   int maxStep;
    1186   int userMaxStep;
    1187   int xNotTransversal = 0;
    1188   int fIrreducible = 1;
    1189   intvec multSeq;
    1190   ideal isES = 1;
    1191 
    1192   error = checkBasering();
    1193   if (error)
    1194   {
    1195     return(0,1);
    1196   }
    1197 
    1198   userMaxStep = getInput(#);
    1199 
    1200   // define a new basering "myRing" with new names for parameters
    1201   // and variables.
    1202   // The new names are 'a(1)', ..., 'a(npars)' for the parameters
    1203   // and 't(1)', ..., 't(nvars-2)', 'x', 'y' for the variables.
    1204   createMyRing(p_F, "ds");
    1205 
    1206   createHNERing();
    1207 
    1208   poly p_f = imap(myRing,p_F);
    1209 
    1210   error = checkPoly(p_f);
    1211   if (error)
    1212   {
    1213     return(0,1);
    1214   }
    1215 
    1216   // compute the multiplicity sequence of p_f.
    1217   multSeq, xNotTransversal, fIrreducible = calcMultSequence(p_f);
    1218 
    1219    if ( ! fIrreducible)
    1220    {
    1221      return(list(0,1));
    1222    }
    1223 
    1224   setring myRing;
    1225 
    1226   if (xNotTransversal)
    1227   {
    1228     p_F = swapXY(p_F);
    1229   }
    1230 
    1231   if (userMaxStep != -1 && userMaxStep <  size(multSeq)-1)
    1232   {
    1233     maxStep = userMaxStep;
    1234   }
    1235   else
    1236   {
    1237     maxStep = size(multSeq)-1;
    1238   }
    1239 
    1240   int isES = equisingTest(p_F, multSeq, maxStep);
    1241 
    1242   return(list(isES, error));
     1896  int check=esComputation (1,p_F,#);
     1897  return(check);
    12431898}
    1244 
    12451899example
    12461900{
    12471901   "EXAMPLE:"; echo=2;
    1248    ring r = 11,(a,b,x,y),ds;
     1902   ring r = 0,(a,b,x,y),ds;
    12491903   poly F = (x2+2xy+y2+x5)+ay3+bx5;
    12501904   isEquising(F);
    1251    isEquising(F,1);
    12521905   ideal I = ideal(a);
    12531906   qring q = std(I);
    12541907   poly F = imap(r,F);
    12551908   isEquising(F);
     1909
     1910   ring rr=0,(A,B,C,x,y),ls;
     1911   poly f=x7+y7+(x-y)^2*x2y2;
     1912   poly F=f+A*y*diff(f,x)+B*x*diff(f,x);
     1913   isEquising(F); 
     1914   isEquising(F,2);    // computation over  Q[a,b] / <a,b>^2
    12561915}
    1257 ///////////////////////////////////////////////////////////////////////////////
    1258 /*
    1259   Weiter Beispiele aus Dipl. von A. Mindnich einfuegen
     1916
     1917
     1918
     1919/*  Examples:
     1920
     1921LIB "equising.lib";
     1922   ring r = 0,(x,y),ds;
     1923   poly p1 = y^2+x^3;
     1924   poly p2 = p1^2+x5y;
     1925   poly p3 = p2^2+x^10*p1;
     1926   poly p=p3^2+x^20*p2;
     1927   p;
     1928   versal(p);
     1929   setring Px;
     1930   poly F=Fs[1,1];
     1931   int t=timer;
     1932   list M=esStratum(F);
     1933   timer-t;  //-> 3
     1934
     1935LIB "equising.lib";
     1936option(prot);
     1937printlevel=2;
     1938ring r=0,(x,y),ds;
     1939poly f=(x-yx+y2)^2-(y+x)^31;
     1940versal(f);
     1941setring Px;
     1942poly F=Fs[1,1];
     1943int t=timer;
     1944list M=esStratum(F);
     1945timer-t;  //-> 233
     1946
     1947
     1948LIB "equising.lib";
     1949printlevel=2;
     1950option(prot);
     1951timer=1;
     1952ring r=32003,(x,y),ds;
     1953poly f=(x4-y4)^2-x10;
     1954versal(f);
     1955setring Px;
     1956poly F=Fs[1,1];
     1957int t=timer;
     1958list M=esStratum(F,3);
     1959timer-t;  //-> 8
     1960
     1961LIB "equising.lib";
     1962printlevel=2;
     1963timer=1;
     1964ring rr=0,(x,y),ls;
     1965poly f=x7+y7+(x-y)^2*x2y2;
     1966list K=esIdeal(f);
     1967// tau_es
     1968vdim(std(K[1])); //-> 22
     1969// tau_es_fix
     1970vdim(std(K[2])); //-> 24
     1971
     1972
     1973LIB "equising.lib";
     1974printlevel=2;
     1975timer=1;
     1976ring rr=0,(x,y),ls;
     1977poly f=x7+y7+(x-y)^2*x2y2+x2y4; // Newton non-deg.
     1978list K=esIdeal(f);
     1979// tau_es
     1980vdim(std(K[1])); //-> 21
     1981// tau_es_fix
     1982vdim(std(K[2])); //-> 23
     1983
     1984LIB "equising.lib";
     1985ring r=0,(w,v),ds;
     1986poly f=w2-v199;
     1987list L=hnexpansion(f);
     1988versal(f);
     1989setring Px;
     1990list L=imap(r,L);
     1991poly F=Fs[1,1];
     1992list M=esStratum(F,2,L);
     1993
     1994LIB "equising.lib";
     1995printlevel=2;
     1996timer=1;
     1997ring rr=0,(A,B,C,x,y),ls;
     1998poly f=x7+y7+(x-y)^2*x2y2;
     1999poly F=f+A*y*diff(f,x)+B*x*diff(f,x)+C*diff(f,y);
     2000list M=esStratum(F,6);
     2001std(M[1]);  // standard basis of equisingularity ideal
     2002
     2003LIB "equising.lib";
     2004printlevel=2;
     2005timer=1;
     2006ring rr=0,(x,y),ls;
     2007poly f=x20+y7+(x-y)^2*x2y2+x2y4; // Newton non-degenerate
     2008list K=esIdeal(f);
     2009
     2010ring rr=0,(x,y),ls;
     2011poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13;
     2012list K=esIdeal(f);
     2013versal(f);
     2014setring Px;
     2015poly F=Fs[1,1];
     2016list M=esStratum(F,2);
     2017
     2018LIB "equising.lib";
     2019printlevel=2;
     2020ring rr=0,(x,y),ls;
     2021poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13;
     2022list K=esIdeal(f);
     2023vdim(std(K[1]));  //-> 51
     2024tau_es(f);        //-> 51
     2025
     2026printlevel=3;
     2027f=f*(y-x2)*(y2-x3)*(x-y5);
     2028int t=timer;
     2029list L=esIdeal(f);
     2030vdim(std(L[1]));  //-> 99
     2031timer-t;   //-> 42
     2032t=timer;
     2033tau_es(f);        //-> 99
     2034timer-t;   //-> 23
     2035
     2036
     2037LIB "equising.lib";
     2038printlevel=3;
     2039ring rr=0,(x,y),ds;
     2040poly f=x4+4x3y+6x2y2+4xy3+y4+2x2y15+4xy16+2y17+xy23+y24+y30+y31;
     2041list K=esIdeal(f);
     2042vdim(std(K[1]));  //-> 68
     2043tau_es(f);        //-> 68
     2044
     2045versal(f);
     2046setring Px;
     2047poly F=Fs[1,1];
     2048int t=timer;
     2049list M=esStratum(F);
     2050timer-t;          //-> 0
     2051
     2052
     2053
    12602054*/
Note: See TracChangeset for help on using the changeset viewer.