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
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="version normal.lib 4.0.0.0 Jun_2013 "; // $Id$
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
11PROCEDURES:
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
27
28SEE ALSO: locnormal_lib;modnormal_lib
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:@*
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\".)@*
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
72         ideal inputJ. This option is only for use in other procedures. Using
73         this option, the result might not be the normalization.@*
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
77         this option, the result might not be the normalization.@*
78         (Option only valid for global algorithm.)@*
79         Options used for computing integral basis (over rings of two
80         variables):@*
81         - \"var1\" -> uses a polynomial in the first variable as
82         universal denominator.@*
83         - \"var2\" -> uses a polynomial in the second variable as universal
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.
178  denomOption = 0;   // The default universal denominator is the smallest
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
249  int isGlobal = attrib(basering,"global");// Checks if the original ring has
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];
382        if(!defined(savebasering)) { def savebasering;}
383        savebasering=basering;
384        setring newR; // must be in a compatible ring to newR
385                      // as ringlist may produce ring-dep. stuff
386        if(!defined(newRList)) { list newRList;}
387        newRList = ringlist(newR);
388        setring savebasering;
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);
410          newRListO = imap(newR, newRList);
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       }
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       }
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   {
2174     def R1=changeord(list(list("dp",1:nvars(basering))));
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
2365   int genus=(d-1)*(d-2) div 2;
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
2437      list nor=normal(I,"withDelta");
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;
2494       ideal I1=imap(newR,radIDsing);
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;
2507          if (inv[2]!=0) { nsing++; }
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;
2515           list pr=minAssGTZ(radIDsing);
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{
2590   intvec save_opt=option(get);
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;
2684      option(set,save_opt);
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;
2693         delt=(tau+1) div 2;
2694         option(set,save_opt);
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;
2702         delt=(tau+2) div 2;
2703         option(set,save_opt);
2704         return(list(d*delt,d*tau,d*(2*delt-tau+1)));
2705      }
2706
2707      int mu=vdim(std(jacob(f)));
2708
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";"";}
2723        option(set,save_opt);
2724        return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
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);
2737         //int st=w_deg(hc,v)+v[1]+v[2];
2738         int st=w_deg(var(1)^NP[size(NP)][1],v)+1;
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];
2756
2757         setring R;
2758         option(set,save_opt);
2759         return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
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;
2773         option(set,save_opt);
2774         return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
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      }
2788      option(set,save_opt);
2789      return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
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);
2796      option(set,save_opt);
2797      return(list(d*delt,d*tau,d));
2798   }
2799   option(set,save_opt);
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///////////////////////////////////////////////////////////////////////////////
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)
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
3546   if ( attrib(basering,"global") != 1)
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{
4021  int i, ii,jj;
4022  intvec V;                          // to be used for variable weights
4023  int y = printlevel-voice+2;
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.
4168      if(attrib(basering,"global") != 1)
4169      {
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
4278  intvec save_opt=option(get);
4279  option(redSB);
4280  option(returnSB);
4281  int step = 0;                       // Number of steps. (for debugging)
4282  int dbg = printlevel - voice + 2;   // dbg = printlevel (default: dbg = 0)
4283  int i;                              // counter
4284  int isGlobal = attrib(basering,"global");
4285
4286  poly c;                     // The common denominator.
4287
4288  def R = basering;
4289
4290//------------------------ Groebner bases and dimension of I-----------------
4291  if(isGlobal == 1)
4292  {
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
4310  option(redSB);
4311  option(returnSB);
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);
4318  if((inputC == 0) && (inputJ == 0))
4319  {
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
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    }
4341    J = groebner(J);
4342  }
4343  else
4344  {
4345    ideal J = fetch(R, inputC);
4346    J = groebner(J);
4347  }
4348
4349  //------------------ We check if the singular locus is empty -------------
4350  if(J[1] == 1)
4351  {
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;
4363    if(withDelta)
4364    {
4365      list output = ideal(1), poly(1), ROut, 0;
4366    }
4367    else
4368    {
4369      list output = ideal(1), poly(1), ROut;
4370    }
4371    option(set,save_opt);
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;
4380  if(J != 0)
4381  {
4382    if(denomOption == 0)
4383    {
4384      poly condu = getSmallest(J);   // Choses the polynomial of smallest degree
4385                                     // of J as universal denominator.
4386    }
4387    else
4388    {
4389      poly condu = getOneVar(J, denomOption);
4390    }
4391    if(dbg >= 1)
4392    {
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.
4401    if((decomp == 1) or (decomp == 2))
4402    {
4403      Id1 = quotient(0, condu);
4404      if(size(Id1) > 0)
4405      {
4406        // We have to split.
4407        if(dbg >= 1)
4408        {
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;
4417        ideal JDefault = 0; // Now it uses the default J;
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;
4421        option(set,save_opt);
4422        return(list(nor1, nor2));
4423      }
4424    }
4425  }
4426  else
4427  {
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.
4437  if(inputJ == 0)
4438  {
4439    if(isGlobal == 1)
4440    {
4441      ideal J = fetch(Q, J);
4442      J = J, I;
4443      if(dbg >= 1)
4444      {
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);
4453      if(locAtZero(J))
4454      {
4455        J = maxideal(1);
4456      }
4457      else
4458      {
4459        J = radical(J);
4460      }
4461    }
4462    else
4463    {
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;
4474      if(dbg >= 1)
4475      {
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    }
4485  }
4486  else
4487  {
4488    ideal J = inputJ;
4489  }
4490
4491  if(dbg >= 1)
4492  {
4493    "The radical of the original singular locus is";
4494    J;
4495    if(dbg>=2){pause();}
4496  }
4497
4498  // ---------------- election of the non zero divisor ---------------------
4499  setring Q;
4500  J = fetch(R, J);
4501  J = interred(J);
4502  if(denomOption == 0)
4503  {
4504    poly D = getSmallest(J);    // Chooses the polynomial of smallest degree as
4505                                // non-zerodivisor.
4506  }
4507  else
4508  {
4509    poly D = getOneVar(J, denomOption);
4510  }
4511  if(dbg >= 1)
4512  {
4513    "The non zero divisor is ", D;
4514    "";
4515  }
4516
4517  // ------- splitting the ideal by the non-zerodivisor (if possible) --------
4518  // If the ideal is equidimensional, but not necessarily prime, we check if D
4519  // is actually a non-zerodivisor of R/I.
4520  // If not, we split I.
4521  if((decomp == 1) or (decomp == 2))
4522  {
4523    // We check if D is actually a non-zerodivisor of R/I.
4524    // If not, we split I.
4525    Id1 = quotient(0, D);
4526    if(size(Id1) > 0)
4527    {
4528      // We have to split.
4529      if(dbg >= 1)
4530      {
4531        "A zerodivisor was found. We split the ideal. The zerodivisor is ", D;
4532      }
4533      setring R;
4534      ideal Id1 = fetch(Q, Id1), I;
4535      Id1 = groebner(Id1);
4536      ideal Id2 = quotient(I, Id1);
4537      // I = Id1 \cap Id2
4538      printlevel = printlevel + 1;
4539
4540      ideal JDefault = 0;  // Now it uses the default J;
4541      list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault)[1];
4542      list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault)[1];
4543      printlevel = printlevel - 1;
4544      option(set,save_opt);
4545      return(list(nor1, nor2));
4546    }
4547  }
4548
4549  // --------------------- normalization ------------------------------------
4550  // We call normalMEqui to compute the normalization.
4551  setring R;
4552  poly D = fetch(Q, D);
4553  poly condu = fetch(Q, condu);
4554  J = fetch(Q, J);
4555  printlevel = printlevel + 1;
4556  list result = normalMEqui(I, J, condu, D, withDelta, denomOption);
4557  printlevel = printlevel - 1;
4558  option(set,save_opt);
4559  return(list(result));
4560}
4561
4562///////////////////////////////////////////////////////////////////////////////
4563
4564static proc normalMEqui(ideal I, ideal origJ, poly condu, poly D, int withDelta)
4565// Here is where the normalization is actually computed.
4566
4567// Computes the normalization of R/I. (basering is R)
4568// I is assumed to be radical and equidimensional.
4569// origJ is the first test ideal.
4570// D is a non-zerodivisor of R/I.
4571// condu is a non-zerodivisor in the conductor or 0 if it was not computed.
4572// If withDelta = 1, computes the delta invariant.
4573{
4574  int step = 0;                       // Number of steps. (for debugging)
4575  int dbg = printlevel - voice + 2;   // dbg = printlevel (default: dbg = 0)
4576  int i;                              // counter
4577  int isNormal = 0;                   // check for exiting the loop
4578  int isGlobal = attrib(basering,"global");
4579  int delt;
4580
4581  def R = basering;
4582  poly c = D;
4583  ideal U;
4584  ideal cJ;
4585  list testOut;                 // Output of proc testIdeal
4586                                // (the test ideal and the ring structure)
4587
4588  qring Q = groebner(I);
4589  intvec save_opt=option(get);
4590  option(redSB);
4591  option(returnSB);
4592  ideal J = imap(R, origJ);
4593  poly c = imap(R, c);
4594  poly D = imap(R, D);
4595  poly condu = imap(R, condu);
4596  ideal cJ;
4597  ideal cJMod;
4598
4599  dbprint(dbg, "Preliminar step begins.");
4600
4601  // --------------------- computation of A1 -------------------------------
4602  dbprint(dbg, "Computing the quotient (DJ : J)...");
4603  ideal U = groebner(quotient(D*J, J));
4604  ideal oldU = 1;
4605
4606  if(dbg >= 2) { "The quotient is"; U; }
4607
4608  // ----------------- Grauer-Remmert criterion check -----------------------
4609  // We check if the equality in Grauert - Remmert criterion is satisfied.
4610  isNormal = checkInclusions(D*oldU, U);
4611  if(isNormal == 0)
4612  {
4613    if(dbg >= 1)
4614    {
4615      "In this step, we have the ring 1/c * U, with c =", c;
4616      "and U = "; U;
4617    }
4618  }
4619  else
4620  {
4621    // The original ring R/I was normal. Nothing to do.
4622    // We define anyway a new ring, equal to R, to be able to return it.
4623    setring R;
4624    list lR = ringlist(R);
4625    def ROut = ring(lR);
4626    setring ROut;
4627    ideal norid = fetch(R, I);
4628    ideal normap = maxideal(1);
4629    export norid;
4630    export normap;
4631    setring R;
4632    if(withDelta)
4633    {
4634      list output = ideal(1), poly(1), ROut, 0;
4635    }
4636    else
4637    {
4638      list output = ideal(1), poly(1), ROut;
4639    }
4640    option(set,save_opt);
4641    return(output);
4642  }
4643
4644  // ----- computation of the chain of ideals A1 c A2 c ... c An ------------
4645  while(isNormal == 0)
4646  {
4647    step++;
4648    if(dbg >= 1) { ""; "Step ", step, " begins."; }
4649    dbprint(dbg, "Computing the test ideal...");
4650
4651    // --------------- computation of the test ideal ------------------------
4652    // Computes a test ideal for the new ring.
4653    // The test ideal will be the radical in the new ring of the original
4654    // test ideal.
4655    setring R;
4656    U = imap(Q, U);
4657    c = imap(Q, c);
4658    testOut = testIdeal(I, U, origJ, c, D);
4659    cJ = testOut[1];
4660
4661    setring Q;
4662    cJ = imap(R, cJ);
4663    cJ = groebner(cJ);
4664
4665    // cJ / c is now the ideal mapped back.
4666    // We have the generators as an ideal in the new ring,
4667    // but we want the generators as an ideal in the original ring.
4668    cJMod = getGenerators(cJ, U, c);
4669
4670    if(dbg >= 2) { "The test ideal in this step is "; cJMod; }
4671
4672    cJ = cJMod;
4673
4674    // ------------- computation of the quotient (DJ : J)--------------------
4675    oldU = U;
4676    dbprint(dbg, "Computing the quotient (c*D*cJ : cJ)...");
4677    U = quotient(c*D*cJ, cJ);
4678    if(dbg >= 2){"The quotient is "; U;}
4679
4680    // ------------- Grauert - Remmert criterion check ----------------------
4681    // We check if the equality in Grauert - Remmert criterion is satisfied.
4682    isNormal = checkInclusions(D*oldU, U);
4683
4684    if(isNormal == 1)
4685    {
4686      // We go one step back. In the last step we didnt get antyhing new,
4687      // we just verified that the ring was already normal.
4688      dbprint(dbg, "The ring in the previous step was already normal.");
4689      dbprint(dbg, "");
4690      U = oldU;
4691    }
4692    else
4693    {
4694      // ------------- preparation for next iteration ----------------------
4695      // We have to go on.
4696      // The new denominator is chosen.
4697      c = D * c;
4698
4699      // If we have a universal denominator of smaller degree than c,
4700      // we replace c by it.
4701      if(condu != 0)
4702      {
4703        if(deg(c) > deg(condu))
4704        {
4705          U = changeDenominatorQ(U, c, condu);
4706          c = condu;
4707        }
4708      }
4709      if(dbg >= 1)
4710      {
4711        "In this step, we have the ring 1/c * U, with c =", c;
4712        "and U = ";
4713        U;
4714        if(dbg>=2){pause();}
4715      }
4716    }
4717  }
4718
4719  // ------------------------- delta computation ----------------------------
4720  if(withDelta)
4721  {
4722    ideal UD = groebner(U);
4723    delt = vdim(std(modulo(UD, c)));
4724  }
4725
4726  // -------------------------- prepare output -----------------------------
4727  setring R;
4728  U = fetch(Q, U);
4729  c = fetch(Q, c);
4730
4731  // Ring structure of the new ring
4732  def ere = testOut[2];
4733  if(withDelta)
4734  {
4735    list output = U, c, ere, delt;
4736  }
4737  else
4738  {
4739    list output = U, c, ere;
4740  }
4741  option(set,save_opt);
4742  return(output);
4743}
4744
4745///////////////////////////////////////////////////////////////////////////////
4746
4747static proc lineUpLast(ideal U, poly c)
4748// Sets c as the last generator of U.
4749{
4750  int i;
4751  ideal newU;
4752  for (i = 1; i <= ncols(U); i++)
4753  {
4754    if(U[i] != c)
4755    {
4756      if(size(newU) == 0)
4757      { newU = U[i]; }
4758      else
4759      { newU = newU, U[i]; }
4760    }
4761  }
4762  if(size(newU) == 0)
4763  { newU = c; }
4764  else
4765  { newU = newU, c; }
4766  return(newU);
4767}
4768
4769///////////////////////////////////////////////////////////////////////////////
4770
4771static proc lineUp(ideal U, poly c)
4772// Sets c as the first generator of U.
4773{
4774  int i;
4775  ideal newU = c;
4776  for (i = 1; i <= ncols(U); i++)
4777  {
4778    if(U[i] != c)
4779    {
4780      newU = newU, U[i];
4781    }
4782  }
4783  return(newU);
4784}
4785
4786///////////////////////////////////////////////////////////////////////////////
4787
4788//WARNING - elim is not working here!! Check!!
4789//It is now replaced by computing an eliminating groebner basis.
4790proc getOneVar(ideal J, int vari)
4791"USAGE:   getOneVar(J, vari); J is a 0-dimensional ideal, vari is an integer.
4792RETURN:  a polynomial of J in the variable indicated by vari of smallest
4793         degree.@*
4794NOTE:    Works only over rings of two variables.@*
4795         It is intended mainly as an auxiliary procedure for computing
4796         integral bases. @*
4797EXAMPLE: example getOneVar; shows an example
4798"
4799{
4800  def R = basering;
4801  list RL = ringlist(R);
4802  // We keep everything from R but we change the ordering to lp, and we
4803  // order the variables as needed.
4804  RL[3] = list(list("lp", 1:2), list("C", 0:1));
4805  RL[2] = list(var(3-vari), var(vari));
4806  RL[4]=ideal(0); // does not work with qrings: Ex.7 of paraplanecurves
4807  def RR = ring(RL);
4808  setring RR;
4809  ideal J = imap(R, J);
4810  J = groebner(J);
4811  poly g = J[1];
4812  setring R;
4813  poly g = imap(RR, g);
4814  return(g);
4815}
4816example
4817{ "EXAMPLE:";
4818  printlevel = printlevel+1;
4819  echo = 2;
4820  ring s = 0,(x,y),dp;
4821  ideal J = x3-y, y3;
4822  getOneVar(J, 1);
4823
4824  echo = 0;
4825  printlevel = printlevel-1;
4826}
4827///////////////////////////////////////////////////////////////////////////////
4828
4829proc getSmallest(ideal J)
4830"USAGE:   getSmallest(J); J is an ideal.
4831RETURN:  the generator of J of smallest degree. If there are more than one, it
4832         chooses the one with smallest number of monomials.@*
4833NOTE:    It looks only at the generator of J, not at all the polynomials in
4834         the ideal.@*
4835         It is intended maninly to compute a good universal denominator in the
4836         normalization algorithms.@*
4837EXAMPLE: example getSmallest; shows an example
4838"
4839{
4840// Computes the polynomial of smallest degree of J.
4841//
4842  int i;
4843  poly p = J[1];
4844  int d = deg(p);
4845  int di;
4846  for(i = 2; i <= ncols(J); i++)
4847  {
4848    if(J[i] != 0)
4849    {
4850      di = deg(J[i]);
4851      if(di < d)
4852      {
4853        p = J[i];
4854        d = di;
4855      }
4856      else
4857      {
4858        if(di == d)
4859        {
4860          if(size(J[i]) < size(p))
4861          {
4862            p = J[i];
4863          }
4864        }
4865      }
4866    }
4867  }
4868  return(p);
4869}
4870example
4871{ "EXAMPLE:";
4872  printlevel = printlevel+1;
4873  echo = 2;
4874  ring s = 0,(x,y),dp;
4875  ideal J = x3-y, y5, x2-y2+1;
4876  getSmallest(J);
4877
4878  echo = 0;
4879  printlevel = printlevel-1;
4880}
4881
4882///////////////////////////////////////////////////////////////////////////////
4883
4884static proc getGenerators(ideal J, ideal U, poly c)
4885{
4886// Computes the generators of J as an ideal in the original ring,
4887// where J is given by generators in the new ring.
4888
4889// The new ring is given by 1/c * U in the total ring of fractions.
4890
4891  int i, j;                             // counters;
4892  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)
4893  poly p;                               // The lifted polynomial
4894  ideal JGr = groebner(J);              // Groebner base of J
4895
4896  if(dbg>1){"Checking for new generators...";}
4897  for(i = 1; i <= ncols(J); i++)
4898  {
4899    for(j = 1; j <= ncols(U); j++)
4900    {
4901      p = lift(c, J[i]*U[j])[1,1];
4902      p = reduce(p, JGr);
4903      if(p != 0)
4904      {
4905        if(dbg>1)
4906        {
4907          "New polynoial added:", p;
4908          if(dbg>4) {pause();}
4909        }
4910        JGr = JGr, p;
4911        JGr = groebner(JGr);
4912        J = J, p;
4913      }
4914    }
4915  }
4916  return(J);
4917}
4918
4919///////////////////////////////////////////////////////////////////////////////
4920
4921static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D)
4922{
4923// Internal procedure, used in normalM.
4924// Computes the test ideal in the new ring.
4925// It takes the original test ideal and computes the radical of it in the
4926// new ring.
4927
4928// The new ring is 1/c * U.
4929// The original test ideal is origJ.
4930// The original ring is R / I, where R is the basering.
4931  int i;                                // counter
4932  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)
4933  def R = basering;                      // We dont work in the quo
4934  ideal J = origJ;
4935
4936  // ---------- computation of the ring structure of 1/c * U ----------------
4937  U = lineUp(U, c);
4938
4939  if(dbg > 1){"Computing the new ring structure...";}
4940  list ele = computeRing(U, I, "noRed");
4941
4942  def origEre = ele[1];
4943  setring origEre;
4944  if(dbg > 1){"The relations are"; norid;}
4945
4946  // ---------------- setting the ring to work in  --------------------------
4947  int isGlobal = attrib(origEre,"global");      // Checks if the original ring has
4948                                         // global ordering.
4949  if(isGlobal != 1)
4950  {
4951    list rl = ringlist(origEre);
4952    list origOrd = rl[3];
4953    list newOrd = list("dp", intvec(1:nvars(origEre))), list("C", 0);
4954    rl[3] = newOrd;
4955    def ere = ring(rl);     // globR is the original ring but
4956                            // with a global ordering.
4957    setring ere;
4958    ideal norid = imap(origEre, norid);
4959  }
4960  else
4961  {
4962    def ere = origEre;
4963  }
4964
4965  ideal I = imap(R, I);
4966  ideal J = imap(R, J);
4967  J = J, norid, I;
4968
4969
4970  // ----- computation of the test ideal using the ring structure of Ai -----
4971  intvec save_opt=option(get);
4972  option(redSB);
4973  option(returnSB);
4974
4975  if(dbg > 1){"Computing the radical of J...";}
4976  J = radical(J);
4977  if(dbg > 1){"Computing the interreduction of the radical...";}
4978  J = groebner(J);
4979  //J = interred(J);
4980  if(dbg > 1)
4981  {
4982    "The radical in the generated ring is";
4983    J;
4984    if(dbg>4){pause();}
4985  }
4986
4987  setring ere;
4988
4989  // -------------- map from Ai to the total ring of fractions ---------------
4990  // Now we must map back this ideal J to U_i / c in the total ring of
4991  // fractions.
4992  // The map sends T_j -> u_j / c.
4993  // The map is built by the following steps:
4994  // 1) We compute the degree of the generators of J with respect to the
4995  //    new variables T_j.
4996  // 2) For each generator, we multiply each term by a power of c, as if
4997  //    taking c^n as a common denominator (considering the new variables as
4998  //    a polynomial in the old variables divided by c).
4999  // 3) We replace the new variables T_j by the corresponding numerator u_j.
5000  // 4) We lift the resulting polynomial to change the denominator
5001  //    from c^n to c.
5002  int nNewVars = nvars(ere) - nvars(R);      // Number of new variables
5003  poly c = imap(R, c);
5004  intvec @v = 1..nNewVars;    // Vector of the new variables.
5005                              // They must be the first ones.
5006  if(dbg > 1){"The indices of the new variables are", @v;}
5007
5008  // ---------------------- step 1 of the mapping ---------------------------
5009  intvec degs;
5010  for(i = 1; i<=ncols(J); i++)
5011  {
5012    degs[i] = degSubring(J[i], @v);
5013  }
5014  if(dbg > 1)
5015  {
5016    "The degrees with respect to the new variables are";
5017    degs;
5018  }
5019
5020  // ---------------------- step 2 of the mapping ---------------------------
5021  ideal mapJ = mapBackIdeal(J, c, @v);
5022
5023  setring R;
5024
5025  // ---------------------- step 3 of the mapping ---------------------------
5026  ideal z;                    // The variables of the original ring in order.
5027  for(i = 1; i<=nvars(R); i++)
5028  {
5029    z[i] = var(i);
5030  }
5031
5032  map f = ere, U[2..ncols(U)], z[1..ncols(z)]; // The map to the original ring.
5033  if(dbg > 1)
5034  {
5035    "The map is ";
5036    f;
5037    if(dbg>4){pause();}
5038  }
5039
5040  if(dbg > 1){ "Computing the map..."; }
5041
5042  J = f(mapJ);
5043  if(dbg > 1)
5044  {
5045    "The ideal J mapped back (before lifting) is";
5046    J;
5047    if(dbg>4){pause();}
5048  }
5049
5050  // ---------------------- step 4 of the mapping ---------------------------
5051  qring Q = groebner(I);
5052  ideal J = imap(R, J);
5053  poly c = imap(R, c);
5054  for(i = 1; i<=ncols(J); i++)
5055  {
5056    if(degs[i]>1)
5057    {
5058      J[i] = lift(c^(degs[i]-1), J[i])[1,1];
5059    }
5060    else
5061    {
5062      if(degs[i]==0) { J[i] = c*J[i]; }
5063    }
5064  }
5065
5066  if(dbg > 1)
5067  {
5068    "The ideal J lifted is";
5069    J;
5070    if(dbg>4){pause();}
5071  }
5072
5073  // --------------------------- prepare output ----------------------------
5074  J = groebner(J);
5075
5076  setring R;
5077  J = imap(Q, J);
5078
5079  option(set,save_opt);
5080  return(list(J, ele[1]));
5081}
5082
5083///////////////////////////////////////////////////////////////////////////////
5084
5085proc changeDenominator(ideal U1, poly c1, poly c2, ideal I)
5086"USAGE:   changeDenominator(U1, c1, c2, I); U1 and I ideals, c1 and c2
5087         polynomials.@*
5088RETURN:  an ideal U2 such that the A-modules 1/c1 * U1 and 1/c2 * U2 are equal,
5089         where A = R/I and R is the basering.@*
5090NOTE:    It assumes that such U2 exists. It is intended maninly as an auxiliary
5091         procedure in the normalization algorithms.@*
5092EXAMPLE: example changeDenominator; shows an example
5093"
5094{
5095// Let A = R / I. Given an A-module in the form 1/c1 * U1 (U1 ideal of A), it
5096// computes a new ideal U2 such that the the A-module is 1/c2 * U2.
5097// The base ring is R, but the computations are to be done in R / I.
5098  int a;      // counter
5099  def R = basering;
5100  qring Q = I;
5101  ideal U1 = fetch(R, U1);
5102  poly c1 = fetch(R, c1);
5103  poly c2 = fetch(R, c2);
5104  ideal U2 = changeDenominatorQ(U1, c1, c2);
5105  setring R;
5106  ideal U2 = fetch(Q, U2);
5107  return(U2);
5108}
5109example
5110{
5111  "EXAMPLE:";
5112  echo = 2;
5113  ring s = 0,(x,y),dp;
5114  ideal I = y5-y4x+4y2x2-x4;
5115  ideal U1 = normal(I)[2][1];
5116  poly c1 = U1[4];
5117  U1;c1;
5118  // 1/c1 * U1 is the normalization of I.
5119  ideal U2 = changeDenominator(U1, c1, x3, I);
5120  U2;
5121  // 1/x3 * U2 is also the normalization of I, but with a different denominator.
5122  echo = 0;
5123}
5124
5125///////////////////////////////////////////////////////////////////////////////
5126
5127static proc changeDenominatorQ(ideal U1, poly c1, poly c2)
5128{
5129// Given a ring in the form 1/c1 * U, it computes a new U2 st the ring
5130// is 1/c2 * U2.
5131// The base ring is already a quotient ring R / I.
5132  int a;      // counter
5133  ideal U2;
5134  poly p;
5135  for(a = 1; a <= ncols(U1); a++)
5136  {
5137    p = lift(c1, c2*U1[a])[1,1];
5138    U2[a] = p;
5139  }
5140  return(U2);
5141}
5142
5143///////////////////////////////////////////////////////////////////////////////
5144
5145static proc checkInclusions(ideal U1, ideal U2)
5146{
5147// Checks if the identity A = Hom(J, J) of Grauert-Remmert criterion is
5148// satisfied.
5149  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)
5150  list reduction1;
5151  list reduction2;
5152
5153  // ---------------------- inclusion Hom(J, J) c A -------------------------
5154  if(dbg > 1){"Checking the inclusion Hom(J, J) c A:";}
5155  // This interred is used only because a bug in groebner!
5156  U1 = groebner(U1);
5157  reduction1 = reduce(U2, U1);
5158  if(dbg > 1){reduction1[1];}
5159
5160  // ---------------------- inclusion A c Hom(J, J) -------------------------
5161  // The following check should always be satisfied.
5162  // This is only used for debugging.
5163  if(dbg > 1)
5164  {
5165    "and the inclusion A c Hom(J, J): (this should always be satisfied)";
5166    // This interred is used only because a bug in groebner!
5167    U2 = groebner(U2);
5168    reduction2 = reduce(U1, groebner(U2));
5169    reduction2[1];
5170    if(size(reduction2[1]) > 0)
5171    {
5172      "Something went wrong... (this inclusion should always be satisfied)";
5173      ~;
5174    }
5175    else
5176    {
5177      if(dbg>4){pause();}
5178    }
5179  }
5180
5181  if(size(reduction1[1]) == 0)
5182  {
5183    // We are done! The ring computed in the last step was normal.
5184    return(1);
5185  }
5186  else
5187  {
5188    return(0);
5189  }
5190}
5191
5192///////////////////////////////////////////////////////////////////////////////
5193
5194static proc degSubring(poly p, intvec @v)
5195{
5196// Computes the degree of a polynomial taking only some variables as variables
5197// and the others as parameters.
5198
5199// The degree is taken w.r.t. the variables indicated in v.
5200  int i;      // Counter
5201  int d = 0;  // The degree
5202  int e;      // Degree (auxiliar variable)
5203  for(i = 1; i <= size(p); i++)
5204  {
5205    e = sum(leadexp(p[i]), @v);
5206    if(e > d){d = e;}
5207  }
5208  return(d);
5209}
5210
5211///////////////////////////////////////////////////////////////////////////////
5212
5213static proc mapBackIdeal(ideal I, poly c, intvec @v)
5214{
5215// Modifies all polynomials in I so that a map x(i) -> y(i)/c can be
5216// carried out.
5217
5218// v indicates wicih variables x(i) of the ring will be mapped to y(i)/c.
5219
5220  int i;  // counter
5221  for(i = 1; i <= ncols(I); i++)
5222  {
5223    I[i] = mapBackPoly(I[i], c, @v);
5224  }
5225  return(I);
5226}
5227
5228///////////////////////////////////////////////////////////////////////////////
5229
5230static proc mapBackPoly(poly p, poly c, intvec @v)
5231{
5232// Multiplies each monomial of p by a power of c so that a map x(i) -> y(i)/c
5233// can be carried out.
5234
5235// v indicates wicih variables x(i) of the ring will be mapped to y(i)/c.
5236  int i;  // counter
5237  int e;  // exponent
5238  int d = degSubring(p, @v);
5239  poly g = 0;
5240  int size_p=size(p);
5241  for(i = 1; i <= size_p; i++)
5242  {
5243    e = sum(leadexp(p[i]), @v);
5244    g = g + p[i] * c^(d-e);
5245  }
5246  return(g);
5247}
5248
5249//                    End procedures for normal
5250///////////////////////////////////////////////////////////////////////////////
5251
5252
5253///////////////////////////////////////////////////////////////////////////////
5254//                  Begin procedures for normalC
5255
5256// We first define resp. copy some attributes to be used in proc normal and
5257// static proc normalizationPrimes, and ..., to speed up computation in
5258// special cases
5259//NOTE:  We use the following attributes:
5260// 1     attrib(id,"isCohenMacaulay");         //--- Cohen Macaulay
5261// 2     attrib(id,"isCompleteIntersection");  //--- complete intersection
5262// 3     attrib(id,"isHypersurface");          //--- hypersurface
5263// 4     attrib(id,"isEquidimensional");       //--- equidimensional ideal
5264// 5     attrib(id,"isPrim");                  //--- prime ideal
5265// 6     attrib(id,"isRegInCodim2");           //--- regular in codimension 2
5266// 7     attrib(id,"isIsolatedSingularity");   //--- isolated singularities
5267// 8     attrib(id,"onlySingularAtZero");      //--- only singular at 0
5268// 9     attrib(id,"isRadical");               //--- radical ideal
5269//Recall: (attrib(id,"xy"),1) sets attrib xy to TRUE and
5270//        (attrib(id,"xy"),0) to FALSE
5271
5272static proc getAttrib (ideal id)
5273"USAGE:   getAttrib(id);  id=ideal
5274COMPUTE: check attributes for id. If the attributes above are defined,
5275         take its value, otherwise define it and set it to 0
5276RETURN:  intvec of size 9, with entries 0 or 1,  values of attributes defined
5277         above (in this order)
5278EXAMPLE: no example
5279"
5280{
5281  int isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad;
5282
5283  if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
5284  {
5285    if( attrib(id,"isCohenMacaulay")==1 )
5286    { isCoM=1; isEq=1; }
5287  }
5288
5289  if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
5290  {
5291    if(attrib(id,"isCompleteIntersection")==1)
5292    { isCoI=1; isCoM=1; isEq=1; }
5293  }
5294
5295  if( typeof(attrib(id,"isHypersurface"))=="int" )
5296  {
5297    if(attrib(id,"isHypersurface")==1)
5298    { isHy=1; isCoI=1; isCoM=1; isEq=1; }
5299  }
5300
5301  if( typeof(attrib(id,"isEquidimensional"))=="int" )
5302  {
5303    if(attrib(id,"isEquidimensional")==1)
5304    { isEq=1; }
5305  }
5306
5307  if( typeof(attrib(id,"isPrim"))=="int" )
5308  {
5309    if(attrib(id,"isPrim")==1)
5310    { isPr=1; }
5311  }
5312
5313  if( typeof(attrib(id,"isRegInCodim2"))=="int" )
5314  {
5315    if(attrib(id,"isRegInCodim2")==1)
5316    { isReg=1; }
5317  }
5318
5319  if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
5320  {
5321    if(attrib(id,"isIsolatedSingularity")==1)
5322    { isIso=1; }
5323  }
5324
5325  if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
5326  {
5327    if(attrib(id,"onlySingularAtZero")==1)
5328    { oSAZ=1; }
5329  }
5330
5331  if( typeof(attrib(id,"isRad"))=="int" )
5332  {
5333    if(attrib(id,"isRad")==1)
5334    { isRad=1; }
5335  }
5336
5337  intvec atr = isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad;
5338  return(atr);
5339}
5340
5341///////////////////////////////////////////////////////////////////////////////
5342
5343static proc setAttrib (ideal id, intvec atr)
5344"USAGE:   setAttrib(id,atr);  id ideal, atr intvec
5345COMPUTE: set attributes to id specified by atr
5346RETURN:  id, with assigned attributes from atr
5347EXAMPLE: no example
5348"
5349{
5350  attrib(id,"isCohenMacaulay",atr[1]);         //--- Cohen Macaulay
5351  attrib(id,"isCompleteIntersection",atr[2]);  //--- complete intersection
5352  attrib(id,"isHypersurface",atr[3]);          //--- hypersurface
5353  attrib(id,"isEquidimensional",atr[4]);       //--- equidimensional ideal
5354  attrib(id,"isPrim",atr[5]);                  //--- prime ideal
5355  attrib(id,"isRegInCodim2",atr[6]);           //--- regular in codimension 2
5356  attrib(id,"isIsolatedSingularity",atr[7]);   //--- isolated singularities
5357  attrib(id,"onlySingularAtZero",atr[8]);      //--- only singular at 0
5358  attrib(id,"isRadical",atr[9]);               //--- radical ideal
5359
5360  return(id);
5361}
5362
5363///////////////////////////////////////////////////////////////////////////////
5364// copyAttribs is not used anywhere so far
5365
5366static proc copyAttribs (ideal id1, ideal id)
5367"USAGE:   copyAttribs(id1,id);  id1, id ideals
5368COMPUTE: copy attributes from id1 to id
5369RETURN:  id, with assigned attributes from id1
5370EXAMPLE: no example
5371"
5372{
5373  if( typeof(attrib(id1,"isCohenMacaulay"))=="int" )
5374  {
5375    if( attrib(id1,"isCohenMacaulay")==1 )
5376    {
5377      attrib(id,"isEquidimensional",1);
5378    }
5379  }
5380  else
5381  {
5382    attrib(id,"isCohenMacaulay",0);
5383  }
5384
5385  if( typeof(attrib(id1,"isCompleteIntersection"))=="int" )
5386  {
5387    if(attrib(id1,"isCompleteIntersection")==1)
5388    {
5389      attrib(id,"isCohenMacaulay",1);
5390      attrib(id,"isEquidimensional",1);
5391    }
5392  }
5393  else
5394  {
5395    attrib(id,"isCompleteIntersection",0);
5396  }
5397
5398  if( typeof(attrib(id1,"isHypersurface"))=="int" )
5399  {
5400    if(attrib(id1,"isHypersurface")==1)
5401    {
5402      attrib(id,"isCompleteIntersection",1);
5403      attrib(id,"isCohenMacaulay",1);
5404      attrib(id,"isEquidimensional",1);
5405    }
5406  }
5407  else
5408  {
5409    attrib(id,"isHypersurface",0);
5410  }
5411
5412  if( (typeof(attrib(id1,"isEquidimensional"))=="int") )
5413  {
5414    if(attrib(id1,"isEquidimensional")==1)
5415    {
5416      attrib(id,"isEquidimensional",1);
5417    }
5418  }
5419  else
5420  {
5421    attrib(id,"isEquidimensional",0);
5422  }
5423
5424  if( typeof(attrib(id1,"isPrim"))=="int" )
5425  {
5426    if(attrib(id1,"isPrim")==1)
5427    {
5428      attrib(id,"isEquidimensional",1);
5429    }
5430  }
5431  else
5432  {
5433    attrib(id,"isPrim",0);
5434  }
5435
5436  if( (typeof(attrib(id1,"isRegInCodim2"))=="int") )
5437  {
5438    if(attrib(id1,"isRegInCodim2")==1)
5439    {
5440      attrib(id,"isRegInCodim2",1);
5441    }
5442  }
5443  else
5444  {
5445    attrib(id,"isRegInCodim2",0);
5446  }
5447
5448  if( (typeof(attrib(id1,"isIsolatedSingularity"))=="int") )
5449  {
5450    if(attrib(id1,"isIsolatedSingularity")==1)
5451    {
5452      attrib(id,"isIsolatedSingularity",1);
5453    }
5454  }
5455  else
5456  {
5457    attrib(id,"isIsolatedSingularity",0);
5458  }
5459
5460  if( typeof(attrib(id1,"onlySingularAtZero"))=="int" )
5461  {
5462    if(attrib(id1,"onlySingularAtZero")==1)
5463    {
5464      attrib(id,"isIsolatedSingularity",1);
5465    }
5466  }
5467  else
5468  {
5469    attrib(id,"onlySingularAtZero",0);
5470  }
5471
5472  if( typeof(attrib(id1,"isRad"))=="int" )
5473  {
5474    if(attrib(id1,"isRad")==1)
5475    {
5476      attrib(id,"isRad",1);
5477    }
5478  }
5479  else
5480  {
5481    attrib(id,"isRad",0);
5482  }
5483  return(id);
5484}
5485///////////////////////////////////////////////////////////////////////////////
5486
5487proc normalC(ideal id, list #)
5488"USAGE:  normalC(id [,choose]);  id = radical ideal, choose = optional list
5489         of string.
5490         Optional parameters in list choose (can be entered in any order):@*
5491         Decomposition:@*
5492         - \"equidim\" -> computes first an equidimensional decomposition,
5493         and then the normalization of each component (default).@*
5494         - \"prim\" -> computes first the minimal associated primes, and then
5495         the normalization of each prime. @*
5496         - \"noDeco\" -> no preliminary decomposition is done. If the ideal is
5497         not equidimensional radical, output might be wrong.@*
5498         - \"isPrim\" -> assumes that the ideal is prime. If the assumption does
5499         not hold, output might be wrong.@*
5500         - \"noFac\" -> factorization is avoided in the computation of the
5501         minimal associated primes;
5502         Other:@*
5503         - \"withGens\" -> the minimal associated primes P_i of id are
5504         computed and for each P_i, algebra generators of the integral closure
5505         of basering/P_i are computed as elements of its quotient field;@*
5506         If choose is not given or empty, the default options are used.@*
5507ASSUME:  The ideal must be radical, for non-radical ideals the output may
5508         be wrong (id=radical(id); makes id radical). However, if option
5509         \"prim\" is set the minimal associated primes are computed first
5510         and hence normalC computes the normalization of the radical of id.
5511         \"isPrim\" should only be used if id is known to be irreducible.
5512RETURN:  a list, say nor, of size 2 (resp. 3 if option \"withGens\" is set).@*
5513         * nor[1] is always a of r rings, where r is the number of associated
5514         primes with option \"prim\" (resp. >= no of equidimenensional
5515         components with option  \"equidim\").@*
5516         Each ring Ri=nor[1][i], i=1..r, contains two ideals with given
5517         names @code{norid} and @code{normap} such that @*
5518         - Ri/norid is the normalization of the i-th component, i.e. the
5519          integral closure in its field of fractions as affine ring, i.e. Ri is
5520          given in the form K[X(1..p),T(1..q)], where K is the ground field;
5521         - normap gives the normalization map from basering/id to
5522           Ri/norid for each i (the j-th element of normap is mapped to the
5523           j-th variable of R).@*
5524         - the direct sum of the rings Ri/norid is the normalization
5525           of basering/id; @*
5526         ** If option \"withGens\" is not set: @*
5527         * nor[2] shows the delta invariants: nor[2] is a list of an intvec of
5528         size r, the delta invariants of the r components, and an integer, the
5529         delta invariant of basering/id. (-1 means infinite, 0 that basering/P_i
5530         resp. basering/input is normal, -2 means that delta resp. delta of one
5531         of the components is not computed (which may happen if \"equidim\" is
5532         given). @*
5533         ** If option \"withGens\" is set:
5534         * nor[2] is a list of ideals Ii=nor[2][i], i=1..r, in the basering,
5535         generating the integral closure of basering/P_i in its quotient field
5536         as K-algebra (K the ground field):@*
5537         If Ii is given by polynomials g_1,...,g_k, then c:=g_k is a non-zero
5538         divisor and the j-th variables of the ring Ri satisfies var(j)=g_j/c,
5539         j=1..k-1, as element in the quotient field of basering/P_i. The
5540         g_j/g_k+1 are K-algebra generators  of the integral closure of
5541         basering/P_i.@*
5542         * nor[3] shows the delta invariant as above.
5543THEORY:  We use the Grauert-Remmert-de Jong algorithm [c.f. G.-M. Greuel,
5544         G. Pfister: A SINGULAR Introduction to Commutative Algebra, 2nd Edition.
5545         Springer Verlag (2007)].
5546         The procedure computes the algebra structure and the delta invariant of
5547         the normalization of R/id:@*
5548         The normalization is an affine algebra over the ground field K
5549         and nor[1] presents it as such: Ri = K[X(1..p),T(1..q)] and Ri/norid
5550         is the integral closure of R/P_i; if option \"withGens\" is set the
5551         X(j) and T(j) are expressed as quotients in the total ring of
5552         fractions. Note that the X(j) and T(j) generate the integral closure
5553         as K-algebra, but not necessarily as R-module (since relations of the
5554         form X(1)=T(1)*T(2) may have been eliminated). Geometrically the
5555         algebra structure is relevant since the variety of the ideal norid in
5556         Ri is the normalization of the variety of the ideal P_i in R.@*
5557         The delta invariant of a reduced ring A is dim_K(normalization(A)/A).
5558         For A=K[x1,...,xn]/id we call this number also the delta invariant of
5559         id. nor[3] returns the delta invariants of the components P_i and of
5560         id.
5561NOTE:    To use the i-th ring type: @code{def R=nor[1][i]; setring R;}.
5562@*       Increasing/decreasing printlevel displays more/less comments
5563         (default: printlevel=0).
5564@*       Not implemented for local or mixed orderings or quotient rings.
5565         For local or mixed orderings use proc 'normal'.
5566@*       If the input ideal id is weighted homogeneous a weighted ordering may
5567         be used (qhweight(id); computes weights).
5568KEYWORDS: normalization; integral closure; delta invariant.
5569SEE ALSO: normal, normalP.
5570EXAMPLE: example normalC; shows an example
5571"
5572{
5573   int i,j;
5574   int withGens, withEqui, withPrim, isPrim, noFac;
5575   int dbg = printlevel-voice+2;
5576   int nvar = nvars(basering);
5577   int chara  = char(basering);
5578   list result, prim, keepresult;
5579
5580  int decomp;   // Preliminar decomposition:
5581                // 0 -> no decomposition (id is assumed to be prime)
5582                // 1 -> no decomposition
5583                //      (id is assumed to be equidimensional radical)
5584                // 2 -> equidimensional decomposition
5585                // 3 -> minimal associated primes
5586
5587   // Default methods:
5588   noFac = 0;         // Use facstd when computing minimal associated primes
5589   decomp = 2;        // Equidimensional decomposition for nvar > 2
5590   if (nvar <= 2)
5591   { decomp = 3; }    // Compute minimal associated primes if nvar <= 2
5592
5593   if ( attrib(basering,"global") != 1 )
5594   {
5595     "";
5596     "// Not implemented for this ordering,";
5597     "// please change to global ordering or use proc normal";
5598     return(result);
5599   }
5600
5601//--------------------------- define the method ---------------------------
5602   string method;                //make all options one string in order to use
5603                                 //all combinations of options simultaneously
5604   for ( i=1; i <= size(#); i++ )
5605   {
5606     if ( typeof(#[i]) == "string" )
5607     {
5608       method = method + #[i];
5609     }
5610   }
5611
5612   //--------------------------- choosen methods -----------------------
5613   // "withGens": computes algebra generators for each irreducible component
5614   // ### the extra code for withGens should be incorporated in the general case
5615
5616   if ( find(method,"withgens") or find(method,"withGens"))
5617   {
5618     withGens = 1;
5619   }
5620
5621   // the general case: either equidim or minAssGTZ or no decomposition
5622
5623   if ( find(method,"isprim") or find(method,"isPrim") )
5624   {decomp = 0; isPrim=1;}
5625
5626   if ( find(method,"nodeco") or find(method,"noDeco") )
5627   {decomp = 1;}
5628
5629   if ( find(method,"equidim") )
5630   { decomp = 2; }
5631
5632   if ( find(method,"prim") )
5633   { decomp = 3; }
5634
5635   if ( find(method,"nofac") or find(method,"noFac") )
5636   { noFac = 1; }
5637
5638   kill #;
5639   list #;
5640
5641//------- Special algorithm with computation of the generators, RETURN -------
5642   //--------------------- method "withGens" ----------------------------------
5643   //the integral closure is computed in proc primeClosure. In the general case
5644   //it is computed in normalizationPrimes. The main difference is that in
5645   //primeClosure the singular locus is only computed in the first iteration,
5646   //that no attributes are used, and that the generators are computed.
5647   //In primeClosure the (algebra) generators for each irreducible component
5648   //are computed in the static proc closureGenerators
5649
5650   if( withGens )
5651   {
5652      if( dbg >= 1 )
5653      {  "";
5654         "// We use method 'withGens'";
5655      }
5656      if ( decomp == 0 or decomp == 1 )
5657      {
5658         prim[1] = id;
5659         if( dbg >= 0 )
5660         {
5661           "";
5662           "// ** WARNING: result is correct if ideal is prime (not checked) **";
5663           "// if procedure is called with string \"prim\", primality is checked";
5664         }
5665      }
5666      else
5667      {
5668         if(dbg >= 1)
5669         {  "// Computing minimal associated primes..."; }
5670
5671         if( noFac )
5672         { prim = minAssGTZ(id,1); }
5673         else
5674         { prim = minAssGTZ(id); }
5675
5676         if(dbg >= 2)
5677         {  prim;""; }
5678         if(dbg >= 1)
5679         {
5680            "// number of irreducible components is", size(prim);
5681         }
5682      }
5683   //----------- compute integral closure for every component -------------
5684      int del;
5685      intvec deli;
5686      list Gens,l,resu,Resu;
5687      ideal gens;
5688      def R = basering;
5689      poly gg;
5690
5691      for(i=1; i<=size(prim); i++)
5692      {
5693         if(dbg>=1)
5694         {
5695            ""; pause(); "";
5696            "// Computing normalization of component",i;
5697            "   ---------------------------------------";
5698         }
5699
5700         if( defined(ker) ) { kill ker; }
5701         ideal ker = prim[i];
5702         export(ker);
5703         l = R;
5704         l = primeClosure(l,1);              //here the work is done
5705         // primeClosure is called with list l consisting of the basering
5706         //### ausprobieren ob primeClosure(l,1) schneller als primeClosure(l)
5707         // 1 bedeutet: kuerzester nzd
5708         // l[size(l)] is the delta invariant
5709
5710         if ( l[size(l)] >= 0 && del >= 0 )
5711         {
5712            del = del + l[size(l)];
5713         }
5714         else
5715         { del = -1; }
5716         deli = l[size(l)],deli;
5717
5718         l = l[1..size(l)-1];
5719         resu = list(l[size(l)]) + resu;
5720         gens = closureGenerators(l);         //computes algebra(!) generators
5721
5722         //NOTE: gens[i]/gens[size(gens)] expresses the ith variable of resu[1]
5723         //(the normalization) as fraction of elements of the basering;
5724         //the variables of resu[1] are algebra generators.
5725         //gens[size(gens)] is a non-zero divisor of basering/i
5726
5727         //divide by the greatest common divisor:
5728         gg = gcd( gens[1],gens[size(gens)] );
5729         for(j=2; j<=size(gens)-1; j++)
5730         {
5731            gg=gcd(gg,gens[j]);
5732         }
5733         for(j=1; j<=size(gens); j++)
5734         {
5735            gens[j]=gens[j]/gg;
5736         }
5737         Gens = list(gens) + Gens;
5738
5739/*       ### Da die gens Algebra-Erzeuger sind, ist reduce nach Bestimmung
5740         der Algebra-Variablen T(i) nicht zulaessig!
5741         for(i=1;i<=size(gens)-1;i++)
5742         {
5743            gens[i]= reduce(gens[i],std(gens[size(gens)]));
5744         }
5745         for(i=size(gens)-1; i>=1; i--)
5746         {
5747            if(gens[i]==0)
5748            { gens = delete(gens,i); }
5749         }
5750*/
5751         if( defined(ker) ) { kill ker; }
5752      }
5753
5754      if ( del >= 0 )
5755      {
5756         int mul = iMult(prim);
5757         del = del + mul;
5758      }
5759      else
5760      { del = -1; }
5761      deli = deli[1..size(deli)-1];
5762      Resu = resu,Gens,list(deli,del);
5763      int sr = size(resu);
5764
5765      if ( dbg >= 0 )
5766      {"";
5767"// 'normalC' created a list, say nor, of three lists:
5768// To see the list type
5769      nor;
5770
5771// * nor[1] is a list of",sr,"ring(s)
5772// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
5773     def R1 = nor[1][1]; setring R1;  norid; normap;
5774// For the other rings type first (if R is the name of your original basering)
5775     setring R;
5776// and then continue as for R1.
5777// Ri/norid is the affine algebra of the normalization of the i-th
5778// component R/P_i (where P_i is an associated prime of the input ideal id)
5779// and normap the normalization map from R to Ri/norid.
5780
5781// * nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of
5782// elements g1..gk of R such that the gj/gk generate the integral
5783// closure of R/P_i as sub-algebra in the quotient field of R/P_i, with
5784// gj/gk being mapped by normap to the j-th variable of Ri;
5785
5786// * nor[3] shows the delta-invariant of each component and of id
5787// (-1 means infinite, and 0 that R/P_i resp. R/id is normal).";
5788      }
5789      return(Resu);
5790   }
5791   //----------------- end method "withGens" --------------------------------
5792
5793//-------- The general case without computation of the generators -----------
5794// (attrib(id,"xy"),1) sets attrib xy to TRUE and (attrib(id,"xy"),0) to FALSE
5795// We use the following attributes:
5796//   attrib(id,"isCohenMacaulay");         //--- Cohen Macaulay
5797//   attrib(id,"isCompleteIntersection");  //--- complete intersection
5798//   attrib(id,"isHypersurface");          //--- hypersurface
5799//   attrib(id,"isEquidimensional",-1);    //--- equidimensional ideal
5800//   attrib(id,"isPrim");                  //--- prime ideal
5801//   attrib(id,"isRegInCodim2");           //--- regular in codimension 2
5802//   attrib(id,"isIsolatedSingularity";    //--- isolated singularities
5803//   attrib(id,"onlySingularAtZero");      //--- only singular at 0
5804
5805 //------------------- first set the attributes ----------------------
5806   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
5807   {
5808      if( attrib(id,"isCohenMacaulay")==1 )
5809      {
5810         attrib(id,"isEquidimensional",1);
5811      }
5812   }
5813   else
5814   {
5815      attrib(id,"isCohenMacaulay",0);
5816   }
5817
5818   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
5819   {
5820      if(attrib(id,"isCompleteIntersection")==1)
5821      {
5822         attrib(id,"isCohenMacaulay",1);
5823         attrib(id,"isEquidimensional",1);
5824      }
5825   }
5826   else
5827   {
5828      attrib(id,"isCompleteIntersection",0);
5829   }
5830
5831   if( typeof(attrib(id,"isHypersurface"))=="int" )
5832   {
5833      if(attrib(id,"isHypersurface")==1)
5834      {
5835         attrib(id,"isCompleteIntersection",1);
5836         attrib(id,"isCohenMacaulay",1);
5837         attrib(id,"isEquidimensional",1);
5838      }
5839   }
5840   else
5841   {
5842      attrib(id,"isHypersurface",0);
5843   }
5844
5845   if( ! (typeof(attrib(id,"isEquidimensional"))=="int") )
5846   {
5847         attrib(id,"isEquidimensional",0);
5848   }
5849
5850   if( typeof(attrib(id,"isPrim"))=="int" )
5851   {
5852      if(attrib(id,"isPrim")==1)
5853      {
5854         attrib(id,"isEquidimensional",1);
5855      }
5856   }
5857   else
5858   {
5859      attrib(id,"isPrim",0);
5860   }
5861
5862   if( ! (typeof(attrib(id,"isRegInCodim2"))=="int") )
5863   {
5864         attrib(id,"isRegInCodim2",0);
5865   }
5866
5867   if( ! (typeof(attrib(id,"isIsolatedSingularity"))=="int") )
5868   {
5869         attrib(id,"isIsolatedSingularity",0);
5870   }
5871
5872   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
5873   {
5874      if(attrib(id,"onlySingularAtZero")==1)
5875      {
5876         attrib(id,"isIsolatedSingularity",1);
5877      }
5878   }
5879   else
5880   {
5881      attrib(id,"onlySingularAtZero",0);
5882   }
5883
5884   //-------------- compute equidimensional decomposition --------------------
5885   //If the method "equidim" is given, compute the equidim decomposition
5886   //and goto the next step (no normalization
5887   //ACHTUNG: equidim berechnet bei nicht reduzierten id die eingebetteten
5888   //Komponenten als niederdim Komponenten, waehrend diese bei primdecGTZ
5889   //nicht auftauchen: ideal(x,y)*xy
5890   //this is default for nvars > 2
5891
5892   if( decomp == 2 )
5893   {
5894      withPrim = 0;                 //this is used to check later that prim
5895                                    //contains equidim but not prime components
5896      if( dbg >= 1 )
5897      {
5898         "// We use method 'equidim'";
5899      }
5900      if( typeof(attrib(id,"isEquidimensional"))=="int" )
5901      {
5902         if(attrib(id,"isEquidimensional")==1)
5903         {
5904            prim[1] = id;
5905         }
5906         else
5907         {
5908            prim = equidim(id);
5909         }
5910      }
5911      else
5912      {
5913         prim = equidim(id);
5914      }
5915      if(dbg>=1)
5916      {  "";
5917         "// number of equidimensional components:", size(prim);
5918      }
5919      if ( !noFac )
5920      {
5921        intvec opt = option(get);
5922        option(redSB);
5923        for(j=1; j<=size(prim); j++)
5924        {
5925           keepresult = keepresult+facstd(prim[j]);
5926        }
5927        prim = keepresult;
5928        if ( size(prim) == 0 )
5929        {
5930          prim=ideal(0);     //Bug in facstd, liefert leere Liste bei 0-Ideal
5931        }
5932
5933        if(dbg>=1)
5934        {  "";
5935         "// number of components after application of facstd:", size(prim);
5936        }
5937        option(set,opt);
5938      }
5939   }
5940
5941   //------------------- compute associated primes -------------------------
5942   //the case where withEqui = 0, here the min. ass. primes are computed
5943   //start with the computation of the minimal associated primes:
5944
5945   else
5946   {
5947    if( isPrim )
5948    {
5949      if( dbg >= 0 )
5950      {
5951         "// ** WARNING: result is correct if ideal is prime";
5952         "// or equidimensional (not checked) **";
5953         "// disable option \"isPrim\" to decompose ideal into prime";
5954         "// or equidimensional components";"";
5955      }
5956      if( dbg >= 1 )
5957      {
5958        "// We use method 'isPrim'";"";
5959      }
5960      prim[1]=id;
5961    }
5962    else
5963    {
5964      withPrim = 1;                 //this is used to check later that prim
5965                                    //contains prime but not equidim components
5966      if( dbg >= 1 )
5967      {
5968         "// We use method 'prim'";
5969      }
5970
5971      if( typeof(attrib(id,"isPrim"))=="int" )
5972      {
5973         if(attrib(id,"isPrim")==1)
5974         {
5975            prim[1]=id;
5976         }
5977         else
5978         {
5979            if( noFac )
5980            { prim=minAssGTZ(id,1); }     //does not use factorizing groebner
5981            else
5982            { prim=minAssGTZ(id); }       //uses factorizing groebner
5983         }
5984      }
5985      else
5986      {
5987            if( noFac )
5988            { prim=minAssGTZ(id,1); }
5989            else
5990            { prim=minAssGTZ(id); }
5991      }
5992      if(dbg>=1)
5993      {  "";
5994         "// number of irreducible components:", size(prim);
5995      }
5996    }
5997   }
5998
5999   //----- for each component (equidim or irred) compute normalization -----
6000   int sr, skr, del;
6001   intvec deli;
6002   int sp = size(prim);     //size of list prim (# irred or equidim comp)
6003
6004   for(i=1; i<=sp; i++)
6005   {
6006      if(dbg>=1)
6007      {  "";
6008         "// computing the normalization of component",i;
6009         "   ----------------------------------------";
6010      }
6011      //-------------- first set attributes for components ------------------
6012      attrib(prim[i],"isEquidimensional",1);
6013      if( withPrim )
6014      {
6015         attrib(prim[i],"isPrim",1);
6016      }
6017      else
6018      { attrib(prim[i],"isPrim",0); }
6019
6020      if(attrib(id,"onlySingularAtZero")==1)
6021      { attrib(prim[i],"onlySingularAtZero",1); }
6022      else
6023      { attrib(prim[i],"onlySingularAtZero",0); }
6024
6025      if(attrib(id,"isIsolatedSingularity")==1)
6026      { attrib(prim[i],"isIsolatedSingularity",1); }
6027      else
6028      { attrib(prim[i],"isIsolatedSingularity",0); }
6029
6030      if( attrib(id,"isHypersurface")==1 )
6031      {
6032         attrib(prim[i],"isHypersurface",1);
6033         attrib(prim[i],"isCompleteIntersection",1);
6034         attrib(prim[i],"isCohenMacaulay",1);
6035      }
6036      else
6037      { attrib(prim[i],"isHypersurface",0); }
6038
6039      if ( sp == 1)         //the case of one component: copy attribs from id
6040      {
6041        if(attrib(id,"isRegInCodim2")==1)
6042        {attrib(prim[i],"isRegInCodim2",1); }
6043        else
6044        {attrib(prim[i],"isRegInCodim2",0); }
6045
6046        if(attrib(id,"isCohenMacaulay")==1)
6047        {attrib(prim[i],"isCohenMacaulay",1); }
6048        else
6049        {attrib(prim[i],"isCohenMacaulay",0); }
6050
6051        if(attrib(id,"isCompleteIntersection")==1)
6052        {attrib(prim[i],"isCompleteIntersection",1); }
6053        else
6054        {attrib(prim[i],"isCompleteIntersection",0); }
6055      }
6056      else
6057      {
6058        attrib(prim[i],"isRegInCodim2",0);
6059        attrib(prim[i],"isCohenMacaulay",0);
6060        attrib(prim[i],"isCompleteIntersection",0);
6061      }
6062
6063      //------ Now compute the normalization of each component ---------
6064      //note: for equidimensional components the "splitting tools" can
6065      //create further decomposition
6066      //We now start normalizationPrimes with
6067      //ihp = partial normalisation map = identity map = maxideal(1)
6068      //del = partial delta invariant = 0
6069      //deli= intvec of partial delta invariants of components
6070      //in normalizationPrimes all the work is done:
6071
6072      keepresult = normalizationPrimes(prim[i],maxideal(1),0,0);
6073
6074      for(j=1; j<=size(keepresult)-1; j++)
6075      {
6076         result=insert(result,keepresult[j]);
6077      }
6078      skr = size(keepresult);
6079
6080      //compute delta:
6081      if( del >= 0 && keepresult[skr][1] >=0 )
6082      {
6083         del = del + keepresult[skr][1];
6084      }
6085      else
6086      {
6087         del = -1;
6088      }
6089      deli = keepresult[skr][2],deli;
6090
6091      if ( dbg>=1 )
6092      {
6093           "// delta of component",i; keepresult[skr][1];
6094      }
6095   }
6096   sr = size(result);
6097
6098   // -------------- Now compute intersection multiplicities -------------
6099   //intersection multiplicities of list prim, sp=size(prim).
6100      if ( dbg>=1 )
6101      {
6102        "// Sum of delta for all components"; del;
6103        if ( sp>1 )
6104        {
6105           "// Compute intersection multiplicities of the components";
6106        }
6107      }
6108
6109      if ( sp > 1 )
6110      {
6111        int mul = iMult(prim);
6112        if ( mul < 0 )
6113        {
6114           del = -1;
6115        }
6116        else
6117        {
6118           del = del + mul;
6119        }
6120      }
6121   deli = deli[1..size(deli)-1];
6122   result = result,list(deli,del);
6123
6124//--------------- Finally print comments and return ------------------
6125   if ( dbg >= 0)
6126   {"";
6127"// 'normalC' created a list, say nor, of two lists:
6128// To see the result, type
6129      nor;
6130
6131// * nor[1] is a list of",sr,"ring(s).
6132// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
6133      def R1 = nor[1][1];  setring R1;  norid;  normap;
6134// and similair for the other rings nor[1][i];
6135// Ri/norid is the affine algebra of the normalization of r/P_i  (where P_i
6136// is an associated prime or an equidimensional part of the input ideal id)
6137// and normap the normalization map from the basering to Ri/norid;
6138
6139// * nor[2] shows the delta-invariant of each component and of id
6140// (-1 means infinite, 0 that r/P_i resp. r/id is normal, and -2 that delta
6141// of a component was not computed).";
6142   }
6143   return(result);
6144}
6145
6146example
6147{ "EXAMPLE:";
6148   printlevel = printlevel+1;
6149   echo = 2;
6150   ring s = 0,(x,y),dp;
6151   ideal i = (x2-y3)*(x2+y2)*x;
6152
6153   list nor = normalC(i);
6154
6155   nor;
6156   // 2 branches have delta = 1, and 1 branch has delta = 0
6157   // the total delta invariant is 13
6158
6159   def R2 = nor[1][2];  setring R2;
6160   norid; normap;
6161
6162   echo = 0;
6163   printlevel = printlevel-1;
6164   pause("   hit return to continue"); echo=2;
6165
6166   ring r = 2,(x,y,z),dp;
6167   ideal i = z3-xy4;
6168   nor = normalC(i);  nor;
6169   // the delta invariant is infinite
6170   // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module
6171   // in its quotient field Quot(r/i)
6172
6173   // the normalization as affine algebra over the ground field:
6174   def R = nor[1][1]; setring R;
6175   norid; normap;
6176
6177   echo = 0;
6178   pause("   hit return to continue");echo = 2;
6179
6180   setring r;
6181   nor = normalC(i, "withGens", "prim");    // a different algorithm
6182   nor;
6183}
6184
6185//////////////////////////////////////////////////////////////////////////////
6186//closureRingtower seems not to be used anywhere
6187static proc closureRingtower(list L)
6188"USAGE:    closureRingtower(list L); L a list of rings
6189CREATE:   rings R(1),...,R(n) such that R(i)=L[i] for all i
6190EXAMPLE:  example closureRingtower; shows an example
6191"
6192{
6193  int n=size(L);
6194  for (int i=1;i<=n;i++)
6195    {
6196      if (defined(R(i)))
6197      {
6198        string s="Fixed name R("+string(i)+") leads to conflict with existing "
6199              +"object having this name";
6200        ERROR(s);
6201      }
6202      def R(i)=L[i];
6203      export R(i);
6204    }
6205
6206  return();
6207}
6208example
6209{
6210  "EXAMPLE:"; echo=2;
6211  ring R=0,(x,y),dp;
6212  ideal I=x4,y4;
6213  list L=primeClosure(ReesAlgebra(I)[1]);
6214  L=delete(L,size(L));
6215  L;
6216  closureRingtower(L);
6217  R(1);
6218  R(4);
6219  kill R(1),R(2),R(3),R(4);
6220}
6221
6222//                Up to here: procedures for normalC
6223///////////////////////////////////////////////////////////////////////////////
6224
6225///////////////////////////////////////////////////////////////////////////////
6226//                From here: miscellaneous procedures
6227
6228// Used for timing and comparing the different normalization procedures.
6229// Option (can be entered in any order)
6230// "normal"   -> uses the new algortihm (normal)
6231// "normalP"  -> uses normalP
6232// "normalC"  -> uses normalC, without "withGens" option
6233// "primCl"   -> uses normalC, with option "withGens".
6234// "111"      -> checks the output of normalM using norTest.
6235// "p"        -> compares the output of norM with the output of normalP
6236//               ("normalP" option must also be set).
6237// "pc"       -> compares the output of norM with the output of normalC with
6238//               option "withGens"
6239//               ("primCl" option must also be set).
6240
6241proc timeNormal(ideal I, list #)
6242{
6243  def r = basering;
6244
6245  //--------------------------- define the method ---------------------------
6246  int isPrim, useRing;
6247  int decomp = -1;
6248  int norM, norC, norP, primCl;
6249  int checkP, check111, checkPC;
6250  int i;
6251  ideal U1, U2, W;
6252  poly c1, c2;
6253  int ch;
6254  string check;
6255  string method;                //make all options one string in order to use
6256                                //all combinations of options simultaneously
6257  for ( i=1; i <= size(#); i++ )
6258  {
6259    if ( typeof(#[i]) == "string" )
6260    {
6261      method = method + #[i];
6262    }
6263  }
6264  if ( find(method, "normal"))
6265  {norM = 1;}
6266  if ( find(method, "normalP") and (char(basering) > 0))
6267  {norP = 1;}
6268  if ( find(method, "normalC"))
6269  {norC = 1;}
6270  if ( find(method, "primCl"))
6271  {primCl = 1;}
6272  if ( find(method, "isprim") or find(method,"isPrim") )
6273  {decomp = 0;}
6274  if ( find(method, "p") )
6275  {checkP = 1;}
6276  if ( find(method, "pc") )
6277  {checkPC = 1;}
6278  if ( find(method, "111") )
6279  {check111 = 1;}
6280
6281  int tt;
6282  if(norM)
6283  {
6284    tt = timer;
6285    if(decomp == 0)
6286    {
6287      "Running normal(useRing, isPrim)...";
6288      list a1 = normal(I, "useRing", "isPrim");
6289      "Time normal(useRing, isPrim): ", timer - tt;
6290    }
6291    else
6292    {
6293      "Running normal(useRing)...";
6294      list a1 = normal(I, "useRing");
6295      "Time normal(useRing): ", timer - tt;
6296    }
6297    "";
6298  }
6299  if(norP)
6300  {
6301    tt = timer;
6302    if(decomp == 0)
6303    {
6304      "Running normalP(isPrim)...";
6305      list a2 = normalP(I, "isPrim");
6306      "Time normalP(isPrim): ", timer - tt;
6307    }
6308    else
6309    {
6310      "Running normalP()...";
6311      list a2 = normalP(I);
6312      "Time normalP(): ", timer - tt;
6313    }
6314    "";
6315  }
6316
6317  if(norC)
6318  {
6319    tt = timer;
6320    if(decomp == 0)
6321    {
6322      "Running normalC(isPrim)...";
6323      list a3 = normalC(I, "isPrim");
6324      "Time normalC(isPrim): ", timer - tt;
6325    }
6326    else
6327    {
6328      "Running normalC()...";
6329      list a3 = normalC(I);
6330      "Time normalC(): ", timer - tt;
6331    }
6332    "";
6333  }
6334
6335  if(primCl)
6336  {
6337    tt = timer;
6338    if(decomp == 0)
6339    {
6340      "Running normalC(withGens, isPrim)...";
6341      list a4 = normalC(I, "isPrim", "withGens");
6342      "Time normalC(withGens, isPrim): ", timer - tt;
6343    }
6344    else
6345    {
6346      "Running normalC(withGens)...";
6347      list a4 = normalC(I, "withGens");
6348      "Time normalC(withGens): ", timer - tt;
6349    }
6350    "";
6351  }
6352
6353  if(check111 and norM)
6354  {
6355    "Checking output with norTest...";
6356    "WARNING: this checking only works if the original ideal was prime.";
6357    norTest(I, a1);
6358    "";
6359  }
6360
6361  if(checkP and norP and norM)
6362  {
6363    "Comparing with normalP output...";
6364    if(size(a2) > 0)
6365    {
6366      "WARNING: this checking only works if the original ideal was prime.";
6367      U1 = a1[2][1];
6368      c1 = U1[size(U1)];
6369      U2 = a2[1][1];
6370      c2 = a2[1][1][size(a2[1][1])];
6371      W = changeDenominator(U1, c1, c2, groebner(I));
6372      qring q = groebner(I);
6373      ideal U2 = fetch(r, U2);
6374      ideal W = fetch(r, W);
6375      ch = 0;
6376      if(size(reduce(U2, groebner(W))) == 0)
6377      {
6378        "U2 c U1";
6379        ch = 1;
6380      }
6381      if(size(reduce(W, groebner(U2))) == 0)
6382      {
6383        "U1 c U2";
6384        ch = ch + 1;
6385      }
6386      if(ch == 2)
6387      {
6388        "Output of normalP is equal.";
6389      }
6390      else
6391      {
6392        "ERROR: Output of normalP is different.";
6393      }
6394      setring r;
6395      kill q;
6396    }
6397    else
6398    {
6399      "normalP returned no output. Comparison is not possible.";
6400    }
6401    "";
6402  }
6403
6404  if(checkPC and norM and primCl)
6405  {
6406    "Comparing with primeClosure output...";
6407    if(size(a4) > 0)
6408    {
6409      "WARNING: this checking only works if the original ideal was prime.";
6410      // primeClosure check
6411      U1 = a1[2][1];
6412      c1 = U1[size(U1)];
6413      U2 = a4[2][1];
6414      c2 = a4[2][1][size(a4[2][1])];
6415      W = changeDenominator(U1, c1, c2, groebner(I));
6416      qring q = groebner(I);
6417      ideal U2 = fetch(r, U2);
6418      ideal W = fetch(r, W);
6419      ch = 0;
6420      if(size(reduce(U2, groebner(W))) == 0)
6421      {
6422        "U2 c U1";
6423        ch = 1;
6424      }
6425      if(size(reduce(W, groebner(U2))) == 0)
6426      {
6427        "U1 c U2";
6428        ch = ch + 1;
6429      }
6430      if(ch == 2)
6431      {
6432        "Output of normalC(withGens) is equal.";
6433      }
6434      else
6435      {
6436        "ERROR: Output of normalC(withGens) is different.";
6437      }
6438      setring r;
6439      kill q;
6440    }
6441    else
6442    {
6443      "normalC(withGens) returned no output. Comparison is not possible.";
6444    }
6445    "";
6446  }
6447}
6448
6449///////////////////////////////////////////////////////////////////////////
6450static proc sqroot(int n);
6451{
6452  int s = 1;
6453  while(s*s < n) { s++; }
6454  return(s);
6455}
6456
6457///////////////////////////////////////////////////////////////////////////
6458proc norTest (ideal i, list nor, list #)
6459"USAGE:   norTest(i,nor,[n]); i=prime ideal, nor=list, n=optional integer
6460ASSUME:  nor is the output of normal(i) (any options) or
6461         normalP(i,"withRing") or normalC(i) (any options).
6462         In particular, the ring nor[1][1] contains the ideal norid
6463         and the map normap: basering/i --> nor[1][1]/norid.
6464RETURN:  an intvec v such that:
6465@format
6466         v[1] = 1 if the normap is injective and 0 otherwise
6467         v[2] = 1 if the normap is finite and 0 otherwise
6468         v[3] = 1 if nor[1][1]/norid is normal and 0 otherwise
6469@end format
6470         If n=1 (resp n=2) only v[1] (resp. v[2]) is computed and returned
6471THEORY:  The procedure can be used to test whether the computation of the
6472         normalization was correct: basering/i --> nor[1][1]/norid is the
6473         normalization of basering/i if and only if v=1,1,0.
6474NOTE:    For big examples it can be hard to fully test correctness; the
6475         partial test norTest(i,nor,2) is usually fast
6476EXAMPLE: example norTest; shows an example
6477"
6478{
6479//### Sollte erweitert werden auf den reduziblen Fall: einen neuen affinen
6480// Ring nor[1][1]+...+nor[1][r] (direkte Summe) erzeugen, map dorthin
6481// definieren und dann testen.
6482
6483    int prl = printlevel - voice + 2;
6484    int a,b,d;
6485    int n,ii;
6486    if (size(#) > 0) {  n = #[1];  }
6487
6488    def BAS = basering;
6489
6490    //### make a copy of nor to have a cpoy of nor[1][1]  (not a reference to)
6491    // in order not to override norid and normap.
6492    // delete nor[2] (if it contains the module generators, which are not used)
6493    // s.t. newnor does not belong to a ring.
6494
6495    list newnor = nor;
6496    if ( size(newnor) == 3 )
6497    {
6498       newnor = delete(newnor,2);
6499    }
6500    def R = newnor[1][1];
6501    qring QAS = std(i);
6502
6503
6504    setring R;
6505    int nva = nvars(R);
6506    string svars = varstr(R);
6507    string svar;
6508
6509    norid = interred(norid);
6510
6511    //--------- create new ring with one dp block keeping weights ------------
6512    list LR = ringlist(R);
6513    list g3 = LR[3];
6514    int n3 = size(g3);
6515    list newg3;
6516    intvec V;
6517
6518    //--------- check first whether variables Z(i),...,A(i) exist -----------
6519    for (ii=90; ii>=65; ii--)
6520    {
6521       if ( find(svars,ASCII(ii)+"(") == 0 )
6522       {
6523          svar = ASCII(ii);  break;
6524       }
6525    }
6526    if ( size(svar) != 0 )
6527    {
6528        for ( ii = 1; ii <= nva; ii++ )
6529        {
6530            LR[2][ii] = svar+"("+string(ii)+")";
6531            V[ii] = 1;
6532        }
6533    }
6534    else
6535    {
6536        for ( ii = 1; ii <= nva; ii++ )
6537        {
6538           LR[2][ii] = "Z("+string(100*nva+ii)+")";
6539           V[ii] = 1;
6540        }
6541    }
6542
6543    if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
6544    {
6545       list gm = g3[n3];       //last blockis module ordering
6546       newg3[1] = list("dp",V);
6547       newg3 = insert(newg3,gm,size(newg3));
6548    }
6549    else
6550    {
6551       list gm = g3[1];              //first block is module ordering
6552       newg3[1] = list("dp",V);
6553       newg3 = insert(newg3,gm);
6554    }
6555    LR[3] = newg3;
6556//LR;"";
6557    def newR = ring(LR);
6558
6559    setring newR;
6560    ideal norid = fetch(R,norid);
6561    ideal normap = fetch(R,normap);
6562    if( defined(lnorid) )  { kill lnorid; }     //um ** redefinig zu beheben
6563    if( defined(snorid) )  { kill snorid; }     //sollte nicht noetig sein
6564
6565    //----------- go to quotient ring for checking injectivity -------------
6566//"mstd";
6567    list lnorid = mstd(norid);
6568    ideal snorid = lnorid[1];
6569//"size mstdnorid:", size(snorid),size(lnorid[2]);
6570//"size string mstdnorid:", size(string(snorid)),size(string(lnorid[2]));
6571    qring QR = snorid;
6572    ideal qnormap = fetch(newR,normap);
6573    //ideal qnormap = imap(newR,normap);
6574    //ideal qnormap = imap(R,normap);
6575    map Qnormap = QAS,qnormap;    //r/id --> R/norid
6576
6577    //------------------------ check injectivity ---------------------------
6578//"injective:";
6579    a = is_injective(Qnormap,QAS);          //a. Test for injectivity of Qnormap
6580    dbprint ( prl, "injective: "+string(a) );
6581    if ( n==1 ) { return (a); }
6582   a;
6583
6584    //------------------------ check finiteness ---------------------------
6585    setring newR;
6586    b = mapIsFinite(normap,BAS,lnorid[2]);  //b. Test for finiteness of normap
6587    dbprint ( prl, "finite: "+string(b) );
6588    if ( n==2 ) { return (intvec(a,b)); }
6589   b;
6590
6591    //------------------------ check normality ---------------------------
6592    list testnor = normal(lnorid[2],"isPrim","noFac", "withDelta");
6593    //### Problem: bei mehrfachem Aufruf von norTest gibt es
6594    // ** redefining norid & ** redefining normap
6595    //Dies produziert Fehler, da alte norid und normap ueberschrieben werden
6596    //norid und normap werden innnerhalb von proc computeRing ueberschrieben
6597    //Die Kopie newR scheint das Problem zu loesen
6598
6599    printlevel=prl;
6600
6601    d = testnor[3][2];             //d = delta
6602    kill testnor;                              //### sollte ueberfluessig sein
6603    int d1 = (d==0);                           //d1=1 if delta=0
6604    dbprint ( prl, "delta: "+string(d) );
6605    return(intvec(a,b,d1));
6606}
6607example
6608{ "EXAMPLE:"; echo = 2;
6609   int prl = printlevel;
6610   printlevel = -1;
6611   ring r = 0,(x,y),dp;
6612   ideal i = (x-y^2)^2 - y*x^3;
6613   list nor = normal(i);
6614   norTest(i,nor);                //1,1,1 means that normal was correct
6615
6616   nor = normalC(i);
6617   norTest(i,nor);                //1,1,1 means that normal was correct
6618
6619   ring s = 2,(x,y),dp;
6620   ideal i = (x-y^2)^2 - y*x^3;
6621   nor = normalP(i,"withRing");
6622   norTest(i,nor);               //1,1,1 means that normalP was correct
6623   printlevel = prl;
6624}
6625
6626///////////////////////////////////////////////////////////////////////////
6627//
6628//                            EXAMPLES
6629//
6630///////////////////////////////////////////////////////////////////////////
6631/*
6632//commands for computing the normalization:
6633// options for normal:  "equidim", "prim"
6634//                      "noDeco", "isPrim", "noFac"
6635//                       (prim by default)
6636// options for normalP: "withRing", "isPrim" or "noFac"
6637// options for normalC: "equidim", "prim", "withGens",
6638//                      "noDeco", "isPrim", "noFac"
6639
6640//Commands for testing 'normal'
6641 list nor = normal(i); nor;
6642 list nor = normal(i,"isPrim");nor;
6643 list nor = normal(i,"equidim");nor;
6644 list nor = normal(i,"prim");nor;
6645 list nor = normal(i,"equidim","noFac");nor;
6646 list nor = normal(i,"prim","noFac");nor;
6647
6648//Commands for testing 'normalP' in positive char
6649 list nor = normalP(i);nor;              //withGens but no ringstructure
6650 list nor = normalP(i,"withRing"); nor;  //compute the ringstructure
6651 list nor = normalP(i,"isPrim"); nor;    //if i is known to be prime
6652
6653//Commands for testing 'normalC'
6654 list nor = normal(i); nor;
6655 list nor = normal(i,"withGens");nor;
6656 list nor = normal(i,"isPrim");nor;
6657 list nor = normal(i,"equidim");nor;
6658 list nor = normal(i,"prim");nor;
6659 list nor = normal(i,"equidim","noFac");nor;
6660 list nor = normal(i,"prim","noFac");nor;
6661
6662//Commands for testing correctness (i must be prime):
6663list nor = normalP(i,"withRing","isPrim");
6664list nor = normal(i,"isPrim");
6665norTest(i,nor);       //full test for not too big examples (1,1,1 => ok)
6666norTest(i,nor,2);     //partial test for big examples (1,1 => ok)
6667factorize(i[1]);      //checks for irreducibility
6668
6669/////////////////////////////////////////////////////////////////////////////
6670
6671//----------------------Examples for normal (new algorithm)------------------
6672// Timings with Computeserver Dual AMD Opteron 242 1.60GHz.
6673// Examples from "Normalization of Rings" paper.
6674
6675// Example 1
6676// char 0 : normal = 0 secs (7 steps) - normalC = 75 secs
6677// char 2 : normal = 0 secs (7 steps) - normalP = 0 secs - normalC = 0 secs
6678// char 5 : normal = 1 secs (7 steps) - normalP = 71 - normalC = 1 secs
6679// char 11 : normal = 2 secs (7 steps) - normalP = 12 secs - normalC doesn't finish
6680// char 32003 : normal = 1 secs (7 steps) - normalP doesn't finish - normalC = 1 sec
6681LIB"normal.lib";
6682ring r = 2, (x, y), dp;
6683ideal i = (x-y)*x*(y+x^2)^3-y^3*(x^3+x*y-y^2);
6684timeNormal(i, "normal", "normalC", "normalP", "isPrim", "p");
6685
6686// Example 2
6687// char 0  : normal = 1 sec (7 steps) - normalC doesn't finish
6688// char 3 : normal = 1 secs (8 steps) - normalP = 0 secs - normalC = 4 secs
6689// char 13 : normal = 1 sec (7 steps) - normalP doesn't finish - normalC = 13 secs
6690// char 32003 : normal = 1 secs (7 steps) - normalP doesn't finish - normalC = 10 sec
6691//Example is reducible in char 5 and 7
6692LIB"normal.lib";
6693ring r = 3, (x, y), dp;
6694ideal i = 55*x^8+66*y^2*x^9+837*x^2*y^6-75*y^4*x^2-70*y^6-97*y^7*x^2;
6695timeNormal(i, "normal", "normalC", "normalP", "p", "isPrim");
6696
6697// Example 3
6698// char 0 : normal = 3 secs (6 steps) - normalC doesn't finish
6699// char 2 : normal = 1 secs (13 steps) - normalP = 0 secs - normalC doesn't finish
6700// char 5 : normal = 0 secs (6 steps) - normalP = 8 secs - normalC doesn't finish
6701LIB"normal.lib";
6702ring r=5,(x, y),dp;
6703ideal i=y9+y8x+y8+y5+y4x+y3x2+y2x3+yx8+x9;
6704timeNormal(i, "normal", "normalC", "normalP", "isPrim");
6705
6706// Example 4
6707// char 0 : normal = 0 secs (1 step) - normalC = 0 secs
6708// char 5 : normal = 0 secs (1 step) - normalP = 3 secs - normalC = 0 secs
6709// char 11 : normal = 0 secs (1 step) - normalP doesn't finish - normalC = 0 secs
6710// char 32003 : normal = 0 secs (1 step) - normalP doesn't finish - normalC = 0 secs
6711LIB"normal.lib";
6712ring r=5,(x,y),dp;   // genus 0 4 nodes and 6 cusps im P2
6713ideal i=(x2+y^2-1)^3 +27x2y2;
6714timeNormal(i, "normal", "normalC", "normalP", "isPrim");
6715
6716// Example 5
6717// char 0 : normal = 0 secs (1 step) - normalC = 0 secs
6718// char 5 : normal = 1 secs (3 step) - normalP doesn't finish - normalC doesn't finish
6719// char 11 : normal = 0 secs (1 step) - normalP 0 secs - normalC = 0 secs
6720// char 32003 : normal = 0 secs (1 step) - normalP doesn't finish - normalC = 0 secs
6721LIB"normal.lib";
6722ring r=11,(x,y),dp;    //24 sing, delta 24
6723ideal i=-x10+x8y2-x6y4-x2y8+2y10-x8+2x6y2+x4y4-x2y6-y8+2x6-x4y2+x2y4+2x4+2x2y2-y4-x2+y2-1;
6724timeNormal(i, "normal", "normalC", "normalP", "isPrim", "p");
6725
6726// Example 6
6727// char 2 : normal = 5 secs (2 steps) - normalP = 25 secs - normalC = 166 secs
6728LIB"normal.lib";
6729ring r=2,(v,u,z,y,x),dp;
6730ideal i = z3+zyx+y3x2+y2x3, uyx+z2,uz+z+y2x+yx2, u2+u+zy+zx, v3+vux+vz2+vzyx+vzx+uz3+uz2y+z3+z2yx2;
6731timeNormal(i, "normal", "normalC", "normalP", "isPrim", "p");
6732
6733// Example 7
6734// char 0 : normal = 11 secs (6 steps) - normalC = 11 secs
6735// char 2 : normal = 11 secs (6 steps) - normalP = 0 secs - normalC = 11 secs
6736// char 5 : normal = 11 secs (6 steps) - normalP = 3 secs - normalC = 11 secs
6737// char 11 : normal = 11 secs (6 steps) - normalP = 43 secs - normalC = 11 secs
6738// char 32003 : normal = 11 secs (6 steps) - normalP doesn't finish - normalC = 11 secs
6739LIB"normal.lib";
6740ring r=11,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
6741ideal i= x2+zw, y3+xwt, xw3+z3t+ywt2, y2w4-xy2z2t-w3t3;
6742timeNormal(i, "normal", "normalC", "normalP", "isPrim");
6743
6744////////////////////////////////////////////////////////////////////////////////
6745
6746// Other examples with new algorithm
6747
6748// Example 1
6749// char 0 : normal = 1 secs (13 steps) - normalC doesn't finish
6750// char 2 : normal = 1 secs (13 steps) - normalP = 0 secs - normalC doesn't finish
6751// char 5 : normal = 1 secs (13 steps) - normalP = 29 secs - normalC doesn't finish
6752ring r=2,(x,y),dp;  //genus 35
6753ideal i=y30+y13x+x4y5+x3*(x+1)^2;
6754timeNormal(i, "normal", "normalC", "normalP");
6755
6756// Example 2
6757// char 0 : normal = 1 secs (13 steps) - normalC doesn't finish
6758// char 3 : normal = 2 secs (13 steps) - normalP = 0 secs - normalC doesn't finish
6759ring r=3,(x,y),dp;  //genus 19, delta 21
6760ideal i=y20+y13x+x4y5+x3*(x+1)^2;
6761timeNormal(i, "normal", "normalC", "normalP");
6762
6763// Example 3
6764// Very fast with all algorithms
6765ring r = 3, (x, y), dp;
6766ideal I = (x-y^2)^2-x*y^3;
6767timeNormal(I, "normal", "normalC", "normalP", "primCl", "111", "p", "pc");
6768
6769
6770
6771//----------------------Test Example for charp -------------------
6772//Zu tun:
6773//### nach minor nur std statt mstd verwenden
6774//***hat bei keinem Beisp etwas gebracht -> wieder zurueck
6775//### wenn interred ok, dann wieder einsetzen (am Schluss)
6776//### bottelnecks bei maps beheben
6777//### minor verbessern
6778//### preimage verbessern (Ist imm Kern map oder imap verwendet?)
6779//### Gleich in Ordnung dp wechseln, ringlist verwenden
6780//### interred ev nur zum Schluss
6781//    (z.B. wenn nacher std; wenn nacher minor: testen )
6782
6783//Zeiten mit normalV5.lib (mstd aktiv, interred inaktiv)
6784
6785//SWANSON EXAMPLES: (Macaulay2, icFracP=normalP, icFractions<->normal)
6786//---------------------------------------------------------------------
6787//1. Series Fp[x,y,u,v]/(x2v-y2u)
6788//-------------------------------
6789//characteristic p   2   3    5    7    11   13   17   37   97
6790//icFracP          0.04 0.03 0.04 0.04 0.04 0.05 0.05 0.13 0.59  Mac
6791//normalP           0   0    0    0     0    0    0    0   1    Sing
6792//icFractions      0.08 0.09 0.09 0.09 0.14 0.15 0.15 0.15 0.15  Mac
6793//normal             0   0    0    0     0    0    0    0    0   Sing
6794
67952. Series Fp[u, v, w, x, y, z]/u2x4+uvy4+v2z4
6796//--------------------------------------------
6797//characteristic p 2    3    5    7   11
6798//icFracP         0.07 0.22 9.67 143 12543
6799//normalP          0    0    5   42  1566
6800//icFractions     1.16   *    *   *    *       *: > 6h
6801//normal            0    0    0   0    0
6802
6803//3. Series Fp[u, v, w, x, y, z]/(u2xp+uvyp+v2zp)
6804//-----------------------------------------------
6805//characteristic p  2    3    5    7    11   13  17 19 23
6806//icFracP          0.06 0.07 0.09 0.27 1.81 4.89 26 56 225
6807//normalP          0     0    0    0    1    2  6  10  27
6808//icFractions      0.16 1.49 75.00 4009 *    *   *  *  *
6809//normal            0     0    2   836
6810//normal(neu)       0     0    1   2    10  155
6811//### p=7 normal braucht 807 sec in:
6812// ideal endid  = phi1(endid);      //### bottelneck'
6813
6814//1.
6815int p = 2;  ring r = p,(u,v,x,y,z),dp; ideal i = x2v-y2u;
6816//2.
6817int p = 7; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2x4+uvy4+v2z4;
6818//3.
6819int p=11; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2*x^p+uv*y^p+v2*z^p;
6820
6821//IRREDUCIBLE EXAMPLES:
6822//---------------------
6823//timing for MacBookPro 2.2GHz Intel Core 2 Duo, 4GB Ram
6824//Sing. ix86Mac-darwin version 3-1-0 (3100-2008101314)  Oct 13 2008 14:46:59
6825//if no time is given: < 1  sec
6826
6827//Apply:
6828list nor = normal(i,"isPrim"); nor;
6829list nor = normalP(i,"withRing","isPrim"); nor;
6830def R=nor[1][1]; setring R; norid; normap;
6831setring r;
6832norTest(i,nor);
6833
6834int tt = timer;
6835list nor = normalP(i,"withRing","isPrim"); nor;
6836timer-tt;
6837int tt = timer;
6838list nor = normal(i,"isPrim");
6839timer-tt;
6840
6841ring r=19,(x,y,u,v),dp;    //delta -1
6842ideal i=x2v-y2u;
6843//norTest 2 sec
6844
6845ring r=2,(y,x2,x1),lp;     //delta -1
6846ideal i=y^4+y^2*x2*x1+x2^3*x1^2+x2^2*x1^3;
6847//### norid hat 1 Element nach interred
6848
6849ring r  = 11,(x,y,z),wp(2,1,2); //alles < 1 sec
6850ideal i=z3 - xy4 + x2;          //not reduced, delta =0 ok
6851ideal i=y4+x5+y2x;              //not reduced, delta -1
6852//interred verkleinert norid
6853
6854ring r=3,(u,v,x,y,z),dp;   //delta -1
6855ideal i=u2x3+uvy3+v2z3;
6856
6857ring r=3,(u,v,x,y,z),dp;   //delta -1
6858ideal i=u2x4+uvy4+v2z4;
6859//norTest(i,nor);  0 sec, norTest(i,nor) haengt!
6860
6861ring r=5,(u,v,x,y,z),dp;   //delta -1
6862ideal i=u2x6+uvy6+v2z6;
6863//normalP 5sec, normalC 1sec
6864//V5: norTest(i,nor); 45 sec bei normalP, V6 12 sec
6865//28 sec bei normal
6866
6867ring r=5,(u,v,x,y,z),dp;   //delta -1
6868ideal i=u2x5+uvy5+v2z5;
6869//normalP 1sec, normalC 1 sec,
6870//norTest lange: minor(jacob(I),h,J) 193 (308)sec, haengt dann bei M = std(M);
6871//norTest(i,nor,2); verwenden!
6872//Sing 3.0-4 orig  >9h! haengt bei Q = mstd(Q)[2];
6873
6874ring r=2,(y,x),wp(12,5);  //delta 3
6875ideal i=y5+y2x4+y2x+yx2+x12;
6876//normalP 0 sec (Test 0 sec), normalC 2 sec (Test 2 sec)
6877//normalC withGens (ohne interred) 0sec
6878
6879ring r=2,(y,x),dp;       //delta= 22
6880ideal i=y9+y8x+y8+y5+y4x+y3x2+y2x3+yx8+x9;
6881//normalP 1sec, interred verkleinert norid betraechtlich
6882//normalC haengt bei minor, ideal im loop wird zu gross ###
6883//interred bei normalC vergroeesert string um Faktor 4000!
6884//withGens haengt bei interred in loop 4 (> 10 h) oder
6885//(nach Ausschalten von interred) bei
6886//int delt=vdim(std(modulo(f,ideal(p)))); (>?h)
6887
6888//Leonard1: (1. Komponente von Leonard),  delta -1
6889ring r=2,(v,u,z,y,x),dp;
6890ideal i = z3+zyx+y3x2+y2x3, uyx+z2,uz+z+y2x+yx2, u2+u+zy+zx,
6891          v3+vux+vz2+vzyx+vzx+uz3+uz2y+z3+z2yx2;
6892//normalP 5 sec (withRing 9 sec), norTest(i,nor,2); 45 sec
6893//normalC 102sec, 99sec
6894//### Zeit wird bei ideal Ann = quotient(SM[2],SL[1]); und bei
6895// f  = quotient(p*J,J); verbraucht
6896//withGens (ohne interred) 131sec, norTest(i,nor,2); 2min25sec
6897//norTest(i,nor,2);  45 sec
6898
6899 ring r=2,(y,x),wp(25,21); //Leonard2, delta 232
6900 ring r=2,(y,x),dp;
6901 ideal i=
6902 y^21+y^20*x +y^18*(x^3+x+1) +y^17*(x^3+1) +y^16*(x^4+x)
6903 +y^15*(x^7+x^6+x^3+x+1) +y^14*x^7 +y^13*(x^8+x^7+x^6+x^4+x^3+1)
6904 +y^12*(x^9+x^8+x^4+1) +y^11*(x^11+x^9+x^8+x^5+x^4+x^3+x^2)
6905 +y^10*(x^12+x^9+x^8+x^7+x^5+x^3+x+1)
6906 +y^9*(x^14+x^13+x^10+x^9+x^8+x^7+x^6+x^3+x^2+1)
6907 +y^8*(x^13+x^9+x^8+x^6+x^4+x^3+x) +y^7*(x^16+x^15+x^13+x^12+x^11+x^7+x^3+x)
6908 +y^6*(x^17+x^16+x^13+x^9+x^8+x) +y^5*(x^17+x^16+x^12+x^7+x^5+x^2+x+1)
6909 +y^4*(x^19+x^16+x^15+x^12+x^6+x^5+x^3+1)
6910 +y^3*(x^18+x^15+x^12+x^10+x^9+x^7+x^4+x)
6911 +y^2*(x^22+x^21+x^20+x^18+x^13+x^12+x^9+x^8+x^7+x^5+x^4+x^3)
6912 +y*(x^23+x^22+x^20+x^17+x^15+x^14+x^12+x^9)
6913 +(x^25+x^23+x^19+x^17+x^15+x^13+x^11+x^5);
6914//normalP: dp 2sec withRing 8sec,
6915//wp 4sec, withRing:51sec Zeit in lin = subst(lin, var(ii), vip); in elimpart ),
6916//norTest(i,nor,2): haengt bei mstd(norid);
6917//### normalC: (m. interred): haengt bei endid = interred(endid);
6918//GEFIXTES INTERRED ABWARTEN. Dann interred aktivieren
6919//interred(norid) haengt u. mst(norid) zu lange
6920//(o.interred): haengt bei  haengt bei list SM = mstd(i);
6921//ideal in der Mitte zu gross
6922//i = Ideal (size 118, 13 var) fuer die neue Normalisierung
6923//normal(neu) haengt bei return(std(i)) (offensichtlich in eineranderen lib)
6924
6925REDUCIBLE EXAMPLES:
6926------------------
6927//Apply:
6928int tt = timer;
6929list nor=normalP(i,"isPrim","withRing");
6930timer-tt;
6931
6932list nor = normal(i); nor;
6933list nor = normalC(i); nor;
6934list nor = normalC(i, "withGens"); nor;
6935list nor = normalP(i,"withRing"); nor;
6936list nor = normalP(i); nor;
6937def R=nor[1][1]; setring R; norid; normap;
6938
6939//Leonhard 4 Komponenten, dim=2, delta: 0,0,0,-1
6940ring r=2,(v,u,z,y,x),dp;      //lp zu lange
6941ideal i=z3+zyx+y3x2+y2x3, uyx+z2, v3+vuyx+vux+vzyx+vzx+uy3x2+uy2x+zy3x+zy2x2;
6942//normalP: 19 sec, withRing: 22 sec
6943//normalC ohne (mit) interred: 112 (113)sec, equidim: 99sec
6944//normalC 1. mal 111 sec, (2.mal) 450sec!! 3.mal 172 sec
6945//(unterschiedlich lange primdec, mit Auswirkungen)
6946//char 19: normalC: 15sec , withGens: 14sec (o.interr.)
6947
6948//----------------------Test Example for special cases -------------------
6949int tt = timer;
6950list nor=normalP(i,"withRing");nor;
6951//list nor=normalP(i,"withRing", "isPrim");nor;
6952timer-tt;
6953def R1 = nor[1][1]; setring R1;  norid; normap; interred(norid);
6954setring r;
6955
6956int tt = timer;
6957list nor=normal(i,"isPrim");nor;
6958timer-tt;
6959
6960ring r = 29,(x,y,z),dp;
6961ideal i = x2y2,x2z2;       //Nicht equidimensional, equidim reduziert nicht, ok
6962ideal i  = xyz*(z3-xy4);   //### interred(norid) verkuerzt
6963//je 0 sec
6964
6965ideal j = x,y;
6966ideal i = j*xy;
6967equidim(i);
6968//hat eingebettete Komponente, equidim rechnet wie in Beschreibung (ok)
6969
6970ring r  = 19,(x,y),dp;
6971   ideal i = x3-y4;                   //delta = 3
6972   ideal i = y*x*(x3-y4);             //delta = 11; 0,0,3
6973   ideal i = (x2-y3)*(x3-y4);         //delta = 13; 1,3
6974   ideal i = (x-y)*(x3+y2)*(x3-y4);   //delta = 23; 0,1,3
6975   ideal i = (x-1)*(x3+y2)*(x2-y3);   //delta = 16; 0,1,1
6976   ideal i = (x-y^2)^2 - y*x^3;       //delta = 3
6977   //singularities at not only at 0, hier rechnet equidim falsch
6978
6979// -------------------------- General Examples  ---------------------------//Huneke, irred., delta=2 (Version 3-0-4: < 1sec)
6980//Version 3-0-6 default: 1sec, mit gens 2sec, mit delta 5 sec
6981//(prim,noFac):ca 7 Min, prim:ca 10 min(wg facstd)
6982//
6983// "equidim" < 1sec irred. 5sec
6984// ring r=31991,(a,b,c,d,e),dp;
6985ring r=2,(a,b,c,d,e),dp;                    //delta=2
6986ideal i=
69875abcde-a5-b5-c5-d5-e5,
6988ab3c+bc3d+a3be+cd3e+ade3,
6989a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
6990abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
6991ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
6992a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
6993a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
6994b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
6995//normalC: char 2, 31991: 0 sec (isPrim); char 2, equidim: 7 sec
6996//norTest(i,nor,2); 1sec
6997//normalP char 2: 1sec (isPrim)
6998//size(norid); size(string(norid));21 1219 interred(norid): 21 1245 (0 sec)
6999
7000int tt = timer;
7001list nor=normalC(i);nor;
7002timer-tt;
7003
7004list nor = normalP(i,"isPrim");
7005
7006//Vasconcelos irred., delta -1 (dauert laenger)
7007//auf macbook pro = 20 sec mit alter Version,
7008//Sing 3-0-6:
7009// Char 32003: "equidim" 30 sec, "noFac": 30sec
7010//gens: nach 9 min abgebr (haengt in Lin = ideal(T*syzf);) !!!! Hans zu tun
7011//Char 2: default (charp) 2 sec, normalC ca 30 sec
7012//ring r=32003,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
7013ring r=2,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
7014ideal i= x2+zw, y3+xwt, xw3+z3t+ywt2, y2w4-xy2z2t-w3t3;
7015//normalC: char 2: 22,  sec (mit und ohne isPrim)
7016//normalP char 2: 0sec (isPrim)      o. interred
7017//char 32003: ### haengt in ideal endid  = phi1(endid);
7018
7019//-------------------------------------------------------
7020//kleine Beispiele:
7021
7022//Theo1, irred, delta=-1
7023//normalC: 1sec, normalP: 3 sec
7024ring r=32003,(x,y,z),wp(2,3,6); //dim 2,dim slocus 1
7025ideal i=zy2-zx3-x6;
7026//normalC: char 2,19,32003: 0  sec (isPrim)
7027//normalP (isPrim) char 2,19: 0sec, char 29: 1sec
7028
7029//Theo1a, CohenMacaulay regular in codim 2, dim slocus=1, delta=0
7030//normalC: 0 sec, normalP: haegt in K=preimage(R,phi,L);
7031ring r=32003,(x,y,z,u),dp;
7032ideal i=zy2-zx3-x6+u2;
7033//normalC: char 2,32003: 0  sec (isPrim)
7034//normalP (isPrim) char 2: 0sec, char 19: haengt in K = preimage(Q,phi,L);
7035
7036//Theo2, irreduzibel, reduziert, < 1sec, delta -1
7037ring r=0,(x,y,z),wp(3,4,12);
7038ideal i=z*(y3-x4)+x8;
7039//normalC: char 2,32003,0: 0  sec (isPrim)
7040//normalP (isPrim) char 2: 0 1sec, char 19: 1sec char 29: 7 sec
7041
7042//Theo2a, reduiziert, 2-dim, dim_slocus=1, alte Version 3 sec,
7043//normalP ca 30 sec, normalC ca 4sec, delta -1
7044//ring r=32003,(T(1..4)),wp(3,4,12,17);
7045//ring r=11,(T(1..4)),dp;
7046ring r=11,(T(1..4)),wp(3,4,12,17);
7047ideal i=
7048T(1)^8-T(1)^4*T(3)+T(2)^3*T(3),
7049T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4),
7050T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4),
7051T(1)^6*T(2)*T(3)+T(1)^2*T(2)^4*T(3)+T(1)^3*T(2)^2*T(4)-T(1)^2*T(2)*T(3)^2+T(4)^2;
7052//normalC: char 2,32003: 0  sec (isPrim)
7053//normalP (isPrim) char 2: 0sec, char 11 2se, char 19: 13sec
7054//norTest 48sec in char11
7055//### interred verkuerzt
7056//char 29: haengt in K = preimage(Q,phi,L);
7057
7058//Theo3, irred, 2-dim, 1-dim sing, < 1sec
7059ring r=11,(x,y,z),wp(3,5,15);
7060ideal i=z*(y3-x5)+x10;
7061//normalC: char 2,0: 0  sec (withRing)
7062//normalP (withRing) char 2,11: 0sec, char 19: 13sec norTest 12sec(char 11)
7063
7064//Theo4 reducible, delta (0,0,0) -1
7065ring r=29,(x,y,z),dp;
7066ideal i=(x-y)*(x-z)*(y-z);
7067//normalC: char 2,32003: 0  sec
7068//normalP char withRing 2, 29: 0sec, 6sec
7069
7070//Theo6
7071ring r=32003,(x,y,z),dp;
7072ideal i=x2y2+x2z2+y2z2;
7073//normalC: char 2,32003: 0  sec
7074//normalP char withRing 2, 29: 0sec, 4sec
7075
7076//Sturmfels, CM, 15 componenten, alle glatt
7077ring r=0,(b,s,t,u,v,w,x,y,z),dp;
7078ideal i= bv+su, bw+tu, sw+tv, by+sx, bz+tx, sz+ty,uy+vx,uz+wx,vz+wy,bvz;
7079//normalC car 11, 0: 1sec, normalP 0 sec
7080
7081//riemenschneider, , dim 3, 5 Komp. delta (0,0,0,0,0), -1
7082ring r=2,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
7083ideal i=xz,vx,ux,su,qu,txy,stx,qtx,uv2z-uwz,uv3-uvw,puv2-puw;
7084//alles 0 sec in char 2
7085
7086//4 Komponenten, alle glatt, 0sec
7087ring r=11,(x,y,z,t),dp;
7088ideal i=x2z+xzt,xyz,xy2-xyt,x2y+xyt;
7089
7090//dim 3, 2 Komponenten delta (-1,0), -1
7091ring r=2,(u,v,w,x,y,z),wp(1,1,1,3,2,1);
7092ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;
7093//alles 0 sec in char 2
7094//---------------------------------------------------------
7095int tt = timer;
7096list nor=normalP(i,"normalC","withRing");nor;
7097timer-tt;
7098
7099//St_S/Y, 3 Komponenten, 2 glatt, 1 normal
7100//charp haengt (in char 20) in K=preimage(R,phi,L);
7101//ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
7102ring r=11,(b,s,t,u,v,w,x,y,z),dp;
7103ideal i=wy-vz,vx-uy,tv-sw,su-bv,tuy-bvz;
7104//normalC: char 2,32003: 0  sec
7105//normalP char withRing 2: 1sec, char 11: 40sec
7106
7107//Horrocks: cahr 0: 17 (8 in char 11) Komponenten alle normal, delta 1
7108//char 11: 8 Komponenten alle normal, delta -1
7109ring r=0,(a,b,c,d,e,f),dp;
7110//ring r=11,(a,b,c,d,e,f),dp; //Charp bis p = 200 ca 3sec
7111ideal i=
7112adef-16000be2f+16001cef2, ad2f+8002bdef+8001cdf2, abdf-16000b2ef+16001bcf2,
7113a2df+8002abef+8001acf2, ad2e-8000bde2-7999cdef, acde-16000bce2+16001c2ef,
7114a2de-8000abe2-7999acef, acd2+8002bcde+8001c2df, abd2-8000b2de-7999bcdf,
7115a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2,
7116abde-8000b2e2-acdf-16001bcef-8001c2f2, abcd-16000b2ce+16001bc2f,
7117a2cd+8002abce+8001ac2f, a2bd-8000ab2e-7999abcf, ab3f-3bdf3,
7118a2b2f-2adf3-16000bef3+16001cf4, a3bf+4aef3, ac3e-10668cde3,
7119a2c2e+10667ade3+16001be4+5334ce3f, a3ce+10669ae3f, bc3d+8001cd3e,
7120ac3d+8000bc3e+16001cd2e2+8001c4f, b2c2d+16001ad4+4000bd3e+12001cd3f,
7121b2c2e-10668bc3f-10667cd2ef, abc2e-cde2f, b3cd-8000bd3f, b3ce-10668b2c2f-10667bd2ef, abc2f-cdef2, a2bce-16000be3f+16001ce2f2,
7122ab3d-8000b4e-8001b3cf+16000bd2f2, ab2cf-bdef2,
7123a2bcf-16000be2f2+16001cef3, a4d-8000a3be+8001a3cf-2ae2f2;
7124//normalC: char 0: 1sec char 11: 0sec
7125//normalP: char 11: 0sec
7126
7127//2sec mit normalC, in char 2 ebenfalls (char 20 mit charp >1 min)
7128//4 Komp. in char 2, delta (0,0,0,0) -1, char 11:delta (-1,0,0,0) -1
7129ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
7130ideal i=
7131wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3,
7132vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2,
7133tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2,
7134sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2,
7135tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3;
7136//normalC: char 2,32003: 1 sec
7137//normalP char withRing 2: 1sec, char 11: 40sec
7138
7139//---------------------------------------------------------
7140//genus:
7141int tt = timer;
7142list nor=normal(i, "noFac");nor;
7143timer-tt;
7144
7145//Yoshihiko Sakai, irred, 0sec, delta = 8
7146ring r=0,(x,y),dp;                    //genus 0 4 nodes and 6 cusps im P2
7147//ring r=7,(x,y),dp;                  //charp haengt in K = preimage(Q,phi,L)
7148ideal i=(x2+y^2-1)^3 +27x2y2;
7149
7150ring r=0,(x,y),dp;   //genus 0
7151ideal i=(x-y^2)^2 - y*x^3;
7152
7153ring r=0,(x,y),dp;  //genus 4
7154ideal i=y3-x6+1;
7155
7156int m=9;           // q=9: genus 0
7157int p=2;
7158int q=9;//2,...,9
7159ring r=0,(x,y),dp;
7160ideal i=y^m - x^p*(x - 1)^q;
7161
7162ring r=0,(x,y),dp;  //genus 19
7163ideal i=55*x^8+66*y^2*x^9+837*x^2*y^6-75*y^4*x^2-70*y^6-97*y^7*x^2;
7164
7165ring r=23,(x,y),dp;  //genus 34, delta 2
7166ideal i=y10+(-2494x2+474)*y8+(84366+2042158x4-660492)*y6
7167        +(128361096x4-47970216x2+6697080-761328152x6)*y4
7168        +(-12024807786x4-506101284x2+15052058268x6+202172841-3212x8)*y2
7169        +34263110700x4-228715574724x6+5431439286x2+201803238
7170        -9127158539954x10-3212722859346x8;
7171//normalC, normalP 0 sec
7172
7173//Rob Koelman
7174//ring r=0,(x,y,z),dp;      //dim sing = 1 (nach ca 15 min abgebrochen)
7175ring r=32003,(x,y,z),dp;
7176ideal i=
7177761328152*x^6*z^4-5431439286*x^2*y^8+2494*x^2*z^8+228715574724*x^6*y^4+
7178 9127158539954*x^10-15052058268*x^6*y^2*z^2+3212722859346*x^8*y^2-
7179 134266087241*x^8*z^2-202172841*y^8*z^2-34263110700*x^4*y^6-6697080*y^6*z^4-
7180 2042158*x^4*z^6-201803238*y^10+12024807786*x^4*y^4*z^2-128361096*x^4*y^2*z^4+
7181 506101284*x^2*z^2*y^6+47970216*x^2*z^4*y^4+660492*x^2*z^6*y^2-
7182 z^10-474*z^8*y^2-84366*z^6*y^4;
7183//normalC char 32003: 10 sec, char 0 :
7184
7185//ring r=0,(x,y),dp;//genus 10  with 26 cusps (nach ca 4 min abgebrochen)
7186ring r=32003,(x,y),dp;    //24 sing, delta 24
7187ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6
7188-5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4
7189+506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6
7190-2042158x4+660492x2y2-84366y4+2494x2-474y2-1;
7191//normalC 32003: 4 sec, char 0: abgebrochen bei pr = facstd(i); ###
7192
7193ring r=0,(x,y),dp;   //irred, genus 1  with 5 cusps, delta 5
7194ideal i=57y5+516x4y-320x4+66y4-340x2y3+73y3+128x2-84x2y2-96x2y;
7195//normalC 0 sec
7196
7197ring r=2,(x,y),dp;  //genus 4, 2 Zweige, delta (13,9) 89
7198ideal i=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
7199//normalC: char 2 : 1sec, char 0: lange
7200//normalP char 2 withRing: 0sec
7201
7202ring r=2,(y,z,w,u),dp; //2 Komp. genus -5
7203ideal i=y2+z2+w2+u2,w4-u4;
7204//normalC: char 2 : 0sec, char 0: 1sec
7205//normalP char 2 withRing: 0sec
7206
7207ring r=0,(y,z,w,u),dp; //irred. genus 9
7208ideal i=y2+z2+w2+u2,z4+w4+u4;
7209//char 0: 0sec
7210
7211ring r=0,(x,y,t),dp;  //irred, delta -1
7212ideal i= 25x8+200x7y+720x6y2+1520x5y3+2064x4y4+1856x3y5+1088x2y6+384xy7+64y8-12x6t2-72x5yt2-184x4y2t2-256x3y3t2-192x2y4t2-64xy5t2-2x4t4-8x3yt4+16xy3t4+16y4t4+4x2t6+8xyt6+8y2t6+t8;
7213//char 0: 0sec
7214
7215ring r=0,(x,y,z,w,u),dp;
7216ideal i=x2+y2+z2+w2+u2,x3+y3+z3,z4+w4+u4;
7217//char 0: 0sec
7218
7219//---------------------------------------------------------
7220//Probleme mit normalC in char 2 und char 0
7221
7222int tt = timer;
7223list nor=normalC(i,"withRing");nor;
7224timer-tt;
7225
7226//Mark van Hoeij
7227ring r=3,(x,y),dp;  //genus 19, delta 21
7228ideal i=y20+y13x+x4y5+x3*(x+1)^2;
7229//normalC: char 2 > 10 min   bei list SM = mstd(i);###
7230//normalP char 2 withRing: 0sec, char 11: haengt bei K = preimage(Q,phi,L);
7231
7232ring r=2,(x,y),dp;  //genus 35
7233ideal i=y30+y13x+x4y5+x3*(x+1)^2;
7234//char 0 abgebrochen bei list SM = mstd(i); ###
7235//char 2 nach ca 30 min
7236//normalC: char 2: abgebr. bei list SM = mstd(i);  //Now the work starts'
7237//normalC, withGens, char 2: abgebrochen bei Q=mstd(Q)[2];
7238//normalP char 2 withRing: 0sec
7239
7240ring r=0,(x,y),dp;   //irred, genus 55, delta 21
7241ideal i=y40+y13x+x4y5+x3*(x+1)^2;
7242//normalC: char 2 lange
7243//normalP char 2 withRing: 0sec
7244
7245ring r=29,(x,y,t),dp; //char 0: genus -5, 4 Komp, delta (-1,-1,0,0), -1
7246ideal i=x8+8x7y+32x6y2+80x5y3+136x4y4+160x3y5+128x2y6+64xy7+16y8+4x6t2+24x5yt2+72x4y2t2+128x3y3t2+144x2y4t2+96xy5t2+32y6t2+14x4t4+56x3yt4+112x2y2t4+112xy3t4+40y4t4+20x2t6+40xyt6+8y2t6+9t8;
7247//normalC: char 29 : 0sec, char 0: 0sec  //char 29 6 Komponenten
7248//normalP char 29 withRing: 1sec
7249
7250//-------------------------- problematic examples ------------------------
7251//ring r=0,(x,y,t),dp;
7252ring r=32003,(x,y,t),dp;
7253ideal i=
725432761x8+786264x7y+8314416x6y2+50590224x5y3+193727376x4y4+478146240x3y5+742996800x2y6+664848000xy7+262440000y8+524176x7t+11007696x6yt+99772992x5y2t+505902240x4y3t+1549819008x3y4t+2868877440x2y5t+2971987200xy6t+1329696000y7t+3674308x6t2+66137544x5yt2+499561128x4y2t2+2026480896x3y3t2+4656222144x2y4t2+5746386240xy5t2+2976652800y6t2+14737840x5t3+221067600x4yt3+1335875904x3y2t3+4064449536x2y3t3+6226336512xy4t3+3842432640y5t3+36997422x4t4+443969064x3yt4+2012198112x2y2t4+4081745520xy3t4+3126751632y4t4+59524208x3t5+535717872x2yt5+1618766208xy2t5+1641991392y3t5+59938996x2t6+359633976xyt6+543382632y2t6+34539344xt7+103618032yt7+8720497t8;
7255//char 0: lange (es liegt an den grossen Zahlen), char 32003: 0 sec
7256
7257//dasselbe Beipiel in char 19: irred
7258ring r=0,(x,y,t),dp;
7259ideal i=
72605x8+6x7y-3x6y2+7x5y3-6x4y4-8x3y5-5x2y6-8y8+4x7t+8x6yt+2x5y2t-6x4y3t+9x3y4t+9x2y5t-xy6t-7x6t2+7x5yt2-x4y2t2+5x3y3t2+7x2y4t2+xy5t2-3y6t2-4x5t3 -3x4yt3+2x3y2t3-7x2y3t3-6xy4t3-3y5t3-5x4t4-3x3yt4-4x2y2t4-8xy3t4 +7y4t4+x3t5+9x2yt5+9xy2t5-8y3t5-2y2t6+4xt7-7yt7-9t8;
7261//normalP: char 2,3: 0sec, norTest 0,2 sec, char 11 haengt bei peimage
7262//normalC: char 3: 0 sec, char 0: 1sec
7263
7264//ring r=0,(x,y),dp;
7265ring r=32003,(x,y),dp;
7266ideal i=
7267x30y21+21x29y20+210x28y19+10x27y19+1330x27y18+190x26y18+5985x26y17
7268+1710x25y17+20349x25y16+45x24y17+9690x24y16+54264x24y15+765x23y16
7269+38760x23y15+116280x23y14+6120x22y15+116280x22y14+120x21y15
7270+203490x22y13+30600x21y14+271320x21y13+1799x20y14+293930x21y12+107100x20y13
7271+503880x20y12+12586x19y13+352716x20y11+278460x19y12+210x18y13+755820x19y11
7272+54509x18y12+352716x19y10+556920x18y11+2723x17y12+923780x18y10+163436x17y11
7273+293930x18y9+875160x17y10+16296x16y11+923780x17y9+359359x16y10+252x15y11
7274+203490x17y8+1093950x16y9+59598x15y10+755820x16y8+598598x15y9+2751x14y10
7275+116280x16y7+1093950x15y8+148610x14y9+503880x15y7+769197x14y8+13650x13y9
7276+54264x15y6+875160x14y7+266805x13y8+210x12y9+271320x14y6+768768x13y7
7277+40635x12y8+20349x14y5+556920x13y6+354816x12y7+1855x11y8+116280x13y5
7278+597597x12y6+80640x11y7+5985x13y4+278460x12y5+353892x11y6+7280x10y7+38760x12y4
7279+358358x11y5+112014x10y6+120x9y7+1330x12y3+107100x11y4+264726x10y5+16660x9y6
7280+9690x11y3+162799x10y4+111132x9y5+805x8y6+210x11y2+30600x10y3+146685x9y4
7281+24500x8y5+1710x10y2+54236x9y3+78750x8y4+2310x7y5+21x10y+6120x9y2+58520x8y3
7282+24010x7y4+45x6y5+190x9y+12509x8y2+39060x7y3+3675x6y4+x9+765x8y+15918x7y2
7283+15680x6y3+204x5y4+10x8+1786x7y+12915x6y2+3500x5y3+45x7+2646x6y+6580x5y2
7284+366x4y3+119x6+2562x5y+1995x4y2+10x3y3+203x5+1610x4y+324x3y2+231x4+630x3y
7285+23x2y2+175x3+141x2y+85x2+16xy+24x+y+4;
7286list nor = normal(i);
7287//normalC: char 0: ### haengt in SB of singular locus JM = mstd(J);
7288//normalC: char 32003,"noFac","equidim": 0sec, "noFac": 1sec
7289// ev neues interred
7290genus(i);         // haengt bei int mu=vdim(std(jacob(f)));
7291                  //### ist das noetig?
7292
7293//Singular rechnet genus richtig, auch im Fall, dass Kurve irreduzibel,
7294//aber nicht absolut irreduzibel ist:
7295ring r = 0,(x,y),dp;
7296ideal i = x2+y2;      //irreduzibel /Q aber reduzibel /C (x-iy)*(x+iy)
7297factorize( x2+y2);    //liefert irreduzibel
7298genus(i);             //sollte 0+0-2+1= -1 sein
7299genus(i,1);           //beides ist korrekt in Singular
7300
7301*/
Note: See TracBrowser for help on using the repository browser.