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

fieker-DuValspielwiese
Last change on this file since e182c8 was e182c8, checked in by Hans Schönemann <hannes@…>, 19 years ago
*lossen/greuel/hannes: new libs git-svn-id: file:///usr/local/Singular/svn/trunk@7847 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 52.2 KB
Line 
1version="$Id: equising.lib,v 1.10 2005-04-19 15:23:39 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.
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(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);
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  // 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
909  list LeadDataFm=determine_coef(Fm);
910  def LeadDataG=coef(G,xy);
911 
912  for (i=1; i<=ncols(LeadDataG); i++)
913  {
914    if (LeadDataG[1,i]==LeadDataFm[2])
915    {
916      poly LeadG = LeadDataG[2,i];  // determine the coefficient of G
917      i=ncols(LeadDataG)+1;
918    }
919  }
920   
921  G=LeadDataFm[1]*G-LeadG*Fm;  // leading terms in y should cancel...
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  dbprint(i_print,"// ");
1092  dbprint(i_print,"// Elimination starts:");
1093  dbprint(i_print,"// -------------------");
1094
1095  poly gg;
1096  int b_left=no_b; 
1097
1098  for (i=1;i<=no_b; i++)
1099  {
1100    attrib(J,"isSB",1); 
1101    gg=reduce(b(i),J);
1102    if (gg==0)
1103    {
1104      b_left = b_left-1;  // another b(i) has to be 0
1105    }                 
1106    J = subst(J, b(i), gg);
1107    attrib(J,"isSB",0);
1108  }
1109  J=simplify(J,10);
1110  if (typ==1) // isEquising
1111  {
1112    if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
1113    if(ideal(nselect(J,1,no_b))<>0)
1114    {
1115      setring old_ring;
1116      option(set,ov);
1117      return(0);
1118    }   
1119  }
1120
1121  ideal J_no_b = nselect(J,1,no_b);
1122  if (size(J) > size(J_no_b))
1123  {
1124    dbprint(i_print,"// std computation started");
1125    // some b(i) didn't appear in linear conditions and have to be eliminated
1126    if (defined(artin_bd))
1127    {
1128      // first we make the ring smaller (removing variables, which are
1129      // forced to 0 by J
1130      list LL=make_ring_small(J);
1131      ideal Shortmap=LL[2];
1132      minPolyStr = "";
1133      if (minpoly !=0)
1134      {
1135        minPolyStr = "minpoly = "+string(minpoly);
1136      }
1137      ordStr = "dp(" + string(b_left) + "),dp";
1138      ideal qId = ideal(basering);
1139   
1140      helpStr = "ring Shortring = ("
1141                + charstr(basering) + "),("+ LL[1] +") , ("+ ordStr  +");";
1142      execute(helpStr);
1143      execute(minPolyStr);
1144      // ring has changed to "Shortring"
1145
1146      ideal MM=maxideal(artin_bd);
1147      MM=subst(MM,x,0);
1148      MM=subst(MM,y,0);
1149      MM=simplify(MM,2);
1150      dbprint(i_print-1,"// maxideal("+string(artin_bd)+") has "
1151                         +string(size(MM))+" elements");
1152      dbprint(i_print-1,"//");
1153
1154      // we change to the qring mod m^artin_bd
1155      // first, we have to check if we were in a qring when starting
1156      ideal qId = imap(myRing, qId); 
1157      if (qId == 0)
1158      {
1159         attrib(MM,"isSB",1);
1160         qring QQ=MM;     
1161      }
1162      else
1163      {   
1164         qId=qId,MM;
1165         qring QQ = std(qId);
1166      }
1167
1168      ideal Shortmap=imap(myRing,Shortmap);
1169      map phiphi=myRing,Shortmap;
1170
1171      ideal J=phiphi(J);
1172      option(redSB);
1173      J=std(J);
1174      J=nselect(J,1,no_b);
1175
1176      setring myRing;
1177      // back to "myRing"
1178
1179      J=nselect(J,1,no_b);
1180      Jnew=imap(QQ,J);
1181
1182      J=J,Jnew;
1183      J=interred(J);
1184      if (defined(artin_bd)){ J=jet(J,artin_bd-1); }
1185    }
1186    else
1187    {
1188      J=std(J);
1189      J=nselect(J,1,no_b);
1190      if (defined(artin_bd)){ J=jet(J,artin_bd-1); }
1191    }
1192  }
1193 
1194  dbprint(i_print,"// finished");
1195  dbprint(i_print,"// ");
1196 
1197  minPolyStr = "";option(set,ov);
1198  if (minpoly !=0)
1199  {
1200   minPolyStr = "minpoly = "+string(minpoly);
1201  }
1202
1203  kill HNEring;
1204 
1205  if (typ==1) // isEquising
1206  {
1207    if (defined(artin_bd)) { J=jet(Jnew,artin_bd-1); }
1208    if(J<>0)
1209    {
1210      setring old_ring;
1211      option(set,ov);
1212      return(0);
1213    }   
1214    else
1215    {
1216      setring old_ring;
1217      option(set,ov);
1218      return(1);
1219    }   
1220  }
1221 
1222  setring old_ring;
1223  // we are back in the original ring
1224
1225  if (npars(myRing)<>0)
1226  {
1227    ideal qIdeal = ideal(basering);
1228    helpStr = "ring ESSring = ("
1229                 + string(char(basering))+ "," + parstr(myRing) +
1230                 ") , ("+ varstr(basering)+") , ("+ ordstr(basering) +");";
1231    execute(helpStr);
1232    execute(minPolyStr); 
1233    // basering has changed to ESSring
1234
1235    ideal qIdeal = fetch(old_ring, qIdeal);
1236    if(qIdeal != 0)
1237    {
1238      def r_base = basering;
1239      kill ESSring;
1240      qring ESSring = std(qIdeal);
1241    }
1242    kill qIdeal;
1243 
1244    ideal SSS;
1245    for (int ii=1;ii<=nvars(basering);ii++)
1246    {
1247      SSS[ii+no_b]=var(ii);
1248    }
1249    map phi=myRing,SSS;   // b(i) variables are mapped to zero
1250
1251    ideal ES=phi(J);
1252    ideal ES_all_triv=phi(Jtriv);
1253    kill phi;
1254 
1255    if (defined(p_F)<=0)
1256    {
1257      poly p_F=fetch(old_ring,p_F);
1258      export(p_F);
1259    }
1260    export(ES);
1261    export(ES_all_triv);
1262    setring old_ring;
1263    dbprint(i_print+2,"
1264// 'esStratum' created a list M of a ring and an integer.
1265// To access the ideal defining the equisingularity stratum, type:
1266        def ESSring = M[1]; setring ESSring;  ES; ");
1267
1268    option(set,ov);
1269    return(list(ESSring,0));
1270  }
1271  else
1272  {
1273    // no new ring definition necessary
1274    ideal SSS;
1275    for (int ii=1;ii<=nvars(basering);ii++)
1276    {
1277      SSS[ii+no_b]=var(ii);
1278    }
1279    map phi=myRing,SSS;  // b(i) variables are mapped to zero
1280
1281    ideal ES=phi(J);
1282    ideal ES_all_triv=phi(Jtriv);
1283    kill phi;
1284
1285    setring old_ring;
1286    dbprint(i_print,"// output of 'esStratum' is a list consisting of:
1287//    _[1][1] = ideal defining the equisingularity stratum
1288//    _[1][2] = ideal defining the part of the equisingularity stratum
1289//              where all equimultiple sections are trivial
1290//    _[2] = 0");
1291
1292    option(set,ov);
1293    return(list(list(ES,ES_all_triv),0));
1294  }
1295 
1296}
1297
1298////////////////////////////////////////////////////////////////////////////////
1299
1300proc tau_es (poly f,list #)
1301"USAGE:   tau_es(f); f poly
1302ASSUME:  f is a reduced bivariate polynomial, the basering has precisely
1303         two variables, is local and no qring.
1304RETURN:  int, the codimension of the mu-const stratum in the semi-universal
1305         deformation base.
1306NOTE:    printlevel>=1 displays additional information.
1307         When called with any additional parameter, the computation of the
1308         Milnor number is avoided (no check for NND).
1309SEE ALSO: esIdeal, tjurina, invariants
1310EXAMPLE: example tau_es; shows an example.
1311"
1312{
1313  int i,j,k,s;
1314  int slope_x, slope_y, upper;
1315  int i_print = printlevel - voice + 3;
1316  string MinPolyStr;
1317
1318  // some checks first
1319  if ( nvars(basering)<>2 )
1320  {
1321    print("// basering has not the correct number (two) of variables !");
1322    print("// computation stopped");
1323    return(0);
1324  }
1325  if ( mult(std(1+var(1)+var(2))) <> 0)
1326  {
1327    print("// basering is not local !");
1328    print("// computation stopped");
1329    return(0);
1330  } 
1331
1332  if (mult(std(f))<=1)
1333  {
1334    // f is rigid
1335    return(0);
1336  }
1337
1338  if ( deg(squarefree(f))!=deg(f) )
1339  {
1340    print("// input polynomial was not reduced");
1341    print("// try    squarefree(f);   first");
1342    return(0);
1343  }
1344
1345  def old_ring=basering;
1346  execute("ring @myRing=("+charstr(basering)+"),("+varstr(basering)+"),ds;");
1347  poly f=imap(old_ring,f);
1348
1349  ideal Jacobi_Id = jacob(f);
1350
1351  // check for A_k singularity
1352  // ----------------------------------------
1353  if (mult(std(f))==2)
1354  {
1355    dbprint(i_print-1,"// ");
1356    dbprint(i_print-1,"// polynomial defined A_k singularity");
1357    dbprint(i_print-1,"// ");
1358    return( vdim(std(Jacobi_Id)) );
1359  }
1360
1361  // check for D_k singularity
1362  // ----------------------------------------
1363  if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3)
1364  {
1365    dbprint(i_print,"// ");
1366    dbprint(i_print,"// polynomial defined D_k singularity");
1367    dbprint(i_print,"// ");   
1368    ideal ES_Id = f, jacob(f);
1369    return( vdim(std(Jacobi_Id)));
1370  }
1371
1372
1373  if (size(#)==0)
1374  {
1375    // check if Newton polygon non-degenerate
1376    // ----------------------------------------
1377    Jacobi_Id=std(Jacobi_Id);
1378    int mu = vdim(Jacobi_Id);
1379    poly f_tilde=f+var(1)^mu+var(2)^mu;  //to obtain convenient Newton-polygon
1380   
1381    list NP=newtonpoly(f_tilde);
1382    dbprint(i_print-1,"// Newton polygon:");
1383    dbprint(i_print-1,NP);
1384    dbprint(i_print-1,"");
1385
1386    if(is_NND(f,mu,NP))          // f is Newton non-degenerate
1387    {
1388      upper=NP[1][2];
1389      ideal ES_Id= x^k*y^upper;
1390      dbprint(i_print-1,"polynomial is Newton non-degenerate");
1391      dbprint(i_print-1,"");
1392      k=0;
1393      for (i=1;i<=size(NP)-1;i++)
1394      {
1395        slope_x=NP[i+1][1]-NP[i][1];
1396        slope_y=NP[i][2]-NP[i+1][2];
1397        for (k=NP[i][1]+1; k<=NP[i+1][1]; k++)
1398        {
1399          while ( slope_x*upper + slope_y*k >=
1400                  slope_x*NP[i][2] + slope_y*NP[i][1])
1401          {
1402            upper=upper-1;
1403          }
1404          upper=upper+1;
1405          ES_Id=ES_Id, x^k*y^upper;
1406        }
1407      }
1408      ES_Id=std(ES_Id);
1409      dbprint(i_print-2,"ideal of monomials above Newton bd. is generated by:");
1410      dbprint(i_print-2,ES_Id);
1411      ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
1412      ES_Id = ES_Id, Jacobi_Id;
1413      ES_Id = std(ES_Id);
1414      dbprint(i_print-1,"// ");
1415      dbprint(i_print-1,"// Equisingularity ideal is computed!");
1416      dbprint(i_print-1,"");
1417      return(vdim(ES_Id));
1418    }
1419    else
1420    {
1421      dbprint(i_print-1,"polynomial is Newton degenerate !");
1422      dbprint(i_print-1,"");
1423    }
1424  }
1425
1426  // for Newton degenerate polynomials, we compute the HN expansion, and
1427  // count the number of free points .....
1428 
1429  dbprint(i_print-1,"// ");
1430  dbprint(i_print-1,"// Compute HN expansion");
1431  dbprint(i_print-1,"// ---------------------");
1432  i=printlevel;
1433  printlevel=printlevel-5;
1434  if (2*size(coeffs(f,x))<size(coeffs(f,y)))
1435  {
1436    f=swapXY(f);
1437  }   
1438  list LLL=hnexpansion(f);
1439  if (size(LLL)==0) { // empty list returned by hnexpansion
1440    setring old_ring;
1441    ERROR("Unable to compute HN expansion !");
1442  }   
1443  else
1444  {
1445    if (typeof(LLL[1])=="ring") {
1446      def HNering = LLL[1];
1447      setring HNering;
1448      def @L=hne;
1449    }
1450    else {
1451      def @L=LLL;
1452    }
1453  }
1454  def HNEring=basering;
1455
1456  printlevel=i;
1457  dbprint(i_print-1,"// finished");
1458  dbprint(i_print-1,"// ");
1459
1460  list M=multsequence(@L);
1461  M=control_Matrix(M);     // this returns the 4 control matrices
1462  intmat Mult=M[1];
1463 
1464  list L1=inf_Tangents(@L,nrows(M[1]));
1465  matrix B=L1[1];
1466 
1467  // determine sum_i m_i(m_i+1)/2 (over inf. near points)
1468  int conditions=0;
1469  for (i=1;i<=nrows(Mult);i++)
1470  {
1471    for (j=1;j<=ncols(Mult);j++)
1472    {
1473      conditions=conditions+(Mult[i,j]*(Mult[i,j]+1)/2);
1474    }
1475  }
1476  int freePts=no_freePoints(M[1],B);
1477  int taues=conditions-freePts-2;
1478 
1479  setring old_ring;
1480  return(taues);
1481}
1482example
1483{
1484   "EXAMPLE:"; echo=2;
1485   ring r=32003,(x,y),ds;
1486   poly f=(x4-y4)^2-x10;
1487   tau_es(f);
1488}
1489
1490
1491////////////////////////////////////////////////////////////////////////////////
1492
1493proc esIdeal (poly f,list #)
1494"USAGE:   esIdeal(f[,any]]); f poly
1495ASSUME:  f is a reduced bivariate polynomial, the basering has precisely
1496         two variables, is local and no qring, and the characteristic of
1497         the ground field does not divide mult(f).
1498RETURN:  if called with only one parameter: list of two ideals,
1499@format
1500          _[1]:  equisingularity ideal of f (in sense of Wahl),
1501          _[2]:  ideal of equisingularity with fixed position of the
1502                 singularity;
1503@end format
1504         if called with more than one parameter: list of three ideals,
1505@format
1506          _[1]:  equisingularity ideal of f (in sense of Wahl)
1507          _[2]:  ideal of equisingularity with fixed position of the
1508                 singularity;
1509          _[3]:  ideal of all g such that the deformation defined by f+eg
1510                 (e^2=0) is isomorphic to an equisingular deformation
1511                 of V(f) with all equimultiple sections being trivial.
1512@end format
1513NOTE:    if some of the above condition is not satisfied then return
1514         value is list(0,0).
1515SEE ALSO: tau_es, esStratum
1516KEYWORDS: equisingularity ideal
1517EXAMPLE: example esIdeal; shows examples.
1518"
1519{
1520
1521  int typ;
1522  if (size(#)>0) { typ=1; }   // I^s is also computed
1523  int i,k,s;
1524  int slope_x, slope_y, upper;
1525  int i_print = printlevel - voice + 3;
1526  string MinPolyStr;
1527
1528  // some checks first
1529  if ( nvars(basering)<>2 )
1530  {
1531    print("// basering has not the correct number (two) of variables !");
1532    print("// computation stopped");
1533    return(list(0,0));
1534  }
1535  if ( mult(std(1+var(1)+var(2))) <> 0)
1536  {
1537    print("// basering is not local !");
1538    print("// computation stopped");
1539    return(list(0,0));
1540  } 
1541
1542  if (mult(std(f))<=1)
1543  {
1544    // f is rigid
1545    if (typ==0)
1546    {
1547      return(list(ideal(1),ideal(1)));
1548    }
1549    else
1550    {
1551      return(list(ideal(1),ideal(1),ideal(1)));
1552    }
1553  }
1554
1555  if ( deg(squarefree(f))!=deg(f) )
1556  {
1557    print("// input polynomial was not squarefree");
1558    print("// try    squarefree(f);   first");
1559    return(list(0,0));
1560  }
1561
1562  if (char(basering)<>0)
1563  {
1564    if (mult(std(f)) mod char(basering)==0)
1565    {
1566      print("// characteristic of ground field divides "
1567            + "multiplicity of polynomial !");
1568      print("// computation stopped");
1569      return(list(0,0));   
1570    }
1571  }
1572 
1573  // check for A_k singularity
1574  // ----------------------------------------
1575  if (mult(std(f))==2)
1576  {
1577    dbprint(i_print,"// ");
1578    dbprint(i_print,"// polynomial defined A_k singularity");
1579    dbprint(i_print,"// ");   
1580    ideal ES_Id = f, jacob(f);
1581    ES_Id = interred(ES_Id);
1582    ideal ESfix_Id = f, maxideal(1)*jacob(f);
1583    ESfix_Id= interred(ESfix_Id);
1584    if (typ==0) // only for computation of I^es and I^es_fix
1585    {
1586      return( list(ES_Id,ESfix_Id) );
1587    }
1588    else
1589    {
1590      return( list(ES_Id,ESfix_Id,ES_Id) );
1591    }
1592  }
1593
1594  // check for D_k singularity
1595  // ----------------------------------------
1596  if (mult(std(f))==3 and size(factorize(jet(f,3))[1])>=3)
1597  {
1598    dbprint(i_print,"// ");
1599    dbprint(i_print,"// polynomial defined D_k singularity");
1600    dbprint(i_print,"// ");   
1601    ideal ES_Id = f, jacob(f);
1602    ES_Id = interred(ES_Id);
1603    ideal ESfix_Id = f, maxideal(1)*jacob(f);
1604    ESfix_Id= interred(ESfix_Id);
1605    if (typ==0) // only for computation of I^es and I^es_fix
1606    {
1607      return( list(ES_Id,ESfix_Id) );
1608    }
1609    else
1610    {
1611      return( list(ES_Id,ESfix_Id,ES_Id) );
1612    }
1613  }
1614
1615  // check if Newton polygon non-degenerate
1616  // ----------------------------------------
1617  int mu = milnor(f);
1618  poly f_tilde=f+var(1)^mu+var(2)^mu;  //to obtain a convenient Newton-polygon
1619 
1620  list NP=newtonpoly(f_tilde);
1621  dbprint(i_print-1,"// Newton polygon:");
1622  dbprint(i_print-1,NP);
1623  dbprint(i_print-1,"");
1624
1625  if(is_NND(f,mu,NP))          // f is Newton non-degenerate
1626  {
1627    upper=NP[1][2];
1628    ideal ES_Id= x^k*y^upper;
1629    dbprint(i_print,"polynomial is Newton non-degenerate");
1630    dbprint(i_print,"");
1631    k=0;
1632    for (i=1;i<=size(NP)-1;i++)
1633    {
1634      slope_x=NP[i+1][1]-NP[i][1];
1635      slope_y=NP[i][2]-NP[i+1][2];
1636      for (k=NP[i][1]+1; k<=NP[i+1][1]; k++)
1637      {
1638        while ( slope_x*upper + slope_y*k >=
1639                slope_x*NP[i][2] + slope_y*NP[i][1])
1640        {
1641          upper=upper-1;
1642        }
1643        upper=upper+1;
1644        ES_Id=ES_Id, x^k*y^upper;
1645      }
1646    }
1647    ES_Id=std(ES_Id);
1648    dbprint(i_print-1,"ideal of monomials above Newton bd. is generated by:");
1649    dbprint(i_print-1,ES_Id);
1650    ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
1651    ES_Id = ES_Id, f, jacob(f);
1652    dbprint(i_print,"// ");
1653    dbprint(i_print,"// equisingularity ideal is computed!");
1654    if (typ==0)
1655    {
1656      return(list(ES_Id,ESfix_Id));
1657    }
1658    else
1659    {
1660      return(list(ES_Id,ESfix_Id,ES_Id));
1661    }
1662  }
1663  else
1664  {
1665    dbprint(i_print,"polynomial is Newton degenerate !");
1666    dbprint(i_print,"");
1667  }
1668 
1669  def old_ring=basering;
1670
1671  dbprint(i_print,"// ");
1672  dbprint(i_print,"// versal deformation with triv. section");
1673  dbprint(i_print,"// =====================================");
1674  dbprint(i_print,"// ");
1675 
1676  ideal JJ=maxideal(1)*jacob(f);
1677  ideal kbase_versal=kbase(std(JJ));
1678  s=size(kbase_versal);
1679  string ring_versal="ring @Px = ("+charstr(basering)+"),(t(1.."+string(s)+"),"
1680                        +varstr(basering)+"),(ds("+string(s)+"),"
1681                        +ordstr(basering)+");";
1682  MinPolyStr = string(minpoly);
1683                           
1684  execute(ring_versal);
1685  if (MinPolyStr<>"0")
1686  {
1687    MinPolyStr = "minpoly="+MinPolyStr;
1688    execute(MinPolyStr);
1689  }
1690  // basering has changed to @Px
1691
1692  poly F=imap(old_ring,f);
1693  ideal kbase_versal=imap(old_ring,kbase_versal);
1694  for (i=1; i<=s; i++)
1695  {
1696    F=F+var(i)*kbase_versal[i];
1697  }
1698  dbprint(i_print-1,F);
1699  dbprint(i_print-1,"");
1700
1701
1702  ideal ES_Id,ES_Id_all_triv;
1703  poly Ftriv=F;
1704 
1705  dbprint(i_print,"// ");
1706  dbprint(i_print,"// Compute equisingularity Stratum over Spec(C[t]/t^2)");
1707  dbprint(i_print,"// ===================================================");
1708  dbprint(i_print,"// ");
1709  list M=esStratum(F,2);
1710  dbprint(i_print,"// finished");
1711  dbprint(i_print,"// ");
1712
1713  if (M[2]==1) // error occured during esStratum computation
1714  {
1715    print("Some error has occured during the computation");
1716    return(list(0,0));
1717  }
1718 
1719  if ( typeof(M[1])=="list" )
1720  {
1721    int defpars = nvars(basering)-2;
1722    poly Fred,Ftrivred;
1723    poly g;
1724    F=reduce(F,std(M[1][1]));
1725    Ftriv=reduce(Ftriv,std(M[1][2]));
1726 
1727    for (i=1; i<=defpars; i++)
1728    {
1729      Fred=reduce(F,std(var(i)));
1730      Ftrivred=reduce(Ftriv,std(var(i)));
1731
1732      g=subst(F-Fred,var(i),1);
1733      ES_Id=ES_Id, g;
1734      F=Fred;
1735
1736      g=subst(Ftriv-Ftrivred,var(i),1);
1737      ES_Id_all_triv=ES_Id_all_triv, g;
1738      Ftriv=Ftrivred;
1739    }
1740
1741    setring old_ring;
1742    // back to original ring
1743
1744    ideal ES_Id = imap(@Px,ES_Id);   
1745    ES_Id = interred(ES_Id);
1746
1747    ideal ES_Id_all_triv = imap(@Px,ES_Id_all_triv);   
1748    ES_Id_all_triv = interred(ES_Id_all_triv);
1749
1750    ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
1751    ES_Id = ES_Id, f, jacob(f);
1752    ES_Id_all_triv = ES_Id_all_triv, f, jacob(f);
1753
1754    if (typ==0)
1755    {
1756      return(list(ES_Id,ESfix_Id));
1757    }
1758    else
1759    {
1760      return(list(ES_Id,ESfix_Id,ES_Id_all_triv));
1761    }
1762  }
1763  else
1764  {
1765    def AuxRing=M[1];
1766 
1767    dbprint(i_print,"// ");
1768    dbprint(i_print,"// change ring to ESSring");
1769 
1770    setring AuxRing;  // contains p_F, ES
1771
1772    int defpars = nvars(basering)-2;
1773    poly Fred,Fredtriv;
1774    poly g;
1775    ideal ES_Id,ES_Id_all_triv;
1776
1777    poly p_Ftriv=p_F
1778
1779    p_F=reduce(p_F,std(ES));
1780    p_Ftriv=reduce(p_Ftriv,std(ES_all_triv));
1781    for (i=1; i<=defpars; i++)
1782    {
1783      Fred=reduce(p_F,std(var(i)));
1784      Fredtriv=reduce(p_Ftriv,std(var(i)));
1785
1786      g=subst(p_F-Fred,var(i),1);
1787      ES_Id=ES_Id, g;
1788      p_F=Fred;
1789
1790      g=subst(p_Ftriv-Fredtriv,var(i),1);
1791      ES_Id_all_triv=ES_Id_all_triv, g;
1792      p_Ftriv=Fredtriv;
1793
1794    }
1795
1796    dbprint(i_print,"// ");
1797    dbprint(i_print,"// back to the original ring");
1798
1799    setring old_ring;
1800    // back to original ring
1801
1802    ideal ES_Id = imap(AuxRing,ES_Id);
1803    ES_Id = interred(ES_Id);
1804
1805    ideal ES_Id_all_triv = imap(AuxRing,ES_Id_all_triv);
1806    ES_Id_all_triv = interred(ES_Id_all_triv);
1807
1808    kill @Px;
1809    kill AuxRing;
1810
1811    ideal ESfix_Id = ES_Id, f, maxideal(1)*jacob(f);
1812    ES_Id = ES_Id, f, jacob(f);
1813    ES_Id_all_triv = ES_Id_all_triv, f, jacob(f);
1814    dbprint(i_print,"// ");
1815    dbprint(i_print,"// equisingularity ideal is computed!");
1816    if (typ==0)
1817    {
1818      return(list(ES_Id,ESfix_Id));
1819    }
1820    else
1821    {
1822      return(list(ES_Id,ESfix_Id,ES_Id_all_triv));
1823    }
1824  }
1825}
1826example
1827{
1828  "EXAMPLE:"; echo=2;
1829  ring r=0,(x,y),ds;
1830  poly f=x7+y7+(x-y)^2*x2y2;
1831  list K=esIdeal(f);
1832  option(redSB);
1833  // Wahl's equisingularity ideal:
1834  std(K[1]);
1835
1836  ring rr=0,(x,y),ds;
1837  poly f=x4+4x3y+6x2y2+4xy3+y4+2x2y15+4xy16+2y17+xy23+y24+y30+y31;
1838  list K=esIdeal(f);
1839  vdim(std(K[1]));
1840  // the latter should be equal to:
1841  tau_es(f);     
1842}
1843
1844///////////////////////////////////////////////////////////////////////////////
1845
1846proc esStratum (poly p_F, list #)
1847"USAGE:   esStratum(F[,m,L]); F poly, m int, L list
1848ASSUME:  F defines a deformation of a reduced bivariate polynomial f
1849         and the characteristic of the basering does not divide mult(f). @*
1850         If nv is the number of variables of the basering, then the first
1851         nv-2 variables are the deformation parameters. @*
1852         If the basering is a qring, ideal(basering) must only depend
1853         on the deformation parameters.
1854COMPUTE: equations for the stratum of equisingular deformations with 
1855         fixed (trivial) section.
1856RETURN:  list l: either consisting of a list and an integer, where
1857@format
1858  l[1][1]=ideal defining the equisingularity stratum
1859  l[1][2]=ideal defining the part of the equisingularity stratum where all
1860          equimultiple sections through the non-nodes of the reduced total
1861          transform are trivial sections
1862  l[2]=1 if some error has occured,  l[2]=0 otherwise;
1863@end format
1864or consisting of a ring and an integer, where
1865@format
1866  l[1]=ESSring is a ring extension of basering containing the ideal ES
1867        (describing the ES-stratum), the ideal ES_all_triv (describing the
1868        part with trival equimultiple sections) and the poly p_F=F,
1869  l[2]=1 if some error has occured,  l[2]=0 otherwise.
1870@end format
1871NOTE:    L is supposed to be the output of hnexpansion (with the given ordering
1872         of the variables appearing in f). @*
1873         If m is given, the ES Stratum over A/maxideal(m) is computed. @*
1874         This procedure uses @code{execute} or calls a procedure using
1875         @code{execute}.
1876         printlevel>=2 displays additional information.
1877SEE ALSO: esIdeal, isEquising
1878KEYWORDS: equisingularity stratum
1879EXAMPLE: example esStratum; shows examples.
1880"
1881{
1882  list l=esComputation (0,p_F,#);
1883  return(l);
1884}
1885example
1886{
1887   "EXAMPLE:"; echo=2;
1888   int p=printlevel;
1889   printlevel=1;
1890   ring r = 0,(a,b,c,d,e,f,g,x,y),ds;
1891   poly F = (x2+2xy+y2+x5)+ax+by+cx2+dxy+ey2+fx3+gx4;
1892   list M = esStratum(F);
1893   M[1][1];
1894   
1895   printlevel=3;     // displays additional information
1896   esStratum(F,2)  ; // ES-stratum over Q[a,b,c,d,e,f,g] / <a,b,c,d,e,f,g>^2
1897
1898   ideal I = f-fa,e+b;
1899   qring q = std(I);
1900   poly F = imap(r,F);
1901   esStratum(F);
1902   printlevel=p;
1903}
1904
1905///////////////////////////////////////////////////////////////////////////////
1906
1907proc isEquising (poly p_F, list #)
1908"USAGE:   isEquising(F[,m,L]); F poly, m int, L list
1909ASSUME:  F defines a deformation of a reduced bivariate polynomial f
1910         and the characteristic of the basering does not divide mult(f). @*
1911         If nv is the number of variables of the basering, then the first
1912         nv-2 variables are the deformation parameters. @*
1913         If the basering is a qring, ideal(basering) must only depend
1914         on the deformation parameters.
1915COMPUTE: tests if the given family is equisingular along the trivial
1916         section.
1917RETURN:  int: 1 if the family is equisingular, 0 otherwise.
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 family is considered over A/maxideal(m). @*
1921         This procedure uses @code{execute} or calls a procedure using
1922         @code{execute}.
1923         printlevel>=2 displays additional information.
1924EXAMPLE: example isEquising; shows examples.
1925"
1926{
1927  int check=esComputation (1,p_F,#);
1928  return(check);
1929}
1930example
1931{
1932   "EXAMPLE:"; echo=2;
1933   ring r = 0,(a,b,x,y),ds;
1934   poly F = (x2+2xy+y2+x5)+ay3+bx5;
1935   isEquising(F);
1936   ideal I = ideal(a);
1937   qring q = std(I);
1938   poly F = imap(r,F);
1939   isEquising(F);
1940
1941   ring rr=0,(A,B,C,x,y),ls;
1942   poly f=x7+y7+(x-y)^2*x2y2;
1943   poly F=f+A*y*diff(f,x)+B*x*diff(f,x);
1944   isEquising(F); 
1945   isEquising(F,2);    // computation over  Q[a,b] / <a,b>^2
1946}
1947
1948
1949
1950/*  Examples:
1951
1952LIB "equising.lib";
1953   ring r = 0,(x,y),ds;
1954   poly p1 = y^2+x^3;
1955   poly p2 = p1^2+x5y;
1956   poly p3 = p2^2+x^10*p1;
1957   poly p=p3^2+x^20*p2;
1958   p;
1959   versal(p);
1960   setring Px;
1961   poly F=Fs[1,1];
1962   int t=timer;
1963   list M=esStratum(F);
1964   timer-t;  //-> 3
1965
1966LIB "equising.lib";
1967option(prot);
1968printlevel=2;
1969ring r=0,(x,y),ds;
1970poly f=(x-yx+y2)^2-(y+x)^31;
1971versal(f);
1972setring Px;
1973poly F=Fs[1,1];
1974int t=timer;
1975list M=esStratum(F);
1976timer-t;  //-> 233
1977
1978
1979LIB "equising.lib";
1980printlevel=2;
1981option(prot);
1982timer=1;
1983ring r=32003,(x,y),ds;
1984poly f=(x4-y4)^2-x10;
1985versal(f);
1986setring Px;
1987poly F=Fs[1,1];
1988int t=timer;
1989list M=esStratum(F,3);
1990timer-t;  //-> 8
1991
1992LIB "equising.lib";
1993printlevel=2;
1994timer=1;
1995ring rr=0,(x,y),ls;
1996poly f=x7+y7+(x-y)^2*x2y2;
1997list K=esIdeal(f);
1998// tau_es
1999vdim(std(K[1])); //-> 22
2000// tau_es_fix
2001vdim(std(K[2])); //-> 24
2002
2003
2004LIB "equising.lib";
2005printlevel=2;
2006timer=1;
2007ring rr=0,(x,y),ls;
2008poly f=x7+y7+(x-y)^2*x2y2+x2y4; // Newton non-deg.
2009list K=esIdeal(f);
2010// tau_es
2011vdim(std(K[1])); //-> 21
2012// tau_es_fix
2013vdim(std(K[2])); //-> 23
2014
2015LIB "equising.lib";
2016ring r=0,(w,v),ds;
2017poly f=w2-v199;
2018list L=hnexpansion(f);
2019versal(f);
2020setring Px;
2021list L=imap(r,L);
2022poly F=Fs[1,1];
2023list M=esStratum(F,2,L);
2024
2025LIB "equising.lib";
2026printlevel=2;
2027timer=1;
2028ring rr=0,(A,B,C,x,y),ls;
2029poly f=x7+y7+(x-y)^2*x2y2;
2030poly F=f+A*y*diff(f,x)+B*x*diff(f,x)+C*diff(f,y);
2031list M=esStratum(F,6);
2032std(M[1][1]);  // standard basis of equisingularity ideal
2033
2034LIB "equising.lib";
2035printlevel=2;
2036timer=1;
2037ring rr=0,(x,y),ls;
2038poly f=x20+y7+(x-y)^2*x2y2+x2y4; // Newton non-degenerate
2039list K=esIdeal(f);
2040K;
2041
2042ring rr=0,(x,y),ls;
2043poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13;
2044list K=esIdeal(f);
2045versal(f);
2046setring Px;
2047poly F=Fs[1,1];
2048list M=esStratum(F,2);
2049
2050LIB "equising.lib";
2051ring R=0,(A,B,C,D,x,y),ds;
2052poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13;
2053poly F=f+Ax9+Bx7y2+Cx9y+Dx8y2;
2054list M=esStratum(F,2);
2055
2056
2057LIB "equising.lib";
2058printlevel=2;
2059ring rr=0,(x,y),ls;
2060poly f=x6y-3x4y4-x4y5+3x2y7-x4y6+2x2y8-y10+2x2y9-y11+x2y10-y12-y13;
2061list K=esIdeal(f);
2062vdim(std(K[1]));  //-> 51
2063tau_es(f);        //-> 51
2064
2065printlevel=3;
2066f=f*(y-x2)*(y2-x3)*(x-y5);
2067int t=timer;
2068list L=esIdeal(f);
2069vdim(std(L[1]));  //-> 99
2070timer-t;   //-> 42
2071t=timer;
2072tau_es(f);        //-> 99
2073timer-t;   //-> 23
2074
2075
2076LIB "equising.lib";
2077printlevel=3;
2078ring rr=0,(x,y),ds;
2079poly f=x4+4x3y+6x2y2+4xy3+y4+2x2y15+4xy16+2y17+xy23+y24+y30+y31;
2080list K=esIdeal(f);
2081vdim(std(K[1]));  //-> 68
2082tau_es(f);        //-> 68
2083
2084versal(f);
2085setring Px;
2086poly F=Fs[1,1];
2087list M=esStratum(F);
2088timer-t;          //-> 0
2089
2090
2091
2092*/
Note: See TracBrowser for help on using the repository browser.