source: git/Singular/LIB/equising.lib @ 7db9488

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