# source:git/Singular/LIB/equising.lib@e7cc147

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