source: git/Singular/LIB/normal.lib @ 9ebd82

jengelh-datetimespielwiese
Last change on this file since 9ebd82 was 9ebd82, checked in by Jakob Kroeker <kroeker@…>, 9 years ago
removed unnesessary changes
  • Property mode set to 100644
File size: 217.4 KB
RevLine 
[380a17b]1//////////////////////////////////////////////////////////////////////////////
[3686937]2version="version normal.lib 4.0.0.0 Jun_2013 "; // $Id$
[1598658]3category="Commutative Algebra";
4info="
5LIBRARY:  normal.lib     Normalization of Affine Rings
6AUTHORS:  G.-M. Greuel,  greuel@mathematik.uni-kl.de,
7@*        S. Laplagne,   slaplagn@dm.uba.ar,
8@*        G. Pfister,    pfister@mathematik.uni-kl.de
9
10
[bf42f0]11PROCEDURES:
[1598658]12 normal(I,[...]);    normalization of an affine ring
13 normalP(I,[...]);   normalization of an affine ring in positive characteristic
14 normalC(I,[...]);   normalization of an affine ring through a chain of rings
15 HomJJ(L);           presentation of End_R(J) as affine ring, J an ideal
16 genus(I);           computes the geometric genus of a projective curve
17 primeClosure(L);    integral closure of R/p, p a prime ideal
18 closureFrac(L);     writes a poly in integral closure as element of Quot(R/p)
19 iMult(L);           intersection multiplicity of the ideals of the list L
20
21 deltaLoc(f,S);      sum of delta invariants at conjugated singular points
22 locAtZero(I);       checks whether the zero set of I is located at 0
23 norTest(I,nor);     checks the output of normal, normalP, normalC
24 getSmallest(J);     computes the polynomial of smallest degree of J
25 getOneVar(J, vari); computes a polynomial of J in the variable vari
26 changeDenominator(U1, c1, c2, I); computes ideal U2 such that 1/c1*U1=1/c2*U2
[1e1ec4]27
28SEE ALSO: locnormal_lib;modnormal_lib
[1598658]29";
30
31LIB "general.lib";
32LIB "poly.lib";
33LIB "sing.lib";
34LIB "primdec.lib";
35LIB "elim.lib";
36LIB "presolve.lib";
37LIB "inout.lib";
38LIB "ring.lib";
39LIB "hnoether.lib";
40LIB "reesclos.lib";
41LIB "algebra.lib";
42
43///////////////////////////////////////////////////////////////////////////////
44
45proc normal(ideal id, list #)
46"USAGE:  normal(id [,choose]); id = radical ideal, choose = list of options. @*
47         Optional parameters in list choose (can be entered in any order):@*
48         Decomposition:@*
[c47547]49         - \"equidim\" -> computes first an equidimensional decomposition of the
50         input ideal, and then the normalization of each component (default).@*
51         - \"prim\" -> computes first the minimal associated primes of the input
52         ideal, and then the normalization of each prime. (When the input ideal
53         is not prime and the minimal associated primes are easy to compute,
54         this method is usually faster than \"equidim\".)@*
[1598658]55         - \"noDeco\" -> no preliminary decomposition is done. If the ideal is
56         not equidimensional radical, output might be wrong.@*
57         - \"isPrim\" -> assumes that the ideal is prime. If this assumption
58         does not hold, the output might be wrong.@*
59         - \"noFac\" -> factorization is avoided in the computation of the
60         minimal associated primes;
61         Other:@*
62         - \"useRing\" -> uses the original ring ordering.@*
63         If this option is set and if the ring ordering is not global, normal
64         will change to a global ordering only for computing radicals and prime
65         or equidimensional decompositions.@*
66         If this option is not set, normal changes to dp ordering and performs
67         all computations with respect to this ordering.@*
68         - \"withDelta\" (or \"wd\") -> returns also the delta invariants.@*
69         If the optional parameter choose is not given or empty, only
70         \"equidim\" but no other option is used.@*
71         - list(\"inputJ\", ideal inputJ) -> takes as initial test ideal the
[ea87a9]72         ideal inputJ. This option is only for use in other procedures. Using
[1598658]73         this option, the result might not be the normalization.@*
[ea87a9]74         (Option only valid for global algorithm.)@*
75         - list(\"inputC\", ideal inputC) -> takes as initial conductor the
76         ideal inputC. This option is only for use in other procedures. Using
[1598658]77         this option, the result might not be the normalization.@*
78         (Option only valid for global algorithm.)@*
[ea87a9]79         Options used for computing integral basis (over rings of two
[1598658]80         variables):@*
[ea87a9]81         - \"var1\" -> uses a polynomial in the first variable as
[1598658]82         universal denominator.@*
[ea87a9]83         - \"var2\" -> uses a polynomial in the second variable as universal
[1598658]84         denominator.@*
85         If the optional parameter choose is not given or empty, only
86         \"equidim\" but no other option is used.@*
87ASSUME:  The ideal must be radical, for non-radical ideals the output may
88         be wrong (id=radical(id); makes id radical). However, when using the
89         \"prim\" option the minimal associated primes of id are computed first
90         and hence normal computes the normalization of the radical of id.@*
91NOTE:    \"isPrim\" should only be used if id is known to be prime.
92RETURN:  a list, say nor, of size 2 (resp. 3 with option \"withDelta\").
93@format  Let R denote the basering and id the input ideal.
94         * nor[1] is a list of r rings, where r is the number of associated
95         primes P_i with option \"prim\" (resp. >= no of equidimenensional
96         components P_i with option \"equidim\").@*
97         Each ring Ri := nor[1][i], i=1..r, contains two ideals with given
98         names @code{norid} and @code{normap} such that: @*
99         - Ri/norid is the normalization of the i-th component, i.e. the
100          integral closure of R/P_i in its field of fractions (as affine ring);
101         - @code{normap} gives the normalization map from R/id to
102           Ri/norid for each i.@*
103         - the direct sum of the rings Ri/norid, i=1,..r, is the normalization
104           of R/id as affine algebra; @*
105         * nor[2] is a list of size r with information on the normalization of
106         the i-th component as module over the basering R:@*
107         nor[2][i] is an ideal, say U, in R such that the integral closure
108         of basering/P_i is generated as module over R by 1/c * U, with c
109         the last element U[size(U)] of U.@*
110         * nor[3] (if option \"withDelta\" is set) is a list of an intvec
111         of size r, the delta invariants of the r components, and an integer,
112         the total delta invariant of basering/id (-1 means infinite, and 0
113         that R/P_i resp. R/id is normal).
114@end format
115THEORY:  We use here a general algorithm described in [G.-M.Greuel, S.Laplagne,
116         F.Seelisch: Normalization of Rings (2009)].@*
117         The procedure computes the R-module structure, the algebra structure
118         and the delta invariant of the normalization of R/id:@*
119         The normalization of R/id is the integral closure of R/id in its total
120         ring of fractions. It is a finitely generated R-module and nor[2]
121         computes R-module generators of it. More precisely: If U:=nor[2][i]
122         and c:=U[size(U)], then c is a non-zero divisor and U/c is an R-module
123         in the total ring of fractions, the integral closure of R/P_i. Since
124         U[size(U)]/c is equal to 1, R/P_i resp. R/id is contained in the
125         integral closure.@*
126         The normalization is also an affine algebra over the ground field
127         and nor[1] presents it as such. For geometric considerations nor[1] is
128         relevant since the variety of the ideal norid in Ri is the
129         normalization of the variety of the ideal P_i in R.@*
130         The delta invariant of a reduced ring A is dim_K(normalization(A)/A).
131         For A=K[x1,...,xn]/id we call this number also the delta invariant of
132         id. nor[3] returns the delta invariants of the components P_i and of
133         id.
134NOTE:    To use the i-th ring type e.g.: @code{def R=nor[1][i]; setring R;}.
135@*       Increasing/decreasing printlevel displays more/less comments
136         (default: printlevel=0).
137@*       Implementation works also for local rings.
138@*       Not implemented for quotient rings.
139@*       If the input ideal id is weighted homogeneous a weighted ordering may
140         be used together with the useRing-option (qhweight(id); computes
141         weights).
142KEYWORDS: normalization; integral closure; delta invariant.
143SEE ALSO: normalC, normalP.
144EXAMPLE: example normal; shows an example
145"
146{
147  intvec opt = option(get);     // Save current options
148
149  int i,j;
150  int decomp;   // Preliminary decomposition:
151                // 0 -> no decomposition (id is assumed to be prime)
152                // 1 -> no decomposition
153                //      (id is assumed to be equidimensional radical)
154                // 2 -> equidimensional decomposition
155                // 3 -> minimal associated primes
156  int noFac, useRing, withDelta;
157  int dbg = printlevel - voice + 2;
158  int nvar = nvars(basering);
159  int chara  = char(basering);
160  int denomOption;   // Method for choosing the conductor
161
162  ideal inputJ = 0;      // Test ideal given in the input (if any).
163  ideal inputC = 0;      // Conductor ideal given in the input (if any).
164
165  list result, resultNew;
166  list keepresult;
167  list ringStruc;
168  ideal U;
169  poly c;
170  int sp;            // Number of components.
171
172  // Default methods:
173  noFac = 0;         // Use facSTD when computing minimal associated primes
174  decomp = 2;        // Equidimensional decomposition
175  useRing = 0;       // Change first to dp ordering, and perform all
176                     // computations there.
177  withDelta = 0;     // Do not compute the delta invariant.
[ea87a9]178  denomOption = 0;   // The default universal denominator is the smallest
[1598658]179                     // degree polynomial.
180
181//--------------------------- define the method ---------------------------
182  for ( i=1; i <= size(#); i++ )
183  {
184    if ( typeof(#[i]) == "string" )
185    {
186//--------------------------- choosen methods -----------------------
187      if ( (#[i]=="isprim") or (#[i]=="isPrim") )
188      {decomp = 0;}
189
190      if ( (#[i]=="nodeco") or (#[i]=="noDeco") )
191      {decomp = 1;}
192
193      if (#[i]=="prim")
194      {decomp = 3;}
195
196      if (#[i]=="equidim")
197      {decomp = 2;}
198
199      if ( (#[i]=="nofac") or (#[i]=="noFac") )
200      {noFac=1;}
201
202      if ( ((#[i]=="useRing") or (#[i]=="usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C"))
203      {useRing = 1;}
204
205      if ( (#[i]=="withDelta") or (#[i]=="wd") or (#[i]=="withdelta"))
206      {
207        if((decomp == 0) or (decomp == 3))
208        {
209          withDelta = 1;
210        }
211        else
212        {
213          decomp = 3;
214          withDelta = 1;
215          //Note: the delta invariants cannot be computed with an equidimensional
216          //decomposition, hence we compute first the minimal primes
217        }
218      }
219      if (#[i]=="var1")
220      {denomOption = 1;}
221      if (#[i]=="var2")
222      {denomOption = 2;}
223    }
224    if(typeof(#[i]) == "list"){
225      if(size(#[i]) == 2){
226        if (#[i][1]=="inputJ"){
227          if(typeof(#[i][2]) == "ideal"){
228            inputJ = #[i][2];
229          }
230        }
231      }
232      if (#[i][1]=="inputC"){
233        if(size(#[i]) == 2){
234          if(typeof(#[i][2]) == "ideal"){
235            inputC = #[i][2];
236          }
237        }
238      }
239    }
240  }
241  kill #;
242
243//------------------------ change ring if required ------------------------
244// If the ordering is not global, we change to dp ordering for computing the
245// min ass primes.
246// If the ordering is global, but not dp, and useRing = 0, we also change to
247// dp ordering.
248
[1e1ec4]249  int isGlobal = attrib(basering,"global");// Checks if the original ring has
[1598658]250                                          // global ordering.
251
252  def origR = basering;   // origR is the original ring
253                          // R is the ring where computations will be done
254
255  if((useRing  == 1) and (isGlobal == 1))
256  {
257    def globR = basering;
258  }
259  else
260  {
261    // We change to dp ordering.
262    list rl = ringlist(origR);
263    list origOrd = rl[3];
264    list newOrd = list("dp", intvec(1:nvars(origR))), list("C", 0);
265    rl[3] = newOrd;
266    def globR = ring(rl);
267    setring globR;
268    ideal id = fetch(origR, id);
269  }
270
271//------------------------ trivial checkings ------------------------
272  id = groebner(id);
273  if((size(id) == 0) or (id[1] == 1))
274  {
275    // The original ring R/I was normal. Nothing to do.
276    // We define anyway a new ring, equal to R, to be able to return it.
277    setring origR;
278    list lR = ringlist(origR);
279    def ROut = ring(lR);
280    setring ROut;
281    ideal norid = fetch(origR, id);
282    ideal normap = maxideal(1);
283    export norid;
284    export normap;
285    setring origR;
286    if(withDelta)
287    {
288      result = list(list(ROut), list(ideal(1)), list(intvec(0), 0));
289    }
290    else
291    {
292      result = list(list(ROut), list(ideal(1)));
293    }
294    sp = 1;      // number of rings in the output
295    option(set, opt);
296    normalOutputText(dbg, withDelta, sp);
297    return(result);
298  }
299//------------------------ preliminary decomposition-----------------------
300  list prim;
301  if(decomp == 2)
302  {
303    dbprint(dbg, "// Computing the equidimensional decomposition...");
304    prim = equidim(id);
305  }
306  if((decomp == 0) or (decomp == 1))
307  {
308    prim = id;
309  }
310  if(decomp == 3)
311  {
312    dbprint(dbg, "// Computing the minimal associated primes...");
313    if( noFac )
314    { prim = minAssGTZ(id,1); }
315    else
316    { prim = minAssGTZ(id); }
317  }
318  sp = size(prim);
319  if(dbg>=1)
320  {
321    prim; "";
322    "// number of components is", sp;
323    "";
324  }
325
326
327//----------------- back to the original ring if required ------------------
328// if ring was not global and useRing is on, we go back to the original ring
329  if((useRing == 1) and (isGlobal != 1))
330  {
331    setring origR;
332    def R = basering;
333    list prim = fetch(globR, prim);
334  }
335  else
336  {
337    def R = basering;
338    ideal inputJ = fetch(origR, inputJ);
339    ideal inputC = fetch(origR, inputC);
340    if(useRing == 0)
341    {
342      ideal U;
343      poly c;
344    }
345  }
346
347// ---------------- normalization of the components-------------------------
348// calls normalM to compute the normalization of each component.
349
350  list norComp;       // The normalization of each component.
351  int delt;
352  int deltI = 0;
353  int totalComps = 0;
354
355  setring origR;
356  def newROrigOrd;
357  list newRListO;
358  setring R;
359  def newR;
360  list newRList;
361
362  for(i=1; i<=size(prim); i++)
363  {
364    if(dbg>=2){pause();}
365    if(dbg>=1)
366    {
367      "// start computation of component",i;
368      "   --------------------------------";
369    }
370    if(groebner(prim[i])[1] != 1)
371    {
372      if(dbg>=2)
373      {
374        "We compute the normalization in the ring"; basering;
375      }
376      printlevel = printlevel + 1;
377      norComp = normalM(prim[i], decomp, withDelta, denomOption, inputJ, inputC);
378      printlevel = printlevel - 1;
379      for(j = 1; j <= size(norComp); j++)
380      {
381        newR = norComp[j][3];
[1e1ec4]382        if(!defined(savebasering)) { def savebasering;}
[aacb5a]383        savebasering=basering;
[6e80ab]384        setring newR; // must be in a compatible ring to newR
385                      // as ringlist may produce ring-dep. stuff
[1e1ec4]386        if(!defined(newRList)) { list newRList;}
[1598658]387        newRList = ringlist(newR);
[6e80ab]388        setring savebasering;
[1598658]389        U = norComp[j][1];
390        c = norComp[j][2];
391        if(withDelta)
392        {
393          delt = norComp[j][4];
394          if((delt >= 0) and (deltI >= 0))
395          {
396            deltI = deltI + delt;
397          }
398          else
399          {
400            deltI = -1;
401          }
402        }
403        // -- incorporate result for this component to the list of results ---
404        if(useRing == 0)
405        {
406          // We go back to the original ring.
407          setring origR;
408          U = fetch(R, U);
409          c = fetch(R, c);
[6e80ab]410          newRListO = imap(newR, newRList);
[1598658]411          // We change the ordering in the new ring.
412          if(nvars(newR) > nvars(origR))
413          {
414            newRListO[3]=insert(origOrd, newRListO[3][1]);
415          }
416          else
417          {
418            newRListO[3] = origOrd;
419          }
420          newROrigOrd = ring(newRListO);
421          setring newROrigOrd;
422          ideal norid = imap(newR, norid);
423          ideal normap = imap(newR, normap);
424          export norid;
425          export normap;
426          setring origR;
427          totalComps++;
428          result[totalComps] = list(U, c, newROrigOrd);
429          if(withDelta)
430          {
431            result[totalComps] = insert(result[totalComps], delt, 3);
432          }
433          setring R;
434        }
435        else
436        {
437          setring R;
438          totalComps++;
439          result[totalComps] = norComp[j];
440        }
441      }
442    }
443  }
444
445// -------------------------- delta computation ----------------------------
446  if(withDelta == 1)
447  {
448    // Intersection multiplicities of list prim, sp=size(prim).
449    if ( dbg >= 1 )
450    {
451      "// Sum of delta for all components: ", deltI;
452    }
453    if(size(prim) > 1)
454    {
455      dbprint(dbg, "// Computing the sum of the intersection multiplicities of the components...");
456      int mul = iMult(prim);
457      if ( mul < 0 )
458      {
459        deltI = -1;
460      }
461      else
462      {
463        deltI = deltI + mul;
464      }
465      if ( dbg >= 1 )
466      {
467        "// Intersection multiplicity is : ", mul;
468      }
469    }
470  }
471
472// -------------------------- prepare output ------------------------------
473  setring origR;
474
475  list RL;      // List of rings
476  list MG;      // Module generators
477  intvec DV;    // Vector of delta's of each component
478  for(i = 1; i <= size(result); i++)
479  {
480    RL[i] = result[i][3];
481    MG[i] = lineUpLast(result[i][1], result[i][2]);
482    if(withDelta)
483    {
484      DV[i] = result[i][4];
485    }
486  }
487  if(withDelta)
488  {
489    resultNew = list(RL, MG, list(DV, deltI));
490  }
491  else
492  {
493    resultNew = list(RL, MG);
494  }
495  sp = size(RL);              //RL = list of rings
496
497  option(set, opt);
498  normalOutputText(dbg, withDelta, sp);
499  return(resultNew);
500}
501
502example
503{ "EXAMPLE:";
504  printlevel = printlevel+1;
505  echo = 2;
506  ring s = 0,(x,y),dp;
507  ideal i = (x2-y3)*(x2+y2)*x;
508  list nor = normal(i, "withDelta", "prim");
509  nor;
510
511  // 2 branches have delta = 1, and 1 branch has delta = 0
512  // the total delta invariant is 13
513
514  def R2 = nor[1][2];  setring R2;
515  norid; normap;
516
517  echo = 0;
518  printlevel = printlevel-1;
519  pause("   hit return to continue"); echo=2;
520
521  ring r = 2,(x,y,z),dp;
522  ideal i = z3-xy4;
523  list nor = normal(i, "withDelta", "prim");  nor;
524  // the delta invariant is infinite
525  // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module
526  // in its quotient field Quot(r/i)
527
528  // the normalization as affine algebra over the ground field:
529  def R = nor[1][1]; setring R;
530  norid; normap;
531}
532
533///////////////////////////////////////////////////////////////////////////////
534// Prints the output text in proc normal.
535//
536static proc normalOutputText(int dbg, int withDelta, int sp)
537// int dbg: printlevel
538// int withDelta: output contains information about the delta invariant
539// int sp: number of output rings.
540{
541  if ( dbg >= 0 )
542  {
543    "";
544    if(!withDelta)
545    {
546      "// 'normal' created a list, say nor, of two elements.";
547    }
548    else
549    {
550      "// 'normal' created a list, say nor, of three elements.";
551    }
552    "// To see the list type";
553    "      nor;";
554    "";
555    "// * nor[1] is a list of", sp, "ring(s).";
556    "// To access the i-th ring nor[1][i], give it a name, say Ri, and type";
557    "     def R1 = nor[1][1]; setring R1; norid; normap;";
558    "// For the other rings type first (if R is the name of your base ring)";
559    "     setring R;";
560    "// and then continue as for R1.";
561    "// Ri/norid is the affine algebra of the normalization of R/P_i where";
562    "// P_i is the i-th component of a decomposition of the input ideal id";
563    "// and normap the normalization map from R to Ri/norid.";
564    "";
565    "// * nor[2] is a list of", sp, "ideal(s). Let ci be the last generator";
566    "// of the ideal nor[2][i]. Then the integral closure of R/P_i is";
567    "// generated as R-submodule of the total ring of fractions by";
568    "// 1/ci * nor[2][i].";
569
570    if(withDelta)
571    { "";
572      "// * nor[3] is a list of an intvec of size", sp, "the delta invariants ";
573      "// of the components, and an integer, the total delta invariant ";
574      "// of R/id (-1 means infinite, and 0 that R/P_i resp. R/id is normal).";
575    }
576  }
577}
578
579
580///////////////////////////////////////////////////////////////////////////////
581
582proc HomJJ (list Li)
583"USAGE:   HomJJ (Li);  Li = list: ideal SBid, ideal id, ideal J, poly p
584ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
585@*       SBid = standard basis of id,
586@*       J    = ideal of P containing the polynomial p,
587@*       p    = nonzero divisor of R
588COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as
589         affine ring, together with the map R --> Hom_R(J,J) of affine rings,
590         where R is the quotient ring of P modulo the standard basis SBid.
591RETURN:  a list l of three objects
592@format
593         l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
594               such that l[1]/endid = Hom_R(J,J) and
595               endphi describes the canonical map R -> Hom_R(J,J)
596         l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
597         l[3] : an integer, = dim_K(Hom_R(J,J)/R) (the contribution to delta)
598                if the dimension is finite, -1 otherwise
599@end format
600NOTE:    printlevel >=1: display comments (default: printlevel=0)
601EXAMPLE: example HomJJ;  shows an example
602"
603{
604//---------- initialisation ---------------------------------------------------
605   int isIso,isPr,isHy,isCo,isRe,isEq,oSAZ,ii,jj,q,y;
606   intvec rw,rw1;
607   list L;
608   y = printlevel-voice+2;  // y=printlevel (default: y=0)
609   def P = basering;
610   ideal SBid, id, J = Li[1], Li[2], Li[3];
611   poly p = Li[4];
612   int noRed = 0;
613   if(size(Li) > 4)
614   {
615     if(Li[5] == 1) { noRed = 1; }
616   }
617
618   attrib(SBid,"isSB",1);
619   int homo = homog(Li[2]);               //is 1 if id is homogeneous, 0 if not
620
621//---- set attributes for special cases where algorithm can be simplified -----
622   if( homo==1 )
623   {
624      rw = ringweights(P);
625   }
626   if( typeof(attrib(id,"isPrim"))=="int" )
627   {
628      if(attrib(id,"isPrim")==1)  { isPr=1; }
629   }
630   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
631   {
632      if(attrib(id,"onlySingularAtZero")==1){oSAZ=1; }
633   }
634   if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
635   {
636      if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; }
637   }
638   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
639   {
640      if(attrib(id,"isCohenMacaulay")==1) { isCo=1; }
641   }
642   if( typeof(attrib(id,"isRegInCodim2"))=="int" )
643   {
644      if(attrib(id,"isRegInCodim2")==1) { isRe=1; }
645   }
646   if( typeof(attrib(id,"isEquidimensional"))=="int" )
647   {
648      if(attrib(id,"isEquidimensional")==1) { isEq=1; }
649   }
650//-------------------------- go to quotient ring ------------------------------
651   qring R  = SBid;
652   ideal id = fetch(P,id);
653   ideal J  = fetch(P,J);
654   poly p   = fetch(P,p);
655   ideal f,rf,f2;
656   module syzf;
657//---------- computation of p*Hom(J,J) as R-ideal -----------------------------
658   if ( y>=1 )
659   {
660     "// compute p*Hom(J,J) = p*J:J";
661     "//   the ideal J:";J;
662   }
663   f  = quotient(p*J,J);
664
665   //### (neu GMG 4.10.08) divide by the greatest common divisor:
666   poly gg = gcd( f[1],p );
667   for(ii=2; ii <=ncols(f); ii++)
668   {
669      gg=gcd(gg,f[ii]);
670   }
671   for(ii=1; ii<=ncols(f); ii++)
672   {
673      f[ii]=f[ii]/gg;
674   }
675   p = p/gg;
676
677   if ( y>=1 )
678   {
679      "//   the non-zerodivisor p:"; p;
680      "//   the module p*Hom(J,J) = p*J:J :"; f;
681      "";
682   }
683   f2 = std(p);
684
685//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------
686
687   //rf = interred(reduce(f,f2));
688   //### interred hier weggelassen, unten zugefuegt
689   rf = reduce(f,f2);       //represents p*Hom(J,J)/p*R = Hom(J,J)/R
690   if ( size(rf) == 0 )
691   {
692      if ( homog(f) && find(ordstr(basering),"s")==0 )
693      {
694         ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp);
695      }
696      else
697      {
698         ring newR1 = char(P),(X(1..nvars(P))),dp;
699      }
700      ideal endphi = maxideal(1);
701      ideal endid = fetch(P,id);
702      endid = simplify(endid,2);
703      L = substpart(endid,endphi,homo,rw);   //## hier substpart
704      def lastRing = L[1];
705      setring lastRing;
706
707      attrib(endid,"onlySingularAtZero",oSAZ);
708      attrib(endid,"isCohenMacaulay",isCo);
709      attrib(endid,"isPrim",isPr);
710      attrib(endid,"isIsolatedSingularity",isIso);
711      attrib(endid,"isRegInCodim2",isRe);
712      attrib(endid,"isEqudimensional",isEq);
713      attrib(endid,"isHypersurface",0);
714      attrib(endid,"isCompleteIntersection",0);
715      attrib(endid,"isRadical",0);
716      L=lastRing;
717      L = insert(L,1,1);
718      dbprint(y,"// case R = Hom(J,J)");
719      if(y>=1)
720      {
721         "//   R=Hom(J,J)";
722         lastRing;
723         "//   the new ideal";
724         endid;
725         "   ";
726         "//   the old ring";
727         P;
728         "//   the old ideal";
729         setring P;
730         id;
731         "   ";
732         setring lastRing;
733         "//   the map to the new ring";
734         endphi;
735         "   ";
736         pause();
737         "";
738      }
739      setring P;
740      L[3]=0;
741      return(L);
742   }
743   if(y>=1)
744   {
745      "// R is not equal to Hom(J,J), we have to try again";
746      pause();
747      "";
748   }
749//---------- Hom(J,J) != R: create new ring and map from old ring -------------
750// the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module
751// f2=p (i.e. ideal generated by p)
752
753   //f = mstd(f)[2];              //### geaendert GMG 04.10.08
754   //ideal ann = quotient(f2,f);  //### f durch rf ersetzt
755   rf = mstd(rf)[2];              //rf = NF(f,p), hence <p,rf> = <p,f>
756   ideal ann = quotient(f2,rf);   //p:f = p:rf
757
758   //------------- compute the contribution to delta ----------
759   //delt=dim_K(Hom(JJ)/R (or -1 if infinite)
760
761   int delt=vdim(std(modulo(f,ideal(p))));
762
763   f = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
764   q = size(f);
765   syzf = syz(f);
766
767   if ( homo==1 )
768   {
769      rw1 = rw,0;
770      for ( ii=2; ii<=q; ii++ )
771      {
772         rw  = rw, deg(f[ii])-deg(f[1]);
773         rw1 = rw1, deg(f[ii])-deg(f[1]);
774      }
775      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp);
776   }
777   else
778   {
779      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp;
780   }
781
782   //map psi1 = P,maxideal(1);          //### psi1 durch fetch ersetzt
783   //ideal SBid = psi1(SBid);
784   ideal SBid = fetch(P,SBid);
785   attrib(SBid,"isSB",1);
786
787   qring newR = std(SBid);
788
789   //map psi = R,ideal(X(1..nvars(R)));  //### psi durch fetch ersetzt
790   //ideal id = psi(id);
791   //ideal f = psi(f);
792   //module syzf = psi(syzf);
793   ideal id = fetch(R,id);
794   ideal f = fetch(R,f);
795   module syzf = fetch(R,syzf);
796   ideal pf,Lin,Quad,Q;
797   matrix T,A;
798   list L1;
799
800//---------- computation of Hom(J,J) as affine ring ---------------------------
801// determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1),
802// Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. This is of course
803// the same as the kernel of R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi.
804// It is a fact, that the kernel is generated by the linear and the quadratic
805// relations
806// f=p,rf, rf=reduce(f,p), generates pJ:J mod(p),
807// i.e. p*Hom(J,J)/p*R as R-module
808
809   pf = f[1]*f;
810   T = matrix(ideal(T(1..q)),1,q);
811   Lin = ideal(T*syzf);
812   if(y>=1)
813   {
814      "// the ring structure of Hom(J,J) as R-algebra";
815      "//   the linear relations:";
816      Lin;
817   }
818
819   poly ff;
820   for (ii=2; ii<=q; ii++ )
821   {
822      for ( jj=2; jj<=ii; jj++ )
823      {
824         ff = NF(f[ii]*f[jj],std(0));       //this makes lift much faster
825         A = lift(pf,ff);                   //ff lin. comb. of elts of pf mod I
826         Quad = Quad, ideal(T(jj)*T(ii) - T*A);  //quadratic relations
827      }
828   }
829
830   if(y>=1)
831   {
832      "//   the quadratic relations";
833      Quad;
834      pause();
835      newline;
836   }
837   Q = Lin,Quad;
838   Q = subst(Q,T(1),1);
839   //Q = mstd(Q)[2];            //### sehr aufwendig, daher weggelassen (GMG)
840   //### ev das neue interred
841   //mstd dient nur zum verkleinern, die SB-Eigenschaft geht spaeter verloren
842   //da in neuen Ring abgebildet und mit id vereinigt
843
844//---------- reduce number of variables by substitution, if possible ----------
845   if (homo==1)
846   {
847      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp);
848   }
849   else
850   {
851      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp;
852   }
853
854   ideal endid  = imap(newR,id),imap(newR,Q);
855   //hier wird Q weiterverwendet, die SB-Eigenschaft wird nicht verwendet.
856   endid = simplify(endid,2);
857   ideal endphi = ideal(X(1..nvars(R)));
858
859
860  if(noRed == 0)
861  {
862    L = substpart(endid,endphi,homo,rw);
863    def lastRing=L[1];
864    setring lastRing;
865    //return(lastRing);
866  }
867  else
868  {
869    list RL = ringlist(newRing);
870    def lastRing = ring(RL);
871    setring lastRing;
872    ideal endid = fetch(newRing, endid);
873    ideal endphi = fetch(newRing, endphi);
874    export(endid);
875    export(endphi);
876    //def lastRing = newRing;
877    //setring R;
878    //return(newR);
879  }
880
881
882//   L = substpart(endid,endphi,homo,rw);
883
884//   def lastRing=L[1];
885//   setring lastRing;
886
887   attrib(endid,"onlySingularAtZero",0);
888   map sigma=R,endphi;
889   ideal an=sigma(ann);
890   export(an);  //noetig?
891   //ideal te=an,endid;
892   //if(isIso && (size(reduce(te,std(maxideal(1))))==0))   //#### ok???
893   // {
894   //    attrib(endid,"onlySingularAtZero",oSAZ);
895   // }
896   //kill te;
897   attrib(endid,"isCohenMacaulay",isCo);                //#### ok???
898   attrib(endid,"isPrim",isPr);
899   attrib(endid,"isIsolatedSingularity",isIso);
900   attrib(endid,"isRegInCodim2",isRe);
901   attrib(endid,"isEquidimensional",isEq);
902   attrib(endid,"isHypersurface",0);
903   attrib(endid,"isCompleteIntersection",0);
904   attrib(endid,"isRadical",0);
905   if(y>=1)
906   {
907      "// the new ring after reduction of the number of variables";
908      lastRing;
909      "//   the new ideal";
910      endid;  "";
911      "// the old ring";
912      P;
913      "//   the old ideal";
914      setring P;
915      id;
916      "   ";
917      setring lastRing;
918      "//   the map to the new ring";
919      endphi;
920      "   ";
921      pause();
922      "";
923   }
924   L = lastRing;
925   L = insert(L,0,1);
926   L[3] = delt;
927   return(L);
928}
929example
930{"EXAMPLE:";  echo = 2;
931  ring r   = 0,(x,y),wp(2,3);
932  ideal id = y^2-x^3;
933  ideal J  = x,y;
934  poly p   = x;
935  list Li  = std(id),id,J,p;
936  list L   = HomJJ(Li);
937  def end = L[1];    // defines ring L[1], containing ideals endid, endphi
938  setring end;       // makes end the basering
939  end;
940  endid;             // end/endid is isomorphic to End(r/id) as ring
941  map psi = r,endphi;// defines the canonical map r/id -> End(r/id)
942  psi;
943  L[3];              // contribution to delta
944}
945
946
947///////////////////////////////////////////////////////////////////////////////
948//compute intersection multiplicities as needed for delta(I) in
949//normalizationPrimes and normalP:
950
951proc iMult (list prim)
952"USAGE:   iMult(L);  L a list of ideals
953RETURN:  int, the intersection multiplicity of the ideals of L;
954         if iMult(L) is infinite, -1 is returned.
955THEORY:  If r=size(L)=2 then iMult(L) = vdim(std(L[1]+L[2])) and in general
956         iMult(L) = sum{ iMult(L[j],Lj) | j=1..r-1 } with Lj the intersection
957         of L[j+1],...,L[r]. If I is the intersection of all ideals in L then
958         we have delta(I) = delta(L[1])+...+delta(L[r]) + iMult(L) where
959         delta(I) = vdim (normalisation(R/I)/(R/I)), R the basering.
960EXAMPLE: example iMult; shows an example
961"
962{    int i,mul,mu;
963     int sp = size(prim);
964     int y = printlevel-voice+2;
965     if ( sp > 1 )
966     {
967        ideal I(sp-1) = prim[sp];
968        mu = vdim(std(I(sp-1)+prim[sp-1]));
969        mul = mu;
970        if ( y>=1 )
971        {
972          "// intersection multiplicity of component",sp,"with",sp-1,":"; mu;
973        }
974        if ( mu >= 0 )
975        {
976           for (i=sp-2; i>=1 ; i--)
977           {
978              ideal I(i) = intersect(I(i+1),prim[i+1]);
979              mu = vdim(std(I(i)+prim[i]));
980              if ( mu < 0 )
981              {
982                break;
983              }
984              mul = mul + mu;
985              if ( y>=1 )
986              {
987                "// intersection multiplicity of components",sp,"...",i+1,"with",i; mu;
988              }
989           }
990        }
991     }
992     return(mul);
993}
994example
995{ "EXAMPLE:"; echo = 2;
996   ring s  = 23,(x,y),dp;
997   list L = (x-y),(x3+y2);
998   iMult(L);
999   L = (x-y),(x3+y2),(x3-y4);
1000   iMult(L);
1001}
1002///////////////////////////////////////////////////////////////////////////////
1003//check if I has a singularity only at zero, as needed in normalizationPrimes
1004
1005proc locAtZero (ideal I)
1006"USAGE:   locAtZero(I);  I = ideal
1007RETURN:  int, 1 if I has only one point which is located at zero, 0 otherwise
1008ASSUME:  I is given as a standard bases in the basering
1009NOTE:    only useful in affine rings, in local rings vdim does the check
1010EXAMPLE: example locAtZero; shows an example
1011"
1012{
1013   int ii,jj, caz;                   //caz: conzentrated at zero
1014   int dbp = printlevel-voice+2;
1015   int nva = nvars(basering);
1016   int vdi = vdim(I);
1017   if ( vdi < 0 )
1018   {
1019      if (dbp >=1)
1020      { "// non-isolated singularitiy";""; }
1021      return(caz);
1022   }
1023
1024   //Now the ideal is 0-dim
1025   //First an easy test
1026   //If I is homogenous and not constant it is concentrated at 0
1027   if( homog(I)==1 && size(jet(I,0))==0)
1028   {
1029      caz=1;
1030      if (dbp >=1)
1031      { "// isolated singularity and homogeneous";""; }
1032      return(caz);
1033   }
1034
1035   //Now the general case with I 0-dim. Choose an appropriate power pot,
1036   //and check each variable x whether x^pot is in I.
1037   int mi1 = mindeg1(lead(I));
1038   int pot = vdi;
1039   if ( (mi1+(mi1==1))^2 < vdi )
1040   {
1041      pot = (mi1+(mi1==1))^2;      //### alternativ: pot = vdi lassen
1042   }
1043
1044   while ( 1 )
1045   {
1046      caz = 1;
1047      for ( ii=1; ii<= nva; ii++ )
1048      {
1049        if ( NF(var(ii)^pot,I) != 0 )
1050        {
1051           caz = 0; break;
1052        }
1053      }
1054      if ( caz == 1 || pot >= vdi )
1055      {
1056        if (dbp >=1)
1057        {
1058          "// mindeg, exponent, vdim used in 'locAtZero':", mi1,pot,vdi; "";
1059        }
1060        return(caz);
1061      }
1062      else
1063      {
1064        if ( pot^2 < vdi )
1065        { pot = pot^2; }
1066        else
1067        { pot = vdi; }
1068      }
1069   }
1070}
1071example
1072{ "EXAMPLE:"; echo = 2;
1073   ring r = 0,(x,y,z),dp;
1074   poly f = z5+y4+x3+xyz;
1075   ideal i = jacob(f),f;
1076   i=std(i);
1077   locAtZero(i);
1078   i= std(i*ideal(x-1,y,z));
1079   locAtZero(i);
1080}
1081
1082///////////////////////////////////////////////////////////////////////////////
1083
1084//The next procedure normalizationPrimes computes the normalization of an
1085//irreducible or an equidimensional ideal i.
1086//- If i is irreducuble, then the returned list, say nor, has size 2
1087//with nor[1] the normalization ring and nor[2] the delta invariant.
1088//- If i is equidimensional, than the "splitting tools" can create a
1089//decomposition of i and nor can have more than 1 ring.
1090
1091static proc normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #)
1092"USAGE:   normalizationPrimes(i,ihp,delt[,si]);  i = equidimensional ideal,
1093         ihp = map (partial normalization), delt = partial delta-invariant,
1094         si = ideal s.t. V(si) contains singular locus (optional)
1095RETURN:   a list of rings, say nor, and an integer, the delta-invariant
1096          at the end of the list.
1097          each ring nor[j], j = 1..size(nor)-1, contains two ideals
1098          with given names norid and normap such that
1099           - the direct sum of the rings nor[j]/norid is
1100             the normalization of basering/i;
1101           - normap gives the normalization map from basering/id
1102             to nor[j]/norid (for each j)
1103          nor[size(nor)] = dim_K(normalisation(P/i) / (P/i)) is the
1104          delta-invariant, where P is the basering.
1105EXAMPLE: example normalizationPrimes; shows an example
1106"
1107{
1108   //Note: this procedure calls itself as long as the test for
1109   //normality, i.e if R==Hom(J,J), is negative.
1110
1111   int printlev = printlevel;   //store printlevel in order to reset it later
1112   int y = printlevel-voice+2;  // y=printlevel (default: y=0)
1113   if(y>=1)
1114   {
1115     "";
1116     "// START a normalization loop with the ideal";
1117     i;  "";
1118     "// in the ring:";
1119     basering;  "";
1120     pause();
1121     "";
1122   }
1123
1124   def BAS=basering;
1125   list result,keepresult1,keepresult2,JM,gnirlist;
1126   ideal J,SB,MB;
1127   int depth,lauf,prdim,osaz;
1128   int ti=timer;
1129
1130   gnirlist = ringlist(BAS);
1131
1132//----------- the trivial case of a zero ideal as input, RETURN ------------
1133   if(size(i)==0)
1134   {
1135      if(y>=1)
1136      {
1137          "// the ideal was the zero-ideal";
1138      }
1139     // execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
1140     //                 +ordstr(basering)+");");
1141         def newR7 = ring(gnirlist);
1142         setring newR7;
1143         ideal norid=ideal(0);
1144         ideal normap=fetch(BAS,ihp);
1145         export norid;
1146         export normap;
1147         result=newR7;
1148         result[size(result)+1]=list(delt,delti);
1149         setring BAS;
1150         return(result);
1151   }
1152
1153//--------------- General NOTATION, compute SB of input -----------------
1154// SM is a list, the result of mstd(i)
1155// SM[1] = SB of input ideal i,
1156// SM[2] = (minimal) generators for i.
1157// We work with SM and will copy the attributes from i to SM[2]
1158// JM will be a list, either JM[1]=maxideal(1),JM[2]=maxideal(1)
1159// in case i has onlySingularAtZero, or JM = mstd(si) where si = #[1],
1160// or JM = mstd(J) where J is the ideal of the singular locus
1161// JM[2] must be (made) radical
1162
1163   if(y>=1)
1164   {
1165     "// SB-computation of the ideal";
1166   }
1167
1168   list SM = mstd(i);              //Now the work starts
1169   int dimSM =  dim(SM[1]);        //dimension of variety to normalize
1170   if(y>=1)
1171   {
1172      "// the dimension is:";  dimSM;
1173   }
1174//----------------- the general case, set attributes ----------------
1175   //Note: onlySingularAtZero is NOT preserved under the ring extension
1176   //basering --> Hom(J,J) (in contrast to isIsolatedSingularity),
1177   //therefore we reset it:
1178
1179   attrib(i,"onlySingularAtZero",0);
1180
1181   if(attrib(i,"isPrim")==1)
1182   {
1183      attrib(SM[2],"isPrim",1);
1184   }
1185   else
1186   {
1187      attrib(SM[2],"isPrim",0);
1188   }
1189   if(attrib(i,"isIsolatedSingularity")==1)
1190   {
1191      attrib(SM[2],"isIsolatedSingularity",1);
1192   }
1193   else
1194   {
1195      attrib(SM[2],"isIsolatedSingularity",0);
1196   }
1197   if(attrib(i,"isCohenMacaulay")==1)
1198   {
1199      attrib(SM[2],"isCohenMacaulay",1);
1200   }
1201   else
1202   {
1203      attrib(SM[2],"isCohenMacaulay",0);
1204   }
1205   if(attrib(i,"isRegInCodim2")==1)
1206   {
1207      attrib(SM[2],"isRegInCodim2",1);
1208   }
1209   else
1210   {
1211      attrib(SM[2],"isRegInCodim2",0);
1212   }
1213   if(attrib(i,"isEquidimensional")==1)
1214   {
1215      attrib(SM[2],"isEquidimensional",1);
1216   }
1217   else
1218   {
1219      attrib(SM[2],"isEquidimensional",0);
1220   }
1221   if(attrib(i,"isCompleteIntersection")==1)
1222   {
1223     attrib(SM[2],"isCompleteIntersection",1);
1224   }
1225   else
1226   {
1227      attrib(SM[2],"isCompleteIntersection",0);
1228   }
1229   if(attrib(i,"isHypersurface")==1)
1230   {
1231     attrib(SM[2],"isHypersurface",1);
1232   }
1233   else
1234   {
1235      attrib(SM[2],"isHypersurface",0);
1236   }
1237
1238   if(attrib(i,"onlySingularAtZero")==1)
1239   {
1240      attrib(SM[2],"onlySingularAtZero",1);
1241   }
1242   else
1243   {
1244      attrib(SM[2],"onlySingularAtZero",0);
1245   }
1246
1247   //------- an easy and cheap test for onlySingularAtZero ---------
1248   if( (attrib(SM[2],"isIsolatedSingularity")==1) && (homog(SM[2])==1) )
1249   {
1250      attrib(SM[2],"onlySingularAtZero",1);
1251   }
1252
1253//-------------------- Trivial cases, in each case RETURN ------------------
1254// input ideal is the ideal of a partial normalization
1255
1256   // ------------ Trivial case: input ideal contains a unit ---------------
1257   if( dimSM == -1)
1258   {  "";
1259      "      // A unit ideal was found.";
1260      "      // Stop with partial result computed so far";"";
1261
1262         MB=SM[2];
1263         intvec rw;
1264         list LL=substpart(MB,ihp,0,rw);
1265         def newR6=LL[1];
1266         setring newR6;
1267         ideal norid=endid;
1268         ideal normap=endphi;
1269         kill endid,endphi;
1270         export norid;
1271         export normap;
1272         result=newR6;
1273         result[size(result)+1]=list(delt,delti);
1274         setring BAS;
1275         return(result);
1276   }
1277
1278   // --- Trivial case: input ideal is zero-dimensional and homog ---
1279   if( (dim(SM[1])==0) && (homog(SM[2])==1) )
1280   {
1281      if(y>=1)
1282      {
1283         "// the ideal was zero-dimensional and homogeneous";
1284      }
1285      MB=maxideal(1);
1286      intvec rw;
1287      list LL=substpart(MB,ihp,0,rw);
1288      def newR5=LL[1];
1289      setring newR5;
1290      ideal norid=endid;
1291      ideal normap=endphi;
1292      kill endid,endphi;
1293      export norid;
1294      export normap;
1295      result=newR5;
1296      result[size(result)+1]=list(delt,delti);
1297      setring BAS;
1298      return(result);
1299   }
1300
1301   // --- Trivial case: input ideal defines a line ---
1302   //the one-dimensional, homogeneous case and degree 1 case
1303   if( (dim(SM[1])==1) && (maxdeg1(SM[2])==1) && (homog(SM[2])==1) )
1304   {
1305      if(y>=1)
1306      {
1307         "// the ideal defines a line";
1308      }
1309      MB=SM[2];
1310      intvec rw;
1311      list LL=substpart(MB,ihp,0,rw);
1312      def newR4=LL[1];
1313      setring newR4;
1314      ideal norid=endid;
1315      ideal normap=endphi;
1316      kill endid,endphi;
1317      export norid;
1318      export normap;
1319      result=newR4;
1320      result[size(result)+1]=list(delt,delti);
1321      setring BAS;
1322      return(result);
1323   }
1324
1325//---------------------- The non-trivial cases start -------------------
1326   //the higher dimensional case
1327   //we test first hypersurface, CohenMacaulay and complete intersection
1328
1329   if( ((size(SM[2])+dim(SM[1])) == nvars(basering)) )
1330   {
1331      //the test for complete intersection
1332      attrib(SM[2],"isCohenMacaulay",1);
1333      attrib(SM[2],"isCompleteIntersection",1);
1334      attrib(SM[2],"isEquidimensional",1);
1335      if(y>=1)
1336      {
1337         "// the ideal is a complete intersection";
1338      }
1339   }
1340   if( size(SM[2]) == 1 )
1341   {
1342      attrib(SM[2],"isHypersurface",1);
1343      if(y>=1)
1344      {
1345         "// the ideal is a hypersurface";
1346      }
1347   }
1348
1349   //------------------- compute the singular locus -------------------
1350   // Computation if singular locus is critical
1351   // Notation: J ideal of singular locus or (if given) containing it
1352   // JM = mstd(J) or maxideal(1),maxideal(1)
1353   // JM[1] SB of singular locus, JM[2] minbasis, dimJ = dim(JM[1])
1354   // SM[1] SB of the input ideal i, SM[2] minbasis
1355   // Computation if singular locus is critical, because it determines the
1356   // size of the ring Hom_R(J,J). We only need a test ideal contained in J.
1357
1358   //----------------------- onlySingularAtZero -------------------------
1359   if( attrib(SM[2],"onlySingularAtZero") )
1360   {
1361       JM = maxideal(1),maxideal(1);
1362       attrib(JM[1],"isSB",1);
1363       attrib(JM[2],"isRadical",1);
1364       if( dim(SM[1]) >=2 )
1365       {
1366         attrib(SM[2],"isRegInCodim2",1);
1367       }
1368   }
1369
1370   //-------------------- not onlySingularAtZero -------------------------
1371   if( attrib(SM[2],"onlySingularAtZero") == 0 )
1372   {
1373      //--- the case where an ideal #[1] is given:
1374      if( size(#)>0 )
1375      {
1376         J = #[1],SM[2];
1377         JM = mstd(J);
1378         if( typeof(attrib(#[1],"isRadical"))!="int" )
1379         {
1380            attrib(JM[2],"isRadical",0);
1381         }
1382      }
1383
1384      //--- the case where an ideal #[1] is not given:
1385      if( (size(#)==0) )
1386      {
1387         if(y >=1 )
1388         {
1389            "// singular locus will be computed";
1390         }
1391
1392         J = SM[1],minor(jacob(SM[2]),nvars(basering)-dim(SM[1]),SM[1]);
1393         if( y >=1 )
1394         {
1395            "// SB of singular locus will be computed";
1396         }
1397         JM = mstd(J);
1398      }
1399
1400      int dimJ = dim(JM[1]);
1401      attrib(JM[1],"isSB",1);
1402      if( y>=1 )
1403      {
1404         "// the dimension of the singular locus is";  dimJ ; "";
1405      }
1406
1407      if(dim(JM[1]) <= dim(SM[1])-2)
1408      {
1409         attrib(SM[2],"isRegInCodim2",1);
1410      }
1411
1412      //------------------ the smooth case, RETURN -------------------
1413      if( dimJ == -1 )
1414      {
1415         if(y>=1)
1416         {
1417            "// the ideal is smooth";
1418         }
1419         MB=SM[2];
1420         intvec rw;
1421         list LL=substpart(MB,ihp,0,rw);
1422         def newR3=LL[1];
1423         setring newR3;
1424         ideal norid=endid;
1425         ideal normap=endphi;
1426         kill endid,endphi;
1427         export norid;
1428         export normap;
1429         result=newR3;
1430         result[size(result)+1]=list(delt,delti);
1431         setring BAS;
1432         return(result);
1433      }
1434
1435      //------- extra check for onlySingularAtZero, relatively cheap ----------
1436      //it uses the procedure 'locAtZero' from for testing
1437      //if an ideal is concentrated at 0
1438       if(y>=1)
1439       {
1440         "// extra test for onlySingularAtZero:";
1441       }
[152996]1442       if ( locAtZero(JM[1]) )
1443       {
1444           attrib(SM[2],"onlySingularAtZero",1);
1445           JM = maxideal(1),maxideal(1);
1446           attrib(JM[1],"isSB",1);
1447           attrib(JM[2],"isRadical",1);
1448       }
1449       else
1450       {
1451            attrib(SM[2],"onlySingularAtZero",0);
1452       }
[1598658]1453   }
1454
1455  //displaying the attributes:
1456   if(y>=2)
1457   {
1458      "// the attributes of the ideal are:";
1459      "// isCohenMacaulay:", attrib(SM[2],"isCohenMacaulay");
1460      "// isCompleteIntersection:", attrib(SM[2],"isCompleteIntersection");
1461      "// isHypersurface:", attrib(SM[2],"isHypersurface");
1462      "// isEquidimensional:", attrib(SM[2],"isEquidimensional");
1463      "// isPrim:", attrib(SM[2],"isPrim");
1464      "// isRegInCodim2:", attrib(SM[2],"isRegInCodim2");
1465      "// isIsolatedSingularity:", attrib(SM[2],"isIsolatedSingularity");
1466      "// onlySingularAtZero:", attrib(SM[2],"onlySingularAtZero");
1467      "// isRad:", attrib(SM[2],"isRad");"";
1468   }
1469
1470   //------------- case: CohenMacaulay in codim 2, RETURN ---------------
1471   if( (attrib(SM[2],"isRegInCodim2")==1) &&
1472       (attrib(SM[2],"isCohenMacaulay")==1) )
1473   {
1474      if(y>=1)
1475      {
1476         "// the ideal was CohenMacaulay and regular in codim 2, hence normal";
1477      }
1478      MB=SM[2];
1479      intvec rw;
1480      list LL=substpart(MB,ihp,0,rw);
1481      def newR6=LL[1];
1482      setring newR6;
1483      ideal norid=endid;
1484      ideal normap=endphi;
1485      kill endid,endphi;
1486      export norid;
1487      export normap;
1488      result=newR6;
1489      result[size(result)+1]=list(delt,delti);
1490      setring BAS;
1491      return(result);
1492   }
1493
1494//---------- case: isolated singularity only at 0, RETURN ------------
1495   // In this case things are easier, we can use the maximal ideal as radical
1496   // of the singular locus;
1497   // JM mstd of ideal of singular locus, SM mstd of input ideal
1498
1499   if( attrib(SM[2],"onlySingularAtZero") )
1500   {
1501   //------ check variables for being a non zero-divizor ------
1502   // SL = ideal of vars not contained in ideal SM[1]:
1503
1504      attrib(SM[2],"isIsolatedSingularity",1);
1505      ideal SL = simplify(reduce(maxideal(1),SM[1]),2);
1506      ideal Ann = quotient(SM[2],SL[1]);
1507      ideal qAnn = simplify(reduce(Ann,SM[1]),2);
1508      //NOTE: qAnn=0 if and only if first var (=SL[1]) not in SM is a nzd of R/SM
1509
1510   //------------- We found a non-zerodivisor of R/SM -----------------------
1511   // here the enlarging of the ring via Hom_R(J,J) starts
1512
1513      if( size(qAnn)==0 )
1514      {
1515         if(y>=1)
1516         {
1517            "";
1518            "// the ideal rad(J):"; maxideal(1);
1519            "";
1520         }
1521
1522      // ------------- test for normality, compute Hom_R(J,J) -------------
1523      // Note:
1524      // HomJJ (ideal SBid, ideal id, ideal J, poly p) with
1525      //        SBid = SB of id, J = radical ideal of basering  P with:
1526      //        nonNormal(R) is in V(J), J contains the nonzero divisor p
1527      //        of R = P/id (J = test ideal)
1528      // returns a list l of three objects
1529      // l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
1530      //        s.t. l[1]/endid = Hom_R(J,J) and endphi= map R -> Hom_R(J,J)
1531      // l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
1532      // l[3] : an integer, = dim_K(Hom_R(J,J)/R) if finite, -1 otherwise
1533
1534         list RR;
1535         RR = SM[1],SM[2],maxideal(1),SL[1];
1536         RR = HomJJ(RR,y);
1537         // --------------------- non-normal case ------------------
1538         //RR[2]==0 means that the test for normality is negative
1539         if( RR[2]==0 )
1540         {
1541            def newR=RR[1];
1542            setring newR;
1543            map psi=BAS,endphi;
1544            list JM = psi(JM); //###
1545            ideal J = JM[2];
1546            if ( delt>=0 && RR[3]>=0 )
1547            {
1548               delt = delt+RR[3];
1549            }
1550            else
1551            { delt = -1; }
1552            delti[size(delti)]=delt;
1553
1554            // ---------- recursive call of normalizationPrimes -----------
1555        //normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #)
1556        //ihp = (partial) normalisation map from basering
1557        //#[1] ideal s.t. V(#[1]) contains singular locus of i (test ideal)
1558
1559            if ( y>=1 )
1560            {
1561            "// case: onlySingularAtZero, non-zerodivisor found";
1562            "// contribution of delta in ringextension R -> Hom_R(J,J):"; delt;
1563            }
1564
1565            //intvec atr=getAttrib(endid);
1566            //"//### case: isolated singularity only at 0, recursive";
1567            //"size endid:", size(endid), size(string(endid));
1568            //"interred:";
1569            //endid = interred(endid);
1570            //endid = setAttrib(endid,atr);
1571            //"size endid:", size(endid), size(string(endid));
1572
1573           printlevel=printlevel+1;
1574           list tluser =
1575                normalizationPrimes(endid,psi(ihp),delt,delti);
1576           //list tluser =
1577           //     normalizationPrimes(endid,psi(ihp),delt,delti,J);
1578           //#### ??? improvement: give also the old ideal of sing locus???
1579
1580           printlevel = printlev;             //reset printlevel
1581           setring BAS;
1582           return(tluser);
1583         }
1584
1585         // ------------------ the normal case, RETURN -----------------
1586         // Now RR[2] must be 1, hence the test for normality was positive
1587         MB=SM[2];
1588         //execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
1589         //             +ordstr(basering)+");");
1590         def newR7 = ring(gnirlist);
1591         setring newR7;
1592         ideal norid=fetch(BAS,MB);
1593         ideal normap=fetch(BAS,ihp);
1594         if ( delt>=0 && RR[3]>=0 )
1595         {
1596               delt = delt+RR[3];
1597         }
1598         else
1599         { delt = -1; }
1600         delti[size(delti)]=delt;
1601
1602         intvec atr = getAttrib(norid);
1603
1604         //"//### case: isolated singularity only at 0, final";
1605         //"size norid:", size(norid), size(string(norid));
1606         //"interred:";
1607         //norid = interred(norid);
1608         //norid = setAttrib(norid,atr);
1609         //"size norid:", size(norid), size(string(norid));
1610
1611         export norid;
1612         export normap;
1613         result=newR7;
1614         result[size(result)+1]=list(delt,delti);
1615         setring BAS;
1616         return(result);
1617      }
1618
1619   //------ zerodivisor of R/SM was found, gives a splitting ------------
1620   //Now the case where qAnn!=0, i.e. SL[1] is a zero divisor of R/SM
1621   //and we have found a splitting: id and id1
1622   //id = Ann defines components of R/SM in the complement of V(SL[1])
1623   //id1 defines components of R/SM in the complement of V(id)
1624
1625      else
1626       {
1627          ideal id = Ann;
1628          attrib(id,"isCohenMacaulay",0);
1629          attrib(id,"isPrim",0);
1630          attrib(id,"isIsolatedSingularity",1);
1631          attrib(id,"isRegInCodim2",0);
1632          attrib(id,"isHypersurface",0);
1633          attrib(id,"isCompleteIntersection",0);
1634          attrib(id,"isEquidimensional",0);
1635          attrib(id,"onlySingularAtZero",1);
1636
1637          ideal id1 = quotient(SM[2],Ann);
1638          attrib(id1,"isCohenMacaulay",0);
1639          attrib(id1,"isPrim",0);
1640          attrib(id1,"isIsolatedSingularity",1);
1641          attrib(id1,"isRegInCodim2",0);
1642          attrib(id1,"isHypersurface",0);
1643          attrib(id1,"isCompleteIntersection",0);
1644          attrib(id1,"isEquidimensional",0);
1645          attrib(id1,"onlySingularAtZero",1);
1646
1647          // ---------- recursive call of normalizationPrimes -----------
1648          if ( y>=1 )
1649          {
1650            "// case: onlySingularAtZero, zerodivisor found, splitting:";
1651            "// total delta before splitting:", delt;
1652            "// splitting in two components:";
1653          }
1654
1655          printlevel = printlevel+1;  //to see comments in normalizationPrimes
1656          keepresult1 = normalizationPrimes(id,ihp,0,0);   //1st split factor
1657          keepresult2 = normalizationPrimes(id1,ihp,0,0);  //2nd split factor
1658          printlevel = printlev;                           //reset printlevel
1659
1660          int delt1 = keepresult1[size(keepresult1)][1];
1661          int delt2 = keepresult2[size(keepresult2)][1];
1662          intvec delti1 = keepresult1[size(keepresult1)][2];
1663          intvec delti2 = keepresult2[size(keepresult2)][2];
1664
1665          if( delt>=0 && delt1>=0 && delt2>=0 )
1666          {  ideal idid1=id,id1;
1667             int mul = vdim(std(idid1));
1668             if ( mul>=0 )
1669             {
1670               delt = delt+mul+delt1+delt2;
1671             }
1672             else
1673             {
1674               delt = -1;
1675             }
1676          }
1677         if ( y>=1 )
1678         {
1679           "// delta of first component:", delt1;
1680           "// delta of second componenet:", delt2;
1681           "// intersection multiplicity of both components:", mul;
1682           "// total delta after splitting:", delt;
1683         }
1684
1685          else
1686          {
1687            delt = -1;
1688          }
1689          for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
1690          {
1691             keepresult1=insert(keepresult1,keepresult2[lauf]);
1692          }
1693          keepresult1[size(keepresult1)]=list(delt,delti);
1694
1695          return(keepresult1);
1696       }
1697   }
1698   // Case "onlySingularAtZero" has finished and returned result
1699
1700//-------------- General case, not onlySingularAtZero, RETURN ---------------
1701   //test for non-normality, i.e. if Hom(I,I)<>R
1702   //we can use Hom(I,I) to continue
1703
1704   //------ check variables for being a non zero-divizor ------
1705   // SL = ideal of vars not contained in ideal SM[1]:
1706
1707   ideal SL = simplify(reduce(JM[2],SM[1]),2);
1708   ideal Ann = quotient(SM[2],SL[1]);
1709   ideal qAnn = simplify(reduce(Ann,SM[1]),2);
1710   //NOTE: qAnn=0 <==> first var (=SL[1]) not contained in SM is a nzd of R/SM
1711
1712   //------------- We found a non-zerodivisor of R/SM -----------------------
1713   //SM = mstd of ideal of variety, JM = mstd of ideal of singular locus
1714
1715   if( size(qAnn)==0 )
1716   {
1717      list RR;
1718      list RS;
1719      // ----------------- Computation of the radical -----------------
1720      if(y>=1)
1721      {
1722         "// radical computation of singular locus";
1723      }
1724      J = radical(JM[2]);   //the radical of singular locus
1725      JM = mstd(J);
1726
1727      if(y>=1)
1728      {
1729        "// radical is equal to:";"";  JM[2];
1730        "";
1731      }
1732      // ------------ choose non-zerodivisor of smaller degree ----------
1733      //### evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen ?
1734      if( deg(SL[1]) > deg(J[1]) )
1735      {
1736         Ann=quotient(SM[2],J[1]);
1737         qAnn=simplify(reduce(Ann,SM[1]),2);
1738         if(size(qAnn)==0)
1739         {
1740           SL[1]=J[1];
1741         }
1742      }
1743
1744      // --------------- computation of Hom(rad(J),rad(J)) --------------
1745      RR=SM[1],SM[2],JM[2],SL[1];
1746
1747     if(y>=1)
1748     {
1749        "// compute Hom(rad(J),rad(J))";
1750     }
1751
1752     RS=HomJJ(RR,y);               //most important subprocedure
1753
1754     // ------------------ the normal case, RETURN -----------------
1755     // RS[2]==1 means that the test for normality was positive
1756     if(RS[2]==1)
1757     {
1758         def lastR=RS[1];
1759         setring lastR;
1760         map psi1=BAS,endphi;
1761         ideal norid=endid;
1762         ideal normap=psi1(ihp);
1763         kill endid,endphi;
1764
1765        intvec atr=getAttrib(norid);
1766
1767        //"//### general case: not isolated singularity only at 0, final";
1768        //"size norid:", size(norid), size(string(norid));
1769        //"interred:";
1770        //norid = interred(norid);
1771        //norid = setAttrib(norid,atr);
1772        //"size norid:", size(norid), size(string(norid));
1773
1774         export norid;
1775         export normap;
1776         result=lastR;
1777         if ( y>=1 )
1778         {
1779            "// case: not onlySingularAtZero, last ring Hom_R(J,J) computed";
1780            "// delta before last ring:", delt;
1781         }
1782
1783         if ( delt>=0 && RS[3]>=0 )
1784         {
1785            delt = delt+RS[3];
1786         }
1787         else
1788         { delt = -1; }
1789
1790        // delti = delti,delt;
1791         delti[size(delti)]=delt;
1792
1793         if ( y>=1 )
1794         {
1795           "// delta of last ring:", delt;
1796         }
1797
1798         result[size(result)+1]=list(delt,delti);
1799         setring BAS;
1800         return(result);
1801     }
1802
1803    // ----- the non-normal case, recursive call of normalizationPrimes -------
1804    // RS=HomJJ(RR,y) was computed above, RS[1] contains endid and endphi
1805    // RS[1] = new ring Hom_R(J,J), RS[2]= 0 or 1, RS[2]=contribution to delta
1806    // now RS[2]must be 0, i.e. the test for normality was negative
1807
1808      int n = nvars(basering);
1809      ideal MJ = JM[2];
1810
1811      def newR=RS[1];
1812      setring newR;
1813      map psi=BAS,endphi;
1814      if ( y>=1 )
1815      {
1816        "// case: not onlySingularAtZero, compute new ring = Hom_R(J,J)";
1817        "// delta of old ring:", delt;
1818      }
1819      if ( delt>=0 && RS[3]>=0 )
1820      {
1821         delt = delt+RS[3];
1822      }
1823      else
1824      { delt = -1; }
1825      if ( y>=1 )
1826      {
1827        "// delta of new ring:", delt;
1828      }
1829
1830      delti[size(delti)]=delt;
1831      intvec atr=getAttrib(endid);
1832
1833      //"//### general case: not isolated singularity only at 0, recursive";
1834      //"size endid:", size(endid), size(string(endid));
1835      //"interred:";
1836      //endid = interred(endid);
1837      //endid = setAttrib(endid,atr);
1838      //"size endid:", size(endid), size(string(endid));
1839
1840      printlevel = printlevel+1;
1841      list tluser=
1842          normalizationPrimes(endid,psi(ihp),delt,delti,psi(MJ));
1843      printlevel = printlev;                //reset printlevel
1844      setring BAS;
1845      return(tluser);
1846   }
1847
1848   //---- A whole singular component was found, RETURN -----
1849   if( Ann == 1)
1850   {
1851      "// Input appeared not to be a radical ideal!";
1852      "// A (everywhere singular) component with ideal";
1853      "// equal to its Jacobian ideal was found";
1854      "// Procedure will stop with partial result computed so far";"";
1855
1856         MB=SM[2];
1857         intvec rw;
1858         list LL=substpart(MB,ihp,0,rw);
1859         def newR6=LL[1];
1860         setring newR6;
1861         ideal norid=endid;
1862         ideal normap=endphi;
1863         kill endid,endphi;
1864         export norid;
1865         export normap;
1866         result=newR6;
1867         result[size(result)+1]=lst(delt,delti);
1868         setring BAS;
1869         return(result);
1870   }
1871
1872   //------ zerodivisor of R/SM was found, gives a splitting ------------
1873   //Now the case where qAnn!=0, i.e. SL[1] is a zero divisor of R/SM
1874   //and we have found a splitting: new1 and new2
1875   //id = Ann defines components of R/SM in the complement of V(SL[1])
1876   //id1 defines components of R/SM in the complement of V(id)
1877
1878   else
1879   {
1880      if(y>=1)
1881      {
1882         "// zero-divisor found";
1883      }
1884      int equi = attrib(SM[2],"isEquidimensional");
1885      int oSAZ = attrib(SM[2],"onlySingularAtZero");
1886      int isIs = attrib(SM[2],"isIsolatedSingularity");
1887
1888      ideal new1 = Ann;
1889      ideal new2 = quotient(SM[2],Ann);
1890      //ideal new2=SL[1],SM[2];
1891
1892      //execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
1893      //                +ordstr(basering)+");");
1894      def newR1 = ring(gnirlist);
1895      setring newR1;
1896
1897      ideal vid = fetch(BAS,new1);
1898      ideal ihp = fetch(BAS,ihp);
1899      attrib(vid,"isCohenMacaulay",0);
1900      attrib(vid,"isPrim",0);
1901      attrib(vid,"isIsolatedSingularity",isIs);
1902      attrib(vid,"isRegInCodim2",0);
1903      attrib(vid,"onlySingularAtZero",oSAZ);
1904      attrib(vid,"isEquidimensional",equi);
1905      attrib(vid,"isHypersurface",0);
1906      attrib(vid,"isCompleteIntersection",0);
1907
1908      // ---------- recursive call of normalizationPrimes -----------
1909      if ( y>=1 )
1910      {
1911        "// total delta before splitting:", delt;
1912        "// splitting in two components:";
1913      }
1914      printlevel = printlevel+1;
1915      keepresult1 =
1916                  normalizationPrimes(vid,ihp,0,0);  //1st split factor
1917
1918      list delta1 = keepresult1[size(keepresult1)];
1919
1920      setring BAS;
1921      //execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
1922      //                +ordstr(basering)+");");
1923      def newR2 = ring(gnirlist);
1924      setring newR2;
1925
1926      ideal vid = fetch(BAS,new2);
1927      ideal ihp = fetch(BAS,ihp);
1928      attrib(vid,"isCohenMacaulay",0);
1929      attrib(vid,"isPrim",0);
1930      attrib(vid,"isIsolatedSingularity",isIs);
1931      attrib(vid,"isRegInCodim2",0);
1932      attrib(vid,"isEquidimensional",equi);
1933      attrib(vid,"isHypersurface",0);
1934      attrib(vid,"isCompleteIntersection",0);
1935      attrib(vid,"onlySingularAtZero",oSAZ);
1936
1937      keepresult2 =
1938                    normalizationPrimes(vid,ihp,0,0);
1939      list delta2 = keepresult2[size(keepresult2)];   //2nd split factor
1940      printlevel = printlev;                          //reset printlevel
1941
1942      setring BAS;
1943
1944      //compute intersection multiplicity of both components:
1945      new1 = new1,new2;
1946      int mul=vdim(std(new1));
1947
1948     // ----- normalizationPrimes finished, add up results, RETURN --------
1949      for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
1950      {
1951         keepresult1 = insert(keepresult1,keepresult2[lauf]);
1952      }
1953      if ( delt >=0 && delta1[1] >=0 && delta2[1] >=0 && mul >=0 )
1954      {
1955         delt = delt+mul+delta1[1]+delta2[1];
1956      }
1957      else
1958      {  delt = -1; }
1959      delti = -2;
1960
1961      if ( y>=1 )
1962      {
1963        "// zero divisor produced a splitting into two components";
1964        "// delta of first component:", delta1;
1965        "// delta of second componenet:", delta2;
1966        "// intersection multiplicity of both components:", mul;
1967        "// total delta after splitting:", delt;
1968      }
1969      keepresult1[size(keepresult1)] = list(delt,delti);
1970      return(keepresult1);
1971   }
1972}
1973example
1974{ "EXAMPLE:";echo = 2;
1975   // Huneke
1976   ring qr=31991,(a,b,c,d,e),dp;
1977   ideal i=
1978   5abcde-a5-b5-c5-d5-e5,
1979   ab3c+bc3d+a3be+cd3e+ade3,
1980   a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1981   abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1982   ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
1983   a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1984   a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
1985   b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1986
1987   list pr=normalizationPrimes(i);
1988   def r1 = pr[1];
1989   setring r1;
1990   norid;
1991   normap;
1992}
1993
1994///////////////////////////////////////////////////////////////////////////////
1995static proc substpart(ideal endid, ideal endphi, int homo, intvec rw)
1996
1997"//Repeated application of elimpart to endid, until no variables can be
1998//directy substituded. homo=1 if input is homogeneous, rw contains
1999//original weights, endphi (partial) normalization map";
2000
2001//NOTE concerning iteration of maps: Let phi: x->f(y,z), y->g(x,z) then
2002//phi: x+y+z->f(y,z)+g(x,z)+z, phi(phi):x+y+z->f(g(x,z),z)+g(f(y,z),z)+z
2003//and so on: none of the x or y will be eliminated
2004//Now subst: first x and then y: x+y+z->f(g(x,z),z)+g(x,z)+z eliminates y
2005//further subst replaces x by y, makes no sense (objects more compicated).
2006//Subst first y and then x eliminates x
2007//In our situation we have triangular form: x->f(y,z), y->g(z).
2008//phi: x+y+z->f(y,z)+g(z)+z, phi(phi):x+y+z->f(g(z),z)+g(z)+z eliminates x,y
2009//subst x,y: x+y+z->f(g(z),z)+g(z)+z, eliminates x,y
2010//subst y,x: x+y+z->f(y,z)+g(z)+z eliminates only x
2011//HENCE: substitute vars depending on most other vars first
2012//However, if the sytem xi-fi is reduced then xi does not appear in any of the
2013//fj and hence the order does'nt matter when substitutinp xi by fi
2014
2015{
2016   def newRing = basering;
2017   int ii,jj;
2018   map phi = newRing,maxideal(1);    //identity map
2019   list Le = elimpart(endid);
2020   //this proc and the next loop try to substitute as many variables as
2021   //possible indices of substituted variables
2022
2023   int q = size(Le[2]);    //q vars, stored in Le[2], have been substitutet
2024   intvec rw1 = 0;         //will become indices of substituted variables
2025   rw1[nvars(basering)] = 0;
2026   rw1 = rw1+1;            //rw1=1,..,1 (as many 1 as nvars(basering))
2027
2028   while( size(Le[2]) != 0 )
2029   {
2030      endid = Le[1];
2031      if ( defined(ps) )
2032      { kill ps; }
2033      map ps = newRing,Le[5];
2034      phi = ps(phi);
2035      for(ii=1;ii<=size(Le[2]);ii++)
2036      {
2037         phi=phi(phi);
2038      }
2039      //eingefuegt wegen x2-y2z2+z3
2040
2041      for( ii=1; ii<=size(rw1); ii++ )
2042      {
2043         if( Le[4][ii]==0 )        //ii = index of var which was substituted
2044         {
2045            rw1[ii]=0;             //substituted vars have entry 0 in rw1
2046         }
2047      }
2048      Le=elimpart(endid);          //repeated application of elimpart
2049      q = q + size(Le[2]);
2050   }
2051   endphi = phi(endphi);
2052//---------- return -----------------------------------------------------------
2053// first the trivial case, where all variable have been eliminated
2054   if( nvars(newRing) == q )
2055   {
2056     ring lastRing = char(basering),T(1),dp;
2057     ideal endid = T(1);
2058     ideal endphi = T(1);
2059     for(ii=2; ii<=q; ii++ )
2060     {
2061        endphi[ii] = 0;
2062     }
2063     export(endid,endphi);
2064     list L = lastRing;
2065     setring newRing;
2066     return(L);
2067   }
2068
2069// in the homogeneous case put weights for the remaining vars correctly, i.e.
2070// delete from rw those weights for which the corresponding entry of rw1 is 0
2071
2072   if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 )
2073   {
2074      jj=1;
2075      for( ii=2; ii<size(rw1); ii++)
2076      {
2077         jj++;
2078         if( rw1[ii]==0 )
2079         {
2080            rw=rw[1..jj-1],rw[jj+1..size(rw)];
2081            jj=jj-1;
2082         }
2083      }
2084      if( rw1[1]==0 ) { rw=rw[2..size(rw)]; }
2085      if( rw1[size(rw1)]==0 ){ rw=rw[1..size(rw)-1]; }
2086
2087      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),(a(rw),dp);
2088   }
2089   else
2090   {
2091      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),dp;
2092   }
2093   ideal lastmap;
2094   jj = 1;
2095
2096   for(ii=1; ii<=size(rw1); ii++ )
2097   {
2098      if ( rw1[ii]==1 ) { lastmap[ii] = T(jj); jj=jj+1; }
2099      if ( rw1[ii]==0 ) { lastmap[ii] = 0; }
2100   }
2101   map phi1 = newRing,lastmap;
2102   ideal endid  = phi1(endid);      //### bottelneck
2103   ideal endphi = phi1(endphi);
2104
2105/*
2106Versuch: subst statt phi
2107   for(ii=1; ii<=size(rw1); ii++ )
2108   {
2109      if ( rw1[ii]==1 ) { endid = subst(endid,var(ii),T(jj)); }
2110      if ( rw1[ii]==0 ) { endid = subst(endid,var(ii),0); }
2111   }
2112*/
2113   export(endid);
2114   export(endphi);
2115   list L = lastRing;
2116   setring newRing;
2117   return(L);
2118}
2119///////////////////////////////////////////////////////////////////////////////
2120static proc deltaP(ideal I)
2121{
2122   def R=basering;
2123   int c,d,i;
2124   int n=nvars(R);
2125   list nor;
2126   if(size(I)>1){ERROR("no hypersurface");}
2127   ideal J=std(slocus(I));
2128   if(dim(J)<=0){return(0);}
2129   poly h;
2130   d=1;
2131   while((d)&&(i<n))
2132   {
2133      i++;
2134      h=var(i);
2135      d=dim(std(J+ideal(h)));
2136   }
2137   i=0;
2138   while(d)
2139   {
2140      i++;
2141      if(i>10){ERROR("delta not found, please inform the authors")};
2142      h=randomLast(100)[n];
2143      d=dim(std(J+ideal(h)));
2144   }
2145   I=I,h-1;
2146   if(char(R)<=19)
2147   {
2148      nor=normalP(I);
2149   }
2150   else
2151   {
2152      nor=normal(I);
2153   }
2154   return(nor[2][2]);
2155}
2156
2157proc genus(ideal I,list #)
2158"USAGE:   genus(i) or genus(i,1); I a 1-dimensional ideal
2159RETURN:  an integer, the geometric genus p_g = p_a - delta of the projective
2160         curve defined by i, where p_a is the arithmetic genus.
2161NOTE:    delta is the sum of all local delta-invariants of the singularities,
2162         i.e. dim(R'/R), R' the normalization of the local ring R of the
2163         singularity. @*
2164         genus(i,1) uses the normalization to compute delta. Usually genus(i,1)
2165         is slower than genus(i) but sometimes not.
2166EXAMPLE: example genus; shows an example
2167"
2168{
2169   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
2170
2171   def R0=basering;
2172   if(char(basering)>0)
2173   {
[1e1ec4]2174     def R1=changeord(list(list("dp",1:nvars(basering))));
[1598658]2175     setring R1;
2176     ideal I=imap(R0,I);
2177     I=radical(I);
2178     I=std(I);
2179     if(dim(I)!=1)
2180     {
2181       if(((homog(I))&&(dim(I)!=2))||(!homog(I)))
2182       {
2183         ERROR("This is not a curve");
2184       }
2185     }
2186     if(homog(I)&&(dim(I)==2))
2187     {
2188       def S=R0;
2189       setring S;
2190       ideal J=I;
2191     }
2192     else
2193     {
2194       def S=changevar(varstr(R0)+",@t");
2195       setring S;
2196       ideal J=imap(R1,I);
2197       J=homog(J,@t);
2198       J=std(J);
2199     }
2200     int pa=1-hilbPoly(J)[1];
2201     pa=pa-deltaP(J);
2202     setring R0;
2203     return(pa);
2204   }
2205   I=std(I);
2206   if(dim(I)!=1)
2207   {
2208      if(((homog(I))&&(dim(I)!=2))||(!homog(I)))
2209      {
2210        // ERROR("This is not a curve");
2211        if(w==1){"** WARNING: Input does not define a curve **"; "";}
2212      }
2213   }
2214   list L=elimpart(I);
2215   if(size(L[2])!=0)
2216   {
2217      map psi=R0,L[5];
2218      I=std(psi(I));
2219   }
2220   if(size(I)==0)
2221   {
2222      return(0);
2223   }
2224   list N=findvars(I,0);
2225   if(size(N[1])==1)
2226   {
2227
2228      poly p=I[1];
2229     // if(deg(squarefree(p))<deg(p)){ERROR("Curve is not reduced");}
2230      return(-deg(p)+1);
2231   }
2232   if(size(N[1]) < nvars(R0))
2233   {
2234      string newvar=string(N[1]);
2235      execute("ring R=("+charstr(R0)+"),("+newvar+"),dp;");
2236      ideal I =imap(R0,I);
2237      attrib(I,"isSB",1);
2238   }
2239   else
2240   {
2241      def R=basering;
2242   }
2243   if(dim(I)==2)
2244   {
2245      def newR=basering;
2246   }
2247   else
2248   {
2249      if(dim(I)==0)
2250      {
2251         execute("ring Rhelp=("+charstr(R0)+"),(@s,@t),dp;");
2252      }
2253      else
2254      {
2255         execute("ring Rhelp=("+charstr(R0)+"),(@s),dp;");
2256      }
2257      def newR=R+Rhelp;
2258      setring newR;
2259      ideal I=imap(R,I);
2260      I=homog(I,@s);
2261      attrib(I,"isSB",1);
2262   }
2263
2264   if((nvars(basering)<=3)&&(size(I)>1))
2265   {
2266       ERROR("This is not equidimensional");
2267   }
2268
2269   intvec hp=hilbPoly(I);
2270   int p_a=1-hp[1];
2271   int d=hp[2];
2272
2273   if(w>=1)
2274   {
2275      "";"The ideal of the projective curve:";"";I;"";
2276      "The coefficients of the Hilbert polynomial";hp;
2277      "arithmetic genus:";p_a;
2278      "degree:";d;"";
2279   }
2280
2281   intvec v = hilb(I,1);
2282   int i,o;
2283   if(nvars(basering)>3)
2284   {
2285      map phi=newR,maxideal(1);
2286      int de;
2287      ideal K,L1;
2288      matrix M;
2289      poly m=var(4);
2290      poly he;
2291      for(i=5;i<=nvars(basering);i++){m=m*var(i);}
2292      K=eliminate(I,m,v);
2293      if(size(K)==1){de=deg(K[1]);}
2294      m=var(1);
2295      for(i=2;i<=nvars(basering)-3;i++){m=m*var(i);}
2296      i=0;
2297      while(d!=de)
2298      {
2299         o=1;
2300         i++;
2301         K=phi(I);
2302         K=eliminate(K,m,v);
2303         if(size(K)==1){de=deg(K[1]);}
2304         if((i==5)&&(d!=de))
2305         {
2306            K=reduce(equidimMax(I),I);
2307            if(size(K)!=0){ERROR("This is not equidimensional");}
2308         }
2309         if(i==10)
2310         {
2311            J;K;
2312            ERROR("genus: did not find a good projection for to
2313                           the plain");
2314         }
2315         if(i<5)
2316         {
2317            M=sparsetriag(nvars(newR),nvars(newR),80-5*i,i);
2318         }
2319         else
2320         {
2321            if(i<8)
2322            {
2323               M=transpose(sparsetriag(nvars(newR),nvars(newR),80-5*i,i));
2324            }
2325            else
2326            {
2327               he=0;
2328               while(he==0)
2329               {
2330                  M=randommat(nvars(newR),nvars(newR),ideal(1),20);
2331                  he=det(M);
2332               }
2333            }
2334         }
2335         L1=M*transpose(maxideal(1));
2336         phi=newR,L1;
2337      }
2338      I=K;
2339   }
2340   poly p=I[1];
2341
2342   execute("ring S=("+charstr(R)+"),(x,y,t),dp;");
2343   ideal L=maxideal(1);
2344   execute("ring C=("+charstr(R)+"),(x,y),ds;");
2345   ideal I;
2346   execute("ring A=("+charstr(R)+"),(x,t),dp;");
2347   map phi=S,1,x,t;
2348   map psi=S,x,1,t;
2349   poly g,h;
2350   ideal I,I1;
2351   execute("ring B=("+charstr(R)+"),(x,t),ds;");
2352
2353   setring S;
2354   if(o)
2355   {
2356     for(i=1;i<=nvars(newR)-3;i++){L[i]=0;}
2357     L=L,maxideal(1);
2358   }
2359   map sigma=newR,L;
2360   poly F=sigma(p);
2361   if(w>=1){"the projected curve:";"";F;"";}
2362
2363   kill newR;
2364
[ed7a55c]2365   int genus=(d-1)*(d-2) div 2;
[1598658]2366   if(w>=1){"the arithmetic genus of the plane curve:";genus;pause();}
2367
2368   int delt,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,l,nsing,
2369       tauloc,tausing,k,rat,nbranchinf,nbranch,nodes,cuspsinf,nodesinf;
2370   list inv;
2371
2372   if(w>=1)
2373     {"";"analyse the singularities at oo";"";"singular locus at (1,x,0):";"";}
2374   setring A;
2375   g=phi(F);
2376   h=psi(F);
2377   I=g,jacob(g),var(2);
2378   I=std(I);
2379   if(deg(I[1])>0)
2380   {
2381      list qr=minAssGTZ(I);
2382      if(w>=1){qr;"";}
2383
2384      for(k=1;k<=size(qr);k++)
2385      {
2386         if(w>=1){ nsing=nsing+vdim(std(qr[k]));}
2387         inv=deltaLoc(g,qr[k]);
2388         deltainf=deltainf+inv[1];
2389         tauinf=tauinf+inv[2];
2390         l=vdim(std(qr[k]));
2391         if(inv[2]==l){nodesinf=nodesinf+l;}
2392         if(inv[2]==2*l){cuspsinf=cuspsinf+l;}
2393         nbranchinf=nbranchinf+inv[3];
2394      }
2395   }
2396   else
2397   {
2398     if(w>=1){"            the curve is smooth at (1,x,0)";"";}
2399   }
2400   if(w>=1){"singular locus at (0,1,0):";"";}
2401   inv=deltaLoc(h,maxideal(1));
2402   if((w>=1)&&(inv[2]!=0)){ nsing++;}
2403   deltainf=deltainf+inv[1];
2404   tauinf=tauinf+inv[2];
2405   if(inv[2]==1){nodesinf++;}
2406   if(inv[2]==2){cuspsinf++;}
2407
2408   if((w>=1)&&(inv[2]==0)){"            the curve is smooth at (0,1,0)";"";}
2409   if(inv[2]>0){nbranchinf=nbranchinf+inv[3];}
2410
2411   if(w>=1)
2412   {
2413      if(tauinf==0)
2414      {
2415        "            the curve is smooth at oo";"";
2416      }
2417      else
2418      {
2419         "number of singularities at oo:";nsing;
2420         "nodes at oo:";nodesinf;
2421         "cusps at oo:";cuspsinf;
2422         "branches at oo:";nbranchinf;
2423         "Tjurina number at oo:";tauinf;
2424         "delta at oo:";deltainf;
2425         "Milnor number at oo:";2*deltainf-nbranchinf+nsing;
2426         pause();
2427      }
2428      "singularities at (x,y,1):";"";
2429   }
2430   execute("ring newR=("+charstr(R)+"),(x,y),dp;");
2431   //the singularities at the affine part
2432   map sigma=S,var(1),var(2),1;
2433   ideal I=sigma(F);
2434
2435   if(size(#)!=0)
2436   {                              //uses the normalization to compute delta
[e1cda9]2437      list nor=normal(I,"withDelta");
[1598658]2438      delt=nor[size(nor)][2];
2439      genus=genus-delt-deltainf;
2440      setring R0;
2441      return(genus);
2442   }
2443
2444   ideal I1=jacob(I);
2445   matrix Hess[2][2]=jacob(I1);
2446   ideal ID=I+I1+ideal(det(Hess));//singular locus of I+I1
2447
2448   ideal radID=std(radical(ID));//the non-nodal locus
2449   if(w>=1){"the non-nodal locus:";"";radID;pause();"";}
2450   if(deg(radID[1])==0)
2451   {
2452     ideal IDsing=1;
2453   }
2454   else
2455   {
2456     ideal IDsing=minor(jacob(ID),2)+radID;//singular locus of ID
2457   }
2458
2459   iglob=vdim(std(IDsing));
2460
2461   if(iglob!=0)//computation of the radical of IDsing
2462   {
2463      ideal radIDsing=reduce(IDsing,radID);
2464      if(size(radIDsing)==0)
2465      {
2466         radIDsing=radID;
2467         attrib(radIDsing,"isSB",1);
2468      }
2469      else
2470      {
2471         radIDsing=std(radical(IDsing));
2472      }
2473      iglob=vdim(radIDsing);
2474      if((w>=1)&&(iglob))
2475          {"the non-nodal-cuspidal locus:";radIDsing;pause();"";}
2476   }
2477   cusps=vdim(radID)-iglob;
2478   nsing=nsing+cusps;
2479
2480   if(iglob==0)
2481   {
2482      if(w>=1){"             there are only cusps and nodes";"";}
2483      tau=vdim(std(I+jacob(I)));
2484      tauinf=tauinf+tau;
2485      nodes=tau-2*cusps;
2486      delt=nodes+cusps;
2487      nbranch=2*tau-3*cusps;
2488      nsing=nsing+nodes;
2489   }
2490   else
2491   {
2492       if(w>=1){"the non-nodal-cuspidal singularities";"";}
2493       setring C;
[03b9ee]2494       ideal I1=imap(newR,radIDsing);
[1598658]2495       iloc=vdim(std(I1));
2496       if(iglob==iloc)
2497       {
2498          if(w>=1){"only cusps and nodes outside (0,0,1)";}
2499          setring newR;
2500          tau=vdim(std(I+jacob(I)));
2501          tauinf=tauinf+tau;
2502          inv=deltaLoc(I[1],maxideal(1));
2503          delt=inv[1];
2504          tauloc=inv[2];
2505          nodes=tau-tauloc-2*cusps;
2506          nsing=nsing+nodes;
[f2e8212]2507          if (inv[2]!=0) { nsing++; }
[1598658]2508          nbranch=inv[3]+ 2*nodes+cusps;
2509          delt=delt+nodes+cusps;
2510          if((w>=1)&&(inv[2]==0)){"smooth at (0,0,1)";}
2511        }
2512        else
2513        {
2514           setring newR;
[03b9ee]2515           list pr=minAssGTZ(radIDsing);
[1598658]2516           if(w>=1){pr;}
2517
2518           for(k=1;k<=size(pr);k++)
2519           {
2520              if(w>=1){nsing=nsing+vdim(std(pr[k]));}
2521              inv=deltaLoc(I[1],pr[k]);
2522              delt=delt+inv[1];
2523              tausing=tausing+inv[2];
2524              nbranch=nbranch+inv[3];
2525           }
2526           tau=vdim(std(I+jacob(I)));
2527           tauinf=tauinf+tau;
2528           nodes=tau-tausing-2*cusps;
2529           nsing=nsing+nodes;
2530           delt=delt+nodes+cusps;
2531           nbranch=nbranch+2*nodes+cusps;
2532        }
2533   }
2534   genus=genus-delt-deltainf;
2535   if(w>=1)
2536   {
2537      "The projected plane curve has locally:";"";
2538      "singularities:";nsing;
2539      "branches:";nbranch+nbranchinf;
2540      "nodes:"; nodes+nodesinf;
2541      "cusps:";cusps+cuspsinf;
2542      "Tjurina number:";tauinf;
2543      "Milnor number:";2*(delt+deltainf)-nbranch-nbranchinf+nsing;
2544      "delta of the projected curve:";delt+deltainf;
2545      "delta of the curve:";p_a-genus;
2546      "genus:";genus;
2547      "====================================================";
2548      "";
2549   }
2550   setring R0;
2551   return(genus);
2552}
2553example
2554{ "EXAMPLE:"; echo = 2;
2555   ring r=0,(x,y),dp;
2556   ideal i=y^9 - x^2*(x - 1)^9;
2557   genus(i);
2558   ring r7=7,(x,y),dp;
2559   ideal i=y^9 - x^2*(x - 1)^9;
2560   genus(i);
2561}
2562
2563///////////////////////////////////////////////////////////////////////////////
2564proc deltaLoc(poly f,ideal singL)
2565"USAGE:  deltaLoc(f,J);  f poly, J ideal
2566ASSUME: f is reduced bivariate polynomial; basering has exactly two variables;
2567        J is irreducible prime component of the singular locus of f (e.g., one
2568        entry of the output of @code{minAssGTZ(I);}, I = <f,jacob(f)>).
2569RETURN:  list L:
2570@texinfo
2571@table @asis
2572@item @code{L[1]}; int:
2573         the sum of (local) delta invariants of f at the (conjugated) singular
2574         points given by J.
2575@item @code{L[2]}; int:
2576         the sum of (local) Tjurina numbers of f at the (conjugated) singular
2577         points given by J.
2578@item @code{L[3]}; int:
2579         the sum of (local) number of branches of f at the (conjugated)
2580         singular points given by J.
2581@end table
2582@end texinfo
2583NOTE:    procedure makes use of @code{execute}; increasing printlevel displays
2584         more comments (default: printlevel=0).
2585SEE ALSO: delta, tjurina
2586KEYWORDS: delta invariant; Tjurina number
2587EXAMPLE: example deltaLoc;  shows an example
2588"
2589{
[651cce]2590   intvec save_opt=option(get);
[1598658]2591   option(redSB);
2592   def R=basering;
2593   execute("ring S=("+charstr(R)+"),(x,y),lp;");
2594   map phi=R,x,y;
2595   ideal singL=phi(singL);
2596   singL=simplify(std(singL),1);
2597   attrib(singL,"isSB",1);
2598   int d=vdim(singL);
2599   poly f=phi(f);
2600   int i;
2601   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
2602   if(d==1)
2603   {
2604      map alpha=S,var(1)-singL[2][2],var(2)-singL[1][2];
2605      f=alpha(f);
2606      execute("ring C=("+charstr(S)+"),("+varstr(S)+"),ds;");
2607      poly f=imap(S,f);
2608      ideal singL=imap(S,singL);
2609      if((w>=1)&&(ord(f)>=2))
2610      {
2611        "local analysis of the singularities";"";
2612        basering;
2613        singL;
2614        f;
2615        pause();
2616      }
2617   }
2618   else
2619   {
2620      poly p;
2621      poly c;
2622      map psi;
2623      number co;
2624
2625      while((deg(lead(singL[1]))>1)&&(deg(lead(singL[2]))>1))
2626      {
2627         psi=S,x,y+random(-100,100)*x;
2628         singL=psi(singL);
2629         singL=std(singL);
2630          f=psi(f);
2631      }
2632
2633      if(deg(lead(singL[2]))==1)
2634      {
2635         p=singL[1];
2636         c=singL[2]-lead(singL[2]);
2637         co=leadcoef(singL[2]);
2638      }
2639      if(deg(lead(singL[1]))==1)
2640      {
2641         psi=S,y,x;
2642         f=psi(f);
2643         singL=psi(singL);
2644         p=singL[2];
2645         c=singL[1]-lead(singL[1]);
2646         co=leadcoef(singL[1]);
2647      }
2648
2649      execute("ring B=("+charstr(S)+"),a,dp;");
2650      map beta=S,a,a;
2651      poly p=beta(p);
2652
2653      execute("ring C=("+charstr(S)+",a),("+varstr(S)+"),ds;");
2654      number p=number(imap(B,p));
2655
2656      minpoly=p;
2657      map iota=S,a,a;
2658      number c=number(iota(c));
2659      number co=iota(co);
2660
2661      map alpha=S,x-c/co,y+a;
2662      poly f=alpha(f);
2663      f=cleardenom(f);
2664      if((w>=1)&&(ord(f)>=2))
2665      {
2666        "local analysis of the singularities";"";
2667        basering;
2668        alpha;
2669        f;
2670        pause();
2671        "";
2672      }
2673   }
2674   option(noredSB);
2675   ideal fstd=std(ideal(f)+jacob(f));
2676   poly hc=highcorner(fstd);
2677   int tau=vdim(fstd);
2678   int o=ord(f);
2679   int delt,nb;
2680
2681   if(tau==0)                 //smooth case
2682   {
2683      setring R;
[651cce]2684      option(set,save_opt);
[1598658]2685      return(list(0,0,1));
2686   }
2687   if((char(basering)>=181)||(char(basering)==0))
2688   {
2689      if(o==2)                //A_k-singularity
2690      {
2691        if(w>=1){"A_k-singularity";"";}
2692         setring R;
[ed7a55c]2693         delt=(tau+1) div 2;
[651cce]2694         option(set,save_opt);
[1598658]2695         return(list(d*delt,d*tau,d*(2*delt-tau+1)));
2696      }
2697      if((lead(f)==var(1)*var(2)^2)||(lead(f)==var(1)^2*var(2)))
2698      {
2699        if(w>=1){"D_k- singularity";"";}
2700
2701         setring R;
[ed7a55c]2702         delt=(tau+2) div 2;
[651cce]2703         option(set,save_opt);
[1598658]2704         return(list(d*delt,d*tau,d*(2*delt-tau+1)));
2705      }
2706
2707      int mu=vdim(std(jacob(f)));
[507427]2708
[1598658]2709      poly g=f+var(1)^mu+var(2)^mu;  //to obtain a convenient Newton-polygon
2710
2711      list NP=newtonpoly(g);
2712      if(w>=1){"Newton-Polygon:";NP;"";}
2713      int s=size(NP);
2714
2715      if(is_NND(f,mu,NP))
2716      { // the Newton-polygon is non-degenerate
2717        // compute nb, the number of branches
2718        for(i=1;i<=s-1;i++)
2719        {
2720          nb=nb+gcd(NP[i][2]-NP[i+1][2],NP[i][1]-NP[i+1][1]);
2721        }
2722        if(w>=1){"Newton-Polygon is non-degenerated";"";}
[651cce]2723        option(set,save_opt);
[ed7a55c]2724        return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
[1598658]2725      }
2726
2727      if(w>=1){"Newton-Polygon is degenerated";"";}
2728
2729      // the following can certainly be made more efficient when replacing
2730      // 'hnexpansion' (used only for computing number of branches) by
2731      // successive blowing-up + test if Newton polygon degenerate:
2732      if(s>2)    //  splitting of f
2733      {
2734         if(w>=1){"Newton polygon can be used for splitting";"";}
2735         intvec v=NP[1][2]-NP[2][2],NP[2][1];
2736         int de=w_deg(g,v);
[507427]2737         //int st=w_deg(hc,v)+v[1]+v[2];
2738         int st=w_deg(var(1)^NP[size(NP)][1],v)+1;
[1598658]2739         poly f1=var(2)^NP[2][2];
2740         poly f2=jet(g,de,v)/var(2)^NP[2][2];
2741         poly h=g-f1*f2;
2742         de=w_deg(h,v);
2743         poly k;
2744         ideal wi=var(2)^NP[2][2],f2;
2745         matrix li;
2746         while(de<st)
2747         {
2748           k=jet(h,de,v);
2749           li=lift(wi,k);
2750           f1=f1+li[2,1];
2751           f2=f2+li[1,1];
2752           h=g-f1*f2;
2753           de=w_deg(h,v);
2754         }
2755         nb=deltaLoc(f1,maxideal(1))[3]+deltaLoc(f2,maxideal(1))[3];
[507427]2756
[1598658]2757         setring R;
[651cce]2758         option(set,save_opt);
[ed7a55c]2759         return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
[1598658]2760      }
2761
2762      f=jet(f,deg(hc)+2);
2763      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
2764      ideal fac=factorize(f,1);
2765      if(size(fac)>1)
2766      {
2767         nb=0;
2768         for(i=1;i<=size(fac);i++)
2769         {
2770            nb=nb+deltaLoc(fac[i],maxideal(1))[3];
2771         }
2772         setring R;
[651cce]2773         option(set,save_opt);
[ed7a55c]2774         return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
[1598658]2775      }
2776      list HNEXP=hnexpansion(f);
2777      if (typeof(HNEXP[1])=="ring") {
2778        def altring = basering;
2779        def HNEring = HNEXP[1]; setring HNEring;
2780        nb=size(hne);
2781        setring R;
2782        kill HNEring;
2783      }
2784      else
2785      {
2786        nb=size(HNEXP);
2787      }
[651cce]2788      option(set,save_opt);
[ed7a55c]2789      return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
[1598658]2790   }
2791   else             //the case of small characteristic
2792   {
2793      f=jet(f,deg(hc)+2);
2794      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
2795      delt=delta(f);
[651cce]2796      option(set,save_opt);
[1598658]2797      return(list(d*delt,d*tau,d));
2798   }
[651cce]2799   option(set,save_opt);
[1598658]2800}
2801example
2802{ "EXAMPLE:"; echo = 2;
2803  ring r=0,(x,y),dp;
2804  poly f=(x2+y^2-1)^3 +27x2y2;
2805  ideal I=f,jacob(f);
2806  I=std(I);
2807  list qr=minAssGTZ(I);
2808  size(qr);
2809  // each component of the singular locus either describes a cusp or a pair
2810  // of conjugated nodes:
2811  deltaLoc(f,qr[1]);
2812  deltaLoc(f,qr[2]);
2813  deltaLoc(f,qr[3]);
2814  deltaLoc(f,qr[4]);
2815  deltaLoc(f,qr[5]);
2816  deltaLoc(f,qr[6]);
2817}
2818///////////////////////////////////////////////////////////////////////////////
[2ab830]2819// compute the weighted degree of p;
2820// this code is an exact copy of the proc in paraplanecurves.lib
2821// (since we do not want to make it non-static)
[1598658]2822static proc w_deg(poly p, intvec v)
2823{
2824   if(p==0){return(-1);}
2825   int d=0;
2826   while(jet(p,d,v)==0){d++;}
2827   d=(transpose(leadexp(jet(p,d,v)))*v)[1];
2828   return(d);
2829}
2830
2831//proc hilbPoly(ideal J)
2832//{
2833//   poly hp;
2834//   int i;
2835//   if(!attrib(J,"isSB")){J=std(J);}
2836//   intvec v = hilb(J,2);
2837//   for(i=1; i<=size(v); i++){ hp=hp+v[i]*(var(1)-i+2);}
2838//   return(hp);
2839//}
2840
2841
2842//////////////////////////////////////////////////////////////////////////////
2843
2844proc primeClosure (list L, list #)
2845"USAGE:    primeClosure(L [,c]); L a list of a ring containing a prime ideal
2846          ker, c an optional integer
2847RETURN:   a list L (of size n+1) consisting of rings L[1],...,L[n] such that
2848          - L[1] is a copy of (not a reference to!) the input ring L[1]
2849          - all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi
2850            such that
2851                    L[1]/ker --> ... --> L[n]/ker
2852            are injections given by the corresponding ideals phi, and L[n]/ker
2853            is the integral closure of L[1]/ker in its quotient field.
2854          - all rings L[i] contain a polynomial nzd such that elements of
2855            L[i]/ker are quotients of elements of L[i-1]/ker with denominator
2856            nzd via the injection phi.
2857            L[n+1] is the delta invariant
2858NOTE:     - L is constructed by recursive calls of primeClosure itself.
2859          - c determines the choice of nzd:
2860               - c not given or equal to 0: first generator of the ideal SL,
2861                 the singular locus of Spec(L[i]/ker)
2862               - c<>0: the generator of SL with least number of monomials.
2863EXAMPLE:  example primeClosure; shows an example
2864"
2865{
2866  //---- Start with a consistency check:
2867
2868  if (!(typeof(L[1])=="ring"))
2869  {
2870      "// Parameter must be a ring or a list containing a ring!";
2871      return(-1);
2872  }
2873
2874  int dblvl = printlevel-voice+2;
2875  list gnirlist = ringlist(basering);
2876
2877  //---- Some auxiliary variables:
2878  int delt;                      //finally the delta invariant
2879  if ( size(L) == 1 )
2880  {
2881      L[2] = delt;              //set delta to 0
2882  }
2883  int n = size(L)-1;            //L without delta invariant
2884
2885  //---- How to choose the non-zerodivisor later on?
2886
2887  int nzdoption=0;
2888  if (size(#)>0)
2889  {
2890      nzdoption=#[1];
2891  }
2892
2893// R0 below is the ring to work with, if we are in step one, make a copy of the
2894// input ring, so that all objects are created in the copy, not in the original
2895// ring (therefore a copy, not a reference is defined).
2896
2897  if (n==1)
2898  {
2899      def R = L[1];
2900      list Rlist = ringlist(R);
2901      def BAS = basering;
2902      setring R;
2903      if (!(typeof(ker)=="ideal"))
2904      {
2905          "// No ideal ker in the input ring!";
2906          return (-1);
2907      }
2908      ker=simplify(interred(ker),15);
2909      //execute ("ring R0="+charstr(R)+",("+varstr(R)+"),("+ordstr(R)+");");
2910      // Rlist may be not defined in this new ring, so we define it again.
2911      list Rlist2 = ringlist(R);
2912      def R0 = ring(Rlist2);
2913      setring R0;
2914      ideal ker=fetch(R,ker);
2915      // check whether we compute the normalization of the blow up of
2916      // an isolated singularity at the origin (checked in normalI)
2917
2918      if (typeof(attrib(L[1],"iso_sing_Rees"))=="int")
2919      {
2920        attrib(R0,"iso_sing_Rees",attrib(L[1],"iso_sing_Rees"));
2921      }
2922      L[1]=R0;
2923  }
2924  else
2925  {
2926      def R0 = L[n];
2927      setring R0;
2928  }
2929
2930// In order to apply HomJJ from normal.lib, we need the radical of the singular
2931// locus of ker, J:=rad(ker):
2932
2933   list SM=mstd(ker);
2934
2935// In the first iteration, we have to compute the singular locus "from
2936// scratch".
2937// In further iterations, we can fetch it from the previous one but
2938// have to compute its radical
2939// the next rings R1 contain already the (fetched) ideal
2940
2941  if (n==1)                              //we are in R0=L[1]
2942  {
2943      if (typeof(attrib(R0,"iso_sing_Rees"))=="int")
2944      {
2945        ideal J;
2946        for (int s=1;s<=attrib(R0,"iso_sing_Rees");s++)
2947        {
2948          J=J,var(s);
2949        }
2950        J = J,SM[2];
2951        list JM = mstd(J);
2952      }
2953      else
2954      {
2955        if ( dblvl >= 1 )
2956        {"";
2957           "// compute the singular locus";
2958        }
2959        //### Berechnung des singulaeren Orts geaendert (ist so schneller)
2960        ideal J = minor(jacob(SM[2]),nvars(basering)-dim(SM[1]),SM[1]);
2961        J = J,SM[2];
2962        list JM = mstd(J);
2963      }
2964
2965      if ( dblvl >= 1 )
2966      {"";
2967         "// dimension of singular locus is", dim(JM[1]);
2968         if (  dblvl >= 2 )
2969         {"";
2970            "// the singular locus is:"; JM[2];
2971         }
2972      }
2973
2974      if ( dblvl >= 1 )
2975      {"";
2976         "// compute radical of singular locus";
2977      }
2978
2979      J = simplify(radical(JM[2]),2);
2980      if ( dblvl >= 1 )
2981      {"";
2982         "// radical of singular locus is:"; J;
2983         pause();
2984      }
2985  }
2986  else
2987  {
2988      if ( dblvl >= 1 )
2989      {"";
2990         "// compute radical of test ideal in ideal of singular locus";
2991      }
2992      J = simplify(radical(J),2);
2993      if ( dblvl >= 1 )
2994      {"";
2995         "// radical of test ideal is:"; J;
2996         pause();
2997      }
2998  }
2999
3000  // having computed the radical J of/in the ideal of the singular locus,
3001  // we now need to pick an element nzd of J;
3002  // NOTE: nzd must be a non-zero divisor mod ker, i.e. not contained in ker
3003
3004  poly nzd = J[1];
3005  poly nzd1 = NF(nzd,SM[1]);
3006  if (nzd1 != 0)
3007  {
3008     if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) )
3009     {
3010        nzd = nzd1;
3011     }
3012  }
3013
3014  if (nzdoption || nzd1==0)
3015  {
3016    for (int ii=2;ii<=ncols(J);ii++)
3017    {
3018      nzd1 = NF(J[ii],SM[1]);
3019      if ( nzd1 != 0 )
3020      {
3021        if ( (deg(nzd)>=deg(J[ii])) && (size(nzd)>size(J[ii])) )
3022        {
3023          nzd=J[ii];
3024        }
3025        if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) )
3026        {
3027          nzd = nzd1;
3028        }
3029      }
3030    }
3031  }
3032
3033  export nzd;
3034  // In this case we do not eliminate variables, so that the maps
3035  // are well defined.
3036  list RR = SM[1],SM[2],J,nzd,1;
3037
3038  if ( dblvl >= 1 )
3039  {"";
3040     "// compute the first ring extension:";
3041     "RR: ";
3042     RR;
3043  }
3044
3045  list RS = HomJJ(RR);
3046  //NOTE: HomJJ creates new ring with variables X(i) and T(j)
3047//-------------------------------------------------------------------------
3048// If we've reached the integral closure (as determined by the result of
3049// HomJJ), then we are done, otherwise we have to prepare the next iteration.
3050
3051  if (RS[2]==1)     // we've reached the integral closure, we are still in R0
3052    {
3053      kill J;
3054      if ( n== 1)
3055      {
3056        def R1 = RS[1];
3057        setring R1;
3058        ideal norid, normap = endid, endphi;
3059        kill endid,  endphi;
3060
3061        //"//### case: primeClosure, final";
3062        //"size norid:", size(norid), size(string(norid));
3063        //"interred:";
3064        //norid = interred(norid);
3065        //"size norid:", size(norid), size(string(norid));
3066
3067        export (norid, normap);
3068        L[1] = R1;
3069      }
3070      return(L);
3071    }
3072  else                        // prepare the next iteration
3073    {
3074      if (n==1)               // In the first iteration: keep only the data
3075      {                       // needed later on.
3076         kill RR,SM;
3077         export(ker);
3078      }
3079      if ( dblvl >= 1 )
3080      {"";
3081         "// computing the next ring extension, we are in loop"; n+1;
3082      }
3083
3084      def R1 = RS[1];         // The data of the next ring R1:
3085      delt = RS[3];           // the delta invariant of the ring extension
3086      setring R1;             // keep only what is necessary and kill
3087      ideal ker=endid;        // everything else.
3088      export(ker);
3089      ideal norid=endid;
3090
3091      //"//### case: primeClosure, loop", n+1;
3092      //"size norid:", size(norid), size(string(norid));
3093      //"interred:";
3094      //norid = interred(norid);        //????
3095      //"size norid:", size(norid), size(string(norid));
3096
3097      export(norid);
3098      kill endid;
3099
3100      map phi = R0,endphi;                        // fetch the singular locus
3101      ideal J = mstd(simplify(phi(J)+ker,4))[2];  // ideal J in R1
3102      export(J);
3103      if(n>1)
3104      {
3105         ideal normap=phi(normap);
3106      }
3107      else
3108      {
3109         ideal normap=endphi;
3110      }
3111      export(normap);
3112      kill phi;              // we save phi as ideal, not as map, so that
3113      ideal phi=endphi;      // we have more flexibility in the ring names
3114      kill endphi;           // later on.
3115      export(phi);
3116      L=insert(L,R1,n);       // Add the new ring R1 and go on with the
3117                              // next iteration
3118      if ( L[size(L)] >= 0 && delt >= 0 )
3119      {
3120         delt = L[size(L)] + delt;
3121      }
3122      else
3123      {
3124         delt = -1;
3125      }
3126      L[size(L)] = delt;
3127
3128      if (size(#)>0)
3129      {
3130          return (primeClosure(L,#));
3131      }
3132      else
3133      {
3134          return(primeClosure(L));         // next iteration.
3135      }
3136    }
3137}
3138example
3139{
3140  "EXAMPLE:"; echo=2;
3141  ring R=0,(x,y),dp;
3142  ideal I=x4,y4;
3143  def K=ReesAlgebra(I)[1];        // K contains ker such that K/ker=R[It]
3144  list L=primeClosure(K);
3145  def R(1)=L[1];                  // L[4] contains ker, L[4]/ker is the
3146  def R(4)=L[4];                  // integral closure of L[1]/ker
3147  setring R(1);
3148  R(1);
3149  ker;
3150  setring R(4);
3151  R(4);
3152  ker;
3153}
3154
3155///////////////////////////////////////////////////////////////////////////////
3156
3157proc closureFrac(list L)
3158"USAGE:    closureFrac (L); L a list of size n+1 as in the result of
3159          primeClosure, L[n] contains an additional polynomial f
3160CREATE:   a list fraction of two elements of L[1], such that
3161          f=fraction[1]/fraction[2] via the injections phi L[i]-->L[i+1].
3162EXAMPLE:  example closureFrac; shows an example
3163"
3164{
3165// Define some auxiliary variables:
3166
3167  int n=size(L)-1;
3168  int i,j,k,l,n2,n3;
3169  intvec V;
3170  string mapstr;
3171  for (i=1; i<=n; i++) { def R(i) = L[i]; }
3172
3173// The quotient representing f is computed as in 'closureGenerators' with
3174// the differences that
3175//   - the loop is done twice: for the numerator and for the denominator;
3176//   - the result is stored in the list fraction and
3177//   - we have to make sure that no more objects of the rings R(i) survive.
3178
3179  for (j=1; j<=2; j++)
3180    {
3181      setring R(n);
3182      if (j==1)
3183      {
3184         poly p=f;
3185      }
3186      else
3187      {
3188         p=1;
3189      }
3190
3191      for (k=n; k>1; k--)
3192      {
3193          if (j==1)
3194          {
3195             map phimap=R(k-1),phi;
3196          }
3197
3198          p=p*phimap(nzd);
3199
3200          if (j==2)
3201          {
3202            kill phimap;
3203          }
3204
3205          if (j==1)
3206          {
3207             //### noch abfragen ob Z(i) definiert ist
3208             list gnirlist = ringlist(R(k));
3209             n2 = size(gnirlist[2]);
3210             n3 = size(gnirlist[3]);
3211             for( i=1; i<=ncols(phi); i++)
3212             {
3213               gnirlist[2][n2+i] = "Z("+string(i)+")";
3214             }
3215             V=0;
3216             V[ncols(phi)]=0; V=V+1;
3217             gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1);
3218             def S(k) = ring(gnirlist);
3219             setring S(k);
3220
3221             //execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
3222             //          Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
3223             //          +"),dp("+string(ncols(phi))+"));");
3224
3225              ideal phi = imap(R(k),phi);
3226              ideal J = imap (R(k),ker);
3227              for (l=1;l<=ncols(phi);l++)
3228              {
3229                  J=J+(Z(l)-phi[l]);
3230              }
3231              J=groebner(J);
3232              poly h=NF(imap(R(k),p),J);
3233          }
3234          else
3235          {
3236              setring S(k);
3237              h=NF(imap(R(k),p),J);
3238              setring R(k);
3239              kill p;
3240          }
3241
3242          setring R(k-1);
3243
3244          if (j==1)
3245          {
3246              ideal maxi;
3247              maxi[nvars(R(k))] = 0;
3248              maxi = maxi,maxideal(1);
3249              map backmap = S(k),maxi;
3250
3251              //mapstr=" map backmap = S(k),";
3252              //for (l=1;l<=nvars(R(k));l++)
3253              //{
3254              //  mapstr=mapstr+"0,";
3255              //}
3256              //execute (mapstr+"maxideal(1);");
3257              poly p;
3258          }
3259          p=NF(backmap(h),std(ker));
3260          if (j==2)
3261          {
3262            kill backmap;
3263          }
3264        }
3265
3266      if (j==1)
3267        {
3268          if (defined(fraction))
3269            {
3270              kill fraction;
3271              list fraction=p;
3272            }
3273          else
3274            {
3275              list fraction=p;
3276            }
3277        }
3278      else
3279        {
3280          fraction=insert(fraction,p,1);
3281        }
3282    }
3283  export(fraction);
3284  return ();
3285}
3286example
3287{
3288  "EXAMPLE:"; echo=2;
3289  ring R=0,(x,y),dp;
3290  ideal ker=x2+y2;
3291  export ker;
3292  list L=primeClosure(R);          // We normalize R/ker
3293  for (int i=1;i<=size(L);i++) { def R(i)=L[i]; }
3294  setring R(2);
3295  kill R;
3296  phi;                             // The map R(1)-->R(2)
3297  poly f=T(2);                     // We will get a representation of f
3298  export f;
3299  L[2]=R(2);
3300  closureFrac(L);
3301  setring R(1);
3302  kill R(2);
3303  fraction;                        // f=fraction[1]/fraction[2] via phi
3304  kill R(1);
3305}
3306
3307///////////////////////////////////////////////////////////////////////////////
3308// closureGenerators is called inside proc normal (option "withGens" )
3309//
3310
3311// INPUT is the output of proc primeClosure (except for the last element, the
3312// delta invariant) : hence input is a list L consisting of rings
3313// L[1],...,L[n] (denoted R(1)...R(n) below) such that
3314// - L[1] is a copy of (not a reference to!) the input ring L[1]
3315// - all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi
3316// such that
3317//                L[1]/ker --> ... --> L[n]/ker
3318// are injections given by the corresponding ideals phi, and L[n]/ker
3319// is the integral closure of L[1]/ker in its quotient field.
3320// - all rings L[i] contain a polynomial nzd such that elements of
3321// L[i]/ker are quotients of elements of L[i-1]/ker with denominator
3322// nzd via the injection phi.
3323
3324// COMPUTE: In the list L of rings R(1),...,R(n), compute representations of
3325// the ring variables of the last ring R(n) as fractions of elements of R(1):
3326// The proc returns an ideal preim s.t. preim[i]/preim[size(preim)] expresses
3327// the ith variable of R(n) as fraction of elements of the basering R(1)
3328// preim[size(preim)] is a non-zero divisor of basering/i.
3329
3330proc closureGenerators(list L);
3331{
3332  def Rees=basering;         // when called inside normalI (in reesclos.lib)
3333                             // the Rees Algebra is the current basering
3334
3335  // ------- First of all we need some variable declarations -----------
3336  int n = size(L);                // the number of rings R(1)-->...-->R(n)
3337  int length = nvars(L[n]);       // the number of variables of the last ring
3338  int j,k,l,n2,n3;
3339  intvec V;
3340  string mapstr;
3341  list preimages;
3342  //Note: the empty list belongs to no ring, hence preimages can be used
3343  //later in R(1)
3344  //this is not possible for ideals (belong always to some ring)
3345
3346  for (int i=1; i<=n; i++)
3347  {
3348     def R(i)=L[i];          //give the rings from L a name
3349  }
3350
3351  // For each variable (counter j) and for each intermediate ring (counter k):
3352  // Find a preimage of var_j*phi(nzd(k-1)) in R(k-1).
3353  // Finally, do the same for nzd.
3354
3355  for (j=1; j <= length+1; j++ )
3356  {
3357      setring R(n);
3358      if (j==1)
3359      {
3360        poly p;
3361      }
3362      if (j <= length )
3363      {
3364        p=var(j);
3365      }
3366      else
3367      {
3368        p=1;
3369      }
3370      //i.e. p=j-th var of R(n) for j<=length and p=1 for j=length+1
3371
3372      for (k=n; k>1; k--)
3373      {
3374
3375        if (j==1)
3376        {
3377          map phimap=R(k-1),phi;   //phimap:R(k-1)-->R(n), k=2..n, is the map
3378                                   //belonging to phi in R(n)
3379        }
3380
3381        p = p*phimap(nzd);
3382
3383          // Compute the preimage of [p mod ker(k)] under phi in R(k-1):
3384          // As p is an element of Image(phi), there is a polynomial h such
3385          // that h is mapped to [p mod ker(k)], and h can be computed as the
3386          // normal form of p w.r.t. a Groebner basis of
3387          // J(k) := <ker(k),Z(l)-phi(k)(l)> in R(k)[Z]=:S(k)
3388
3389        if (j==1)   // In the first iteration: Create S(k), fetch phi and
3390                    // ker(k) and construct the ideal J(k).
3391        {
3392         //### noch abfragen ob Z(i) definiert ist
3393         list gnirlist = ringlist(R(k));
3394         n2 = size(gnirlist[2]);
3395         n3 = size(gnirlist[3]);
3396         for( i=1; i<=ncols(phi); i++)
3397         {
3398            gnirlist[2][n2+i] = "Z("+string(i)+")";
3399         }
3400         V=0;
3401         V[ncols(phi)]=0;
3402         V=V+1;
3403         gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1);
3404         def S(k) = ring(gnirlist);
3405         setring S(k);
3406
3407        // execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
3408        //           Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
3409        //           +"),dp("+string(ncols(phi))+"));");
3410
3411          ideal phi = imap(R(k),phi);
3412          ideal J = imap (R(k),ker);
3413          for ( l=1; l<=ncols(phi); l++ )
3414          {
3415             J=J+(Z(l)-phi[l]);
3416          }
3417          J = groebner(J);
3418          poly h = NF(imap(R(k),p),J);
3419        }
3420        else
3421        {
3422           setring S(k);
3423           h = NF(imap(R(k),p),J);
3424        }
3425
3426        setring R(k-1);
3427
3428        if (j==1)  // In the first iteration: Compute backmap:S(k)-->R(k-1)
3429        {
3430           ideal maxi;
3431           maxi[nvars(R(k))] = 0;
3432           maxi = maxi,maxideal(1);
3433           map backmap = S(k),maxi;
3434
3435           //mapstr=" map backmap = S(k),";
3436           //for (l=1;l<=nvars(R(k));l++)
3437           //{
3438           //  mapstr=mapstr+"0,";
3439           //}
3440           //execute (mapstr+"maxideal(1);");
3441
3442           poly p;
3443        }
3444        p = NF(backmap(h),std(ker));
3445     }
3446     // Whe are down to R(1), store here the result in the list preimages
3447     preimages = insert(preimages,p,j-1);
3448  }
3449  ideal preim;                  //make the list preimages to an ideal preim
3450  for ( i=1; i<=size(preimages); i++ )
3451  {
3452     preim[i] = preimages[i];
3453  }
3454  // R(1) was a copy of Rees, so we have to get back to the basering Rees from
3455  // the beginning and fetch the result (the ideal preim) to this ring.
3456  setring Rees;
3457  return (fetch(R(1),preim));
3458}
3459
3460///////////////////////////////////////////////////////////////////////////////
3461//                From here: procedures for char p with Frobenius
3462///////////////////////////////////////////////////////////////////////////////
3463
3464proc normalP(ideal id,list #)
3465"USAGE:  normalP(id [,choose]); id = radical ideal, choose = optional list of
3466         strings.
3467         Optional parameters in list choose (can be entered in any order):@*
3468         \"withRing\", \"isPrim\", \"noFac\", \"noRed\", where@*
3469         - \"noFac\" -> factorization is avoided during the computation
3470         of the minimal associated primes.@*
3471         - \"isPrim\" -> assumes that the ideal is prime. If the assumption
3472         does not hold, output might be wrong.@*
3473         - \"withRing\" -> the ring structure of the normalization is
3474         computed. The number of variables in the new ring is reduced as much
3475         as possible.@*
3476         - \"noRed\" -> when computing the ring structure, no reduction on the
3477         number of variables is done, it creates one new variable for every
3478         new module generator of the integral closure in the quotient field.@*
3479ASSUME:  The characteristic of the ground field must be positive. If the
3480         option \"isPrim\" is not set, the minimal associated primes of id
3481         are computed first and hence normalP computes the normalization of
3482         the radical of id. If option \"isPrim\" is set, the ideal must be
3483         a prime ideal otherwise the result may be wrong.
3484RETURN:  a list, say 'nor' of size 2 (resp. 3 if \"withRing\" is set).@*
3485         ** If option \"withRing\" is not set: @*
3486         Only the module structure is computed: @*
3487         * nor[1] is a list of ideals Ii, i=1..r, in the basering R where r
3488         is the number of minimal associated prime ideals P_i of the input
3489         ideal id, describing the module structure:@*
3490         If Ii is given by polynomials g_1,...,g_k in R, then c:=g_k is
3491         non-zero in the ring R/P_i and g_1/c,...,g_k/c generate the integral
3492         closure of R/P_i as R-module in the quotient field of R/P_i.@*
3493         * nor[2] shows the delta invariants: it is a list of an intvec
3494         of size r, the delta invariants of the r components, and an integer,
3495         the total delta invariant of R/id
3496         (-1 means infinite, and 0 that R/P_i resp. R/id is normal). @*
3497         ** If option \"withRing\" is set: @*
3498         The ring structure is also computed, and in this case:@*
3499         * nor[1] is a list of r rings.@*
3500         Each ring Ri = nor[1][i], i=1..r, contains two ideals with given
3501         names @code{norid} and @code{normap} such that @*
3502         - Ri/norid is the normalization of R/P_i, i.e. isomorphic as
3503           K-algebra (K the ground field) to the integral closure of R/P_i in
3504           the field of fractions of R/P_i; @*
3505         - the direct sum of the rings Ri/norid is the normalization
3506           of R/id; @*
3507         - @code{normap} gives the normalization map from R to Ri/norid.@*
3508         * nor[2] gives the module generators of the normalization of R/P_i,
3509         it is the same as nor[1] if \"withRing\" is not set.@*
3510         * nor[3] shows the delta invariants, it is the same as nor[2] if
3511         \"withRing\" is not set.
3512THEORY:  normalP uses the Leonard-Pellikaan-Singh-Swanson algorithm (using the
3513         Frobenius) cf. [A. K. Singh, I. Swanson: An algorithm for computing
3514         the integral closure, arXiv:0901.0871].
3515         The delta invariant of a reduced ring A is dim_K(normalization(A)/A).
3516         For A=K[x1,...,xn]/id we call this number also the delta invariant of
3517         id. The procedure returns the delta invariants of the components P_i
3518         and of id.
3519NOTE:    To use the i-th ring type: @code{def R=nor[1][i]; setring R;}.
3520@*       Increasing/decreasing printlevel displays more/less comments
3521         (default: printlevel = 0).
3522@*       Not implemented for local or mixed orderings or quotient rings.
3523         For local or mixed orderings use proc 'normal'.
3524@*       If the input ideal id is weighted homogeneous a weighted ordering may
3525         be used (qhweight(id); computes weights).
3526@*       Works only in characteristic p > 0; use proc normal in char 0.
3527KEYWORDS: normalization; integral closure; delta invariant.
3528SEE ALSO: normal, normalC
3529EXAMPLE: example normalP; shows an example
3530"
3531{
3532   int i,j,y, sr, del, co;
3533   intvec deli;
3534   list resu, Resu, prim, Gens, mstdid;
3535   ideal gens;
3536
3537   // Default options
3538   int wring = 0;           // The ring structure is not computed.
3539   int noRed = 0;           // No reduction is done in the ring structure
3540   int isPrim = 0;          // Ideal is not assumed to be prime
3541   int noFac = 0;           // Use facstd when computing the decomposition
3542
3543
3544   y = printlevel-voice+2;
3545
[1e1ec4]3546   if ( attrib(basering,"global") != 1)
[1598658]3547   {
3548     "";
3549     "// Not implemented for this ordering,";
3550     "// please change to global ordering!";
3551     return(resu);
3552   }
3553   if ( char(basering) <= 0)
3554   {
3555     "";
3556     "// Algorithm works only in positive characteristic,";
3557     "// use procedure 'normal' if the characteristic is 0";
3558     return(resu);
3559   }
3560
3561//--------------------------- define the method ---------------------------
3562   string method;                //make all options one string in order to use
3563                                 //all combinations of options simultaneously
3564   for ( i=1; i<= size(#); i++ )
3565   {
3566     if ( typeof(#[i]) == "string" )
3567     {
3568       method = method + #[i];
3569     }
3570   }
3571
3572   if ( find(method,"withring") or find(method,"withRing") )
3573   {
3574     wring=1;
3575   }
3576   if ( find(method,"noRed") or find(method,"nored") )
3577   {
3578     noRed=1;
3579   }
3580   if ( find(method,"isPrim") or find(method,"isprim") )
3581   {
3582     isPrim=1;
3583   }
3584   if ( find(method,"noFac") or find(method,"nofac"))
3585   {
3586     noFac=1;
3587   }
3588
3589   kill #;
3590   list #;
3591//--------------------------- start computation ---------------------------
3592   ideal II,K1,K2;
3593
3594   //----------- check first (or ignore) if input id is prime -------------
3595
3596   if ( isPrim )
3597   {
3598      prim[1] = id;
3599      if( y >= 0 )
3600      { "";
3601    "// ** WARNING: result is correct if ideal is prime (not checked) **";
3602    "// disable option \"isPrim\" to decompose ideal into prime components";"";
3603      }
3604   }
3605   else
3606   {
3607      if(y>=1)
3608      {  "// compute minimal associated primes"; }
3609
3610      if( noFac )
3611      { prim = minAssGTZ(id,1); }
3612      else
3613      { prim = minAssGTZ(id); }
3614
3615      if(y>=1)
3616      {
3617         prim;"";
3618         "// number of irreducible components is", size(prim);
3619      }
3620   }
3621
3622   //----------- compute integral closure for every component -------------
3623
3624      for(i=1; i<=size(prim); i++)
3625      {
3626         if(y>=1)
3627         {
3628            ""; pause(); "";
3629            "// start computation of component",i;
3630            "   --------------------------------";
3631         }
3632         if(y>=1)
3633         {  "// compute SB of ideal";
3634         }
3635         mstdid = mstd(prim[i]);
3636         if(y>=1)
3637         {  "// dimension of component is", dim(mstdid[1]);"";}
3638
3639      //------- 1-st main subprocedure: compute module generators ----------
3640         printlevel = printlevel+1;
3641         II = normalityTest(mstdid);
3642
3643      //------ compute also the ringstructure if "withRing" is given -------
3644         if ( wring )
3645         {
3646         //------ 2-nd main subprocedure: compute ring structure -----------
3647           if(noRed == 0){
3648             resu = list(computeRing(II,prim[i])) + resu;
3649           }
3650           else
3651           {
3652             resu = list(computeRing(II,prim[i], "noRed")) + resu;
3653           }
3654         }
3655         printlevel = printlevel-1;
3656
3657      //----- rearrange module generators s.t. denominator comes last ------
3658         gens=0;
3659         for( j=2; j<=size(II); j++ )
3660         {
3661            gens[j-1]=II[j];
3662         }
3663         gens[size(gens)+1]=II[1];
3664         Gens = list(gens) + Gens;
3665      //------------------------------ compute delta -----------------------
3666         K1 = mstdid[1]+II;
3667         K1 = std(K1);
3668         K2 = mstdid[1]+II[1];
3669         K2 = std(K2);
3670         // K1 = std(mstdid[1],II);      //### besser
3671         // K2 = std(mstdid[1],II[1]);   //### besser: Hannes, fixen!
3672         co = codim(K1,K2);
3673         deli = co,deli;
3674         if ( co >= 0 && del >= 0 )
3675         {
3676            del = del + co;
3677         }
3678         else
3679         { del = -1; }
3680      }
3681
3682      if ( del >= 0 )
3683      {
3684         int mul = iMult(prim);
3685         del = del + mul;
3686      }
3687      else
3688      { del = -1; }
3689
3690      deli = deli[1..size(deli)-1];
3691      if ( wring )
3692      { Resu = resu,Gens,list(deli,del); }
3693      else
3694      { Resu = Gens,list(deli,del); }
3695
3696   sr = size(prim);
3697
3698//-------------------- Finally print comments and return --------------------
3699   if(y >= 0)
3700   {"";
3701     if ( wring )
3702     {
3703"// 'normalP' created a list, say nor, of three lists:
3704// To see the result, type
3705     nor;
3706
3707// * nor[1] is a list of",sr,"ring(s):
3708// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
3709     def R1 = nor[1][1]; setring R1;  norid; normap;
3710// for the other rings type first setring R; (if R is the name of your
3711// original basering) and then continue as for R1;
3712// Ri/norid is the affine algebra of the normalization of the i-th
3713// component R/P_i (where P_i is a min. associated prime of the input ideal)
3714// and normap the normalization map from R to Ri/norid;
3715
3716// * nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of
3717// elements g1..gk of r such that the gj/gk generate the integral
3718// closure of R/P_i as R-module in the quotient field of R/P_i.
3719
3720// * nor[3] shows the delta-invariant of each component and of the input
3721// ideal (-1 means infinite, and 0 that r/P_i is normal).";
3722     }
3723     else
3724     {
3725"// 'normalP' computed a list, say nor, of two lists:
3726// To see the result, type
3727     nor;
3728
3729// * nor[1] is a list of",sr,"ideal(s), where each ideal nor[1][i] consists
3730// of elements g1..gk of the basering R such that gj/gk generate the integral
3731// closure of R/P_i (where P_i is a min. associated prime of the input ideal)
3732// as R-module in the quotient field of R/P_i;
3733
3734// * nor[2] shows the delta-invariant of each component and of the input ideal
3735// (-1 means infinite, and 0 that R/P_i is normal).";
3736     }
3737   }
3738
3739   return(Resu);
3740}
3741example
3742{ "EXAMPLE:"; echo = 2;
3743   ring r  = 11,(x,y,z),wp(2,1,2);
3744   ideal i = x*(z3 - xy4 + x2);
3745   list nor= normalP(i); nor;
3746   //the result says that both components of i are normal, but i itself
3747   //has infinite delta
3748   pause("hit return to continue");
3749
3750   ring s = 2,(x,y),dp;
3751   ideal i = y*((x-y^2)^2 - x^3);
3752   list nor = normalP(i,"withRing"); nor;
3753
3754   def R2  = nor[1][2]; setring R2;
3755   norid; normap;
3756}
3757
3758///////////////////////////////////////////////////////////////////////////////
3759// Assume: mstdid is the result of mstd(prim[i]), prim[i] a prime component of
3760// the input ideal id of normalP.
3761// Output is an ideal U s.t. U[i]/U[1] are module generators.
3762
3763static proc normalityTest(list mstdid)
3764{
3765   int y = printlevel-voice+2;
3766   intvec op = option(get);
3767   option(redSB);
3768   def R = basering;
3769   int n, p = nvars(R), char(R);
3770   int ii;
3771
3772   ideal J = mstdid[1];         //J is the SB of I
3773   ideal I = mstdid[2];
3774   int h = n-dim(J);            //codimension of V(I), I is a prime ideal
3775
3776   //-------------------------- compute singular locus ----------------------
3777   qring Q = J;                 //pass to quotient ring
3778   ideal I = imap(R,I);
3779   ideal J = imap(R,J);
3780   attrib(J,"isSB",1);
3781   if ( y >= 1)
3782   { "start normality test";  "compute singular locus";}
3783
3784   ideal M = minor(jacob(I),h,J); //use the command minor modulo J (hence J=0)
3785   M = std(M);                    //this makes M much smaller
3786   //keep only minors which are not 0 mod I (!) this is important since we
3787   //need a nzd mod I
3788
3789   //---------------- choose nzd from ideal of singular locus --------------
3790   ideal D = M[1];
3791   for( ii=2; ii<=size(M); ii++ )            //look for the shortest one
3792   {
3793      if( size(M[ii]) < size(D[1]) )
3794      {
3795          D = M[ii];
3796      }
3797   }
3798
3799   //--------------- start p-th power algorithm and return ----------------
3800   ideal F = var(1)^p;
3801   for(ii=2; ii<=n; ii++)
3802   {
3803      F=F,var(ii)^p;
3804   }
3805
3806   ideal Dp=D^(p-1);
3807   ideal U=1;
3808   ideal K,L;
3809   map phi=Q,F;
3810   if ( y >= 1)
3811   {  "compute module generators of integral closure";
3812      "denominator D is:";  D;
3813      pause();
3814   }
3815
3816   ii=0;
3817   list LK;
3818   while(1)
3819   {
3820      ii=ii+1;
3821      if ( y >= 1)
3822      { "iteration", ii; }
3823      L = U*Dp + I;
3824      //### L=interred(L) oder mstd(L)[2]?
3825      //Wird dadurch kleiner aber string(L) wird groesser
3826      K = preimage(Q,phi,L);    //### Improvement by block ordering?
3827      option(returnSB);
3828      K = intersect(U,K);          //K is the new U, it is a SB
3829      LK = mstd(K);
3830      K = LK[2];
3831
3832   //---------------------------- simplify output --------------------------
3833      if(size(reduce(U,LK[1]))==0)  //previous U coincides with new U
3834      {                             //i.e. we reached the integral closure
3835         U=simplify(reduce(U,groebner(D)),2);
3836         U = D,U;
3837         poly gg = gcd(U[1],U[size(U)]);
3838         for(ii=2; ii<=size(U)-1 ;ii++)
3839         {
3840            gg = gcd(gg,U[ii]);
3841         }
3842         for(ii=1; ii<=size(U); ii++)
3843         {
3844            U[ii]=U[ii]/gg;
3845         }
3846         U = simplify(U,6);
3847         //if ( y >= 1)
3848         //{ "module generators are U[i]/U[1], with U:"; U;
3849         //  ""; pause(); }
3850         option(set,op);
3851         setring R;
3852         ideal U = imap(Q,U);
3853         return(U);
3854      }
3855      U=K;
3856   }
3857}
3858
3859///////////////////////////////////////////////////////////////////////////////
3860
3861static proc substpartSpecial(ideal endid, ideal endphi)
3862{
3863   //Note: newRing is of the form (R the original basering):
3864   //char(R),(T(1..N),X(1..nvars(R))),(dp(N),...);
3865
3866   int ii,jj,kk;
3867   def BAS = basering;
3868   int n = nvars(basering);
3869
3870   list Le = elimpart(endid);
3871   int q = size(Le[2]);                   //q variables have been substituted
3872//Le;"";
3873   if ( q == 0 )
3874   {
3875      ideal normap = endphi;
3876      ideal norid = endid;
3877      export(norid);
3878      export(normap);
3879      list L = BAS;
3880      return(L);
3881   }
3882
3883      list gnirlist = ringlist(basering);
3884      endid = Le[1];
3885//endphi;"";
3886      for( ii=1; ii<=n; ii++)
3887      {
3888         if( Le[4][ii] == 0 )            //ii=index of substituted var
3889         {
3890            endphi = subst(endphi,var(ii),Le[5][ii]);
3891         }
3892      }
3893//endphi;"";
3894      list g2 = gnirlist[2];             //the varlist
3895      list g3 = gnirlist[3];             //contains blocks of orderings
3896      int n3 = size(g3);
3897
3898   //----------------- first identify module ordering ------------------
3899      if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
3900      {
3901         list gm = g3[n3];              //last blockis module ordering
3902         g3 = delete(g3,n3);
3903         int m = 0;
3904      }
3905      else
3906      {
3907         list gm = g3[1];              //first block is module ordering
3908         g3 = delete(g3,1);
3909         int m = 1;
3910      }
3911   //---- delete variables which were substituted and weights  --------
3912      intvec V;
3913      int n(0);
3914      list newg2;
3915      list newg3;
3916      for ( ii=1; ii<=n3-1; ii++ )
3917      {
3918        // If order is a matrix ordering, it is replaced by dp ordering.
3919        // TODO: replace it only when some of the original
3920        //       variables are eliminated.
3921        if(g3[ii][1] == "M"){
3922          g3[ii][1] = "dp";
3923          g3[ii][2] = (1..sqroot(size(g3[ii][2])))*0+1;
3924        }
3925        V = V,g3[ii][2];           //copy weights for ordering in each block
3926        if ( ii==1 )               //into one intvector
3927        {
3928           V = V[2..size(V)];
3929        }
3930        // int n(ii) = size(g3[ii][2]);
3931        int n(ii) = size(V);
3932        intvec V(ii);
3933
3934        for ( jj = n(ii-1)+1; jj<=n(ii); jj++)
3935        {
3936          if( Le[4][jj] !=0 )     //jj=index of var which was not substituted
3937          {
3938            kk=kk+1;
3939            newg2[kk] = g2[jj];   //not substituted var from varlist
3940            V(ii)=V(ii),V[jj];    //weight of not substituted variable
3941          }
3942        }
3943        if ( size(V(ii)) >= 2 )
3944        {
3945           V(ii) = V(ii)[2..size(V(ii))];
3946           list g3(ii)=g3[ii][1],V(ii);
3947           newg3 = insert(newg3,g3(ii),size(newg3));
3948//"newg3"; newg3;
3949        }
3950      }
3951//"newg3"; newg3;
3952      //newg3 = delete(newg3,1);    //delete empty list
3953
3954/*
3955//### neue Ordnung, 1 Block fuer alle vars, aber Gewichte erhalten;
3956//vorerst nicht realisiert, da bei leonhard1 alte Version (neue Variable T(i)
3957//ein neuer Block) ein kuerzeres Ergebnis liefert
3958      kill g3;
3959      list g3;
3960      V=0;
3961      for ( ii= 1; ii<=n3-1; ii++ )
3962      {
3963        V=V,V(ii);
3964      }
3965      V = V[2..size(V)];
3966
3967      if ( V==1 )
3968      {
3969         g3[1] = list("dp",V);
3970      }
3971      else
3972      {
3973         g3[1] = lis("wp",V);
3974      }
3975      newg3 = g3;
3976
3977//"newg3";newg3;"";
3978//### Ende neue Ordnung
3979*/
3980
3981      if ( m == 0 )
3982      {
3983         newg3 = insert(newg3,gm,size(newg3));
3984      }
3985      else
3986      {
3987         newg3 = insert(newg3,gm);
3988      }
3989      gnirlist[2] = newg2;
3990      gnirlist[3] = newg3;
3991
3992//gnirlist;
3993      def newBAS = ring(gnirlist);            //change of ring to less vars
3994      setring newBAS;
3995      ideal normap = imap(BAS,endphi);
3996      //normap = simplify(normap,2);
3997      ideal norid =  imap(BAS,endid);
3998      export(norid);
3999      export(normap);
4000      list L = newBAS;
4001      return(L);
4002
4003   //Hier scheint interred gut zu sein, da es Ergebnis relativ schnell
4004   //verkleinert. Hier wird z.B. bei leonard1 size(norid) verkleinert aber
4005   //size(string(norid)) stark vergroessert, aber es hat keine Auswirkungen
4006   //da keine map mehr folgt.
4007   //### Bei Leonard2 haengt interred (BUG)
4008   //mstd[2] verkleinert norid nocheinmal um die Haelfte, dauert aber 3.71 sec
4009   //### Ev. Hinweis auf mstd in der Hilfe?
4010
4011}
4012
4013///////////////////////////////////////////////////////////////////////////////
4014// Computes the ring structure of a ring given by module generators.
4015// Assume: J[i]/J[1] are the module generators in the quotient field
4016// with J[1] as universal denominator.
4017// If option "noRed" is not given, a reduction in the number of variables is
4018// attempted.
4019static proc computeRing(ideal J, ideal I, list #)
4020{
[9ebd82]4021  int i, ii,jj;
[1598658]4022  intvec V;                          // to be used for variable weights
[9ebd82]4023  int y = printlevel-voice+2;
[1598658]4024  def R = basering;
4025  poly c = J[1];                     // the denominator
4026  list gnirlist = ringlist(basering);
4027  string svars = varstr(basering);
4028  int nva = nvars(basering);
4029  string svar;
4030  ideal maxid = maxideal(1);
4031
4032  int noRed = 0;     // By default, we try to reduce the number of generators.
4033  if(size(#) > 0){
4034    if ( typeof(#[1]) == "string" )
4035    {
4036      if (#[1] == "noRed"){noRed = 1;}
4037    }
4038  }
4039
4040  if ( y >= 1){"// computing the ring structure...";}
4041
4042  if(c==1)
4043  {
4044/*    if( defined(norid) )  { kill norid; }
4045      if( defined(normap) ) { kill normap; }
4046      ideal norid = I;
4047      ideal normap =  maxid;  */
4048
4049    list gnirlist = ringlist(R);
4050    def R1 = ring(gnirlist);
4051    setring R1;
4052    ideal norid = imap(R, I);
4053    ideal normap = imap(R, maxid);
4054    export norid;
4055    export normap;
4056
4057    if(noRed == 1){
4058      setring R;
4059      return(R1);
4060    }
4061    else
4062    {
4063      list L = substpartSpecial(norid,normap);
4064      def lastRing = L[1];
4065      setring R;
4066      return(lastRing);
4067    }
4068  }
4069
4070
4071  //-------------- Enlarge ring by creating new variables ------------------
4072  //check first whether variables T(i) and then whether Z(i),...,A(i) exist
4073  //old variable names are not touched
4074
4075  if ( find(svars,"T(") == 0 )
4076  {
4077    svar = "T";
4078  }
4079  else
4080  {
4081    for (ii=90; ii>=65; ii--)
4082    {
4083      if ( find(svars,ASCII(ii)+"(") == 0 )
4084      {
4085        svar = ASCII(ii);  break;
4086      }
4087    }
4088  }
4089
4090  int q = size(J)-1;
4091  if ( size(svar) != 0 )
4092  {
4093    for ( ii=q; ii>=1; ii-- )
4094    {
4095      gnirlist[2] = insert(gnirlist[2],svar+"("+string(ii)+")");
4096    }
4097  }
4098  else
4099  {
4100    for ( ii=q; ii>=1; ii-- )
4101    {
4102      gnirlist[2] = insert(gnirlist[2],"T("+string(100*nva+ii)+")");
4103    }
4104  }
4105
4106  V[q]=0;                        //create intvec of variable weights
4107  V=V+1;
4108  gnirlist[3] = insert(gnirlist[3],list("dp",V));
4109
4110  //this is a block ordering with one dp-block (1st block) for new vars
4111  //the remaining weights and blocks for old vars are kept
4112  //### perhaps better to make only one block, keeping weights ?
4113  //this might effect syz below
4114  //alt: ring newR = char(R),(X(1..nvars(R)),T(1..q)),dp;
4115  //Reihenfolge geaendert:neue Variablen kommen zuerst, Namen ev. nicht T(i)
4116
4117  def newR = ring(gnirlist);
4118  setring newR;                //new extended ring
4119  ideal I = imap(R,I);
4120
4121  //------------- Compute linear and quadratic relations ---------------
4122  if(y>=1)
4123  {
4124     "// compute linear relations:";
4125  }
4126  qring newQ = std(I);
4127
4128  ideal f = imap(R,J);
4129  module syzf = syz(f);
4130  ideal pf = f[1]*f;
4131  //f[1] is the denominator D from normalityTest, a non zero divisor of R/I
4132
4133  ideal newT = maxideal(1);
4134  newT = 1,newT[1..q];
4135  //matrix T = matrix(ideal(1,T(1..q)),1,q+1);   //alt
4136  matrix T = matrix(newT,1,q+1);
4137  ideal Lin = ideal(T*syzf);
4138  //Lin=interred(Lin);
4139  //### interred reduziert ev size aber size(string(LIN)) wird groesser
4140
4141  if(y>=1)
4142  {
4143    if(y>=3)
4144    {
4145      "//   the linear relations:";  Lin; pause();"";
4146    }
4147      "// the ring structure of the normalization as affine algebra";
4148      "//   number of linear relations:", size(Lin);
4149  }
4150
4151  if(y>=1)
4152  {
4153    "// compute quadratic relations:";
4154  }
4155  matrix A;
4156  ideal Quad;
4157  poly ff;
4158  newT = newT[2..size(newT)];
4159  matrix u;  // The units for non-global orderings.
4160
4161  // Quadratic relations
4162  for (ii=2; ii<=q+1; ii++ )
4163  {
4164    for ( jj=2; jj<=ii; jj++ )
4165    {
4166      ff = NF(f[ii]*f[jj],std(0));     // this makes lift much faster
4167      // For non-global orderings, we have to take care of the units.
[1e1ec4]4168      if(attrib(basering,"global") != 1)
4169      {
[1598658]4170        A = lift(pf, ff, u);
4171        Quad = Quad,ideal(newT[jj-1]*newT[ii-1] * u[1, 1]- T*A);
4172      }
4173      else
4174      {
4175        A = lift(pf,ff);              // ff lin. comb. of elts of pf mod I
4176        Quad = Quad,ideal(newT[jj-1]*newT[ii-1] - T*A);
4177      }
4178      //A = lift(pf, f[ii]*f[jj]);
4179      //Quad = Quad, ideal(T(jj-1)*T(ii-1) - T*A);
4180    }
4181  }
4182  Quad = Quad[2..ncols(Quad)];
4183
4184  if(y>=1)
4185  {
4186    if(y>=3)
4187    {
4188      "//   the quadratic relations"; Quad; pause();"";
4189    }
4190      "//   number of quadratic relations:", size(Quad);
4191  }
4192  ideal Q1 = Lin,Quad;     //elements of Q1 are in NF w.r.t. I
4193
4194  //Q1 = mstd(Q1)[2];
4195  //### weglassen, ist sehr zeitaufwendig.
4196  //Ebenso interred, z.B. bei Leonard1 (1. Komponente von Leonard):
4197  //"size Q1:", size(Q1), size(string(Q1));   //75 60083
4198  //Q1 = interred(Q1);
4199  //"size Q1:", size(Q1), size(string(Q1));   //73 231956 (!)
4200  //### Speicherueberlauf bei substpartSpecial bei 'ideal norid  = phi1(endid)'
4201  //Beispiel fuer Hans um map zu testen!
4202
4203  setring newR;
4204  ideal endid  = imap(newQ,Q1),I;
4205  ideal endphi = imap(R,maxid);
4206
4207  if(noRed == 0){
4208    list L=substpartSpecial(endid,endphi);
4209    def lastRing=L[1];
4210    if(y>=1)
4211    {
4212      "//   number of substituted variables:", nvars(newR)-nvars(lastRing);
4213      pause();"";
4214    }
4215    return(lastRing);
4216  }
4217  else
4218  {
4219    ideal norid = endid;
4220    ideal normap = endphi;
4221    export(norid);
4222    export(normap);
4223    setring R;
4224    return(newR);
4225  }
4226}
4227
4228//                Up to here: procedures for char p with Frobenius
4229///////////////////////////////////////////////////////////////////////////////
4230
4231///////////////////////////////////////////////////////////////////////////////
4232//                From here: subprocedures for normal
4233
4234// inputJ is used in parametrization of rational curves algorithms, to specify
4235// a different test ideal.
4236
4237static proc normalM(ideal I, int decomp, int withDelta, int denomOption, ideal inputJ, ideal inputC){
4238// Computes the normalization of a ring R / I using the module structure as far
4239// as possible.
4240// The ring R is the basering.
4241// Input: ideal I
4242// Output: a list of 3 elements (resp 4 if withDelta = 1), say nor.
4243// - nor[1] = U, an ideal of R.
4244// - nor[2] = c, an element of R.
4245// U and c satisfy that 1/c * U is the normalization of R / I in the
4246// quotient field Q(R/I).
4247// - nor[3] = ring say T, containing two ideals norid and normap such that
4248// normap gives the normalization map from R / I to T / norid.
4249// - nor[4] = the delta invariant, if withDelta = 1.
4250
4251// Details:
4252// --------
4253// Computes the ideal of the minors in the first step and then reuses it in all
4254// steps.
4255// In step s, the denominator is D^s, where D is a nzd of the original quotient
4256// ring, contained in the radical of the singular locus.
4257// This denominator is used except when the degree of D^i is greater than the
4258// degree of a universal denominator.
4259// The nzd is taken as the smallest degree polynomial in the radical of the
4260// singular locus.
4261
4262// It assumes that the ideal I is equidimensional radical. This is not checked
4263// in the procedure!
4264// If decomp = 0 or decomp = 3 it assumes further that I is prime. Therefore
4265// any non-zero element in the jacobian ideal is assumed to be a
4266// non-zerodivisor.
4267
4268// It works over the given basering.
4269// If it has a non-global ordering, it changes it to dp global only for
4270// computing radical.
4271
4272// The delta invariant is only computed if withDelta = 1, and decomp = 0 or
4273// decomp = 3 (meaning that the ideal is prime).
4274
4275// denomOption = 0      -> Uses the smallest degree polynomial
4276// denomOption = i > 0  -> Uses a polynomial in the i-th variable
4277
[651cce]4278  intvec save_opt=option(get);
4279  option(redSB);
4280  option(returnSB);
[1598658]4281  int step = 0;                       // Number of steps. (for debugging)
4282  int dbg = printlevel - voice + 2;   // dbg = printlevel (default: dbg = 0)
4283  int i;                              // counter
[1e1ec4]4284  int isGlobal = attrib(basering,"global");
[1598658]4285
4286  poly c;                     // The common denominator.
4287
4288  def R = basering;
4289
4290//------------------------ Groebner bases and dimension of I-----------------
[651cce]4291  if(isGlobal == 1)
4292  {
[1598658]4293    list IM = mstd(I);
4294    I = IM[1];
4295    ideal IMin = IM[2];   // A minimal set of generators in the groebner basis.
4296  }
4297  else
4298  {
4299    // The minimal set of generators is not computed by mstd for
4300    // non-global orderings.
4301    I = groebner(I);
4302    ideal IMin = I;
4303  }
4304  int d = dim(I);
4305
4306  // ---------------- computation of the singular locus ---------------------
4307  // We compute the radical of the ideal of minors modulo the original ideal.
4308  // This is done only in the first step.
4309  qring Q = I;   // We work in the quotient by the groebner base of the ideal I
[651cce]4310  option(redSB);
4311  option(returnSB);
[1598658]4312
4313  // If a conductor ideal was given as input, we use it instead of the
4314  // singular locus. If a test ideal was given as input, we do not compute the
4315  // singular locus.
4316  ideal inputC = fetch(R, inputC);
4317  ideal inputJ = fetch(R, inputJ);
[651cce]4318  if((inputC == 0) && (inputJ == 0))
4319  {
[1598658]4320    // We compute the radical of the ideal of minors modulo the original ideal.
4321    // This is done only in the first step.
4322    ideal I = fetch(R, I);
4323    attrib(I, "isSB", 1);
4324    ideal IMin = fetch(R, IMin);
4325
4326    dbprint(dbg, "Computing the jacobian ideal...");
4327
4328    // If a given conductor ideal is given, we use it.
4329    // If a given test ideal is given, we don't need to compute the jacobian
[1e1ec4]4330
4331    // reduction mod I in 'minor' is not working for local orderings!
4332    if(attrib(basering,"global"))
4333    {
4334      ideal J = minor(jacob(IMin), nvars(basering) - d, I);
4335    }
4336    else
4337    {
4338      ideal J = minor(jacob(IMin), nvars(basering) - d);
4339      J = reduce(J, groebner(I));
4340    }
[1598658]4341    J = groebner(J);
[651cce]4342  }
4343  else
4344  {
[1598658]4345    ideal J = fetch(R, inputC);
4346    J = groebner(J);
4347  }
4348
4349  //------------------ We check if the singular locus is empty -------------
[651cce]4350  if(J[1] == 1)
4351  {
[1598658]4352    // The original ring R/I was normal. Nothing to do.
4353    // We define anyway a new ring, equal to R, to be able to return it.
4354    setring R;
4355    list lR = ringlist(R);
4356    def ROut = ring(lR);
4357    setring ROut;
4358    ideal norid = fetch(R, I);
4359    ideal normap = maxideal(1);
4360    export norid;
4361    export normap;
4362    setring R;
[651cce]4363    if(withDelta)
4364    {
[1598658]4365      list output = ideal(1), poly(1), ROut, 0;
4366    }
4367    else
4368    {
4369      list output = ideal(1), poly(1), ROut;
4370    }
[651cce]4371    option(set,save_opt);
[1598658]4372    return(list(output));
4373  }
4374
4375
4376  // -------------------- election of the universal denominator----------------
4377  // We first check if a conductor ideal was computed. If not, we don't
4378  // compute a universal denominator.
4379  ideal Id1;
[651cce]4380  if(J != 0)
4381  {
4382    if(denomOption == 0)
4383    {
[1598658]4384      poly condu = getSmallest(J);   // Choses the polynomial of smallest degree
4385                                     // of J as universal denominator.
[651cce]4386    }
4387    else
4388    {
[1598658]4389      poly condu = getOneVar(J, denomOption);
4390    }
[651cce]4391    if(dbg >= 1)
4392    {
[1598658]4393      "";
4394      "The universal denominator is ", condu;
4395    }
4396
4397    // ----- splitting the ideal by the universal denominator (if possible) -----
4398    // If the ideal is equidimensional, but not necessarily prime, we check if
4399    // the universal denominator is a non-zerodivisor of R/I.
4400    // If not, we split I.
[651cce]4401    if((decomp == 1) or (decomp == 2))
4402    {
[1598658]4403      Id1 = quotient(0, condu);
[651cce]4404      if(size(Id1) > 0)
4405      {
[1598658]4406        // We have to split.
[651cce]4407        if(dbg >= 1)
4408        {
[1598658]4409          "A zerodivisor was found. We split the ideal. The zerodivisor is ", condu;
4410        }
4411        setring R;
4412        ideal Id1 = fetch(Q, Id1), I;
4413        Id1 = groebner(Id1);
4414        ideal Id2 = quotient(I, Id1);
4415        // I = I1 \cap I2
4416        printlevel = printlevel + 1;
[c623f27]4417        ideal JDefault = 0; // Now it uses the default J;
[1598658]4418        list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault)[1];
4419        list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault)[1];
4420        printlevel = printlevel - 1;
[651cce]4421        option(set,save_opt);
[1598658]4422        return(list(nor1, nor2));
4423      }
4424    }
[651cce]4425  }
4426  else
4427  {
[1598658]4428    poly condu = 0;
4429  }
4430
4431  // --------------- computation of the first test ideal ---------------------
4432  // To compute the radical we go back to the original ring.
4433  // If we are using a non-global ordering, we must change to the global
4434  // ordering.
4435  setring R;
4436  // If a test ideal is given at the input, we use it.
[651cce]4437  if(inputJ == 0)
4438  {
4439    if(isGlobal == 1)
4440    {
[1598658]4441      ideal J = fetch(Q, J);
4442      J = J, I;
[651cce]4443      if(dbg >= 1)
4444      {
[1598658]4445        "The original singular locus is";
4446        groebner(J);
4447        if(dbg >= 2){pause();}
4448        "";
4449      }
4450      // We check if the only singular point is the origin.
4451      // If so, the radical is the maximal ideal at the origin.
4452      J = groebner(J);
[651cce]4453      if(locAtZero(J))
4454      {
[1598658]4455        J = maxideal(1);
[651cce]4456      }
4457      else
4458      {
[1598658]4459        J = radical(J);
4460      }
[651cce]4461    }
4462    else
4463    {
[1598658]4464      // We change to global dp ordering.
4465      list rl = ringlist(R);
4466      list origOrd = rl[3];
4467      list newOrd = list("dp", intvec(1:nvars(R))), list("C", 0);
4468      rl[3] = newOrd;
4469      def globR = ring(rl);
4470      setring globR;
4471      ideal J = fetch(Q, J);
4472      ideal I = fetch(R, I);
4473      J = J, I;
[651cce]4474      if(dbg >= 1)
4475      {
[1598658]4476        "The original singular locus is";
4477        groebner(J);
4478        if(dbg>=2){pause();}
4479        "";
4480      }
4481      J = radical(J);
4482      setring R;
4483      ideal J = fetch(globR, J);
4484    }
[651cce]4485  }
4486  else
4487  {
[1598658]4488    ideal J = inputJ;
4489  }
4490
[651cce]4491  if(dbg >= 1)
4492  {
[1598658]4493    "The radical of the original singular locus is";