source: git/Singular/LIB/equising.lib @ f6e355

spielwiese
Last change on this file since f6e355 was f6e355, checked in by Hans Schönemann <hannes@…>, 19 years ago
*lossen: new lib git-svn-id: file:///usr/local/Singular/svn/trunk@7836 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 51.0 KB
RevLine 
[f6e355]1version="$Id: equising.lib,v 1.9 2005-04-15 15:20:12 Singular Exp $";
[fd3fb7]2category="Singularities";
[e7cc147]3info="
[558eb2]4LIBRARY:  equising.lib  Equisingularity Stratum of a Family of Plane Curves
[f6e355]5AUTHOR:   Christoph Lossen, lossen@mathematik.uni-kl.de
6          Andrea Mindnich, mindnich@mathematik.uni-kl.de
[4ac997]7
[f6e355]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
[5d9ec3]16";
[e7cc147]17
[f6e355]18LIB "hnoether.lib";
[e7cc147]19LIB "poly.lib";
20LIB "elim.lib";
[f6e355]21LIB "deform.lib";
22LIB "sing.lib";
23
24////////////////////////////////////////////////////////////////////////////////
25//
26//  The following (static) procedures are used by esComputation
27//
28////////////////////////////////////////////////////////////////////////////////
[e7cc147]29// COMPUTES  a weight vector. x and y get weight 1 and all other
30//           variables get weight 0.
31static proc xyVector()
32{
33  intvec iv ;
34  iv[nvars(basering)]=0 ;
35  iv[rvar(x)] =1;
36  iv[rvar(y)] =1;
37  return (iv);
38}
[f6e355]39
40////////////////////////////////////////////////////////////////////////////////
41// exchanges the variables x and y in the polynomial f
[e7cc147]42static proc swapXY(poly f)
43{
44  def r_base = basering;
45  ideal MI = maxideal(1);
46  MI[rvar(x)]=y;
47  MI[rvar(y)]=x;
48  map phi = r_base, MI;
49  f=phi(f);
50  return (f);
51}
[f6e355]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);
[e7cc147]56{
[f6e355]57  intvec w=xyVector(); 
58  poly Fd=jet(F,m,w); 
59  return(Fd);
[e7cc147]60}
[f6e355]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"
[e7cc147]93{
[f6e355]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];
[e7cc147]100  }
[f6e355]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++)
[e7cc147]107  {
[f6e355]108    dummy=1;
109    for (j=1;j<=ncols(M[2]);j++)
[e7cc147]110    {
[f6e355]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;
[e7cc147]121    }
122  }
123
[f6e355]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++)
[e7cc147]138  {
[f6e355]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    }
[e7cc147]148  }
[f6e355]149  // adding up multiplicities
150  for (i=nrows(S);i>=2;i--) 
[e7cc147]151  {
[f6e355]152    for (j=1;j<=ncols(S);j++)
153    {
154      maxDeg[i-1,j]=maxDeg[i-1,j]+maxDeg[i,j];
155    }
[e7cc147]156  }
[f6e355]157 
158  list L=S,T,U,maxDeg;
159  return(L);
[e7cc147]160}
161
162
[f6e355]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,
168{
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++)
[e7cc147]177  {
[f6e355]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++)
[e7cc147]186    {
[f6e355]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      }
[e7cc147]241    }
[f6e355]242  }
243  L=B,V;
244  return(L);
245}
[e7cc147]246
[f6e355]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
253{
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++)
[e7cc147]259    {
[f6e355]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     
[e7cc147]268    }
269  }
[f6e355]270  counter=counter+ncols(U);
271  return(counter);
[e7cc147]272}
273
[f6e355]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
282{
283  int i,j,k,counter;
284  for (j=1;j<=ncols(Mult);j++)
[e7cc147]285  {
[f6e355]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    }
[e7cc147]296  }
[f6e355]297  return(counter);
[e7cc147]298}
299
300
301///////////////////////////////////////////////////////////////////////////////
302// COMPUTES string(minpoly) and substitutes the parameter by newParName
303static proc makeMinPolyString (string newParName)
304{
305  int i;
306  string parName = parstr(basering);
307  int parNameSize = size(parName);
308
309  string oldMinPolyStr = string (minpoly);
310  int minPolySize = size(oldMinPolyStr);
311
312  string newMinPolyStr = "";
313
314  for (i=1;i <= minPolySize; i++)
315  {
316    if (oldMinPolyStr[i,parNameSize] == parName)
317    {
318      newMinPolyStr = newMinPolyStr + newParName;
319      i = i + parNameSize-1;
320    }
321    else
322    {
323      newMinPolyStr = newMinPolyStr + oldMinPolyStr[i];
324    }
325  }
326
327  return(newMinPolyStr);
328}
329
330
331///////////////////////////////////////////////////////////////////////////////
[f6e355]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)
[e7cc147]341{
[f6e355]342  def r_old = basering;
[e7cc147]343
[f6e355]344  int chara = char(basering);
345  string charaStr;
346  int i;
347  string helpStr;
348  int nDefParams = nvars(r_old)-2;
[e7cc147]349
[f6e355]350  ideal qIdeal = ideal(basering);
[e7cc147]351
[f6e355]352  if ((npars(basering)==0) and (minPolyStr==""))
[e7cc147]353  {
[f6e355]354    helpStr = "ring myRing1 ="
355              + string(chara)+ ", (t(1..nDefParams), x, y),("+ ordStr +");";
356    execute(helpStr);
[e7cc147]357  }
[f6e355]358  else
[e7cc147]359  {
[f6e355]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);
[e7cc147]369
[f6e355]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
[e7cc147]381    {
[f6e355]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);
[e7cc147]387    }
[f6e355]388  }
[e7cc147]389
[f6e355]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;
[e7cc147]394
[f6e355]395  // Extension by no_b auxiliary variables
396  if (no_b>0)
[e7cc147]397  {
[f6e355]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    }
[e7cc147]451
[f6e355]452    poly p_F = imap(myRing1, p_F);
453    ideal mIdeal = imap(myRing1, mIdeal);
454    export p_F,mIdeal;
455    kill myRing1;
[e7cc147]456  }
[f6e355]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;   
[e7cc147]474  }
[f6e355]475 
476  setring r_old;
477  return(myRing);
[e7cc147]478}
[f6e355]479
480////////////////////////////////////////////////////////////////////////////////
481// returns list of coef, leadmonomial 
482//
483static proc determine_coef (poly Fm)
[e7cc147]484{
[f6e355]485  def r_base = basering; // is assumed to be the result of createMyRing
[e7cc147]486
487  int chara = char(basering);
488  string charaStr;
489  int i;
490  string minPolyStr = "";
[f6e355]491  string helpStr = "";
[e7cc147]492
493  if (npars(basering) == 0)
494  {
[f6e355]495    helpStr = "ring myRing1 ="
496              + string(chara)+ ", (y,x),ds;";
[e7cc147]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");
[f6e355]507        helpStr = "ring myRing1 = (" + string(chara) + ",a), (y,x),ds;";
[e7cc147]508        execute(helpStr);
509
510        helpStr = "minpoly =" + minPolyStr + ";";
511        execute (helpStr);
512      }
[f6e355]513      else // no minpoly given
[e7cc147]514      {
[f6e355]515        helpStr = "ring myRing1 =
516                  (" + string(chara) + ",a(1..npars(basering)) ), (y,x),ds;";
[e7cc147]517        execute(helpStr);
518      }
519    }
520    else
521    {
522      i = find (charaStr,",");
523
[f6e355]524      helpStr = " ring myRing1 = (" + charaStr[1,i] + "a), (y,x),ds;";
[e7cc147]525      execute (helpStr);
526    }
527  }
[f6e355]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++)
[e7cc147]536  {
[f6e355]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    }
[e7cc147]542  }
[f6e355]543  return(list(h,g));
[e7cc147]544}
545
546///////////////////////////////////////////////////////////////////////////////
[f6e355]547// RETURNS: 1, if p_f = 0 or char(basering) divides the order of p_f
548//             or p_f is not squarefree.
549//          0, otherwise
550static proc checkPoly (poly p_f)
[e7cc147]551{
[f6e355]552  int i_print = printlevel - voice + 3;
553  int i_ord;
[e7cc147]554
[f6e355]555  if (p_f == 0)
[e7cc147]556  {
[f6e355]557    print("Input is a 'deformation'  of the zero polynomial!");
558    return(1);
[e7cc147]559  }
560
[f6e355]561  i_ord = mindeg1(p_f);
[e7cc147]562
[f6e355]563  if (number(i_ord) == 0)
[e7cc147]564  {
[f6e355]565    print("Characteristic of coefficient field "
566            +"divides order of zero-fiber !");
567    return(1);
[e7cc147]568  }
569
[f6e355]570  if (squarefree(p_f) != p_f)
[e7cc147]571  {
[f6e355]572    print("Original polynomial (= zero-fiber) is not reduced!");
573    return(1);
[e7cc147]574  }
575
[f6e355]576  return(0);
[e7cc147]577}
578
[f6e355]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++)
[e7cc147]588  {
[f6e355]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    } 
[e7cc147]608  }
[f6e355]609  list L=newvar,newmap,counter;
610  attrib(J,"isSB",0);
611  return(L);
[e7cc147]612}
613
614///////////////////////////////////////////////////////////////////////////////
[f6e355]615//  The following procedure is called by esStratum (typ=0), resp. by
616//  isEquising (typ=1)
617///////////////////////////////////////////////////////////////////////////////
[e7cc147]618
[f6e355]619static proc esComputation (int typ, poly p_F, list #)
620{
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)
[e7cc147]638  {
[f6e355]639    print("family is trivial (no deformation parameters)!");
640    if (typ==1) //isEquising
[e7cc147]641    {
[f6e355]642      return(1); 
[e7cc147]643    }
644    else
645    {
[f6e355]646      return(list(ideal(0),0));
[e7cc147]647    }
648  }
649
[f6e355]650  if (size(#)>0)
[e7cc147]651  {
[f6e355]652    if (typeof(#[1])=="int")
[e7cc147]653    {
[f6e355]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
[e7cc147]665        {
[f6e355]666          return(list(ideal(0),int(1)));
[e7cc147]667        }
668      }
[f6e355]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);
[e7cc147]746    }
747    else
748    {
[f6e355]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      }
[e7cc147]757    }
[f6e355]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
[e7cc147]777  }
778
[f6e355]779  // Determine maximal number of needed auxiliary parameters (free tangents):
780  no_b=Determine_no_b(M[3],B); 
[e7cc147]781
[f6e355]782  // test whether HNexpansion needed field extension....
783  string minPolyStr = "";
784  if (minpoly !=0)
785  {
786    minPolyStr = makeMinPolyString("a");
787    minPolyStr = "minpoly =" + minPolyStr + ";";
788  }
[e7cc147]789
[f6e355]790  setring old_ring;
[e7cc147]791
[f6e355]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
[e7cc147]797
[f6e355]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); 
[e7cc147]808  }
[f6e355]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
[e7cc147]829  {
[f6e355]830    J=jet(J,artin_bd-1);
[e7cc147]831  }
832
[f6e355]833  J=interred(J);
[e7cc147]834
[f6e355]835// J=std(J);
[e7cc147]836
[f6e355]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];
[e7cc147]890
[f6e355]891  }
[e7cc147]892
[f6e355]893  list LeadDataFm=determine_coef(Fm);
894  def LeadDataG=coef(G,xy);
895 
896  for (i=1; i<=ncols(LeadDataG); i++)
[e7cc147]897  {
[f6e355]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    }
[e7cc147]903  }
[f6e355]904   
905  G=LeadDataFm[1]*G-LeadG*Fm;  // leading terms in y should cancel...
[e7cc147]906
[f6e355]907  coef_Mat = coef(G,xy);
908  Jnew=coef_Mat[2,1..ncols(coef_Mat)];
[e7cc147]909
[f6e355]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  }
[e7cc147]915
[f6e355]916  // simplification of Jnew
917  Jnew=interred(Jnew);
[e7cc147]918
[f6e355]919  J=J,Jnew;
920  if (typ==1) // isEquising
921  {
922    if(ideal(nselect(J,1,no_b))<>0)
[e7cc147]923    {
[f6e355]924      setring old_ring;
925      return(0);
926    }   
927  }
[e7cc147]928
929
[f6e355]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)
[e7cc147]937    {
[f6e355]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)
[e7cc147]942      {
[f6e355]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        }
[e7cc147]950      }
[f6e355]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        }
[e7cc147]1033
[f6e355]1034        // simplification of J
1035        Jnew=interred(Jnew);
[e7cc147]1036
[f6e355]1037        J=J,Jnew;
1038        if (typ==1) // isEquising
[e7cc147]1039        {
[f6e355]1040          if(ideal(nselect(J,1,no_b))<>0)
1041          {
1042            setring old_ring;
1043            return(0);
1044          }   
[e7cc147]1045        }
[f6e355]1046
[e7cc147]1047      }
1048    }
[f6e355]1049    if (number_of_branches>=2)
[e7cc147]1050    {
[f6e355]1051      J=interred(J);
1052      if (typ==1) // isEquising
[e7cc147]1053      {
[f6e355]1054        if(ideal(nselect(J,1,no_b))<>0)
1055        {
1056          setring old_ring;
1057          return(0);
1058        }   
[e7cc147]1059      }
1060    }
[f6e355]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);
[e7cc147]1070
[f6e355]1071  dbprint(i_print,"// ");
1072  dbprint(i_print,"// Elimination starts:");
1073  dbprint(i_print,"// -------------------");
[e7cc147]1074
[f6e355]1075  poly gg;
1076  int b_left=no_b; 
[e7cc147]1077
[f6e355]1078  for (i=1;i<=no_b; i++)
[e7cc147]1079  {
[f6e355]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);
[e7cc147]1088  }
[f6e355]1089  J=simplify(J,10);
1090  if (typ==1) // isEquising
[e7cc147]1091  {
[f6e355]1092    if(ideal(nselect(J,1,no_b))<>0)
1093    {
1094      setring old_ring;
1095      return(0);
1096    }   
[e7cc147]1097  }
1098
[f6e355]1099  ideal J_no_b = nselect(J,1,no_b);
1100  if (size(J) > size(J_no_b))
[e7cc147]1101  {
[f6e355]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 = "";
1111      if (minpoly !=0)
1112      {
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)
1136      {
1137         attrib(MM,"isSB",1);
1138         qring QQ=MM;     
1139      }
1140        else
1141        {   
1142           qId=qId,MM;
1143           qring QQ = std(qId);
1144        }
[e7cc147]1145
[f6e355]1146      ideal Shortmap=imap(myRing,Shortmap);
1147      map phiphi=myRing,Shortmap;
[e7cc147]1148
[f6e355]1149      ideal J=phiphi(J);
1150      J=std(J);
1151      J=nselect(J,1,no_b);
[e7cc147]1152
[f6e355]1153      setring myRing;
1154      // back to "myRing"
[e7cc147]1155
[f6e355]1156      J=nselect(J,1,no_b);
1157      Jnew=imap(QQ,J);
[e7cc147]1158
[f6e355]1159      J=J,Jnew;
1160      J=interred(J);
1161    }
1162    else
1163    {
1164      J=std(J);
1165      J=nselect(J,1,no_b);
1166    }
[e7cc147]1167  }
[f6e355]1168 
1169  dbprint(i_print,"// finished");
1170  dbprint(i_print,"// ");
1171 
1172  minPolyStr = "";
1173  if (minpoly !=0)
[e7cc147]1174  {
[f6e355]1175   minPolyStr = "minpoly = "+string(minpoly);
[e7cc147]1176  }
1177
[f6e355]1178  kill HNEring;
1179 
1180  if (typ==1) // isEquising
[e7cc147]1181  {
[f6e355]1182    if(J<>0)
1183    {
1184      setring old_ring;
1185      return(0);
1186    }   
1187    else
1188    {
1189      setring old_ring;
1190      return(1);
1191    }   
[e7cc147]1192  }
[f6e355]1193 
1194  setring old_ring;
1195  // we are back in the original ring
[e7cc147]1196
[f6e355]1197  if (npars(myRing)<>0)
[e7cc147]1198  {
[f6e355]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));
[e7cc147]1241  }
1242  else
[f6e355]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));
[e7cc147]1263  }
[f6e355]1264 
[e7cc147]1265}
1266
[f6e355]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.
1280"
[e7cc147]1281{
[f6e355]1282  int i,j,k,s;
1283  int slope_x, slope_y, upper;
1284  int i_print = printlevel - voice + 3;
1285  string MinPolyStr;
[e7cc147]1286
[f6e355]1287  // some checks first
1288  if ( nvars(basering)<>2 )
[e7cc147]1289  {
[f6e355]1290    print("// basering has not the correct number (two) of variables !");
1291    print("// computation stopped");
1292    return(0);
[e7cc147]1293  }
[f6e355]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  } 
[e7cc147]1300
[f6e355]1301  if (mult(std(f))<=1)
[e7cc147]1302  {
[f6e355]1303    // f is rigid
1304    return(0);
[e7cc147]1305  }
1306
[f6e355]1307  if ( deg(squarefree(f))!=deg(f) )
[e7cc147]1308  {
[f6e355]1309    print("// input polynomial was not reduced");
1310    print("// try    squarefree(f);   first");
1311    return(0);
[e7cc147]1312  }
1313
[f6e355]1314  def old_ring=basering;
1315  execute("ring @myRing=("+charstr(basering)+"),("+varstr(basering)+"),ds;");
1316  poly f=imap(old_ring,f);
[e7cc147]1317
[f6e355]1318  ideal Jacobi_Id = jacob(f);
[e7cc147]1319
[f6e355]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  }
[e7cc147]1329
[f6e355]1330  // check for D_k singularity
1331  // ----------------------------------------
1332  if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3)
[e7cc147]1333  {
[f6e355]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)));
[e7cc147]1339  }
1340
1341
[f6e355]1342  if (size(#)==0)
[e7cc147]1343  {
[f6e355]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    }
[e7cc147]1393  }
1394
[f6e355]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  }   
1412  else
1413  {
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);
[e7cc147]1450}
[f6e355]1451example
[e7cc147]1452{
[f6e355]1453   "EXAMPLE:"; echo=2;
1454   ring r=32003,(x,y),ds;
1455   poly f=(x4-y4)^2-x10;
1456   tau_es(f);
1457}
[e7cc147]1458
1459
[f6e355]1460////////////////////////////////////////////////////////////////////////////////
[e7cc147]1461
[f6e355]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{
[e7cc147]1489
[f6e355]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 )
[e7cc147]1499  {
[f6e355]1500    print("// basering has not the correct number (two) of variables !");
1501    print("// computation stopped");
1502    return(list(0,0));
[e7cc147]1503  }
[f6e355]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  } 
[e7cc147]1510
[f6e355]1511  if (mult(std(f))<=1)
[e7cc147]1512  {
[f6e355]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    }
[e7cc147]1522  }
1523
[f6e355]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  }
[e7cc147]1530
[f6e355]1531  if (char(basering)<>0)
[e7cc147]1532  {
[f6e355]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    }
[e7cc147]1540  }
[f6e355]1541 
1542  // check for A_k singularity
1543  // ----------------------------------------
1544  if (mult(std(f))==2)
[e7cc147]1545  {
[f6e355]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
[e7cc147]1554    {
[f6e355]1555      return( list(ES_Id,ESfix_Id) );
[e7cc147]1556    }
1557    else
1558    {
[f6e355]1559      return( list(ES_Id,ESfix_Id,ES_Id) );
[e7cc147]1560    }
1561  }
1562
[f6e355]1563  // check for D_k singularity
1564  // ----------------------------------------
1565  if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3)
[e7cc147]1566  {
[f6e355]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    }
[e7cc147]1582  }
1583
[f6e355]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
[e7cc147]1595  {
[f6e355]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++)
[e7cc147]1602    {
[f6e355]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));
[e7cc147]1630    }
1631  }
[f6e355]1632  else
[e7cc147]1633  {
[f6e355]1634    dbprint(i_print,"polynomial is Newton degenerate !");
1635    dbprint(i_print,"");
[e7cc147]1636  }
[f6e355]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
[e7cc147]1660
[f6e355]1661  poly F=imap(old_ring,f);
1662  ideal kbase_versal=imap(old_ring,kbase_versal);
1663  for (i=1; i<=s; i++)
[e7cc147]1664  {
[f6e355]1665    F=F+var(i)*kbase_versal[i];
[e7cc147]1666  }
[f6e355]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)));
[e7cc147]1700
[f6e355]1701      g=subst(F-Fred,var(i),1);
1702      ES_Id=ES_Id, g;
1703      F=Fred;
[e7cc147]1704
[f6e355]1705      g=subst(Ftriv-Ftrivred,var(i),1);
1706      ES_Id_all_triv=ES_Id_all_triv, g;
1707      Ftriv=Ftrivred;
1708    }
[e7cc147]1709
[f6e355]1710    setring old_ring;
1711    // back to original ring
[e7cc147]1712
[f6e355]1713    ideal ES_Id = imap(@Px,ES_Id);   
1714    ES_Id = interred(ES_Id);
[e7cc147]1715
[f6e355]1716    ideal ES_Id_all_triv = imap(@Px,ES_Id_all_triv);   
1717    ES_Id_all_triv = interred(ES_Id_all_triv);
[e7cc147]1718
[f6e355]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);
[e7cc147]1722
[f6e355]1723    if (typ==0)
[e7cc147]1724    {
[f6e355]1725      return(list(ES_Id,ESfix_Id));
[e7cc147]1726    }
[f6e355]1727    else
[e7cc147]1728    {
[f6e355]1729      return(list(ES_Id,ESfix_Id,ES_Id_all_triv));
[e7cc147]1730    }
[f6e355]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++)
[e7cc147]1751    {
[f6e355]1752      Fred=reduce(p_F,std(var(i)));
1753      Fredtriv=reduce(p_Ftriv,std(var(i)));
[e7cc147]1754
[f6e355]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;
[e7cc147]1762
1763    }
1764
[f6e355]1765    dbprint(i_print,"// ");
1766    dbprint(i_print,"// back to the original ring");
[e7cc147]1767
[f6e355]1768    setring old_ring;
1769    // back to original ring
[e7cc147]1770
[f6e355]1771    ideal ES_Id = imap(AuxRing,ES_Id);
1772    ES_Id = interred(ES_Id);
[e7cc147]1773
[f6e355]1774    ideal ES_Id_all_triv = imap(AuxRing,ES_Id_all_triv);
1775    ES_Id_all_triv = interred(ES_Id_all_triv);
[e7cc147]1776
[f6e355]1777    kill @Px;
1778    kill AuxRing;
[e7cc147]1779
[f6e355]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));
[e7cc147]1792    }
1793  }
1794}
[f6e355]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
[e7cc147]1813///////////////////////////////////////////////////////////////////////////////
1814
[f6e355]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
[7b3971]1818         and the characteristic of the basering does not divide mult(f). @*
[f6e355]1819         If nv is the number of variables of the basering, then the first
1820         nv-2 variables are the deformation parameters. @*
[e7cc147]1821         If the basering is a qring, ideal(basering) must only depend
1822         on the deformation parameters.
[f6e355]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
[50cbdc]1826@format
[f6e355]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;
[7b3971]1832@end format
[f6e355]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
[7b3971]1844         @code{execute}.
[f6e355]1845         printlevel>=2 displays additional information.
1846SEE ALSO: esIdeal, isEquising
1847KEYWORDS: equisingularity stratum
1848EXAMPLE: example esStratum; shows examples.
[e7cc147]1849"
1850{
[f6e355]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;
1860   poly F = (x2+2xy+y2+x5)+ax+by+cx2+dxy+ey2+fx3+gx4;
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
[e7cc147]1866
[f6e355]1867   ideal I = f-fa,e+b;
1868   qring q = std(I);
1869   poly F = imap(r,F);
1870   esStratum(F);
1871   printlevel=p;
1872}
[e7cc147]1873
[f6e355]1874///////////////////////////////////////////////////////////////////////////////
[e7cc147]1875
[f6e355]1876proc isEquising (poly p_F, list #)
1877"USAGE:   isEquising(F[,m,L]); F poly, m int, L list
1878ASSUME:  F defines a deformation of a reduced bivariate polynomial f
1879         and the characteristic of the basering does not divide mult(f). @*
1880         If nv is the number of variables of the basering, then the first
1881         nv-2 variables are the deformation parameters. @*
1882         If the basering is a qring, ideal(basering) must only depend
1883         on the deformation parameters.
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
1891         @code{execute}.
1892         printlevel>=2 displays additional information.
1893EXAMPLE: example isEquising; shows examples.
1894"
1895{
1896  int check=esComputation (1,p_F,#);
1897  return(check);
[e7cc147]1898}
1899example
1900{
1901   "EXAMPLE:"; echo=2;
[f6e355]1902   ring r = 0,(a,b,x,y),ds;
[7b3971]1903   poly F = (x2+2xy+y2+x5)+ay3+bx5;
[e7cc147]1904   isEquising(F);
[7b3971]1905   ideal I = ideal(a);
[e7cc147]1906   qring q = std(I);
1907   poly F = imap(r,F);
[7b3971]1908   isEquising(F);
[f6e355]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
[e7cc147]1915}
[f6e355]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
[558eb2]2054*/
Note: See TracBrowser for help on using the repository browser.