# source:git/Singular/LIB/equising.lib@50cbdc

jengelh-datetimespielwiese
Last change on this file since 50cbdc was 50cbdc, checked in by Hans Schönemann <hannes@…>, 22 years ago
• Property mode set to 100644
File size: 28.1 KB
Line
1version="\$Id: equising.lib,v 1.8 2001-08-27 14:47:49 Singular Exp \$";
2category="Singularities";
3info="
4LIBRARY:  equising.lib  Equisingularity Stratum of a Family of Plane Curves
5AUTHOR:   Andrea Mindnich, mindnich@mathematik.uni-kl.de
6
7PROCEDURES:
8 esStratum(F[,m]);   computes the equisingularity stratum of a deformation
9 isEquising(F[,m]);  tests if a given deformation is equisingular
10";
11
12LIB "poly.lib";
13LIB "elim.lib";
14LIB "hnoether.lib";
15///////////////////////////////////////////////////////////////////////////////
16// COMPUTES  a weight vector. x and y get weight 1 and all other
17//           variables get weight 0.
18static proc xyVector()
19{
20  intvec iv ;
21  iv[nvars(basering)]=0 ;
22  iv[rvar(x)] =1;
23  iv[rvar(y)] =1;
24  return (iv);
25}
26///////////////////////////////////////////////////////////////////////////////
27// exchanges the variables x and y in the polynomial p_f
28static proc swapXY(poly f)
29{
30  def r_base = basering;
31  ideal MI = maxideal(1);
32  MI[rvar(x)]=y;
33  MI[rvar(y)]=x;
34  map phi = r_base, MI;
35  f=phi(f);
36  return (f);
37}
38///////////////////////////////////////////////////////////////////////////////
39// ASSUME:   p_mon is a monomial
40//            and p_g is the product of the variables in p_mon.
41// COMPUTES the coefficient of p_mon in p_h.
42static proc coefficient(poly p_h, poly p_mon, poly p_g)
43{
44  matrix coefMatrix = coef(p_h,p_g);
45  int nc = ncols(coefMatrix);
46  int ii=1;
47  poly p_c=0;
48  while(ii<=nc)
49  {
50    if (coefMatrix[1,ii] == p_mon)
51    {
52      p_c = coefMatrix[2,ii];
53      break;
54    }
55    ii++;
56  }
57  return (p_c);
58}
59///////////////////////////////////////////////////////////////////////////////
60// in p_F the variable p_vari is substituted by the polynomial p_g.
61static proc eSubst(poly p_F, poly p_vari, poly p_g)
62{
63  def r_base = basering;
64  ideal MI;
65  map phi;
66  int ii = rvar(p_vari);
67  if (ii != 0)
68  {
69    MI = maxideal(1);
70    MI[ii] = p_g;
71    phi = r_base, MI;
72    p_F = phi(p_F);
73  }
74  return (p_F);
75}
76///////////////////////////////////////////////////////////////////////////////
77// All ring variables of p_F which occur in (the set of generators of) the
78// ideal  Id, are substituted by 0
79static proc SimplifyF(poly p_F, ideal Id)
80{
81  int i=1;
82  int si = size(Id);
83  for (i=1; i <= si; i++)
84  {
85    if (rvar(Id[i]))
86    {
87      p_F = subst(p_F, Id[i], 0);
88    }
89  }
90  return(p_F);
91}
92
93
94///////////////////////////////////////////////////////////////////////////////
95
96
97// Checks, if the basering is admissible.
98static proc checkBasering ()
99{
100  int error = 0;
101  if(find(charstr(basering),"real"))
102  {
103    ERROR ("cannot compute esStratum with 'real' as coefficient field");
104  }
105  if (nvars(basering) <= 2)
106  {
107   ERROR ("there are to less ringvariables to compute esStratum")
108  }
109  error = checkQIdeal(ideal(basering));
110  return(error);
111}
112///////////////////////////////////////////////////////////////////////////////
113static  proc getInput (list #)
114{
115  def r_base = basering;
116
117  int maxStep = -1;
118
119  if (size(#) >= 1)
120  {
121    if (typeof(#[1]) == "int")
122    {
123      maxStep = #[1];
124    }
125
126    else
127    {
128      ERROR("expected esStratum('poly', 'int') ");
129    }
130  }
131
132  return(maxStep);
133}
134//////////////////////////////////////////////////////////////////////////////
135// RETURNS: 0, if the ideal cond only depends on the deformation parameters
136//          1, otherwise.
137static proc checkQIdeal (ideal cond)
138{
139  def r_base = basering;
140
141  int i_print = printlevel-voice + 4;
142  int i_nvars = nvars(basering);
143
144  ideal id_help = subst(cond,var(i_nvars),0,var(i_nvars-1),0) - cond;
145
146  // cond depends only on the first i_nvars-2 variables <=>
147  //  id_help == <0>
148  if ( simplify(id_help, 2) != 0)
149  {
150    dbprint(i_print,
151    "ideal(basering) must only depend on the deformation parameters");
152    return(1);
153  }
154
155  return(0);
156}
157
158
159///////////////////////////////////////////////////////////////////////////////
160// COMPUTES string(minpoly) and substitutes the parameter by newParName
161static proc makeMinPolyString (string newParName)
162{
163  int i;
164  string parName = parstr(basering);
165  int parNameSize = size(parName);
166
167  string oldMinPolyStr = string (minpoly);
168  int minPolySize = size(oldMinPolyStr);
169
170  string newMinPolyStr = "";
171
172  for (i=1;i <= minPolySize; i++)
173  {
174    if (oldMinPolyStr[i,parNameSize] == parName)
175    {
176      newMinPolyStr = newMinPolyStr + newParName;
177      i = i + parNameSize-1;
178    }
179    else
180    {
181      newMinPolyStr = newMinPolyStr + oldMinPolyStr[i];
182    }
183  }
184
185  return(newMinPolyStr);
186}
187
188
189///////////////////////////////////////////////////////////////////////////////
190// Defines a new ring without deformation-parameters.
191static proc createHNERing()
192{
193  string str;
194  string minpolyStr = string(minpoly);
195
196  str = " ring HNERing = (" + charstr(basering) + "), (x,y), ls;";
197  execute (str);
198
199  str = "minpoly ="+ minpolyStr+";";
200  execute(str);
201
202  keepring(HNERing);
203}
204///////////////////////////////////////////////////////////////////////////////
205// RETURNS: 1, if p_f = 0 or char(basering) divides the order of p_f
206//             or p_f is not squarefree.
207//          0, otherwise
208static proc checkPoly (poly p_f)
209{
210  int i_print = printlevel - voice + 3;
211  int i_ord;
212
213  if (p_f == 0)
214  {
215    dbprint(i_print,"The Input is a 'deformation'  of the zero polynomial");
216    return(1);
217  }
218
219  i_ord = mindeg1(p_f);
220
221  if (number(i_ord) == 0)
222  {
223    dbprint(i_print,
224            "The characteristic of the coefficient field
225             divides the order of the equation");
226    return(1);
227  }
228
229  if (squarefree(p_f) != p_f)
230  {
231    dbprint(i_print, "The curve is reducible");
232    return(1);
233  }
234
235  return(0);
236}
237///////////////////////////////////////////////////////////////////////////////
238// COMPUTES  the multiplicity sequence of p_f
239static proc calcMultSequence (poly p_f)
240{
241  int i_print = printlevel-voice + 3;
242  intvec multSeq=0;
243  list hneList;
244  int xNotTransversal;
245  int fIrreducible = 1;
246
247  // if basering =  (p,a) or (p,a(1..s)),
248  //            p prime, a algebraic, a(1..s) transcendent use reddevelop
249  // otherwise use develop
250  if (char(basering) != 0
251      && npars(basering) !=0
252      && charstr(basering) == string(char(basering)) + "," + parstr(basering))
253  {
254    hneList = reddevelop(p_f, -1);
255    if ( size(hneList)>=2)
256    {
257      fIrreducible = 0;
258      dbprint(i_print, "The curve is reducible");
259      return(multSeq, xNotTransversal, fIrreducible);
260    }
261
262    hneList = hneList[1];
263
264    xNotTransversal= hneList[3];
265  }
266  else
267  {
268    hneList = develop(p_f, -1);
269
270    xNotTransversal= hneList[3];
271    fIrreducible = hneList[5];
272  }
273
274  if ( ! fIrreducible)
275  {
276    dbprint(i_print, "The curve is reducible");
277    return(multSeq, xNotTransversal, fIrreducible);
278  }
279
280  multSeq = multsequence (hneList);
281  return(multSeq, xNotTransversal, fIrreducible);
282}
283///////////////////////////////////////////////////////////////////////////////
284// ASSUME:  The basering is no qring and has at least 3 variables
285// DEFINES: A new basering, "myRing",
286//          with new names for the parameters and variables.
287//          The new names for the parameters are a(1..k),
288//          and t(1..s),x,y for the variables
289//          The ring ordering is ordStr.
290// NOTE:    This proc uses 'execute'.
291static proc createMyRing(poly p_F, string ordStr )
292{
293  def r_old = basering;
294
295  int chara = char(basering);
296  string charaStr;
297  int i;
298  string minPolyStr = "";
299  string helpStr;
300  int nDefParams = nvars(r_old)-2;
301
302
303  ideal qIdeal = ideal(basering);
304
305  if (npars(basering) == 0)
306  {
307    helpStr = "ring myRing ="
308              + string(chara)+ ", (t(1..nDefParams), x, y),"+ ordStr +";";
309    execute(helpStr);
310  }
311
312  else
313  {
314    charaStr = charstr(basering);
315    if (charaStr == string(chara) + "," + parstr(basering))
316    {
317      if (minpoly !=0)
318      {
319        minPolyStr = makeMinPolyString("a");
320
321        helpStr = "ring myRing =
322                 (" + string(chara) + ",a),
323                 (t(1..nDefParams), x, y)," + ordStr + ";";
324        execute(helpStr);
325
326        helpStr = "minpoly =" + minPolyStr + ";";
327        execute (helpStr);
328      }
329      else
330      {
331        helpStr = "ring myRing =
332                  (" + string(chara) + ",a(1..npars(basering)) ),
333                  (t(1..nDefParams), x, y)," + ordStr + ";";
334        execute(helpStr);
335      }
336    }
337    else
338    {
339      i = find (charaStr,",");
340
341      helpStr = " ring myRing =
342              (" + charaStr[1,i] + "a),
343              (t(1..nDefParams), x, y)," + ordStr + ";";
344
345      execute (helpStr);
346    }
347  }
348
349  ideal qIdeal = fetch(r_old, qIdeal);
350
351  if(qIdeal != 0)
352  {
353    def r_base = basering;
354    kill myRing;
355    qring myRing = std(qIdeal);
356  }
357
358  poly p_F = fetch(r_old, p_F);
359  ideal ES;
360
361  keepring(myRing);
362}
363
364///////////////////////////////////////////////////////////////////////////////
365/////////// procedures to compute the equisingularity stratum /////////////////
366///////////////////////////////////////////////////////////////////////////////
367// DEFINES a new basering, myRing,  which has one variable
368//         more than the old ring.
369//         The name for the new variable is "H(nhelpV)".
370//         p_F and ES are "imaped" into the new ring.
371static proc extendRing (poly p_F,  ideal ES, int nHelpV, ideal HCond)
372{
373  def r_old = basering;
374
375  string helpStr;
376  string minPolyStr = "";
377
378  ideal qIdeal = ideal(basering);
379
380  if (minpoly != 0)
381  {
382    if (charstr(basering) == string(char(basering)) + "," + parstr(basering))
383    {
384      minPolyStr = string(minpoly);
385    }
386  }
387
388  string str =
389 "ring myRing = (" + charstr(r_old) + "),
390                (H(" + string(nHelpV)+ ")," + string(maxideal(1)) + "),
391                (dp(" + string(nHelpV) + "),dp);";
392  execute (str);
393
394  if (minPolyStr != "")
395  {
396    helpStr = "minpoly =" + minPolyStr + ";";
397    execute(helpStr);
398  }
399
400  ideal qIdeal = imap(r_old, qIdeal);
401  if(qIdeal != 0)
402  {
403    def r_base = basering;
404
405    kill myRing;
406
407    attrib(qIdeal,"isSB",1);
408    qring myRing = qIdeal;
409  }
410
411
412  poly p_F = imap(r_old, p_F);
413  ideal ES = imap(r_old, ES);
414  ideal HCond = imap(r_old, HCond);
415
416  keepring(myRing);
417}
418///////////////////////////////////////////////////////////////////////////////
419// COMPUTES an ideal equimultCond, such that F_(n-1) mod equimultCond =0,
420//          where F_(n-1) is the (nNew-1)-jet of p_F with respect to x,y.
421static proc calcEquimultCond(poly p_F, int nNew)
422{
423  ideal equimultCond = 0;
424  poly p_FnMinus1;
425  matrix coefMatrix;
426  int nc;
427  int ii = 1;
428
429  p_FnMinus1 = jet(p_F, nNew-1, xyVector());
430
431  coefMatrix = coef(p_FnMinus1, xy);
432
433  nc = ncols(coefMatrix);
434
435  for (ii=1; ii<=nc; ii++)
436  {
437    equimultCond[ii] = NF(coefMatrix[2,ii],std(0));
438  }
439
440  p_F = p_F - p_FnMinus1;
441  p_F = SimplifyF(p_F, equimultCond);
442
443  return(equimultCond, p_F);
444}
445///////////////////////////////////////////////////////////////////////////////
446// COMPUTES smallest integer >= nNew/nOld -1
447static proc calcNZeroSteps (int nOld,int nNew)
448{
449 int nZeroSteps;
450
451 if (nOld mod nNew == 0)
452 {
453   nZeroSteps = nOld div nNew -1;
454 }
455
456 else
457 {
458  nZeroSteps  = nOld div nNew;
459 }
460
461 return(nZeroSteps);
462}
463///////////////////////////////////////////////////////////////////////////////
464// ASSUME: ord_(X,Y)(F)=nNew
465// COMPUTES an ideal I such that (p_F mod I)_nNew  = p_c*y^nNew.
466static proc purePowerOfY (poly p_F, int nNew)
467{
468  ideal id_help = 0;
469  poly  p_Fn;
470  matrix coefMatrix;
471  int nc;
472  poly p_c;
473  int ii=1;
474
475  p_Fn = jet(p_F, nNew, xyVector());
476
477  coefMatrix = coef(p_Fn, xy);
478  nc = ncols(coefMatrix);
479
480  p_c = coefMatrix[2,nc];
481
482  for (ii=1; ii <= nc-1; ii++)
483  {
484    id_help[ii] = NF(coefMatrix[2,ii],std(0));
485  }
486
487  p_F = p_F - p_Fn + p_c*y^nNew;
488
489  p_F = SimplifyF(p_F, id_help);
490
491  return(id_help, p_F, p_c);
492}
493///////////////////////////////////////////////////////////////////////////////
494// ASSUME: ord_(X,Y)(F)=nNew
495// COMPUTES an ideal I such that p_Fn mod I = p_c*(y-p_a*x)^nNew,
496//          where p_Fn is the homogeneous part of p_F of order nNew.
497static proc purePowerOfLin (poly p_F, ideal HCond, int nNew, int nHelpV)
498{
499  ideal id_help = 0;
500  poly p_Fn;
501  matrix coefMatrix;
502  poly p_c;
503  poly p_ca;
504  poly p_help;
505  poly p_a;
506  int ii;
507  int jj;
508  int bico;
509
510  p_Fn = jet(p_F, nNew, xyVector());
511
512  coefMatrix = coeffs(subst(p_Fn,x,1),y);
513
514  p_c = coefMatrix[nNew+1,1];
515  p_ca = coefMatrix[nNew,1]/(-nNew);
516
517  if (npars(basering)==1 && charstr(basering) != string(char(basering)) + "," + parstr(basering))
518  {
519    p_a = H(nHelpV);
520    HCond = HCond + ideal(p_ca - p_a*p_c);
521  }
522  else
523  {
524    p_help = p_ca/p_c;
525    if (p_help * p_c == p_ca)
526    {
527      p_a = p_help;
528    }
529    else
530    {
531      p_a = H(nHelpV);
532      HCond = HCond + ideal(p_ca - p_a*p_c);
533    }
534  }
535
536  bico = (nNew*(nNew-1))/2;
537
538  for (ii = 2; ii <= nNew ; ii++)
539  {
540
541    if (coefMatrix[nNew+1-ii,1] == 0)
542    {
543      if (number(bico) != 0)
544      // Then a^i=0 since c is a unit
545      {
546        id_help = id_help + ideal(NF(p_a^(ii),std(0)));
547        for (jj = ii+1; jj <= nNew; jj++)
548        // the remaining coefficients (of y^(nnew-k)*x^k with k>i)
549        // are also zero.
550        {
551          id_help = id_help
552                    + ideal(NF(coefMatrix[nNew+1-jj,1],std(0)));
553        }
554        break;
555      }
556    }
557
558    else
559    {
560      id_help = id_help +
561            ideal(NF(coefMatrix[nNew+1-ii,1] - bico*p_c*(-p_a)^ii,std(0)));
562    }
563
564    bico = (bico*(nNew-ii))/(ii+1);
565  }
566
567  p_F = SimplifyF(p_F, id_help);
568
569  return(id_help, HCond, p_F, p_c, p_a);
570}
571///////////////////////////////////////////////////////////////////////////////
572// eliminates the variables H(1),..,H(nHelpV) from the ideal ES + HCond
573static proc helpVarElim(ideal ES, ideal HCond, int nHelpV)
574{
575  ES = ES + HCond;
576  ES = std(ES);
577  ES = nselect(ES,1,nHelpV);
578
579  return(ES);
580}
581///////////////////////////////////////////////////////////////////////////////
582// ASSUME:   ord(F)=nNew and p_c(y-p_a*x)^n is the nNew-jet of F with respect
583//           to X,Y
584// COMPUTES  F(x,yx+a*x)/x^n
585static proc formalBlowUp(poly p_F, poly p_c, poly p_a, int nNew)
586{
587
588  p_F = p_F - jet(p_F, nNew, xyVector());
589
590  if (p_a != 0)
591  {
592    p_F = eSubst(p_F, y , yx + p_a*x);
593  }
594  else
595  {
596    p_F = subst(p_F, y, xy);
597  }
598
599  p_F = p_F/(x^nNew);
600
601  p_F = p_F + p_c * y^nNew;
602
603  return (p_F);
604}
605///////////////////////////////////////////////////////////////////////////////
606// ASSUME:  p_F in K[t(1)..t(s),x,y]
607// COMPUTES the minimal ideal ES, such that the deformation p_F mod ES is
608//          equisingular.
609//          The computation is done up to iteration step 'maxstep'.
610// RETURNS: list l, such that
611//          l[1]=1 if some error has occured,
612//          l[1]=0 otherwise and then l[2] = es_cond.
613static proc calcEsCond(poly p_F, intvec multSeq, int maxStep)
614{
615  def r_old = basering;
616
617  ideal ES = 0;
618
619  int ii;
620  int step = 1;
621  int nNew = multSeq[step];
622  int nOld = nNew;
623  int nZeroSteps;
624  int nHelpV = 1;
625  ideal HCond = 0;
626  int maxDeg = 0;
627  int z = printlevel - voice + 3;
628  string str;
629
630  extendRing(p_F, ES, nHelpV, HCond);
631
632  poly p_c;
633  poly p_a;
634  ideal I;
635
636  for (ii = 1; ii <= size(multSeq); ii++)
637  {
638    maxDeg = maxDeg + multSeq[ii];
639  }
640
641  while (step <= maxStep)
642  {
643
644    nOld = nNew;
645    nNew = multSeq[step];
646
647    p_F = jet(p_F, maxDeg, xyVector());
648    maxDeg = maxDeg - nNew;
649
650    if (nNew<nOld)
651    //start a new line in the HNE of F
652    //               _            _
653    // for the next | nold/nnew -1 | iteration steps the coefficient 'a'
654    // in the leading form Fn = c(y-ax)^n should  be zero.
655    {
656      p_F = swapXY(p_F);
657      nZeroSteps = calcNZeroSteps (nOld, nNew);
658    }
659
660    I, p_F = calcEquimultCond (p_F, nNew);
661    ES = ES + I;
662
663    if (z>1)
664    {
665      "// conditions for equimultiplicity in step:", step;
666      I;
667      if (nHelpV >1)
668      {
669        str = string(nHelpV);
670        "// conditions for help variables H(1),..,H("+str+"):";
671        HCond;
672      }
673      pause("press <return> to continue");
674    }
675
676    if (nZeroSteps > 0)
677    {
678      nZeroSteps--;
679
680      // compute conditions, s.th. Fn = c*y^nnew ?
681      I, p_F, p_c = purePowerOfY (p_F, nNew);
682      ES = ES + I;
683
684      if (z>1)
685      {
686        "// conditions for pure power in step:", step;
687        I;
688        if (nHelpV > 1)
689        {
690          str = string(nHelpV);
691          "// conditions for help variables H(1),..,H("+str+"):";
692          HCond;
693        }
694        pause("press <return> to continue");
695      }
696      p_a =0;
697    }
698
699    else
700    {
701      I, HCond, p_F, p_c, p_a = purePowerOfLin (p_F, HCond, nNew, nHelpV);
702
703      ES = ES + I;
704
705      if (z>1)
706      {
707        "// conditions for pure power in step:", step;
708        I;
709        str = string(nHelpV);
710        "// conditions for help variables H(1),..,H("+str+"):";
711        HCond;
712        pause("press <return> to continue");
713      }
714    }
715
716    p_F = formalBlowUp (p_F, p_c, p_a, nNew);
717
718    if (p_a == H(nHelpV))
719    {
720      nHelpV++;
721
722      def r_base = basering;
723      kill myRing;
724
725      extendRing(p_F, ES, nHelpV, HCond);
726
727      kill r_base;
728
729      poly p_a;
730      poly p_c;
731      ideal I;
732    }
733    step++;
734  }
735  if (nHelpV > 1)
736  {
737    ES = helpVarElim(ES, HCond, nHelpV);
738  }
739
740  if (nameof(basering)=="myRing")
741  {
742    setring r_old;
743    ES = imap(myRing, ES);
744  }
745
746  return(ES);
747}
748
749///////////////////////////////////////////////////////////////////////////////
750//                           main procedure
751///////////////////////////////////////////////////////////////////////////////
752
753proc esStratum (poly p_F, list #)
754"USAGE:   esStratum(F[,m]); F poly, m int
755ASSUME:  F defines a deformation of an irreducible bivariate polynomial f
756         and the characteristic of the basering does not divide mult(f). @*
757         If nv is the number of variables of the basering, then the first nv-2
758         variables are the deformation parameters. @*
759         If the basering is a qring, ideal(basering) must only depend
760         on the deformation parameters.
761RETURN:  list l of an ideal and an integer, where
762@format
763  l[1] is the ideal in the deformation parameters, defining the ES-stratum of F,
764  l[2]=1 if some error has occured,  l[2]=0 otherwise.
765@end format
766NOTE:    If m is given, the computation stops after m steps of the iteration. @*
767         printlevel > 0 displays comments and pauses after intermediate
768         computations (default: printlevel=0) @*
769         This procedure uses @code{execute} or calls a procedure using
770         @code{execute}.
771EXAMPLE: example esStratum; shows an example
772"
773{
774  def r_user = basering;
775
776  int ii = 1;
777  int i_nvars = nvars(basering);
778  int error = 0;
779  int xNotTransversal;
780  int fIrreducible;
781  int maxStep;
782  int userMaxStep;
783  ideal cond;
784  intvec multSeq;
785  ideal ES = 0;
786
787  error = checkBasering();
788  if (error)
789  {
790    return(list(ES,error));
791  }
792
793  userMaxStep = getInput(#);
794
795  // define a new basering "myRing" with new names for parameters
796  // and variables.
797  // The new names are 'a(1)', ..., 'a(npars)' for the parameters
798  // and 't(1)', ..., 't(nvars-2)', 'x', 'y' for the variables.
799  createMyRing(p_F, "dp");
800
801  // define a ring without deformation parameters, to compute the HNE
802  // of F mod <t_1,..,t_s>
803  createHNERing();
804
805  poly p_f = imap(myRing,p_F);
806
807  error = checkPoly(p_f);
808  if (error)
809  {
810    setring r_user;
811    return(list( ideal(0),error));
812  }
813
814  // compute the multiplicitysequence of p_f.
815  multSeq, xNotTransversal, fIrreducible = calcMultSequence(p_f);
816
817  if ( ! fIrreducible)
818  {
819    setring r_user;
820    return(list(ideal(0),1));
821  }
822
823  setring myRing;
824
825  if (xNotTransversal)
826  {
827    p_F = swapXY(p_F);
828  }
829
830  if (userMaxStep != -1 && userMaxStep <  size(multSeq)-1)
831  {
832    maxStep = userMaxStep;
833  }
834  else
835  {
836    maxStep = size(multSeq)-1;
837  }
838
839  ES = calcEsCond(p_F, multSeq, maxStep);
840
841  setring r_user;
842  ES = fetch(myRing, ES);
843
844  return(list(ES, error));
845}
846
847example
848{
849   "EXAMPLE:"; echo=2;
850   ring r = 11,(a,b,c,d,e,f,g,x,y),ds;
851   poly F = (x2+2xy+y2+x5)+ax+by+cx2+dxy+ey2+fx3+gx4;
852   esStratum(F);
853   esStratum(F,2);
854   ideal I = f-fa,e+b;
855   qring q = std(I);
856   poly F = imap(r,F);
857   esStratum(F);
858}
859
860///////////////////////////////////////////////////////////////////////////////
861//                  procedures for equisingularity test
862///////////////////////////////////////////////////////////////////////////////
863
864// DEFINES a new basering, myRing,  which has one variable
865//         more than the old ring.
866//         The name for the new variable is "H(nhelpV)".
867static proc T_extendRing(poly p_F, int nHelpV, ideal HCond)
868{
869  def r_old = basering;
870
871  ideal qIdeal = ideal(basering);
872
873  string helpStr;
874  string minPolyStr = "";
875
876  if(minpoly != 0)
877  {
878    if (charstr(basering) == string(char(basering)) + "," + parstr(basering))
879    {
880      minPolyStr = string(minpoly);
881    }
882  }
883
884  string str = "ring myRing =
885               (" + charstr(r_old) + "),
886               (H(" + string( nHelpV)+ ")," + string(maxideal(1)) + "),
887               (dp(" + string( nHelpV) + "), ds);";
888  execute (str);
889
890  if (minPolyStr != "")
891  {
892    helpStr = "minpoly =" + minPolyStr + ";";
893    execute(helpStr);
894  }
895
896  ideal qIdeal = imap(r_old, qIdeal);
897  if(qIdeal != 0)
898  {
899    def r_base = basering;
900    kill myRing;
901    qring myRing = std(qIdeal);
902  }
903
904  poly p_F =imap(r_old, p_F);
905  ideal HCond = imap(r_old, HCond);
906
907  keepring(myRing);
908}
909///////////////////////////////////////////////////////////////////////////////
910// tests, if ord p_F = nNew.
911static proc equimultTest (poly p_F, int nHelpV, int nNew, ideal HCond)
912{
913  poly p_FnMinus1;
914  ideal id_help;
915  matrix coefMatrix;
916  int i;
917  int nc;
918
919  p_FnMinus1 = jet(p_F, nNew-1, xyVector());
920
921  coefMatrix = coef(p_FnMinus1, xy);
922  nc = ncols(coefMatrix);
923
924  for (i=1; i<=nc; i++)
925  {
926    id_help[i] = coefMatrix[2,i];
927  }
928
929  id_help = T_helpVarElim(id_help, HCond, nHelpV);
930
931  if (reduce(id_help, std(0)) !=0 )
932  {
933    return(0, p_F);
934  }
935
936  p_F = p_F - p_FnMinus1;
937
938  return(1, p_F);
939}
940///////////////////////////////////////////////////////////////////////////////
941// ASSUME: ord(p_F)=nNew
942// tests, if p_F =  p_c*y^nNew for some p_c.
943static proc pPOfYTest (poly p_F, int nHelpV, int nNew, ideal HCond)
944{
945  poly  p_Fn;
946  poly p_c;
947  ideal id_help;
948  int nc;
949  int i=1;
950  matrix coefMatrix;
951
952  p_Fn = jet(p_F, nNew, xyVector());
953
954  coefMatrix = coef(p_Fn, xy);
955  nc = ncols(coefMatrix);
956
957  p_c = coefMatrix[2,1];
958
959  for (i = 2; i <= nc; i++)
960  {
961    id_help[i] = coefMatrix[2,i];
962  }
963
964  id_help = T_helpVarElim(id_help, HCond, nHelpV);
965
966  if (reduce(id_help, std(0)) !=0 )
967  {
968    return(0, p_c);
969  }
970
971  return(1, p_c);
972}
973///////////////////////////////////////////////////////////////////////////////
974// ASSUME: ord(p_F)=nNew
975// tests, if p_F =  p_c*(y - p_a*x)^nNew for some p_a, p_c.
976static proc pPOfLinTest(poly p_F, int nNew, int nHelpV, ideal HCond)
977{
978  poly p_Fn;
979  poly p_c;
980  poly p_ca;
981  poly p_help;
982  poly p_a;
983  ideal id_help;
984
985  p_Fn = jet(p_F, nNew, xyVector());
986
987  p_c  = coefficient(p_Fn,y^nNew,y);
988  p_ca = coefficient(p_Fn,y^(nNew-1)*x,xy)/-nNew;
989
990  if (npars(basering)==1
991  && charstr(basering) != string(char(basering)) + "," + parstr(basering))
992  {
993    p_a = H(nHelpV);
994    HCond = HCond + ideal(p_ca - p_a*p_c);
995  }
996  else
997  {
998    p_help = p_ca/p_c;
999    if (p_help * p_c == p_ca)
1000    {
1001      p_a = p_help;
1002    }
1003    else
1004    {
1005      p_a = H(nHelpV);
1006      HCond = HCond + ideal(p_ca - p_a*p_c);
1007    }
1008  }
1009
1010  id_help = ideal(p_Fn - p_c *(y - p_a *x)^nNew);
1011  id_help = T_helpVarElim(id_help, HCond, nHelpV);
1012
1013  if (reduce(id_help, std(0)) != 0 )
1014  {
1015    return(0, p_F, p_c, p_a, HCond);
1016  }
1017
1018  return(1, p_F, p_c, p_a, HCond);
1019}
1020//////////////////////////////////////////////////////////////////////////////
1021// eliminates the variables H(1),..,H(nHelpV) from the ideal ES + HCond
1022static proc T_helpVarElim(ideal ES, ideal HCond, int nHelpV)
1023{
1024
1025  def r_old = basering;
1026
1027  ideal qIdeal = ideal(basering);
1028
1029  string helpStr;
1030  string minPolyStr = "";
1031
1032  if(minpoly != 0)
1033  {
1034    if (charstr(basering) == string(char(basering)) + "," + parstr(basering))
1035    {
1036      minPolyStr = string(minpoly);
1037    }
1038  }
1039
1040  string str = "ring myRing =
1041               (" + charstr(r_old) + "),(" + string(maxideal(1)) + "),
1042               (dp(" + string( nHelpV) + "), dp);";
1043  execute (str);
1044
1045  if (minPolyStr != "")
1046  {
1047    helpStr = "minpoly =" + minPolyStr + ";";
1048    execute(helpStr);
1049  }
1050
1051  ideal qIdeal = imap(r_old, qIdeal);
1052  if(qIdeal != 0)
1053  {
1054    def r_base = basering;
1055    kill myRing;
1056    qring myRing = std(qIdeal);
1057  }
1058
1059  ideal ES = imap(r_old, ES);
1060  ideal HCond = imap(r_old, HCond);
1061
1062  ES = ES + HCond;
1063  ES = std(ES);
1064  ES = nselect(ES,1,nHelpV);
1065
1066  setring r_old;
1067  ES = imap (myRing,ES);
1068
1069  return(ES);
1070}
1071///////////////////////////////////////////////////////////////////////////////
1072// ASSUME:  F in K[t(1)..t(s),x,y],
1073//          the ringordering is ds
1074// RETURNS: list l, such that
1075//          l[1]=1 if some error has occured,
1076//          l[1]=0 otherwise and then
1077//          l[2] = 1, if the deformation is equisingular and
1078//          l[2] = 0 otherwise.
1079static proc equisingTest (poly p_F, intvec multSeq, int maxStep)
1080{
1081  def r_old = basering;
1082
1083  ideal id_Es = 0;
1084
1085  int isES = 1;
1086  int step = 1;
1087  int nNew = multSeq[step];
1088  int nOld = nNew;
1089  int zeroSteps;
1090  ideal  HCond = 0;
1091  int  nHelpV = 1;
1092
1093  T_extendRing (p_F, nHelpV, HCond);
1094
1095  poly p_c;
1096  poly p_a;
1097
1098  while (step <= maxStep)
1099  {
1100    nOld = nNew;
1101    nNew = multSeq[step];
1102
1103    if (nNew < nOld)
1104    //start a new line in the HNE of F
1105    //               _            _
1106    // for the next | nold/nnew -1 | iteration steps the coefficient 'a'
1107    // in the leading form Fn = c(y-ax) should  be zero
1108    {
1109      p_F = swapXY(p_F);
1110      zeroSteps = calcNZeroSteps (nOld, nNew);
1111    }
1112
1113    isES, p_F = equimultTest (p_F, nHelpV, nNew, HCond);
1114
1115    if (! isES)
1116    {
1117      return(0);
1118    }
1119
1120    if (zeroSteps > 0)
1121    {
1122      zeroSteps--;
1123
1124      isES, p_c = pPOfYTest (p_F,  nHelpV, nNew, HCond);
1125      p_a = 0;
1126    }
1127    else
1128    {
1129      isES, p_F, p_c, p_a, HCond = pPOfLinTest (p_F, nNew, nHelpV, HCond);
1130    }
1131
1132    if (! isES)
1133    {
1134      return(0);
1135    }
1136
1137    p_F = formalBlowUp (p_F, p_c, p_a, nNew);
1138
1139    if (p_a == H(nHelpV))
1140    {
1141      nHelpV++;
1142
1143      def r_base = basering;
1144      kill myRing;
1145
1146      T_extendRing(p_F, nHelpV, HCond);
1147
1148      kill r_base;
1149
1150      poly p_a;
1151      poly p_c;
1152    }
1153
1154    step++;
1155  }
1156
1157  return(1);
1158}
1159///////////////////////////////////////////////////////////////////////////////
1160
1161proc isEquising (poly p_F, list #)
1162"USAGE:   isEquising(F[,m]); F poly, m int
1163ASSUME:   F defines a deformation of an irreducible bivariate polynomial f
1164         and the characteristic of the basering does not divide mult(f). @*
1165         If nv is the number of variables of the basering, then the first nv-2
1166         variables are the deformation parameters. @*
1167         If the basering is a qring, ideal(basering) must only depend
1168         on the deformation parameters.
1169RETURN:  list l of two integers, where
1170@format
1171   l[1]=1 if F is an equisingular deformation,  l[1]=0 otherwise.
1172   l[2]=1 if some error has occured,  l[2]=0 otherwise.
1173@end format
1174NOTE:    If m is given, the computation stops after m steps of the iteration. @*
1175         This procedure uses @code{execute} or calls a procedure using
1176         @code{execute}.
1177EXAMPLE: example isEquising; shows an example
1178"
1179{
1180  def r_user = basering;
1181
1182  int ii = 1;
1183  int i_nvars = nvars(basering);
1184  int error = 0;
1185  int maxStep;
1186  int userMaxStep;
1187  int xNotTransversal = 0;
1188  int fIrreducible = 1;
1189  intvec multSeq;
1190  ideal isES = 1;
1191
1192  error = checkBasering();
1193  if (error)
1194  {
1195    return(0,1);
1196  }
1197
1198  userMaxStep = getInput(#);
1199
1200  // define a new basering "myRing" with new names for parameters
1201  // and variables.
1202  // The new names are 'a(1)', ..., 'a(npars)' for the parameters
1203  // and 't(1)', ..., 't(nvars-2)', 'x', 'y' for the variables.
1204  createMyRing(p_F, "ds");
1205
1206  createHNERing();
1207
1208  poly p_f = imap(myRing,p_F);
1209
1210  error = checkPoly(p_f);
1211  if (error)
1212  {
1213    return(0,1);
1214  }
1215
1216  // compute the multiplicity sequence of p_f.
1217  multSeq, xNotTransversal, fIrreducible = calcMultSequence(p_f);
1218
1219   if ( ! fIrreducible)
1220   {
1221     return(list(0,1));
1222   }
1223
1224  setring myRing;
1225
1226  if (xNotTransversal)
1227  {
1228    p_F = swapXY(p_F);
1229  }
1230
1231  if (userMaxStep != -1 && userMaxStep <  size(multSeq)-1)
1232  {
1233    maxStep = userMaxStep;
1234  }
1235  else
1236  {
1237    maxStep = size(multSeq)-1;
1238  }
1239
1240  int isES = equisingTest(p_F, multSeq, maxStep);
1241
1242  return(list(isES, error));
1243}
1244
1245example
1246{
1247   "EXAMPLE:"; echo=2;
1248   ring r = 11,(a,b,x,y),ds;
1249   poly F = (x2+2xy+y2+x5)+ay3+bx5;
1250   isEquising(F);
1251   isEquising(F,1);
1252   ideal I = ideal(a);
1253   qring q = std(I);
1254   poly F = imap(r,F);
1255   isEquising(F);
1256}
1257///////////////////////////////////////////////////////////////////////////////
1258/*
1259  Weiter Beispiele aus Dipl. von A. Mindnich einfuegen
1260*/
Note: See TracBrowser for help on using the repository browser.