source: git/Singular/LIB/equising.lib @ 6fe9a5

spielwiese
Last change on this file since 6fe9a5 was bf42f0, checked in by Hans Schoenemann <hannes@…>, 13 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@13483 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 53.2 KB
Line 
1//
2version="$Id$";
3category="Singularities";
4info="
5LIBRARY:  equising.lib  Equisingularity Stratum of a Family of Plane Curves
6AUTHOR:   Christoph Lossen, lossen@mathematik.uni-kl.de
7          Andrea Mindnich, mindnich@mathematik.uni-kl.de
8
9PROCEDURES:
10 tau_es(f);             codim of mu-const stratum in semi-universal def. base
11 esIdeal(f);            (Wahl's) equisingularity ideal of f
12 esStratum(F[,m,L]);    equisingularity stratum of a family F
13 isEquising(F[,m,L]);   tests if a given deformation is equisingular
14
15 control_Matrix(M);     computes list of blowing-up data
16";
17
18LIB "hnoether.lib";
19LIB "poly.lib";
20LIB "elim.lib";
21LIB "deform.lib";
22LIB "sing.lib";
23
24////////////////////////////////////////////////////////////////////////////////
25//
26//  The following (static) procedures are used by esComputation
27//
28////////////////////////////////////////////////////////////////////////////////
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}
39
40////////////////////////////////////////////////////////////////////////////////
41// exchanges the variables x and y in the polynomial f
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}
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);
56{
57  intvec w=xyVector();
58  poly Fd=jet(F,m,w);
59  return(Fd);
60}
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"
93{
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);
160}
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,
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++)
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);
245}
246
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++)
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);
272}
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
282{
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);
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///////////////////////////////////////////////////////////////////////////////
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)
341{
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(size(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(size(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);
478}
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
546///////////////////////////////////////////////////////////////////////////////
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)
551{
552  int i_print = printlevel - voice + 3;
553  int i_ord;
554
555  if (p_f == 0)
556  {
557    print("Input is a 'deformation'  of the zero polynomial!");
558    return(1);
559  }
560
561  i_ord = mindeg1(p_f);
562
563  if (number(i_ord) == 0)
564  {
565    print("Characteristic of coefficient field "
566            +"divides order of zero-fiber !");
567    return(1);
568  }
569
570  if (squarefree(p_f) != p_f)
571  {
572    print("Original polynomial (= zero-fiber) is not reduced!");
573    return(1);
574  }
575
576  return(0);
577}
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
614///////////////////////////////////////////////////////////////////////////////
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 #)
620{
621  intvec ov=option(get);  // store options set at beginning
622  option(redSB);
623  // Initialize variables
624  int branch=1;
625  int blowup=1;
626  int auxVar=1;
627  int nVars;
628
629  intvec upper_bound, upper_bound_old, fertig, soll;
630  list blowup_string;
631  int i_print= printlevel-voice+2;
632
633  int no_b, number_of_branches, swapped;
634  int i,j,k,m, counter, dummy;
635  string helpStr = "";
636  string ordStr = "";
637  string MinPolyStr = "";
638
639  if (nvars(basering)<=2)
640  {
641    print("family is trivial (no deformation parameters)!");
642    if (typ==1) //isEquising
643    {
644      option(set,ov);
645      return(1);
646    }
647    else
648    {
649      option(set,ov);
650      return(list(ideal(0),0));
651    }
652  }
653
654  if (size(#)>0)
655  {
656    if (typeof(#[1])=="int")
657    {
658      def artin_bd=#[1];  // compute modulo maxideal(artin_bd)
659      if (artin_bd <= 1)
660      {
661        print("Do you really want to compute over Basering/maxideal("
662              +string(artin_bd)+") ?");
663        print("No computation performed !");
664        if (typ==1) //isEquising
665        {
666          option(set,ov);
667          return(1);
668        }
669        else
670        {
671          option(set,ov);
672          return(list(ideal(0),int(1)));
673        }
674      }
675      if (size(#)>1)
676      {
677        if (typeof(#[2])=="list")
678        {
679          def @L=#[2];  // is assumed to be the Hamburger-Noether matrix
680        }
681      }
682    }
683    else
684    {
685      if (typeof(#)=="list")
686      {
687        def @L=#;  // is assumed to be the Hamburger-Noether matrix
688      }
689    }
690  }
691  int ring_is_changed;
692  def old_ring=basering;
693  if(defined(@L)<=0)
694  {
695    // define a new ring without deformation-parameters and change to it:
696    string str;
697    string minpolyStr = string(minpoly);
698    str = " ring HNERing = (" + charstr(basering) + "), (x,y), ls;";
699    execute (str);
700    str = "minpoly ="+ minpolyStr+";";
701    execute(str);
702    ring_is_changed=1;
703    // Basering changed to HNERing (variables x,y, with ls ordering)
704
705    k=nvars(old_ring);
706    matrix Map_Phi[1][k];
707    Map_Phi[1,k-1]=x;
708    Map_Phi[1,k]=y;
709    map phi=old_ring,Map_Phi;
710    poly f=phi(p_F);
711
712    // Heuristics: if x,y are transversal parameters then computation of HNE
713    // can be much faster when exchanging variables...!
714    if (2*size(coeffs(f,x))<size(coeffs(f,y)))
715    {
716      swapped=1;
717      f=swapXY(f);
718    }
719
720    int error=checkPoly(f);
721    if (error)
722    {
723      setring old_ring;
724      if (typ==1) //isEquising
725      {
726        print("Return value (=0) has no meaning!");
727        option(set,ov);
728        return(0);
729      }
730      else
731      {
732        option(set,ov);
733        return(list( ideal(0),error));
734      }
735    }
736
737    dbprint(i_print,"// ");
738    dbprint(i_print,"// Compute HN expansion");
739    dbprint(i_print,"// ---------------------");
740    i=printlevel;
741    printlevel=printlevel-5;
742    list LLL=hnexpansion(f);
743
744    if (size(LLL)==0) { // empty list returned by hnexpansion
745      setring old_ring;
746      print(i_print,"Unable to compute HN expansion !");
747      if (typ==1) //isEquising
748      {
749        print("Return value (=0) has no meaning!");
750        option(set,ov);
751        return(0);
752      }
753      else
754      {
755        option(set,ov);
756        return(list(ideal(0),int(1)));
757      }
758      option(set,ov);
759      return(0);
760    }
761    else
762    {
763      if (typeof(LLL[1])=="ring") {
764        def HNering = LLL[1];
765        setring HNering;
766        def @L=stripHNE(hne);
767      }
768      else {
769        def @L=stripHNE(LLL);
770      }
771    }
772    printlevel=i;
773    dbprint(i_print,"// finished");
774    dbprint(i_print,"// ");
775  }
776  def HNEring=basering;
777  list M=multsequence(@L);
778  M=control_Matrix(M);     // this returns the 4 control matrices
779  def maxDeg=M[4];
780
781  list L1=inf_Tangents(@L,nrows(M[1]));
782  matrix B=L1[1];
783  intvec V=L1[2];
784  kill L1;
785
786  // if we have computed the HNE for f after swapping x and y, we have
787  // to reinterprete the (swap) matrix V:
788  if (swapped==1)
789  {
790    for (i=1;i<=size(V);i++) { V[i]=V[i]-1; } // turns 0 into -1, 1 into 0
791  }
792
793  // Determine maximal number of needed auxiliary parameters (free tangents):
794  no_b=Determine_no_b(M[3],B);
795
796  // test whether HNexpansion needed field extension....
797  string minPolyStr = "";
798  if (minpoly !=0)
799  {
800    minPolyStr = makeMinPolyString("a");
801    minPolyStr = "minpoly =" + minPolyStr + ";";
802  }
803
804  setring old_ring;
805
806  def myRing=createMyRing_new(p_F,"dp",minPolyStr,no_b);
807  setring myRing;  // comes with mIdeal
808  map hole=HNEring,mIdeal;
809  // basering has changed to myRing, in particular, the "old"
810  // variable names, e.g., A,B,C,z,y are replaced by t(1),t(2),t(3),x,y
811
812  ideal bNodes;
813
814  // Initialize some variables:
815  map phi;
816  poly G, F_save;
817  poly b_dummy;
818  ideal J,Jnew,final_Map;
819  number_of_branches=ncols(M[1]);
820  for (i=1;i<=number_of_branches;i++)
821  {
822    poly F(i);
823    ideal bl_Map(i);
824  }
825  upper_bound[number_of_branches]=0;
826  upper_bound[1]=number_of_branches;
827  upper_bound_old=upper_bound;
828  fertig[number_of_branches]=0;
829  for (i=1;i<=number_of_branches;i++){ soll[i]=1; }
830
831  // Hole:  B = matrix of blowup points
832  if (ring_is_changed==0) { matrix B=hole(B); }
833  else                    { matrix B=imap(HNEring,B); }
834  m=M[1][blowup,branch];    // multiplicity at 0
835
836  // now, we start by checking equimultiplicity along trivial section
837  poly Fm=m_Jet(p_F,m-1);
838
839  matrix coef_Mat = coef(Fm,xy);
840  Jnew=coef_Mat[2,1..ncols(coef_Mat)];
841  J=J,Jnew;
842
843  if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
844                         // deformation parameters can be cutted off
845  {
846    J=jet(J,artin_bd-1);
847  }
848
849  J=interred(J);
850  if (defined(artin_bd)) { J=jet(J,artin_bd-1); }
851
852// J=std(J);
853
854  if (typ==1) // isEquising
855  {
856    if(ideal(nselect(J,1..no_b))<>0)
857    {
858      setring old_ring;
859      option(set,ov);
860      return(0);
861    }
862  }
863
864  F(1)=p_F;
865
866  // and reduce the remaining terms in F(1):
867  bl_Map(1)=maxideal(1);
868
869  attrib(J,"isSB",1);
870  bl_Map(1)=reduce(bl_Map(1),J);
871  attrib(J,"isSB",0);
872
873  phi=myRing,bl_Map(1);
874  F(1)=phi(F(1));
875
876  // simplify F(1)
877  attrib(J,"isSB",1);
878  F(1)=reduce(F(1),J);
879  attrib(J,"isSB",0);
880
881  // now we compute the m-jet:
882  Fm=m_Jet(F(1),m);
883
884  G=1;
885  counter=branch;
886  k=upper_bound[branch];
887
888  F_save=F(1);  // is truncated differently in the following loop
889
890  while(counter<=k)
891  {
892    F(counter)=m_Jet(F_save,maxDeg[blowup,counter]);
893    if (V[counter]==0) // 2nd ring variable is tangent to this branch
894    {
895      G=G*(y-(b(auxVar)+B[blowup,counter])*x)^(M[3][blowup,counter]);
896    }
897    else // 1st ring variable is tangent to this branch
898    {
899      G=G*(x-(b(auxVar)+B[blowup,counter])*y)^(M[3][blowup,counter]);
900      F(counter)=swapXY(F(counter));
901    }
902    bl_Map(counter)=maxideal(1);
903    bl_Map(counter)[nvars(basering)]=xy+(b(auxVar)+B[blowup,counter])*x;
904
905    bNodes[counter]=b(auxVar);
906
907    auxVar=auxVar+1;
908    upper_bound[counter]=counter+M[2][blowup+1,counter]-1;
909    counter=counter+M[2][blowup+1,counter];
910
911  }
912
913  list LeadDataFm=determine_coef(Fm);
914  def LeadDataG=coef(G,xy);
915
916  for (i=1; i<=ncols(LeadDataG); i++)
917  {
918    if (LeadDataG[1,i]==LeadDataFm[2])
919    {
920      poly LeadG = LeadDataG[2,i];  // determine the coefficient of G
921      i=ncols(LeadDataG)+1;
922    }
923  }
924
925  G=LeadDataFm[1]*G-LeadG*Fm;  // leading terms in y should cancel...
926
927  coef_Mat = coef(G,xy);
928  Jnew=coef_Mat[2,1..ncols(coef_Mat)];
929
930  // simplification of Jnew
931
932  if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
933                         // deformation parameters can be cutted off
934  {
935    Jnew=jet(Jnew,artin_bd-1);
936  }
937  Jnew=interred(Jnew);
938  if (defined(artin_bd)) { Jnew=jet(Jnew,artin_bd-1); }
939  J=J,Jnew;
940
941  if (typ==1) // isEquising
942  {
943    if(ideal(nselect(J,1..no_b))<>0)
944    {
945      setring old_ring;
946      option(set,ov);
947      return(0);
948    }
949  }
950
951  while (fertig<>soll and blowup<nrows(M[3]))
952  {
953    upper_bound_old=upper_bound;
954    dbprint(i_print,"// Blowup Step "+string(blowup)+" completed");
955    blowup=blowup+1;
956
957    for (branch=1;branch<=number_of_branches;branch=branch+1)
958    {
959      Jnew=0;
960
961      // First we check if the branch still has to be considered:
962      if (branch==upper_bound_old[branch] and fertig[branch]<>1)
963      {
964        if (M[3][blowup-1,branch]==1 and
965               ((B[blowup,branch]<>x and B[blowup,branch]<>y)
966            or (blowup==nrows(M[3])) ))
967        {
968          fertig[branch]=1;
969          dbprint(i_print,"// 1 branch finished");
970        }
971      }
972
973      if (branch<=upper_bound_old[branch] and fertig[branch]<>1)
974      {
975        for (i=branch;i>=1;i--)
976        {
977          if (M[1][blowup-1,i]<>0)
978          {
979            m=M[1][blowup-1,i]; // multiplicity before blowup
980            i=0;
981          }
982        }
983
984        // we blow up the branch and take the strict transform:
985        attrib(J,"isSB",1);
986        bl_Map(branch)=reduce(bl_Map(branch),J);
987        attrib(J,"isSB",0);
988
989        phi=myRing,bl_Map(branch);
990        F(branch)=phi(F(branch))/x^m;
991
992        // simplify F
993        attrib(Jnew,"isSB",1);
994
995        F(branch)=reduce(F(branch),Jnew);
996        attrib(Jnew,"isSB",0);
997
998        m=M[1][blowup,branch]; // multiplicity after blowup
999        Fm=m_Jet(F(branch),m); // homogeneous part of lowest degree
1000
1001
1002        // we check for Fm=F[k]*...*F[k+s] where
1003        //
1004        //    F[j]=(y-b'(j)*x)^m(j), respectively F[j]=(-b'(j)*y+x)^m(j)
1005        //
1006        // according to the entries m(j)= M[3][blowup,j] and
1007        //                          b'(j) mod m_A = B[blowup,j]
1008        // computed from the HNE of the special fibre of the family:
1009        G=1;
1010        counter=branch;
1011        k=upper_bound[branch];
1012
1013        F_save=F(branch);
1014
1015        while(counter<=k)
1016        {
1017          F(counter)=m_Jet(F_save,maxDeg[blowup,counter]);
1018
1019          if (B[blowup,counter]<>x and B[blowup,counter]<>y)
1020          {
1021            G=G*(y-(b(auxVar)+B[blowup,counter])*x)^(M[3][blowup,counter]);
1022            bl_Map(counter)=maxideal(1);
1023            bl_Map(counter)[nvars(basering)]=
1024                                        xy+(b(auxVar)+B[blowup,counter])*x;
1025            bNodes[counter]=b(auxVar);
1026            auxVar=auxVar+1;
1027          }
1028          else
1029          {
1030            if (B[blowup,counter]==x)
1031            {
1032              G=G*x^(M[3][blowup,counter]);  // branch has tangent x !!
1033              F(counter)=swapXY(F(counter)); // will turn x to y for blow up
1034              bl_Map(counter)=maxideal(1);
1035              bl_Map(counter)[nvars(basering)]=xy;
1036            }
1037            else
1038            {
1039              G=G*y^(M[3][blowup,counter]); // tangent has to be y
1040              bl_Map(counter)=maxideal(1);
1041              bl_Map(counter)[nvars(basering)]=xy;
1042            }
1043            bNodes[counter]=0;
1044          }
1045          upper_bound[counter]=counter+M[2][blowup+1,counter]-1;
1046          counter=counter+M[2][blowup+1,counter];
1047        }
1048        G=determine_coef(Fm)[1]*G-Fm;  // leading terms in y should cancel
1049        coef_Mat = coef(G,xy);
1050        Jnew=coef_Mat[2,1..ncols(coef_Mat)];
1051        if (defined(artin_bd)) // the artin_bd-th power of the maxideal of
1052                              // deformation parameters can be cutted off
1053        {
1054          Jnew=jet(Jnew,artin_bd-1);
1055        }
1056
1057        // simplification of J
1058        Jnew=interred(Jnew);
1059
1060        J=J,Jnew;
1061        if (typ==1) // isEquising
1062        {
1063          if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
1064          if(ideal(nselect(J,1..no_b))<>0)
1065          {
1066            setring old_ring;
1067            option(set,ov);
1068            return(0);
1069          }
1070        }
1071      }
1072    }
1073    if (number_of_branches>=2)
1074    {
1075      J=interred(J);
1076      if (typ==1) // isEquising
1077      {
1078        if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
1079        if(ideal(nselect(J,1..no_b))<>0)
1080        {
1081          setring old_ring;
1082          option(set,ov);
1083          return(0);
1084        }
1085      }
1086    }
1087  }
1088
1089  // Computation for all equimultiple sections being trivial (I^s(f))
1090  ideal Jtriv=J;
1091  for (i=1;i<=no_b; i++)
1092  {
1093    if (reduce(b(i),std(bNodes))!=0){
1094      Jtriv=subst(Jtriv,b(i),0);
1095    }
1096  }
1097  Jtriv=std(Jtriv);
1098
1099
1100
1101  dbprint(i_print,"// ");
1102  dbprint(i_print,"// Elimination starts:");
1103  dbprint(i_print,"// -------------------");
1104
1105  poly gg;
1106  int b_left=no_b;
1107
1108  for (i=1;i<=no_b; i++)
1109  {
1110    attrib(J,"isSB",1);
1111    gg=reduce(b(i),J);
1112    if (gg==0)
1113    {
1114      b_left = b_left-1;  // another b(i) has to be 0
1115    }
1116    J = subst(J, b(i), gg);
1117    attrib(J,"isSB",0);
1118  }
1119  J=simplify(J,10);
1120  if (typ==1) // isEquising
1121  {
1122    if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
1123    if(ideal(nselect(J,1..no_b))<>0)
1124    {
1125      setring old_ring;
1126      option(set,ov);
1127      return(0);
1128    }
1129  }
1130
1131  //new CL 11/06:  check in which equations b(k) appears and remove those b(k)
1132  //               which appear in exactly one of the equations (by removing this
1133  //               equation)
1134  dbprint(i_print,"// ");
1135  dbprint(i_print,"// Remove superfluous equations:");
1136  dbprint(i_print,"// -----------------------------");
1137  int Z,App_in;
1138  ideal J_Tmp;
1139  int ncJ=ncols(J);
1140
1141  intmat Mdet[ncJ][1];
1142  for (Z=1;Z<=ncJ;Z++){ Mdet[Z,1]=Z; }
1143
1144  for (i=1;i<=no_b; i++)
1145  {
1146    ideal b_appears_in(i);            // Eintraege sind spaeter 1 oder 0
1147    intmat b_app_in(i)[1][ncJ];       // Eintraege sind spaeter 1 oder 0
1148    b_appears_in(i)[ncJ]=0;
1149    J_Tmp = matrix(J)-subst(J,b(i),0);
1150    for (Z=1; Z<=ncJ; Z++) {
1151      if (J_Tmp[Z]<>0) {      // b(i) appear in J_Tmp[Z]
1152        b_appears_in(i)[Z]=1;
1153        b_app_in(i)[1,Z]=1;
1154      }
1155    }
1156    if (size(b_appears_in(i))==1) { //b(i) appears only in one J_Tmp[Z]
1157      App_in = (b_app_in(i)*Mdet)[1,1];  // determines Z
1158      J[App_in]=0;
1159      b_appears_in(i)[App_in]=0;
1160      b_app_in(i)[1,App_in]=0;
1161    }
1162  }
1163
1164  for (i=1;i<=no_b; i++)
1165  {
1166    if (size(b_appears_in(i))==1) { //b(i) appears only in one J_Tmp[Z]
1167      App_in = (b_app_in(i)*Mdet)[1,1];  // determines Z
1168      J[App_in]=0;
1169      b_appears_in(i)[App_in]=0;
1170      b_app_in(i)[1,Z]=1;
1171      i=0;
1172    }
1173  }
1174
1175  Jtriv = nselect(Jtriv,1..no_b);
1176  ideal J_no_b = nselect(J,1..no_b);
1177  if (size(J) > size(J_no_b))
1178  {
1179    dbprint(i_print,"// std computation started");
1180    // some b(i) didn't appear in linear conditions and have to be eliminated
1181    if (defined(artin_bd))
1182    {
1183      // first we make the ring smaller (removing variables, which are
1184      // forced to 0 by J
1185      list LL=make_ring_small(J);
1186      ideal Shortmap=LL[2];
1187      minPolyStr = "";
1188      if (minpoly !=0)
1189      {
1190        minPolyStr = "minpoly = "+string(minpoly);
1191      }
1192      ordStr = "dp(" + string(b_left) + "),dp";
1193      ideal qId = ideal(basering);
1194
1195      helpStr = "ring Shortring = ("
1196                + charstr(basering) + "),("+ LL[1] +") , ("+ ordStr  +");";
1197      execute(helpStr);
1198      execute(minPolyStr);
1199      // ring has changed to "Shortring"
1200
1201      ideal MM=maxideal(artin_bd);
1202      MM=subst(MM,x,0);
1203      MM=subst(MM,y,0);
1204      MM=simplify(MM,2);
1205      dbprint(i_print-1,"// maxideal("+string(artin_bd)+") has "
1206                         +string(size(MM))+" elements");
1207      dbprint(i_print-1,"//");
1208
1209      // we change to the qring mod m^artin_bd
1210      // first, we have to check if we were in a qring when starting
1211      ideal qId = imap(myRing, qId);
1212      if (qId == 0)
1213      {
1214         attrib(MM,"isSB",1);
1215         qring QQ=MM;
1216      }
1217      else
1218      {
1219         qId=qId,MM;
1220         qring QQ = std(qId);
1221      }
1222
1223      ideal Shortmap=imap(myRing,Shortmap);
1224      map phiphi=myRing,Shortmap;
1225
1226      ideal J=phiphi(J);
1227      option(redSB);
1228      J=std(J);
1229      J=nselect(J,1..no_b);
1230
1231      setring myRing;
1232      // back to "myRing"
1233
1234      J=nselect(J,1..no_b);
1235      Jnew=imap(QQ,J);
1236
1237      J=J,Jnew;
1238      J=interred(J);
1239      if (defined(artin_bd)){ J=jet(J,artin_bd-1); }
1240    }
1241    else
1242    {
1243      J=std(J);
1244      J=nselect(J,1..no_b);
1245      if (defined(artin_bd)){ J=jet(J,artin_bd-1); }
1246    }
1247  }
1248
1249  dbprint(i_print,"// finished");
1250  dbprint(i_print,"// ");
1251
1252  minPolyStr = "";option(set,ov);
1253  if (minpoly !=0)
1254  {
1255   minPolyStr = "minpoly = "+string(minpoly);
1256  }
1257
1258  kill HNEring;
1259
1260  if (typ==1) // isEquising
1261  {
1262    if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
1263    if(J<>0)
1264    {
1265      setring old_ring;
1266      option(set,ov);
1267      return(0);
1268    }
1269    else
1270    {
1271      setring old_ring;
1272      option(set,ov);
1273      return(1);
1274    }
1275  }
1276
1277  setring old_ring;
1278  // we are back in the original ring
1279
1280  if (npars(myRing)<>0)
1281  {
1282    ideal qIdeal = ideal(basering);
1283    helpStr = "ring ESSring = ("
1284                 + string(char(basering))+ "," + parstr(myRing) +
1285                 ") , ("+ varstr(basering)+") , ("+ ordstr(basering) +");";
1286    execute(helpStr);
1287    execute(minPolyStr);
1288    // basering has changed to ESSring
1289
1290    ideal qIdeal = fetch(old_ring, qIdeal);
1291    if(qIdeal != 0)
1292    {
1293      def r_base = basering;
1294      kill ESSring;
1295      qring ESSring = std(qIdeal);
1296    }
1297    kill qIdeal;
1298
1299    ideal SSS;
1300    for (int ii=1;ii<=nvars(basering);ii++)
1301    {
1302      SSS[ii+no_b]=var(ii);
1303    }
1304    map phi=myRing,SSS;   // b(i) variables are mapped to zero
1305
1306    ideal ES=phi(J);
1307    ideal ES_all_triv=phi(Jtriv);
1308    kill phi;
1309
1310    if (defined(p_F)<=0)
1311    {
1312      poly p_F=fetch(old_ring,p_F);
1313      export(p_F);
1314    }
1315    export(ES);
1316    export(ES_all_triv);
1317    setring old_ring;
1318    dbprint(i_print+2,"
1319// 'esStratum' created a list M of a ring and an integer.
1320// To access the ideal defining the equisingularity stratum, type:
1321        def ESSring = M[1]; setring ESSring;  ES; ");
1322
1323    option(set,ov);
1324    return(list(ESSring,0));
1325  }
1326  else
1327  {
1328    // no new ring definition necessary
1329    ideal SSS;
1330    for (int ii=1;ii<=nvars(basering);ii++)
1331    {
1332      SSS[ii+no_b]=var(ii);
1333    }
1334    map phi=myRing,SSS;  // b(i) variables are mapped to zero
1335
1336    ideal ES=phi(J);
1337    ideal ES_all_triv=phi(Jtriv);
1338    kill phi;
1339
1340    setring old_ring;
1341    dbprint(i_print,"// output of 'esStratum' is a list consisting of:
1342//    _[1][1] = ideal defining the equisingularity stratum
1343//    _[1][2] = ideal defining the part of the equisingularity stratum
1344//              where all equimultiple sections are trivial
1345//    _[2] = 0");
1346
1347    option(set,ov);
1348    return(list(list(ES,ES_all_triv),0));
1349  }
1350
1351}
1352
1353////////////////////////////////////////////////////////////////////////////////
1354
1355proc tau_es (poly f,list #)
1356"USAGE:   tau_es(f); f poly
1357ASSUME:  f is a reduced bivariate polynomial, the basering has precisely
1358         two variables, is local and no qring.
1359RETURN:  int, the codimension of the mu-const stratum in the semi-universal
1360         deformation base.
1361NOTE:    printlevel>=1 displays additional information.
1362         When called with any additional parameter, the computation of the
1363         Milnor number is avoided (no check for NND).
1364SEE ALSO: esIdeal, tjurina, invariants
1365EXAMPLE: example tau_es; shows an example.
1366"
1367{
1368  int i,j,k,s;
1369  int slope_x, slope_y, upper;
1370  int i_print = printlevel - voice + 3;
1371  string MinPolyStr;
1372
1373  // some checks first
1374  if ( nvars(basering)<>2 )
1375  {
1376    print("// basering has not the correct number (two) of variables !");
1377    print("// computation stopped");
1378    return(0);
1379  }
1380  if ( mult(std(1+var(1)+var(2))) <> 0)
1381  {
1382    print("// basering is not local !");
1383    print("// computation stopped");
1384    return(0);
1385  }
1386
1387  if (mult(std(f))<=1)
1388  {
1389    // f is rigid
1390    return(0);
1391  }
1392
1393  if ( deg(squarefree(f))!=deg(f) )
1394  {
1395    print("// input polynomial was not reduced");
1396    print("// try    squarefree(f);   first");
1397    return(0);
1398  }
1399
1400  def old_ring=basering;
1401  execute("ring @myRing=("+charstr(basering)+"),("+varstr(basering)+"),ds;");
1402  poly f=imap(old_ring,f);
1403
1404  ideal Jacobi_Id = jacob(f);
1405
1406  // check for A_k singularity
1407  // ----------------------------------------
1408  if (mult(std(f))==2)
1409  {
1410    dbprint(i_print-1,"// ");
1411    dbprint(i_print-1,"// polynomial defined A_k singularity");
1412    dbprint(i_print-1,"// ");
1413    return( vdim(std(Jacobi_Id)) );
1414  }
1415
1416  // check for D_k singularity
1417  // ----------------------------------------
1418  if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3)
1419  {
1420    dbprint(i_print,"// ");
1421    dbprint(i_print,"// polynomial defined D_k singularity");
1422    dbprint(i_print,"// ");
1423    ideal ES_Id = f, jacob(f);
1424    return( vdim(std(Jacobi_Id)));
1425  }
1426
1427
1428  if (size(#)==0)
1429  {
1430    // check if Newton polygon non-degenerate
1431    // ----------------------------------------
1432    Jacobi_Id=std(Jacobi_Id);
1433    int mu = vdim(Jacobi_Id);
1434    poly f_tilde=f+var(1)^mu+var(2)^mu;  //to obtain convenient Newton-polygon
1435
1436    list NP=newtonpoly(f_tilde);
1437    dbprint(i_print-1,"// Newton polygon:");
1438    dbprint(i_print-1,NP);
1439    dbprint(i_print-1,"");
1440
1441    if(is_NND(f,mu,NP))          // f is Newton non-degenerate
1442    {
1443      upper=NP[1][2];
1444      ideal ES_Id= x^k*y^upper;
1445      dbprint(i_print-1,"polynomial is Newton non-degenerate");
1446      dbprint(i_print-1,"");
1447      k=0;
1448      for (i=1;i<=size(NP)-1;i++)
1449      {
1450        slope_x=NP[i+1][1]-NP[i][1];
1451        slope_y=NP[i][2]-NP[i+1][2];
1452        for (k=NP[i][1]+1; k<=NP[i+1][1]; k++)
1453        {
1454          while ( slope_x*upper + slope_y*k >=
1455                  slope_x*NP[i][2] + slope_y*NP[i][1])
1456          {
1457            upper=upper-1;
1458          }
1459          upper=upper+1;
1460          ES_Id=ES_Id, x^k*y^upper;
1461        }
1462      }
1463      ES_Id=std(ES_Id);
1464      dbprint(i_print-2,"ideal of monomials above Newton bd. is generated by:");
1465      dbprint(i_print-2,ES_Id);
1466      ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
1467      ES_Id = ES_Id, Jacobi_Id;
1468      ES_Id = std(ES_Id);
1469      dbprint(i_print-1,"// ");
1470      dbprint(i_print-1,"// Equisingularity ideal is computed!");
1471      dbprint(i_print-1,"");
1472      return(vdim(ES_Id));
1473    }
1474    else
1475    {
1476      dbprint(i_print-1,"polynomial is Newton degenerate !");
1477      dbprint(i_print-1,"");
1478    }
1479  }
1480
1481  // for Newton degenerate polynomials, we compute the HN expansion, and
1482  // count the number of free points .....
1483
1484  dbprint(i_print-1,"// ");
1485  dbprint(i_print-1,"// Compute HN expansion");
1486  dbprint(i_print-1,"// ---------------------");
1487  i=printlevel;
1488  printlevel=printlevel-5;
1489  if (2*size(coeffs(f,x))<size(coeffs(f,y)))
1490  {
1491    f=swapXY(f);
1492  }
1493  list LLL=hnexpansion(f);
1494  if (size(LLL)==0) { // empty list returned by hnexpansion
1495    setring old_ring;
1496    ERROR("Unable to compute HN expansion !");
1497  }
1498  else
1499  {
1500    if (typeof(LLL[1])=="ring") {
1501      def HNering = LLL[1];
1502      setring HNering;
1503      def @L=hne;
1504    }
1505    else {
1506      def @L=LLL;
1507    }
1508  }
1509  def HNEring=basering;
1510
1511  printlevel=i;
1512  dbprint(i_print-1,"// finished");
1513  dbprint(i_print-1,"// ");
1514
1515  list M=multsequence(@L);
1516  M=control_Matrix(M);     // this returns the 4 control matrices
1517  intmat Mult=M[1];
1518
1519  list L1=inf_Tangents(@L,nrows(M[1]));
1520  matrix B=L1[1];
1521
1522  // determine sum_i m_i(m_i+1)/2 (over inf. near points)
1523  int conditions=0;
1524  for (i=1;i<=nrows(Mult);i++)
1525  {
1526    for (j=1;j<=ncols(Mult);j++)
1527    {
1528      conditions=conditions+(Mult[i,j]*(Mult[i,j]+1)/2);
1529    }
1530  }
1531  int freePts=no_freePoints(M[1],B);
1532  int taues=conditions-freePts-2;
1533
1534  setring old_ring;
1535  return(taues);
1536}
1537example
1538{
1539   "EXAMPLE:"; echo=2;
1540   ring r=32003,(x,y),ds;
1541   poly f=(x4-y4)^2-x10;
1542   tau_es(f);
1543}
1544
1545
1546////////////////////////////////////////////////////////////////////////////////
1547
1548proc esIdeal (poly f,list #)
1549"USAGE:   esIdeal(f[,any]]); f poly
1550ASSUME:  f is a reduced bivariate polynomial, the basering has precisely
1551         two variables, is local and no qring, and the characteristic of
1552         the ground field does not divide mult(f).
1553RETURN:  if called with only one parameter: list of two ideals,
1554@format
1555          _[1]:  equisingularity ideal of f (in sense of Wahl),
1556          _[2]:  ideal of equisingularity with fixed position of the
1557                 singularity;
1558@end format
1559         if called with more than one parameter: list of three ideals,
1560@format
1561          _[1]:  equisingularity ideal of f (in sense of Wahl)
1562          _[2]:  ideal of equisingularity with fixed position of the
1563                 singularity;
1564          _[3]:  ideal of all g such that the deformation defined by f+eg
1565                 (e^2=0) is isomorphic to an equisingular deformation
1566                 of V(f) with all equimultiple sections being trivial.
1567@end format
1568NOTE:    if some of the above condition is not satisfied then return
1569         value is list(0,0).
1570SEE ALSO: tau_es, esStratum
1571KEYWORDS: equisingularity ideal
1572EXAMPLE: example esIdeal; shows examples.
1573"
1574{
1575
1576  int typ;
1577  if (size(#)>0) { typ=1; }   // I^s is also computed
1578  int i,k,s;
1579  int slope_x, slope_y, upper;
1580  int i_print = printlevel - voice + 3;
1581  string MinPolyStr;
1582
1583  // some checks first
1584  if ( nvars(basering)<>2 )
1585  {
1586    print("// basering has not the correct number (two) of variables !");
1587    print("// computation stopped");
1588    return(list(0,0));
1589  }
1590  if ( mult(std(1+var(1)+var(2))) <> 0)
1591  {
1592    print("// basering is not local !");
1593    print("// computation stopped");
1594    return(list(0,0));
1595  }
1596
1597  if (mult(std(f))<=1)
1598  {
1599    // f is rigid
1600    if (typ==0)
1601    {
1602      return(list(ideal(1),ideal(1)));
1603    }
1604    else
1605    {
1606      return(list(ideal(1),ideal(1),ideal(1)));
1607    }
1608  }
1609
1610  if ( deg(squarefree(f))!=deg(f) )
1611  {
1612    print("// input polynomial was not squarefree");
1613    print("// try    squarefree(f);   first");
1614    return(list(0,0));
1615  }
1616
1617  if (char(basering)<>0)
1618  {
1619    if (mult(std(f)) mod char(basering)==0)
1620    {
1621      print("// characteristic of ground field divides "
1622            + "multiplicity of polynomial !");
1623      print("// computation stopped");
1624      return(list(0,0));
1625    }
1626  }
1627
1628  // check for A_k singularity
1629  // ----------------------------------------
1630  if (mult(std(f))==2)
1631  {
1632    dbprint(i_print,"// ");
1633    dbprint(i_print,"// polynomial defined A_k singularity");
1634    dbprint(i_print,"// ");
1635    ideal ES_Id = f, jacob(f);
1636    ES_Id = interred(ES_Id);
1637    ideal ESfix_Id = f, maxideal(1)*jacob(f);
1638    ESfix_Id= interred(ESfix_Id);
1639    if (typ==0) // only for computation of I^es and I^es_fix
1640    {
1641      return( list(ES_Id,ESfix_Id) );
1642    }
1643    else
1644    {
1645      return( list(ES_Id,ESfix_Id,ES_Id) );
1646    }
1647  }
1648
1649  // check for D_k singularity
1650  // ----------------------------------------
1651  if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3)
1652  {
1653    dbprint(i_print,"// ");
1654    dbprint(i_print,"// polynomial defined D_k singularity");
1655    dbprint(i_print,"// ");
1656    ideal ES_Id = f, jacob(f);
1657    ES_Id = interred(ES_Id);
1658    ideal ESfix_Id = f, maxideal(1)*jacob(f);
1659    ESfix_Id= interred(ESfix_Id);
1660    if (typ==0) // only for computation of I^es and I^es_fix
1661    {
1662      return( list(ES_Id,ESfix_Id) );
1663    }
1664    else
1665    {
1666      return( list(ES_Id,ESfix_Id,ES_Id) );
1667    }
1668  }
1669
1670  // check if Newton polygon non-degenerate
1671  // ----------------------------------------
1672  int mu = milnor(f);
1673  poly f_tilde=f+var(1)^mu+var(2)^mu;  //to obtain a convenient Newton-polygon
1674
1675  list NP=newtonpoly(f_tilde);
1676  dbprint(i_print-1,"// Newton polygon:");
1677  dbprint(i_print-1,NP);
1678  dbprint(i_print-1,"");
1679
1680  if(is_NND(f,mu,NP))          // f is Newton non-degenerate
1681  {
1682    upper=NP[1][2];
1683    ideal ES_Id= x^k*y^upper;
1684    dbprint(i_print,"polynomial is Newton non-degenerate");
1685    dbprint(i_print,"");
1686    k=0;
1687    for (i=1;i<=size(NP)-1;i++)
1688    {
1689      slope_x=NP[i+1][1]-NP[i][1];
1690      slope_y=NP[i][2]-NP[i+1][2];
1691      for (k=NP[i][1]+1; k<=NP[i+1][1]; k++)
1692      {
1693        while ( slope_x*upper + slope_y*k >=
1694                slope_x*NP[i][2] + slope_y*NP[i][1])
1695        {
1696          upper=upper-1;
1697        }
1698        upper=upper+1;
1699        ES_Id=ES_Id, x^k*y^upper;
1700      }
1701    }
1702    ES_Id=std(ES_Id);
1703    dbprint(i_print-1,"ideal of monomials above Newton bd. is generated by:");
1704    dbprint(i_print-1,ES_Id);
1705    ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
1706    ES_Id = ES_Id, f, jacob(f);
1707    dbprint(i_print,"// ");
1708    dbprint(i_print,"// equisingularity ideal is computed!");
1709    if (typ==0)
1710    {
1711      return(list(ES_Id,ESfix_Id));
1712    }
1713    else
1714    {
1715      return(list(ES_Id,ESfix_Id,ES_Id));
1716    }
1717  }
1718  else
1719  {
1720    dbprint(i_print,"polynomial is Newton degenerate !");
1721    dbprint(i_print,"");
1722  }
1723
1724  def old_ring=basering;
1725
1726  dbprint(i_print,"// ");
1727  dbprint(i_print,"// versal deformation with triv. section");
1728  dbprint(i_print,"// =====================================");
1729  dbprint(i_print,"// ");
1730
1731  ideal JJ=maxideal(1)*jacob(f);
1732  ideal kbase_versal=kbase(std(JJ));
1733  s=size(kbase_versal);
1734  string ring_versal="ring @Px = ("+charstr(basering)+"),(t(1.."+string(s)+"),"
1735                        +varstr(basering)+"),(ds("+string(s)+"),"
1736                        +ordstr(basering)+");";
1737  MinPolyStr = string(minpoly);
1738
1739  execute(ring_versal);
1740  if (MinPolyStr<>"0")
1741  {
1742    MinPolyStr = "minpoly="+MinPolyStr;
1743    execute(MinPolyStr);
1744  }
1745  // basering has changed to @Px
1746
1747  poly F=imap(old_ring,f);
1748  ideal kbase_versal=imap(old_ring,kbase_versal);
1749  for (i=1; i<=s; i++)
1750  {
1751    F=F+var(i)*kbase_versal[i];
1752  }
1753  dbprint(i_print-1,F);
1754  dbprint(i_print-1,"");
1755
1756
1757  ideal ES_Id,ES_Id_all_triv;
1758  poly Ftriv=F;
1759
1760  dbprint(i_print,"// ");
1761  dbprint(i_print,"// Compute equisingularity Stratum over Spec(C[t]/t^2)");
1762  dbprint(i_print,"// ===================================================");
1763  dbprint(i_print,"// ");
1764  list M=esStratum(F,2);
1765  dbprint(i_print,"// finished");
1766  dbprint(i_print,"// ");
1767
1768  if (M[2]==1) // error occured during esStratum computation
1769  {
1770    print("Some error has occured during the computation");
1771    return(list(0,0));
1772  }
1773
1774  if ( typeof(M[1])=="list" )
1775  {
1776    int defpars = nvars(basering)-2;
1777    poly Fred,Ftrivred;
1778    poly g;
1779    F=reduce(F,std(M[1][1]));
1780    Ftriv=reduce(Ftriv,std(M[1][2]));
1781
1782    for (i=1; i<=defpars; i++)
1783    {
1784      Fred=reduce(F,std(var(i)));
1785      Ftrivred=reduce(Ftriv,std(var(i)));
1786
1787      g=subst(F-Fred,var(i),1);
1788      ES_Id=ES_Id, g;
1789      F=Fred;
1790
1791      g=subst(Ftriv-Ftrivred,var(i),1);
1792      ES_Id_all_triv=ES_Id_all_triv, g;
1793      Ftriv=Ftrivred;
1794    }
1795
1796    setring old_ring;
1797    // back to original ring
1798
1799    ideal ES_Id = imap(@Px,ES_Id);
1800    ES_Id = interred(ES_Id);
1801
1802    ideal ES_Id_all_triv = imap(@Px,ES_Id_all_triv);
1803    ES_Id_all_triv = interred(ES_Id_all_triv);
1804
1805    ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
1806    ES_Id = ES_Id, f, jacob(f);
1807    ES_Id_all_triv = ES_Id_all_triv, f, jacob(f);
1808
1809    if (typ==0)
1810    {
1811      return(list(ES_Id,ESfix_Id));
1812    }
1813    else
1814    {
1815      return(list(ES_Id,ESfix_Id,ES_Id_all_triv));
1816    }
1817  }
1818  else
1819  {
1820    def AuxRing=M[1];
1821
1822    dbprint(i_print,"// ");
1823    dbprint(i_print,"// change ring to ESSring");
1824
1825    setring AuxRing;  // contains p_F, ES
1826
1827    int defpars = nvars(basering)-2;
1828    poly Fred,Fredtriv;
1829    poly g;
1830    ideal ES_Id,ES_Id_all_triv;
1831
1832    poly p_Ftriv=p_F;
1833
1834    p_F=reduce(p_F,std(ES));
1835    p_Ftriv=reduce(p_Ftriv,std(ES_all_triv));
1836    for (i=1; i<=defpars; i++)
1837    {
1838      Fred=reduce(p_F,std(var(i)));
1839      Fredtriv=reduce(p_Ftriv,std(var(i)));
1840
1841      g=subst(p_F-Fred,var(i),1);
1842      ES_Id=ES_Id, g;
1843      p_F=Fred;
1844
1845      g=subst(p_Ftriv-Fredtriv,var(i),1);
1846      ES_Id_all_triv=ES_Id_all_triv, g;
1847      p_Ftriv=Fredtriv;
1848
1849    }
1850
1851    dbprint(i_print,"// ");
1852    dbprint(i_print,"// back to the original ring");
1853
1854    setring old_ring;
1855    // back to original ring
1856
1857    ideal ES_Id = imap(AuxRing,ES_Id);
1858    ES_Id = interred(ES_Id);
1859
1860    ideal ES_Id_all_triv = imap(AuxRing,ES_Id_all_triv);
1861    ES_Id_all_triv = interred(ES_Id_all_triv);
1862
1863    kill @Px;
1864    kill AuxRing;
1865
1866    ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
1867    ES_Id = ES_Id, f, jacob(f);
1868    ES_Id_all_triv = ES_Id_all_triv, f, jacob(f);
1869    dbprint(i_print,"// ");
1870    dbprint(i_print,"// equisingularity ideal is computed!");
1871    if (typ==0)
1872    {
1873      return(list(ES_Id,ESfix_Id));
1874    }
1875    else
1876    {
1877      return(list(ES_Id,ESfix_Id,ES_Id_all_triv));
1878    }
1879  }
1880}
1881example
1882{
1883  "EXAMPLE:"; echo=2;
1884  ring r=0,(x,y),ds;
1885  poly f=x7+y7+(x-y)^2*x2y2;
1886  list K=esIdeal(f);
1887  option(redSB);
1888  // Wahl's equisingularity ideal:
1889  std(K[1]);
1890
1891  ring rr=0,(x,y),ds;
1892  poly f=x4+4x3y+6x2y2+4xy3+y4+2x2y15+4xy16+2y17+xy23+y24+y30+y31;
1893  list K=esIdeal(f);
1894  vdim(std(K[1]));
1895  // the latter should be equal to:
1896  tau_es(f);
1897}
1898
1899///////////////////////////////////////////////////////////////////////////////
1900
1901proc esStratum (poly p_F, list #)
1902"USAGE:   esStratum(F[,m,L]); F poly, m int, L list
1903ASSUME:  F defines a deformation of a reduced bivariate polynomial f
1904         and the characteristic of the basering does not divide mult(f). @*
1905         If nv is the number of variables of the basering, then the first
1906         nv-2 variables are the deformation parameters. @*
1907         If the basering is a qring, ideal(basering) must only depend
1908         on the deformation parameters.
1909COMPUTE: equations for the stratum of equisingular deformations with
1910         fixed (trivial) section.
1911RETURN:  list l: either consisting of a list and an integer, where
1912@format
1913  l[1][1]=ideal defining the equisingularity stratum
1914  l[1][2]=ideal defining the part of the equisingularity stratum where all
1915          equimultiple sections through the non-nodes of the reduced total
1916          transform are trivial sections
1917  l[2]=1 if some error has occured,  l[2]=0 otherwise;
1918@end format
1919or consisting of a ring and an integer, where
1920@format
1921  l[1]=ESSring is a ring extension of basering containing the ideal ES
1922        (describing the ES-stratum), the ideal ES_all_triv (describing the
1923        part with trival equimultiple sections) and the polynomial p_F=F,
1924  l[2]=1 if some error has occured,  l[2]=0 otherwise.
1925@end format
1926NOTE:    L is supposed to be the output of hnexpansion (with the given ordering
1927         of the variables appearing in f). @*
1928         If m is given, the ES Stratum over A/maxideal(m) is computed. @*
1929         This procedure uses @code{execute} or calls a procedure using
1930         @code{execute}.
1931         printlevel>=2 displays additional information.
1932SEE ALSO: esIdeal, isEquising
1933KEYWORDS: equisingularity stratum
1934EXAMPLE: example esStratum; shows examples.
1935"
1936{
1937  list l=esComputation (0,p_F,#);
1938  return(l);
1939}
1940example
1941{
1942   "EXAMPLE:"; echo=2;
1943   int p=printlevel;
1944   printlevel=1;
1945   ring r = 0,(a,b,c,d,e,f,g,x,y),ds;
1946   poly F = (x2+2xy+y2+x5)+ax+by+cx2+dxy+ey2+fx3+gx4;
1947   list M = esStratum(F);
1948   M[1][1];
1949
1950   printlevel=3;     // displays additional information
1951   esStratum(F,2)  ; // ES-stratum over Q[a,b,c,d,e,f,g] / <a,b,c,d,e,f,g>^2
1952
1953   ideal I = f-fa,e+b;
1954   qring q = std(I);
1955   poly F = imap(r,F);
1956   esStratum(F);
1957   printlevel=p;
1958}
1959
1960///////////////////////////////////////////////////////////////////////////////
1961
1962proc isEquising (poly p_F, list #)
1963"USAGE:   isEquising(F[,m,L]); F poly, m int, L list
1964ASSUME:  F defines a deformation of a reduced bivariate polynomial f
1965         and the characteristic of the basering does not divide mult(f). @*
1966         If nv is the number of variables of the basering, then the first
1967         nv-2 variables are the deformation parameters. @*
1968         If the basering is a qring, ideal(basering) must only depend
1969         on the deformation parameters.
1970COMPUTE: tests if the given family is equisingular along the trivial
1971         section.
1972RETURN:  int: 1 if the family is equisingular, 0 otherwise.
1973NOTE:    L is supposed to be the output of hnexpansion (with the given ordering
1974         of the variables appearing in f). @*
1975         If m is given, the family is considered over A/maxideal(m). @*
1976         This procedure uses @code{execute} or calls a procedure using
1977         @code{execute}.
1978         printlevel>=2 displays additional information.
1979EXAMPLE: example isEquising; shows examples.
1980"
1981{
1982  int check=esComputation (1,p_F,#);
1983  return(check);
1984}
1985example
1986{
1987   "EXAMPLE:"; echo=2;
1988   ring r = 0,(a,b,x,y),ds;
1989   poly F = (x2+2xy+y2+x5)+ay3+bx5;
1990   isEquising(F);
1991   ideal I = ideal(a);
1992   qring q = std(I);
1993   poly F = imap(r,F);
1994   isEquising(F);
1995
1996   ring rr=0,(A,B,C,x,y),ls;
1997   poly f=x7+y7+(x-y)^2*x2y2;
1998   poly F=f+A*y*diff(f,x)+B*x*diff(f,x);
1999   isEquising(F);
2000   isEquising(F,2);    // computation over  Q[a,b] / <a,b>^2
2001}
2002
2003
2004
2005/*  Examples:
2006
2007LIB "equising.lib";
2008ring r = 0,(x,y),ds;
2009poly p1 = y^2+x^3;
2010poly p2 = p1^2+x5y;
2011poly p3 = p2^2+x^10*p1;
2012poly p=p3^2+x^20*p2;
2013p;
2014list L=versal(p);
2015def Px=L[1];
2016setring Px;
2017poly F=Fs[1,1];
2018int t=timer;
2019list M=esStratum(F);
2020timer-t;  //-> 3
2021
2022LIB "equising.lib";
2023option(prot);
2024printlevel=2;
2025ring r=0,(x,y),ds;
2026poly f=(x-yx+y2)^2-(y+x)^31;
2027list L=versal(f);
2028def Px=L[1];
2029setring Px;
2030poly F=Fs[1,1];
2031int t=timer;
2032list M=esStratum(F);
2033timer-t;  //-> 233
2034
2035LIB "equising.lib";
2036printlevel=2;
2037option(prot);
2038timer=1;
2039ring r=0,(x,y),ls;
2040poly f=(x4-y4)^2-x10;
2041list L=versal(f);
2042def Px=L[1];
2043setring Px;
2044poly F=Fs[1,1];
2045int t=timer;
2046list M=esStratum(F,3);
2047timer-t;  //-> 8
2048
2049LIB "equising.lib";
2050printlevel=2;
2051timer=1;
2052ring rr=0,(x,y),ls;
2053poly f=x7+y7+(x-y)^2*x2y2;
2054list K=esIdeal(f);
2055// tau_es
2056vdim(std(K[1])); //-> 22
2057// tau_es_fix
2058vdim(std(K[2])); //-> 24
2059
2060LIB "equising.lib";
2061printlevel=2;
2062timer=1;
2063ring rr=0,(x,y),ls;
2064poly f=x7+y7+(x-y)^2*x2y2+x2y4; // Newton non-deg.
2065list K=esIdeal(f);
2066// tau_es
2067vdim(std(K[1])); //-> 21
2068// tau_es_fix
2069vdim(std(K[2])); //-> 23
2070
2071LIB "equising.lib";
2072ring r=0,(w,v),ds;
2073poly f=w2-v199;
2074list L=hnexpansion(f);
2075list LL=versal(f);
2076def Px=LL[1];
2077setring Px;
2078list L=imap(r,L);
2079poly F=Fs[1,1];
2080list M=esStratum(F,2,L);
2081
2082LIB "equising.lib";
2083printlevel=2;
2084timer=1;
2085ring rr=0,(A,B,C,x,y),ls;
2086poly f=x7+y7+(x-y)^2*x2y2;
2087poly F=f+A*y*diff(f,x)+B*x*diff(f,x)+C*diff(f,y);
2088list M=esStratum(F,6);
2089std(M[1][1]);  // standard basis of equisingularity ideal
2090
2091LIB "equising.lib";
2092printlevel=2;
2093timer=1;
2094ring rr=0,(x,y),ls;
2095poly f=x20+y7+(x-y)^2*x2y2+x2y4; // Newton non-degenerate
2096list K=esIdeal(f);
2097K;
2098
2099ring rr=0,(x,y),ls;
2100poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13;
2101list K=esIdeal(f);
2102list L=versal(f);
2103def Px=L[1];
2104setring Px;
2105poly F=Fs[1,1];
2106list M=esStratum(F,2);
2107
2108LIB "equising.lib";
2109ring R=0,(A,B,C,D,x,y),ds;
2110poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13;
2111poly F=f+Ax9+Bx7y2+Cx9y+Dx8y2;
2112list M=esStratum(F,2);
2113
2114LIB "equising.lib";
2115printlevel=2;
2116ring rr=0,(x,y),ls;
2117poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13;
2118list K=esIdeal(f);
2119vdim(std(K[1]));  //-> 51
2120tau_es(f);        //-> 51
2121
2122printlevel=3;
2123f=f*(y-x2)*(y2-x3)*(x-y5);
2124int t=timer;
2125list L=esIdeal(f);
2126vdim(std(L[1]));  //-> 99
2127timer-t;   //-> 42
2128t=timer;
2129tau_es(f);        //-> 99
2130timer-t;   //-> 23
2131
2132LIB "equising.lib";
2133printlevel=3;
2134ring rr=0,(x,y),ds;
2135poly f=x4+4x3y+6x2y2+4xy3+y4+2x2y15+4xy16+2y17+xy23+y24+y30+y31;
2136list K=esIdeal(f);
2137vdim(std(K[1]));  //-> 68
2138tau_es(f);        //-> 68
2139
2140list L=versal(f);
2141def Px=L[1];
2142setring Px;
2143poly F=Fs[1,1];
2144list M=esStratum(F);
2145timer-t;          //-> 0
2146
2147*/
Note: See TracBrowser for help on using the repository browser.