# source:git/Singular/LIB/equising.lib@120c517

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