source: git/Singular/LIB/normal.lib @ 19190c3

fieker-DuValspielwiese
Last change on this file since 19190c3 was 19190c3, checked in by Hans Schoenemann <hannes@…>, 10 years ago
format
  • Property mode set to 100644
File size: 219.6 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  ASSUME(0, not isQuotientRing(basering) );
148
149  intvec opt = option(get);     // Save current options
150
151  int i,j;
152  int decomp;   // Preliminary decomposition:
153                // 0 -> no decomposition (id is assumed to be prime)
154                // 1 -> no decomposition
155                //      (id is assumed to be equidimensional radical)
156                // 2 -> equidimensional decomposition
157                // 3 -> minimal associated primes
158  int noFac, useRing, withDelta;
159  int dbg = printlevel - voice + 2;
160  int nvar = nvars(basering);
161  int chara  = char(basering);
162  int denomOption;   // Method for choosing the conductor
163
164  ideal inputJ = 0;      // Test ideal given in the input (if any).
165  ideal inputC = 0;      // Conductor ideal given in the input (if any).
166
167  list result, resultNew;
168  list keepresult;
169  list ringStruc;
170  ideal U;
171  poly c;
172  int sp;            // Number of components.
173
174  // Default methods:
175  noFac = 0;         // Use facSTD when computing minimal associated primes
176  decomp = 2;        // Equidimensional decomposition
177  useRing = 0;       // Change first to dp ordering, and perform all
178                     // computations there.
179  withDelta = 0;     // Do not compute the delta invariant.
180  denomOption = 0;   // The default universal denominator is the smallest
181                     // degree polynomial.
182
183//--------------------------- define the method ---------------------------
184  for ( i=1; i <= size(#); i++ )
185  {
186    if ( typeof(#[i]) == "string" )
187    {
188//--------------------------- choosen methods -----------------------
189      if ( (#[i]=="isprim") or (#[i]=="isPrim") )
190      {decomp = 0;}
191
192      if ( (#[i]=="nodeco") or (#[i]=="noDeco") )
193      {decomp = 1;}
194
195      if (#[i]=="prim")
196      {decomp = 3;}
197
198      if (#[i]=="equidim")
199      {decomp = 2;}
200
201      if ( (#[i]=="nofac") or (#[i]=="noFac") )
202      {noFac=1;}
203
204      if ( ((#[i]=="useRing") or (#[i]=="usering")) and (ordstr(basering) != "dp("+string(nvars(basering))+"),C"))
205      {useRing = 1;}
206
207      if ( (#[i]=="withDelta") or (#[i]=="wd") or (#[i]=="withdelta"))
208      {
209        if((decomp == 0) or (decomp == 3))
210        {
211          withDelta = 1;
212        }
213        else
214        {
215          decomp = 3;
216          withDelta = 1;
217          //Note: the delta invariants cannot be computed with an equidimensional
218          //decomposition, hence we compute first the minimal primes
219        }
220      }
221      if (#[i]=="var1")
222      {denomOption = 1;}
223      if (#[i]=="var2")
224      {denomOption = 2;}
225    }
226    if(typeof(#[i]) == "list"){
227      if(size(#[i]) == 2){
228        if (#[i][1]=="inputJ"){
229          if(typeof(#[i][2]) == "ideal"){
230            inputJ = #[i][2];
231          }
232        }
233      }
234      if (#[i][1]=="inputC"){
235        if(size(#[i]) == 2){
236          if(typeof(#[i][2]) == "ideal"){
237            inputC = #[i][2];
238          }
239        }
240      }
241    }
242  }
243  kill #;
244
245//------------------------ change ring if required ------------------------
246// If the ordering is not global, we change to dp ordering for computing the
247// min ass primes.
248// If the ordering is global, but not dp, and useRing = 0, we also change to
249// dp ordering.
250
251  int isGlobal = attrib(basering,"global");// Checks if the original ring has
252                                          // global ordering.
253
254  def origR = basering;   // origR is the original ring
255                          // R is the ring where computations will be done
256
257  if((useRing  == 1) and (isGlobal == 1))
258  {
259    def globR = basering;
260  }
261  else
262  {
263    // We change to dp ordering.
264    list rl = ringlist(origR);
265    list origOrd = rl[3];
266    list newOrd = list("dp", intvec(1:nvars(origR))), list("C", 0);
267    rl[3] = newOrd;
268    def globR = ring(rl);
269    setring globR;
270    ideal id = fetch(origR, id);
271  }
272
273//------------------------ trivial checkings ------------------------
274  id = groebner(id);
275  if((size(id) == 0) or (id[1] == 1))
276  {
277    // The original ring R/I was normal. Nothing to do.
278    // We define anyway a new ring, equal to R, to be able to return it.
279    setring origR;
280    list lR = ringlist(origR);
281    def ROut = ring(lR);
282    setring ROut;
283    ideal norid = fetch(origR, id);
284    ideal normap = maxideal(1);
285    export norid;
286    export normap;
287    setring origR;
288    if(withDelta)
289    {
290      result = list(list(ROut), list(ideal(1)), list(intvec(0), 0));
291    }
292    else
293    {
294      result = list(list(ROut), list(ideal(1)));
295    }
296    sp = 1;      // number of rings in the output
297    option(set, opt);
298    normalOutputText(dbg, withDelta, sp);
299    return(result);
300  }
301//------------------------ preliminary decomposition-----------------------
302  list prim;
303  if(decomp == 2)
304  {
305    dbprint(dbg, "// Computing the equidimensional decomposition...");
306    prim = equidim(id);
307  }
308  if((decomp == 0) or (decomp == 1))
309  {
310    prim = id;
311  }
312  if(decomp == 3)
313  {
314    dbprint(dbg, "// Computing the minimal associated primes...");
315    if( noFac )
316    { prim = minAssGTZ(id,1); }
317    else
318    { prim = minAssGTZ(id); }
319  }
320  sp = size(prim);
321  if(dbg>=1)
322  {
323    prim; "";
324    "// number of components is", sp;
325    "";
326  }
327
328
329//----------------- back to the original ring if required ------------------
330// if ring was not global and useRing is on, we go back to the original ring
331  if((useRing == 1) and (isGlobal != 1))
332  {
333    setring origR;
334    def R = basering;
335    list prim = fetch(globR, prim);
336  }
337  else
338  {
339    def R = basering;
340    ideal inputJ = fetch(origR, inputJ);
341    ideal inputC = fetch(origR, inputC);
342    if(useRing == 0)
343    {
344      ideal U;
345      poly c;
346    }
347  }
348
349// ---------------- normalization of the components-------------------------
350// calls normalM to compute the normalization of each component.
351
352  list norComp;       // The normalization of each component.
353  int delt;
354  int deltI = 0;
355  int totalComps = 0;
356
357  setring origR;
358  def newROrigOrd;
359  list newRListO;
360  setring R;
361  def newR;
362  list newRList;
363
364  for(i=1; i<=size(prim); i++)
365  {
366    if(dbg>=2){pause();}
367    if(dbg>=1)
368    {
369      "// start computation of component",i;
370      "   --------------------------------";
371    }
372    if(groebner(prim[i])[1] != 1)
373    {
374      if(dbg>=2)
375      {
376        "We compute the normalization in the ring"; basering;
377      }
378      printlevel = printlevel + 1;
379      norComp = normalM(prim[i], decomp, withDelta, denomOption, inputJ, inputC);
380      printlevel = printlevel - 1;
381      for(j = 1; j <= size(norComp); j++)
382      {
383        newR = norComp[j][3];
384        if(!defined(savebasering)) { def savebasering;}
385        savebasering=basering;
386        setring newR; // must be in a compatible ring to newR
387                      // as ringlist may produce ring-dep. stuff
388        if(!defined(newRList)) { list newRList;}
389        newRList = ringlist(newR);
390        setring savebasering;
391        U = norComp[j][1];
392        c = norComp[j][2];
393        if(withDelta)
394        {
395          delt = norComp[j][4];
396          if((delt >= 0) and (deltI >= 0))
397          {
398            deltI = deltI + delt;
399          }
400          else
401          {
402            deltI = -1;
403          }
404        }
405        // -- incorporate result for this component to the list of results ---
406        if(useRing == 0)
407        {
408          // We go back to the original ring.
409          setring origR;
410          U = fetch(R, U);
411          c = fetch(R, c);
412          newRListO = imap(newR, newRList);
413          // We change the ordering in the new ring.
414          if(nvars(newR) > nvars(origR))
415          {
416            newRListO[3]=insert(origOrd, newRListO[3][1]);
417          }
418          else
419          {
420            newRListO[3] = origOrd;
421          }
422          newROrigOrd = ring(newRListO);
423          setring newROrigOrd;
424          ideal norid = imap(newR, norid);
425          ideal normap = imap(newR, normap);
426          export norid;
427          export normap;
428          setring origR;
429          totalComps++;
430          result[totalComps] = list(U, c, newROrigOrd);
431          if(withDelta)
432          {
433            result[totalComps] = insert(result[totalComps], delt, 3);
434          }
435          setring R;
436        }
437        else
438        {
439          setring R;
440          totalComps++;
441          result[totalComps] = norComp[j];
442        }
443      }
444    }
445  }
446
447// -------------------------- delta computation ----------------------------
448  if(withDelta == 1)
449  {
450    // Intersection multiplicities of list prim, sp=size(prim).
451    if ( dbg >= 1 )
452    {
453      "// Sum of delta for all components: ", deltI;
454    }
455    if(size(prim) > 1)
456    {
457      dbprint(dbg, "// Computing the sum of the intersection multiplicities of the components...");
458      int mul = iMult(prim);
459      if ( mul < 0 )
460      {
461        deltI = -1;
462      }
463      else
464      {
465        deltI = deltI + mul;
466      }
467      if ( dbg >= 1 )
468      {
469        "// Intersection multiplicity is : ", mul;
470      }
471    }
472  }
473
474// -------------------------- prepare output ------------------------------
475  setring origR;
476
477  list RL;      // List of rings
478  list MG;      // Module generators
479  intvec DV;    // Vector of delta's of each component
480  for(i = 1; i <= size(result); i++)
481  {
482    RL[i] = result[i][3];
483    MG[i] = lineUpLast(result[i][1], result[i][2]);
484    if(withDelta)
485    {
486      DV[i] = result[i][4];
487    }
488  }
489  if(withDelta)
490  {
491    resultNew = list(RL, MG, list(DV, deltI));
492  }
493  else
494  {
495    resultNew = list(RL, MG);
496  }
497  sp = size(RL);              //RL = list of rings
498
499  option(set, opt);
500  normalOutputText(dbg, withDelta, sp);
501  return(resultNew);
502}
503
504example
505{ "EXAMPLE:";
506  printlevel = printlevel+1;
507  echo = 2;
508  ring s = 0,(x,y),dp;
509  ideal i = (x2-y3)*(x2+y2)*x;
510  list nor = normal(i, "withDelta", "prim");
511  nor;
512
513  // 2 branches have delta = 1, and 1 branch has delta = 0
514  // the total delta invariant is 13
515
516  def R2 = nor[1][2];  setring R2;
517  norid; normap;
518
519  echo = 0;
520  printlevel = printlevel-1;
521  pause("   hit return to continue"); echo=2;
522
523  ring r = 2,(x,y,z),dp;
524  ideal i = z3-xy4;
525  list nor = normal(i, "withDelta", "prim");  nor;
526  // the delta invariant is infinite
527  // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module
528  // in its quotient field Quot(r/i)
529
530  // the normalization as affine algebra over the ground field:
531  def R = nor[1][1]; setring R;
532  norid; normap;
533}
534
535///////////////////////////////////////////////////////////////////////////////
536// Prints the output text in proc normal.
537//
538static proc normalOutputText(int dbg, int withDelta, int sp)
539// int dbg: printlevel
540// int withDelta: output contains information about the delta invariant
541// int sp: number of output rings.
542{
543  if ( dbg >= 0 )
544  {
545    "";
546    if(!withDelta)
547    {
548      "// 'normal' created a list, say nor, of two elements.";
549    }
550    else
551    {
552      "// 'normal' created a list, say nor, of three elements.";
553    }
554    "// To see the list type";
555    "      nor;";
556    "";
557    "// * nor[1] is a list of", sp, "ring(s).";
558    "// To access the i-th ring nor[1][i], give it a name, say Ri, and type";
559    "     def R1 = nor[1][1]; setring R1; norid; normap;";
560    "// For the other rings type first (if R is the name of your base ring)";
561    "     setring R;";
562    "// and then continue as for R1.";
563    "// Ri/norid is the affine algebra of the normalization of R/P_i where";
564    "// P_i is the i-th component of a decomposition of the input ideal id";
565    "// and normap the normalization map from R to Ri/norid.";
566    "";
567    "// * nor[2] is a list of", sp, "ideal(s). Let ci be the last generator";
568    "// of the ideal nor[2][i]. Then the integral closure of R/P_i is";
569    "// generated as R-submodule of the total ring of fractions by";
570    "// 1/ci * nor[2][i].";
571
572    if(withDelta)
573    { "";
574      "// * nor[3] is a list of an intvec of size", sp, "the delta invariants ";
575      "// of the components, and an integer, the total delta invariant ";
576      "// of R/id (-1 means infinite, and 0 that R/P_i resp. R/id is normal).";
577    }
578  }
579}
580
581
582///////////////////////////////////////////////////////////////////////////////
583
584proc HomJJ (list Li)
585"USAGE:   HomJJ (Li);  Li = list: ideal SBid, ideal id, ideal J, poly p
586ASSUME:  R    = P/id,  P = basering, a polynomial ring, id an ideal of P,
587@*       SBid = standard basis of id,
588@*       J    = ideal of P containing the polynomial p,
589@*       p    = nonzero divisor of R
590COMPUTE: Endomorphism ring End_R(J)=Hom_R(J,J) with its ring structure as
591         affine ring, together with the map R --> Hom_R(J,J) of affine rings,
592         where R is the quotient ring of P modulo the standard basis SBid.
593RETURN:  a list l of three objects
594@format
595         l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
596               such that l[1]/endid = Hom_R(J,J) and
597               endphi describes the canonical map R -> Hom_R(J,J)
598         l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
599         l[3] : an integer, = dim_K(Hom_R(J,J)/R) (the contribution to delta)
600                if the dimension is finite, -1 otherwise
601@end format
602NOTE:    printlevel >=1: display comments (default: printlevel=0)
603EXAMPLE: example HomJJ;  shows an example
604"
605{
606   ASSUME(0, not isQuotientRing(basering) );
607
608//---------- initialisation ---------------------------------------------------
609   int isIso,isPr,isHy,isCo,isRe,isEq,oSAZ,ii,jj,q,y;
610   intvec rw,rw1;
611   list L;
612   y = printlevel-voice+2;  // y=printlevel (default: y=0)
613   def P = basering;
614   ideal SBid, id, J = Li[1], Li[2], Li[3];
615   poly p = Li[4];
616   int noRed = 0;
617   if(size(Li) > 4)
618   {
619     if(Li[5] == 1) { noRed = 1; }
620   }
621
622   attrib(SBid,"isSB",1);
623   int homo = homog(Li[2]);               //is 1 if id is homogeneous, 0 if not
624
625//---- set attributes for special cases where algorithm can be simplified -----
626   if( homo==1 )
627   {
628      rw = ringweights(P);
629   }
630   if( typeof(attrib(id,"isPrim"))=="int" )
631   {
632      if(attrib(id,"isPrim")==1)  { isPr=1; }
633   }
634   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
635   {
636      if(attrib(id,"onlySingularAtZero")==1){oSAZ=1; }
637   }
638   if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
639   {
640      if(attrib(id,"isIsolatedSingularity")==1) { isIso=1; }
641   }
642   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
643   {
644      if(attrib(id,"isCohenMacaulay")==1) { isCo=1; }
645   }
646   if( typeof(attrib(id,"isRegInCodim2"))=="int" )
647   {
648      if(attrib(id,"isRegInCodim2")==1) { isRe=1; }
649   }
650   if( typeof(attrib(id,"isEquidimensional"))=="int" )
651   {
652      if(attrib(id,"isEquidimensional")==1) { isEq=1; }
653   }
654//-------------------------- go to quotient ring ------------------------------
655   qring R  = SBid;
656   ideal id = fetch(P,id);
657   ideal J  = fetch(P,J);
658   poly p   = fetch(P,p);
659   ideal f,rf,f2;
660   module syzf;
661//---------- computation of p*Hom(J,J) as R-ideal -----------------------------
662   if ( y>=1 )
663   {
664     "// compute p*Hom(J,J) = p*J:J";
665     "//   the ideal J:";J;
666   }
667   f  = quotient(p*J,J);
668
669   //### (neu GMG 4.10.08) divide by the greatest common divisor:
670   poly gg = gcd( f[1],p );
671   for(ii=2; ii <=ncols(f); ii++)
672   {
673      gg=gcd(gg,f[ii]);
674   }
675   for(ii=1; ii<=ncols(f); ii++)
676   {
677      f[ii]=f[ii]/gg;
678   }
679   p = p/gg;
680
681   if ( y>=1 )
682   {
683      "//   the non-zerodivisor p:"; p;
684      "//   the module p*Hom(J,J) = p*J:J :"; f;
685      "";
686   }
687   f2 = std(p);
688
689//---------- Test: Hom(J,J) == R ?, if yes, go home ---------------------------
690
691   //rf = interred(reduce(f,f2));
692   //### interred hier weggelassen, unten zugefuegt
693   rf = reduce(f,f2);       //represents p*Hom(J,J)/p*R = Hom(J,J)/R
694   if ( size(rf) == 0 )
695   {
696      if ( homog(f) && find(ordstr(basering),"s")==0 )
697      {
698         ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp);
699      }
700      else
701      {
702         ring newR1 = char(P),(X(1..nvars(P))),dp;
703      }
704      ideal endphi = maxideal(1);
705      ideal endid = fetch(P,id);
706      endid = simplify(endid,2);
707      L = substpart(endid,endphi,homo,rw);   //## hier substpart
708      def lastRing = L[1];
709      setring lastRing;
710
711      attrib(endid,"onlySingularAtZero",oSAZ);
712      attrib(endid,"isCohenMacaulay",isCo);
713      attrib(endid,"isPrim",isPr);
714      attrib(endid,"isIsolatedSingularity",isIso);
715      attrib(endid,"isRegInCodim2",isRe);
716      attrib(endid,"isEqudimensional",isEq);
717      attrib(endid,"isHypersurface",0);
718      attrib(endid,"isCompleteIntersection",0);
719      attrib(endid,"isRadical",0);
720      L=lastRing;
721      L = insert(L,1,1);
722      dbprint(y,"// case R = Hom(J,J)");
723      if(y>=1)
724      {
725         "//   R=Hom(J,J)";
726         lastRing;
727         "//   the new ideal";
728         endid;
729         "   ";
730         "//   the old ring";
731         P;
732         "//   the old ideal";
733         setring P;
734         id;
735         "   ";
736         setring lastRing;
737         "//   the map to the new ring";
738         endphi;
739         "   ";
740         pause();
741         "";
742      }
743      setring P;
744      L[3]=0;
745      return(L);
746   }
747   if(y>=1)
748   {
749      "// R is not equal to Hom(J,J), we have to try again";
750      pause();
751      "";
752   }
753//---------- Hom(J,J) != R: create new ring and map from old ring -------------
754// the ring newR1/SBid+syzf will be isomorphic to Hom(J,J) as R-module
755// f2=p (i.e. ideal generated by p)
756
757   //f = mstd(f)[2];              //### geaendert GMG 04.10.08
758   //ideal ann = quotient(f2,f);  //### f durch rf ersetzt
759   rf = mstd(rf)[2];              //rf = NF(f,p), hence <p,rf> = <p,f>
760   ideal ann = quotient(f2,rf);   //p:f = p:rf
761
762   //------------- compute the contribution to delta ----------
763   //delt=dim_K(Hom(JJ)/R (or -1 if infinite)
764
765   int delt=vdim(std(modulo(f,ideal(p))));
766
767   f = p,rf;          // generates pJ:J mod(p), i.e. p*Hom(J,J)/p*R as R-module
768   q = size(f);
769   syzf = syz(f);
770
771   if ( homo==1 )
772   {
773      rw1 = rw,0;
774      for ( ii=2; ii<=q; ii++ )
775      {
776         rw  = rw, deg(f[ii])-deg(f[1]);
777         rw1 = rw1, deg(f[ii])-deg(f[1]);
778      }
779      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp);
780   }
781   else
782   {
783      ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp;
784   }
785
786   //map psi1 = P,maxideal(1);          //### psi1 durch fetch ersetzt
787   //ideal SBid = psi1(SBid);
788   ideal SBid = fetch(P,SBid);
789   attrib(SBid,"isSB",1);
790
791   qring newR = std(SBid);
792
793   //map psi = R,ideal(X(1..nvars(R)));  //### psi durch fetch ersetzt
794   //ideal id = psi(id);
795   //ideal f = psi(f);
796   //module syzf = psi(syzf);
797   ideal id = fetch(R,id);
798   ideal f = fetch(R,f);
799   module syzf = fetch(R,syzf);
800   ideal pf,Lin,Quad,Q;
801   matrix T,A;
802   list L1;
803
804//---------- computation of Hom(J,J) as affine ring ---------------------------
805// determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1),
806// Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. This is of course
807// the same as the kernel of R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi.
808// It is a fact, that the kernel is generated by the linear and the quadratic
809// relations
810// f=p,rf, rf=reduce(f,p), generates pJ:J mod(p),
811// i.e. p*Hom(J,J)/p*R as R-module
812
813   pf = f[1]*f;
814   T = matrix(ideal(T(1..q)),1,q);
815   Lin = ideal(T*syzf);
816   if(y>=1)
817   {
818      "// the ring structure of Hom(J,J) as R-algebra";
819      "//   the linear relations:";
820      Lin;
821   }
822
823   poly ff;
824   for (ii=2; ii<=q; ii++ )
825   {
826      for ( jj=2; jj<=ii; jj++ )
827      {
828         ff = NF(f[ii]*f[jj],std(0));       //this makes lift much faster
829         A = lift(pf,ff);                   //ff lin. comb. of elts of pf mod I
830         Quad = Quad, ideal(T(jj)*T(ii) - T*A);  //quadratic relations
831      }
832   }
833
834   if(y>=1)
835   {
836      "//   the quadratic relations";
837      Quad;
838      pause();
839      newline;
840   }
841   Q = Lin,Quad;
842   Q = subst(Q,T(1),1);
843   //Q = mstd(Q)[2];            //### sehr aufwendig, daher weggelassen (GMG)
844   //### ev das neue interred
845   //mstd dient nur zum verkleinern, die SB-Eigenschaft geht spaeter verloren
846   //da in neuen Ring abgebildet und mit id vereinigt
847
848//---------- reduce number of variables by substitution, if possible ----------
849   if (homo==1)
850   {
851      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),(a(rw),dp);
852   }
853   else
854   {
855      ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp;
856   }
857
858   ideal endid  = imap(newR,id),imap(newR,Q);
859   //hier wird Q weiterverwendet, die SB-Eigenschaft wird nicht verwendet.
860   endid = simplify(endid,2);
861   ideal endphi = ideal(X(1..nvars(R)));
862
863
864  if(noRed == 0)
865  {
866    L = substpart(endid,endphi,homo,rw);
867    def lastRing=L[1];
868    setring lastRing;
869    //return(lastRing);
870  }
871  else
872  {
873    list RL = ringlist(newRing);
874    def lastRing = ring(RL);
875    setring lastRing;
876    ideal endid = fetch(newRing, endid);
877    ideal endphi = fetch(newRing, endphi);
878    export(endid);
879    export(endphi);
880    //def lastRing = newRing;
881    //setring R;
882    //return(newR);
883  }
884
885
886//   L = substpart(endid,endphi,homo,rw);
887
888//   def lastRing=L[1];
889//   setring lastRing;
890
891   attrib(endid,"onlySingularAtZero",0);
892   map sigma=R,endphi;
893   ideal an=sigma(ann);
894   export(an);  //noetig?
895   //ideal te=an,endid;
896   //if(isIso && (size(reduce(te,std(maxideal(1))))==0))   //#### ok???
897   // {
898   //    attrib(endid,"onlySingularAtZero",oSAZ);
899   // }
900   //kill te;
901   attrib(endid,"isCohenMacaulay",isCo);                //#### ok???
902   attrib(endid,"isPrim",isPr);
903   attrib(endid,"isIsolatedSingularity",isIso);
904   attrib(endid,"isRegInCodim2",isRe);
905   attrib(endid,"isEquidimensional",isEq);
906   attrib(endid,"isHypersurface",0);
907   attrib(endid,"isCompleteIntersection",0);
908   attrib(endid,"isRadical",0);
909   if(y>=1)
910   {
911      "// the new ring after reduction of the number of variables";
912      lastRing;
913      "//   the new ideal";
914      endid;  "";
915      "// the old ring";
916      P;
917      "//   the old ideal";
918      setring P;
919      id;
920      "   ";
921      setring lastRing;
922      "//   the map to the new ring";
923      endphi;
924      "   ";
925      pause();
926      "";
927   }
928   L = lastRing;
929   L = insert(L,0,1);
930   L[3] = delt;
931   setring(P);
932   return(L);
933}
934example
935{"EXAMPLE:";  echo = 2;
936  ring r   = 0,(x,y),wp(2,3);
937  ideal id = y^2-x^3;
938  ideal J  = x,y;
939  poly p   = x;
940  list Li  = std(id),id,J,p;
941  list L   = HomJJ(Li);
942  def end = L[1];    // defines ring L[1], containing ideals endid, endphi
943  setring end;       // makes end the basering
944  end;
945  endid;             // end/endid is isomorphic to End(r/id) as ring
946  map psi = r,endphi;// defines the canonical map r/id -> End(r/id)
947  psi;
948  L[3];              // contribution to delta
949}
950
951
952///////////////////////////////////////////////////////////////////////////////
953//compute intersection multiplicities as needed for delta(I) in
954//normalizationPrimes and normalP:
955
956proc iMult (list prim)
957"USAGE:   iMult(L);  L a list of ideals
958RETURN:  int, the intersection multiplicity of the ideals of L;
959         if iMult(L) is infinite, -1 is returned.
960THEORY:  If r=size(L)=2 then iMult(L) = vdim(std(L[1]+L[2])) and in general
961         iMult(L) = sum{ iMult(L[j],Lj) | j=1..r-1 } with Lj the intersection
962         of L[j+1],...,L[r]. If I is the intersection of all ideals in L then
963         we have delta(I) = delta(L[1])+...+delta(L[r]) + iMult(L) where
964         delta(I) = vdim (normalisation(R/I)/(R/I)), R the basering.
965EXAMPLE: example iMult; shows an example
966"
967{
968     ASSUME(0, not isQuotientRing(basering) );
969
970     int i,mul,mu;
971     int sp = size(prim);
972     int y = printlevel-voice+2;
973     if ( sp > 1 )
974     {
975        ideal I(sp-1) = prim[sp];
976        mu = vdim(std(I(sp-1)+prim[sp-1]));
977        mul = mu;
978        if ( y>=1 )
979        {
980          "// intersection multiplicity of component",sp,"with",sp-1,":"; mu;
981        }
982        if ( mu >= 0 )
983        {
984           for (i=sp-2; i>=1 ; i--)
985           {
986              ideal I(i) = intersect(I(i+1),prim[i+1]);
987              mu = vdim(std(I(i)+prim[i]));
988              if ( mu < 0 )
989              {
990                break;
991              }
992              mul = mul + mu;
993              if ( y>=1 )
994              {
995                "// intersection multiplicity of components",sp,"...",i+1,"with",i; mu;
996              }
997           }
998        }
999     }
1000     return(mul);
1001}
1002example
1003{ "EXAMPLE:"; echo = 2;
1004   ring s  = 23,(x,y),dp;
1005   list L = (x-y),(x3+y2);
1006   iMult(L);
1007   L = (x-y),(x3+y2),(x3-y4);
1008   iMult(L);
1009}
1010///////////////////////////////////////////////////////////////////////////////
1011//check if I has a singularity only at zero, as needed in normalizationPrimes
1012
1013proc locAtZero (ideal I)
1014"USAGE:   locAtZero(I);  I = ideal
1015RETURN:  int, 1 if I has only one point which is located at zero, 0 otherwise
1016ASSUME:  I is given as a standard bases in the basering
1017NOTE:    only useful in affine rings, in local rings vdim does the check
1018EXAMPLE: example locAtZero; shows an example
1019"
1020{
1021   ASSUME(0, not isQuotientRing(basering) );
1022
1023   int ii,jj, caz;                   //caz: conzentrated at zero
1024   int dbp = printlevel-voice+2;
1025   int nva = nvars(basering);
1026   int vdi = vdim(I);
1027   if ( vdi < 0 )
1028   {
1029      if (dbp >=1)
1030      { "// non-isolated singularitiy";""; }
1031      return(caz);
1032   }
1033
1034   //Now the ideal is 0-dim
1035   //First an easy test
1036   //If I is homogenous and not constant it is concentrated at 0
1037   if( homog(I)==1 && size(jet(I,0))==0)
1038   {
1039      caz=1;
1040      if (dbp >=1)
1041      { "// isolated singularity and homogeneous";""; }
1042      return(caz);
1043   }
1044
1045   //Now the general case with I 0-dim. Choose an appropriate power pot,
1046   //and check each variable x whether x^pot is in I.
1047   int mi1 = mindeg1(lead(I));
1048   int pot = vdi;
1049   if ( (mi1+(mi1==1))^2 < vdi )
1050   {
1051      pot = (mi1+(mi1==1))^2;      //### alternativ: pot = vdi lassen
1052   }
1053
1054   while ( 1 )
1055   {
1056      caz = 1;
1057      for ( ii=1; ii<= nva; ii++ )
1058      {
1059        if ( NF(var(ii)^pot,I) != 0 )
1060        {
1061           caz = 0; break;
1062        }
1063      }
1064      if ( caz == 1 || pot >= vdi )
1065      {
1066        if (dbp >=1)
1067        {
1068          "// mindeg, exponent, vdim used in 'locAtZero':", mi1,pot,vdi; "";
1069        }
1070        return(caz);
1071      }
1072      else
1073      {
1074        if ( pot^2 < vdi )
1075        { pot = pot^2; }
1076        else
1077        { pot = vdi; }
1078      }
1079   }
1080}
1081example
1082{ "EXAMPLE:"; echo = 2;
1083   ring r = 0,(x,y,z),dp;
1084   poly f = z5+y4+x3+xyz;
1085   ideal i = jacob(f),f;
1086   i=std(i);
1087   locAtZero(i);
1088   i= std(i*ideal(x-1,y,z));
1089   locAtZero(i);
1090}
1091
1092///////////////////////////////////////////////////////////////////////////////
1093
1094//The next procedure normalizationPrimes computes the normalization of an
1095//irreducible or an equidimensional ideal i.
1096//- If i is irreducuble, then the returned list, say nor, has size 2
1097//with nor[1] the normalization ring and nor[2] the delta invariant.
1098//- If i is equidimensional, than the "splitting tools" can create a
1099//decomposition of i and nor can have more than 1 ring.
1100
1101static proc normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #)
1102"USAGE:   normalizationPrimes(i,ihp,delt[,si]);  i = equidimensional ideal,
1103         ihp = map (partial normalization), delt = partial delta-invariant,
1104         si = ideal s.t. V(si) contains singular locus (optional)
1105RETURN:   a list of rings, say nor, and an integer, the delta-invariant
1106          at the end of the list.
1107          each ring nor[j], j = 1..size(nor)-1, contains two ideals
1108          with given names norid and normap such that
1109           - the direct sum of the rings nor[j]/norid is
1110             the normalization of basering/i;
1111           - normap gives the normalization map from basering/id
1112             to nor[j]/norid (for each j)
1113          nor[size(nor)] = dim_K(normalisation(P/i) / (P/i)) is the
1114          delta-invariant, where P is the basering.
1115EXAMPLE: example normalizationPrimes; shows an example
1116"
1117{
1118   ASSUME(1, not isQuotientRing(basering) );
1119   //Note: this procedure calls itself as long as the test for
1120   //normality, i.e if R==Hom(J,J), is negative.
1121
1122   int printlev = printlevel;   //store printlevel in order to reset it later
1123   int y = printlevel-voice+2;  // y=printlevel (default: y=0)
1124   if(y>=1)
1125   {
1126     "";
1127     "// START a normalization loop with the ideal";
1128     i;  "";
1129     "// in the ring:";
1130     basering;  "";
1131     pause();
1132     "";
1133   }
1134
1135   def BAS=basering;
1136   list result,keepresult1,keepresult2,JM,gnirlist;
1137   ideal J,SB,MB;
1138   int depth,lauf,prdim,osaz;
1139   int ti=timer;
1140
1141   gnirlist = ringlist(BAS);
1142
1143//----------- the trivial case of a zero ideal as input, RETURN ------------
1144   if(size(i)==0)
1145   {
1146      if(y>=1)
1147      {
1148          "// the ideal was the zero-ideal";
1149      }
1150     // execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
1151     //                 +ordstr(basering)+");");
1152         def newR7 = ring(gnirlist);
1153         setring newR7;
1154         ideal norid=ideal(0);
1155         ideal normap=fetch(BAS,ihp);
1156         export norid;
1157         export normap;
1158         result=newR7;
1159         result[size(result)+1]=list(delt,delti);
1160         setring BAS;
1161         return(result);
1162   }
1163
1164//--------------- General NOTATION, compute SB of input -----------------
1165// SM is a list, the result of mstd(i)
1166// SM[1] = SB of input ideal i,
1167// SM[2] = (minimal) generators for i.
1168// We work with SM and will copy the attributes from i to SM[2]
1169// JM will be a list, either JM[1]=maxideal(1),JM[2]=maxideal(1)
1170// in case i has onlySingularAtZero, or JM = mstd(si) where si = #[1],
1171// or JM = mstd(J) where J is the ideal of the singular locus
1172// JM[2] must be (made) radical
1173
1174   if(y>=1)
1175   {
1176     "// SB-computation of the ideal";
1177   }
1178
1179   list SM = mstd(i);              //Now the work starts
1180   int dimSM =  dim(SM[1]);        //dimension of variety to normalize
1181   if(y>=1)
1182   {
1183      "// the dimension is:";  dimSM;
1184   }
1185//----------------- the general case, set attributes ----------------
1186   //Note: onlySingularAtZero is NOT preserved under the ring extension
1187   //basering --> Hom(J,J) (in contrast to isIsolatedSingularity),
1188   //therefore we reset it:
1189
1190   attrib(i,"onlySingularAtZero",0);
1191
1192   if(attrib(i,"isPrim")==1)
1193   {
1194      attrib(SM[2],"isPrim",1);
1195   }
1196   else
1197   {
1198      attrib(SM[2],"isPrim",0);
1199   }
1200   if(attrib(i,"isIsolatedSingularity")==1)
1201   {
1202      attrib(SM[2],"isIsolatedSingularity",1);
1203   }
1204   else
1205   {
1206      attrib(SM[2],"isIsolatedSingularity",0);
1207   }
1208   if(attrib(i,"isCohenMacaulay")==1)
1209   {
1210      attrib(SM[2],"isCohenMacaulay",1);
1211   }
1212   else
1213   {
1214      attrib(SM[2],"isCohenMacaulay",0);
1215   }
1216   if(attrib(i,"isRegInCodim2")==1)
1217   {
1218      attrib(SM[2],"isRegInCodim2",1);
1219   }
1220   else
1221   {
1222      attrib(SM[2],"isRegInCodim2",0);
1223   }
1224   if(attrib(i,"isEquidimensional")==1)
1225   {
1226      attrib(SM[2],"isEquidimensional",1);
1227   }
1228   else
1229   {
1230      attrib(SM[2],"isEquidimensional",0);
1231   }
1232   if(attrib(i,"isCompleteIntersection")==1)
1233   {
1234     attrib(SM[2],"isCompleteIntersection",1);
1235   }
1236   else
1237   {
1238      attrib(SM[2],"isCompleteIntersection",0);
1239   }
1240   if(attrib(i,"isHypersurface")==1)
1241   {
1242     attrib(SM[2],"isHypersurface",1);
1243   }
1244   else
1245   {
1246      attrib(SM[2],"isHypersurface",0);
1247   }
1248
1249   if(attrib(i,"onlySingularAtZero")==1)
1250   {
1251      attrib(SM[2],"onlySingularAtZero",1);
1252   }
1253   else
1254   {
1255      attrib(SM[2],"onlySingularAtZero",0);
1256   }
1257
1258   //------- an easy and cheap test for onlySingularAtZero ---------
1259   if( (attrib(SM[2],"isIsolatedSingularity")==1) && (homog(SM[2])==1) )
1260   {
1261      attrib(SM[2],"onlySingularAtZero",1);
1262   }
1263
1264//-------------------- Trivial cases, in each case RETURN ------------------
1265// input ideal is the ideal of a partial normalization
1266
1267   // ------------ Trivial case: input ideal contains a unit ---------------
1268   if( dimSM == -1)
1269   {  "";
1270      "      // A unit ideal was found.";
1271      "      // Stop with partial result computed so far";"";
1272
1273         MB=SM[2];
1274         intvec rw;
1275         list LL=substpart(MB,ihp,0,rw);
1276         def newR6=LL[1];
1277         setring newR6;
1278         ideal norid=endid;
1279         ideal normap=endphi;
1280         kill endid,endphi;
1281         export norid;
1282         export normap;
1283         result=newR6;
1284         result[size(result)+1]=list(delt,delti);
1285         setring BAS;
1286         return(result);
1287   }
1288
1289   // --- Trivial case: input ideal is zero-dimensional and homog ---
1290   if( (dim(SM[1])==0) && (homog(SM[2])==1) )
1291   {
1292      if(y>=1)
1293      {
1294         "// the ideal was zero-dimensional and homogeneous";
1295      }
1296      MB=maxideal(1);
1297      intvec rw;
1298      list LL=substpart(MB,ihp,0,rw);
1299      def newR5=LL[1];
1300      setring newR5;
1301      ideal norid=endid;
1302      ideal normap=endphi;
1303      kill endid,endphi;
1304      export norid;
1305      export normap;
1306      result=newR5;
1307      result[size(result)+1]=list(delt,delti);
1308      setring BAS;
1309      return(result);
1310   }
1311
1312   // --- Trivial case: input ideal defines a line ---
1313   //the one-dimensional, homogeneous case and degree 1 case
1314   if( (dim(SM[1])==1) && (maxdeg1(SM[2])==1) && (homog(SM[2])==1) )
1315   {
1316      if(y>=1)
1317      {
1318         "// the ideal defines a line";
1319      }
1320      MB=SM[2];
1321      intvec rw;
1322      list LL=substpart(MB,ihp,0,rw);
1323      def newR4=LL[1];
1324      setring newR4;
1325      ideal norid=endid;
1326      ideal normap=endphi;
1327      kill endid,endphi;
1328      export norid;
1329      export normap;
1330      result=newR4;
1331      result[size(result)+1]=list(delt,delti);
1332      setring BAS;
1333      return(result);
1334   }
1335
1336//---------------------- The non-trivial cases start -------------------
1337   //the higher dimensional case
1338   //we test first hypersurface, CohenMacaulay and complete intersection
1339
1340   if( ((size(SM[2])+dim(SM[1])) == nvars(basering)) )
1341   {
1342      //the test for complete intersection
1343      attrib(SM[2],"isCohenMacaulay",1);
1344      attrib(SM[2],"isCompleteIntersection",1);
1345      attrib(SM[2],"isEquidimensional",1);
1346      if(y>=1)
1347      {
1348         "// the ideal is a complete intersection";
1349      }
1350   }
1351   if( size(SM[2]) == 1 )
1352   {
1353      attrib(SM[2],"isHypersurface",1);
1354      if(y>=1)
1355      {
1356         "// the ideal is a hypersurface";
1357      }
1358   }
1359
1360   //------------------- compute the singular locus -------------------
1361   // Computation if singular locus is critical
1362   // Notation: J ideal of singular locus or (if given) containing it
1363   // JM = mstd(J) or maxideal(1),maxideal(1)
1364   // JM[1] SB of singular locus, JM[2] minbasis, dimJ = dim(JM[1])
1365   // SM[1] SB of the input ideal i, SM[2] minbasis
1366   // Computation if singular locus is critical, because it determines the
1367   // size of the ring Hom_R(J,J). We only need a test ideal contained in J.
1368
1369   //----------------------- onlySingularAtZero -------------------------
1370   if( attrib(SM[2],"onlySingularAtZero") )
1371   {
1372       JM = maxideal(1),maxideal(1);
1373       attrib(JM[1],"isSB",1);
1374       attrib(JM[2],"isRadical",1);
1375       if( dim(SM[1]) >=2 )
1376       {
1377         attrib(SM[2],"isRegInCodim2",1);
1378       }
1379   }
1380
1381   //-------------------- not onlySingularAtZero -------------------------
1382   if( attrib(SM[2],"onlySingularAtZero") == 0 )
1383   {
1384      //--- the case where an ideal #[1] is given:
1385      if( size(#)>0 )
1386      {
1387         J = #[1],SM[2];
1388         JM = mstd(J);
1389         if( typeof(attrib(#[1],"isRadical"))!="int" )
1390         {
1391            attrib(JM[2],"isRadical",0);
1392         }
1393      }
1394
1395      //--- the case where an ideal #[1] is not given:
1396      if( (size(#)==0) )
1397      {
1398         if(y >=1 )
1399         {
1400            "// singular locus will be computed";
1401         }
1402
1403         J = SM[1],minor(jacob(SM[2]),nvars(basering)-dim(SM[1]),SM[1]);
1404         if( y >=1 )
1405         {
1406            "// SB of singular locus will be computed";
1407         }
1408         JM = mstd(J);
1409      }
1410
1411      int dimJ = dim(JM[1]);
1412      attrib(JM[1],"isSB",1);
1413      if( y>=1 )
1414      {
1415         "// the dimension of the singular locus is";  dimJ ; "";
1416      }
1417
1418      if(dim(JM[1]) <= dim(SM[1])-2)
1419      {
1420         attrib(SM[2],"isRegInCodim2",1);
1421      }
1422
1423      //------------------ the smooth case, RETURN -------------------
1424      if( dimJ == -1 )
1425      {
1426         if(y>=1)
1427         {
1428            "// the ideal is smooth";
1429         }
1430         MB=SM[2];
1431         intvec rw;
1432         list LL=substpart(MB,ihp,0,rw);
1433         def newR3=LL[1];
1434         setring newR3;
1435         ideal norid=endid;
1436         ideal normap=endphi;
1437         kill endid,endphi;
1438         export norid;
1439         export normap;
1440         result=newR3;
1441         result[size(result)+1]=list(delt,delti);
1442         setring BAS;
1443         return(result);
1444      }
1445
1446      //------- extra check for onlySingularAtZero, relatively cheap ----------
1447      //it uses the procedure 'locAtZero' from for testing
1448      //if an ideal is concentrated at 0
1449       if(y>=1)
1450       {
1451         "// extra test for onlySingularAtZero:";
1452       }
1453       if ( locAtZero(JM[1]) )
1454       {
1455           attrib(SM[2],"onlySingularAtZero",1);
1456           JM = maxideal(1),maxideal(1);
1457           attrib(JM[1],"isSB",1);
1458           attrib(JM[2],"isRadical",1);
1459       }
1460       else
1461       {
1462            attrib(SM[2],"onlySingularAtZero",0);
1463       }
1464   }
1465
1466  //displaying the attributes:
1467   if(y>=2)
1468   {
1469      "// the attributes of the ideal are:";
1470      "// isCohenMacaulay:", attrib(SM[2],"isCohenMacaulay");
1471      "// isCompleteIntersection:", attrib(SM[2],"isCompleteIntersection");
1472      "// isHypersurface:", attrib(SM[2],"isHypersurface");
1473      "// isEquidimensional:", attrib(SM[2],"isEquidimensional");
1474      "// isPrim:", attrib(SM[2],"isPrim");
1475      "// isRegInCodim2:", attrib(SM[2],"isRegInCodim2");
1476      "// isIsolatedSingularity:", attrib(SM[2],"isIsolatedSingularity");
1477      "// onlySingularAtZero:", attrib(SM[2],"onlySingularAtZero");
1478      "// isRad:", attrib(SM[2],"isRad");"";
1479   }
1480
1481   //------------- case: CohenMacaulay in codim 2, RETURN ---------------
1482   if( (attrib(SM[2],"isRegInCodim2")==1) &&
1483       (attrib(SM[2],"isCohenMacaulay")==1) )
1484   {
1485      if(y>=1)
1486      {
1487         "// the ideal was CohenMacaulay and regular in codim 2, hence normal";
1488      }
1489      MB=SM[2];
1490      intvec rw;
1491      list LL=substpart(MB,ihp,0,rw);
1492      def newR6=LL[1];
1493      setring newR6;
1494      ideal norid=endid;
1495      ideal normap=endphi;
1496      kill endid,endphi;
1497      export norid;
1498      export normap;
1499      result=newR6;
1500      result[size(result)+1]=list(delt,delti);
1501      setring BAS;
1502      return(result);
1503   }
1504
1505//---------- case: isolated singularity only at 0, RETURN ------------
1506   // In this case things are easier, we can use the maximal ideal as radical
1507   // of the singular locus;
1508   // JM mstd of ideal of singular locus, SM mstd of input ideal
1509
1510   if( attrib(SM[2],"onlySingularAtZero") )
1511   {
1512   //------ check variables for being a non zero-divizor ------
1513   // SL = ideal of vars not contained in ideal SM[1]:
1514
1515      attrib(SM[2],"isIsolatedSingularity",1);
1516      ideal SL = simplify(reduce(maxideal(1),SM[1]),2);
1517      ideal Ann = quotient(SM[2],SL[1]);
1518      ideal qAnn = simplify(reduce(Ann,SM[1]),2);
1519      //NOTE: qAnn=0 if and only if first var (=SL[1]) not in SM is a nzd of R/SM
1520
1521   //------------- We found a non-zerodivisor of R/SM -----------------------
1522   // here the enlarging of the ring via Hom_R(J,J) starts
1523
1524      if( size(qAnn)==0 )
1525      {
1526         if(y>=1)
1527         {
1528            "";
1529            "// the ideal rad(J):"; maxideal(1);
1530            "";
1531         }
1532
1533      // ------------- test for normality, compute Hom_R(J,J) -------------
1534      // Note:
1535      // HomJJ (ideal SBid, ideal id, ideal J, poly p) with
1536      //        SBid = SB of id, J = radical ideal of basering  P with:
1537      //        nonNormal(R) is in V(J), J contains the nonzero divisor p
1538      //        of R = P/id (J = test ideal)
1539      // returns a list l of three objects
1540      // l[1] : a polynomial ring, containing two ideals, 'endid' and 'endphi'
1541      //        s.t. l[1]/endid = Hom_R(J,J) and endphi= map R -> Hom_R(J,J)
1542      // l[2] : an integer which is 1 if phi is an isomorphism, 0 if not
1543      // l[3] : an integer, = dim_K(Hom_R(J,J)/R) if finite, -1 otherwise
1544
1545         list RR;
1546         RR = SM[1],SM[2],maxideal(1),SL[1];
1547         RR = HomJJ(RR,y);
1548         // --------------------- non-normal case ------------------
1549         //RR[2]==0 means that the test for normality is negative
1550         if( RR[2]==0 )
1551         {
1552            def newR=RR[1];
1553            setring newR;
1554            map psi=BAS,endphi;
1555            list JM = psi(JM); //###
1556            ideal J = JM[2];
1557            if ( delt>=0 && RR[3]>=0 )
1558            {
1559               delt = delt+RR[3];
1560            }
1561            else
1562            { delt = -1; }
1563            delti[size(delti)]=delt;
1564
1565            // ---------- recursive call of normalizationPrimes -----------
1566        //normalizationPrimes(ideal i,ideal ihp,int delt,intvec delti,list #)
1567        //ihp = (partial) normalisation map from basering
1568        //#[1] ideal s.t. V(#[1]) contains singular locus of i (test ideal)
1569
1570            if ( y>=1 )
1571            {
1572            "// case: onlySingularAtZero, non-zerodivisor found";
1573            "// contribution of delta in ringextension R -> Hom_R(J,J):"; delt;
1574            }
1575
1576            //intvec atr=getAttrib(endid);
1577            //"//### case: isolated singularity only at 0, recursive";
1578            //"size endid:", size(endid), size(string(endid));
1579            //"interred:";
1580            //endid = interred(endid);
1581            //endid = setAttrib(endid,atr);
1582            //"size endid:", size(endid), size(string(endid));
1583
1584           printlevel=printlevel+1;
1585           list tluser =
1586                normalizationPrimes(endid,psi(ihp),delt,delti);
1587           //list tluser =
1588           //     normalizationPrimes(endid,psi(ihp),delt,delti,J);
1589           //#### ??? improvement: give also the old ideal of sing locus???
1590
1591           printlevel = printlev;             //reset printlevel
1592           setring BAS;
1593           return(tluser);
1594         }
1595
1596         // ------------------ the normal case, RETURN -----------------
1597         // Now RR[2] must be 1, hence the test for normality was positive
1598         MB=SM[2];
1599         //execute("ring newR7="+charstr(basering)+",("+varstr(basering)+"),("
1600         //             +ordstr(basering)+");");
1601         def newR7 = ring(gnirlist);
1602         setring newR7;
1603         ideal norid=fetch(BAS,MB);
1604         ideal normap=fetch(BAS,ihp);
1605         if ( delt>=0 && RR[3]>=0 )
1606         {
1607               delt = delt+RR[3];
1608         }
1609         else
1610         { delt = -1; }
1611         delti[size(delti)]=delt;
1612
1613         intvec atr = getAttrib(norid);
1614
1615         //"//### case: isolated singularity only at 0, final";
1616         //"size norid:", size(norid), size(string(norid));
1617         //"interred:";
1618         //norid = interred(norid);
1619         //norid = setAttrib(norid,atr);
1620         //"size norid:", size(norid), size(string(norid));
1621
1622         export norid;
1623         export normap;
1624         result=newR7;
1625         result[size(result)+1]=list(delt,delti);
1626         setring BAS;
1627         return(result);
1628      }
1629
1630   //------ zerodivisor of R/SM was found, gives a splitting ------------
1631   //Now the case where qAnn!=0, i.e. SL[1] is a zero divisor of R/SM
1632   //and we have found a splitting: id and id1
1633   //id = Ann defines components of R/SM in the complement of V(SL[1])
1634   //id1 defines components of R/SM in the complement of V(id)
1635
1636      else
1637       {
1638          ideal id = Ann;
1639          attrib(id,"isCohenMacaulay",0);
1640          attrib(id,"isPrim",0);
1641          attrib(id,"isIsolatedSingularity",1);
1642          attrib(id,"isRegInCodim2",0);
1643          attrib(id,"isHypersurface",0);
1644          attrib(id,"isCompleteIntersection",0);
1645          attrib(id,"isEquidimensional",0);
1646          attrib(id,"onlySingularAtZero",1);
1647
1648          ideal id1 = quotient(SM[2],Ann);
1649          attrib(id1,"isCohenMacaulay",0);
1650          attrib(id1,"isPrim",0);
1651          attrib(id1,"isIsolatedSingularity",1);
1652          attrib(id1,"isRegInCodim2",0);
1653          attrib(id1,"isHypersurface",0);
1654          attrib(id1,"isCompleteIntersection",0);
1655          attrib(id1,"isEquidimensional",0);
1656          attrib(id1,"onlySingularAtZero",1);
1657
1658          // ---------- recursive call of normalizationPrimes -----------
1659          if ( y>=1 )
1660          {
1661            "// case: onlySingularAtZero, zerodivisor found, splitting:";
1662            "// total delta before splitting:", delt;
1663            "// splitting in two components:";
1664          }
1665
1666          printlevel = printlevel+1;  //to see comments in normalizationPrimes
1667          keepresult1 = normalizationPrimes(id,ihp,0,0);   //1st split factor
1668          keepresult2 = normalizationPrimes(id1,ihp,0,0);  //2nd split factor
1669          printlevel = printlev;                           //reset printlevel
1670
1671          int delt1 = keepresult1[size(keepresult1)][1];
1672          int delt2 = keepresult2[size(keepresult2)][1];
1673          intvec delti1 = keepresult1[size(keepresult1)][2];
1674          intvec delti2 = keepresult2[size(keepresult2)][2];
1675
1676          if( delt>=0 && delt1>=0 && delt2>=0 )
1677          {  ideal idid1=id,id1;
1678             int mul = vdim(std(idid1));
1679             if ( mul>=0 )
1680             {
1681               delt = delt+mul+delt1+delt2;
1682             }
1683             else
1684             {
1685               delt = -1;
1686             }
1687          }
1688         if ( y>=1 )
1689         {
1690           "// delta of first component:", delt1;
1691           "// delta of second componenet:", delt2;
1692           "// intersection multiplicity of both components:", mul;
1693           "// total delta after splitting:", delt;
1694         }
1695
1696          else
1697          {
1698            delt = -1;
1699          }
1700          for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
1701          {
1702             keepresult1=insert(keepresult1,keepresult2[lauf]);
1703          }
1704          keepresult1[size(keepresult1)]=list(delt,delti);
1705
1706          return(keepresult1);
1707       }
1708   }
1709   // Case "onlySingularAtZero" has finished and returned result
1710
1711//-------------- General case, not onlySingularAtZero, RETURN ---------------
1712   //test for non-normality, i.e. if Hom(I,I)<>R
1713   //we can use Hom(I,I) to continue
1714
1715   //------ check variables for being a non zero-divizor ------
1716   // SL = ideal of vars not contained in ideal SM[1]:
1717
1718   ideal SL = simplify(reduce(JM[2],SM[1]),2);
1719   ideal Ann = quotient(SM[2],SL[1]);
1720   ideal qAnn = simplify(reduce(Ann,SM[1]),2);
1721   //NOTE: qAnn=0 <==> first var (=SL[1]) not contained in SM is a nzd of R/SM
1722
1723   //------------- We found a non-zerodivisor of R/SM -----------------------
1724   //SM = mstd of ideal of variety, JM = mstd of ideal of singular locus
1725
1726   if( size(qAnn)==0 )
1727   {
1728      list RR;
1729      list RS;
1730      // ----------------- Computation of the radical -----------------
1731      if(y>=1)
1732      {
1733         "// radical computation of singular locus";
1734      }
1735      J = radical(JM[2]);   //the radical of singular locus
1736      JM = mstd(J);
1737
1738      if(y>=1)
1739      {
1740        "// radical is equal to:";"";  JM[2];
1741        "";
1742      }
1743      // ------------ choose non-zerodivisor of smaller degree ----------
1744      //### evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen ?
1745      if( deg(SL[1]) > deg(J[1]) )
1746      {
1747         Ann=quotient(SM[2],J[1]);
1748         qAnn=simplify(reduce(Ann,SM[1]),2);
1749         if(size(qAnn)==0)
1750         {
1751           SL[1]=J[1];
1752         }
1753      }
1754
1755      // --------------- computation of Hom(rad(J),rad(J)) --------------
1756      RR=SM[1],SM[2],JM[2],SL[1];
1757
1758     if(y>=1)
1759     {
1760        "// compute Hom(rad(J),rad(J))";
1761     }
1762
1763     RS=HomJJ(RR,y);               //most important subprocedure
1764
1765     // ------------------ the normal case, RETURN -----------------
1766     // RS[2]==1 means that the test for normality was positive
1767     if(RS[2]==1)
1768     {
1769         def lastR=RS[1];
1770         setring lastR;
1771         map psi1=BAS,endphi;
1772         ideal norid=endid;
1773         ideal normap=psi1(ihp);
1774         kill endid,endphi;
1775
1776        intvec atr=getAttrib(norid);
1777
1778        //"//### general case: not isolated singularity only at 0, final";
1779        //"size norid:", size(norid), size(string(norid));
1780        //"interred:";
1781        //norid = interred(norid);
1782        //norid = setAttrib(norid,atr);
1783        //"size norid:", size(norid), size(string(norid));
1784
1785         export norid;
1786         export normap;
1787         result=lastR;
1788         if ( y>=1 )
1789         {
1790            "// case: not onlySingularAtZero, last ring Hom_R(J,J) computed";
1791            "// delta before last ring:", delt;
1792         }
1793
1794         if ( delt>=0 && RS[3]>=0 )
1795         {
1796            delt = delt+RS[3];
1797         }
1798         else
1799         { delt = -1; }
1800
1801        // delti = delti,delt;
1802         delti[size(delti)]=delt;
1803
1804         if ( y>=1 )
1805         {
1806           "// delta of last ring:", delt;
1807         }
1808
1809         result[size(result)+1]=list(delt,delti);
1810         setring BAS;
1811         return(result);
1812     }
1813
1814    // ----- the non-normal case, recursive call of normalizationPrimes -------
1815    // RS=HomJJ(RR,y) was computed above, RS[1] contains endid and endphi
1816    // RS[1] = new ring Hom_R(J,J), RS[2]= 0 or 1, RS[2]=contribution to delta
1817    // now RS[2]must be 0, i.e. the test for normality was negative
1818
1819      int n = nvars(basering);
1820      ideal MJ = JM[2];
1821
1822      def newR=RS[1];
1823      setring newR;
1824      map psi=BAS,endphi;
1825      if ( y>=1 )
1826      {
1827        "// case: not onlySingularAtZero, compute new ring = Hom_R(J,J)";
1828        "// delta of old ring:", delt;
1829      }
1830      if ( delt>=0 && RS[3]>=0 )
1831      {
1832         delt = delt+RS[3];
1833      }
1834      else
1835      { delt = -1; }
1836      if ( y>=1 )
1837      {
1838        "// delta of new ring:", delt;
1839      }
1840
1841      delti[size(delti)]=delt;
1842      intvec atr=getAttrib(endid);
1843
1844      //"//### general case: not isolated singularity only at 0, recursive";
1845      //"size endid:", size(endid), size(string(endid));
1846      //"interred:";
1847      //endid = interred(endid);
1848      //endid = setAttrib(endid,atr);
1849      //"size endid:", size(endid), size(string(endid));
1850
1851      printlevel = printlevel+1;
1852      list tluser=
1853          normalizationPrimes(endid,psi(ihp),delt,delti,psi(MJ));
1854      printlevel = printlev;                //reset printlevel
1855      setring BAS;
1856      return(tluser);
1857   }
1858
1859   //---- A whole singular component was found, RETURN -----
1860   if( Ann == 1)
1861   {
1862      "// Input appeared not to be a radical ideal!";
1863      "// A (everywhere singular) component with ideal";
1864      "// equal to its Jacobian ideal was found";
1865      "// Procedure will stop with partial result computed so far";"";
1866
1867         MB=SM[2];
1868         intvec rw;
1869         list LL=substpart(MB,ihp,0,rw);
1870         def newR6=LL[1];
1871         setring newR6;
1872         ideal norid=endid;
1873         ideal normap=endphi;
1874         kill endid,endphi;
1875         export norid;
1876         export normap;
1877         result=newR6;
1878         result[size(result)+1]=lst(delt,delti);
1879         setring BAS;
1880         return(result);
1881   }
1882
1883   //------ zerodivisor of R/SM was found, gives a splitting ------------
1884   //Now the case where qAnn!=0, i.e. SL[1] is a zero divisor of R/SM
1885   //and we have found a splitting: new1 and new2
1886   //id = Ann defines components of R/SM in the complement of V(SL[1])
1887   //id1 defines components of R/SM in the complement of V(id)
1888
1889   else
1890   {
1891      if(y>=1)
1892      {
1893         "// zero-divisor found";
1894      }
1895      int equi = attrib(SM[2],"isEquidimensional");
1896      int oSAZ = attrib(SM[2],"onlySingularAtZero");
1897      int isIs = attrib(SM[2],"isIsolatedSingularity");
1898
1899      ideal new1 = Ann;
1900      ideal new2 = quotient(SM[2],Ann);
1901      //ideal new2=SL[1],SM[2];
1902
1903      //execute("ring newR1="+charstr(basering)+",("+varstr(basering)+"),("
1904      //                +ordstr(basering)+");");
1905      def newR1 = ring(gnirlist);
1906      setring newR1;
1907
1908      ideal vid = fetch(BAS,new1);
1909      ideal ihp = fetch(BAS,ihp);
1910      attrib(vid,"isCohenMacaulay",0);
1911      attrib(vid,"isPrim",0);
1912      attrib(vid,"isIsolatedSingularity",isIs);
1913      attrib(vid,"isRegInCodim2",0);
1914      attrib(vid,"onlySingularAtZero",oSAZ);
1915      attrib(vid,"isEquidimensional",equi);
1916      attrib(vid,"isHypersurface",0);
1917      attrib(vid,"isCompleteIntersection",0);
1918
1919      // ---------- recursive call of normalizationPrimes -----------
1920      if ( y>=1 )
1921      {
1922        "// total delta before splitting:", delt;
1923        "// splitting in two components:";
1924      }
1925      printlevel = printlevel+1;
1926      keepresult1 =
1927                  normalizationPrimes(vid,ihp,0,0);  //1st split factor
1928
1929      list delta1 = keepresult1[size(keepresult1)];
1930
1931      setring BAS;
1932      //execute("ring newR2="+charstr(basering)+",("+varstr(basering)+"),("
1933      //                +ordstr(basering)+");");
1934      def newR2 = ring(gnirlist);
1935      setring newR2;
1936
1937      ideal vid = fetch(BAS,new2);
1938      ideal ihp = fetch(BAS,ihp);
1939      attrib(vid,"isCohenMacaulay",0);
1940      attrib(vid,"isPrim",0);
1941      attrib(vid,"isIsolatedSingularity",isIs);
1942      attrib(vid,"isRegInCodim2",0);
1943      attrib(vid,"isEquidimensional",equi);
1944      attrib(vid,"isHypersurface",0);
1945      attrib(vid,"isCompleteIntersection",0);
1946      attrib(vid,"onlySingularAtZero",oSAZ);
1947
1948      keepresult2 =
1949                    normalizationPrimes(vid,ihp,0,0);
1950      list delta2 = keepresult2[size(keepresult2)];   //2nd split factor
1951      printlevel = printlev;                          //reset printlevel
1952
1953      setring BAS;
1954
1955      //compute intersection multiplicity of both components:
1956      new1 = new1,new2;
1957      int mul=vdim(std(new1));
1958
1959     // ----- normalizationPrimes finished, add up results, RETURN --------
1960      for(lauf=1;lauf<=size(keepresult2)-1;lauf++)
1961      {
1962         keepresult1 = insert(keepresult1,keepresult2[lauf]);
1963      }
1964      if ( delt >=0 && delta1[1] >=0 && delta2[1] >=0 && mul >=0 )
1965      {
1966         delt = delt+mul+delta1[1]+delta2[1];
1967      }
1968      else
1969      {  delt = -1; }
1970      delti = -2;
1971
1972      if ( y>=1 )
1973      {
1974        "// zero divisor produced a splitting into two components";
1975        "// delta of first component:", delta1;
1976        "// delta of second componenet:", delta2;
1977        "// intersection multiplicity of both components:", mul;
1978        "// total delta after splitting:", delt;
1979      }
1980      keepresult1[size(keepresult1)] = list(delt,delti);
1981      return(keepresult1);
1982   }
1983}
1984example
1985{ "EXAMPLE:";echo = 2;
1986   // Huneke
1987   ring qr=31991,(a,b,c,d,e),dp;
1988   ideal i=
1989   5abcde-a5-b5-c5-d5-e5,
1990   ab3c+bc3d+a3be+cd3e+ade3,
1991   a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
1992   abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
1993   ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
1994   a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
1995   a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
1996   b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
1997
1998   list pr=normalizationPrimes(i);
1999   def r1 = pr[1];
2000   setring r1;
2001   norid;
2002   normap;
2003}
2004
2005///////////////////////////////////////////////////////////////////////////////
2006static proc substpart(ideal endid, ideal endphi, int homo, intvec rw)
2007
2008"//Repeated application of elimpart to endid, until no variables can be
2009//directy substituded. homo=1 if input is homogeneous, rw contains
2010//original weights, endphi (partial) normalization map";
2011
2012//NOTE concerning iteration of maps: Let phi: x->f(y,z), y->g(x,z) then
2013//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
2014//and so on: none of the x or y will be eliminated
2015//Now subst: first x and then y: x+y+z->f(g(x,z),z)+g(x,z)+z eliminates y
2016//further subst replaces x by y, makes no sense (objects more compicated).
2017//Subst first y and then x eliminates x
2018//In our situation we have triangular form: x->f(y,z), y->g(z).
2019//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
2020//subst x,y: x+y+z->f(g(z),z)+g(z)+z, eliminates x,y
2021//subst y,x: x+y+z->f(y,z)+g(z)+z eliminates only x
2022//HENCE: substitute vars depending on most other vars first
2023//However, if the sytem xi-fi is reduced then xi does not appear in any of the
2024//fj and hence the order does'nt matter when substitutinp xi by fi
2025
2026{
2027   ASSUME(1, not isQuotientRing(basering) );
2028
2029   def newRing = basering;
2030   int ii,jj;
2031   map phi = newRing,maxideal(1);    //identity map
2032   list Le = elimpart(endid);
2033   //this proc and the next loop try to substitute as many variables as
2034   //possible indices of substituted variables
2035
2036   int q = size(Le[2]);    //q vars, stored in Le[2], have been substitutet
2037   intvec rw1 = 0;         //will become indices of substituted variables
2038   rw1[nvars(basering)] = 0;
2039   rw1 = rw1+1;            //rw1=1,..,1 (as many 1 as nvars(basering))
2040
2041   while( size(Le[2]) != 0 )
2042   {
2043      endid = Le[1];
2044      if ( defined(ps) )
2045      { kill ps; }
2046      map ps = newRing,Le[5];
2047      phi = ps(phi);
2048      for(ii=1;ii<=size(Le[2]);ii++)
2049      {
2050         phi=phi(phi);
2051      }
2052      //eingefuegt wegen x2-y2z2+z3
2053
2054      for( ii=1; ii<=size(rw1); ii++ )
2055      {
2056         if( Le[4][ii]==0 )        //ii = index of var which was substituted
2057         {
2058            rw1[ii]=0;             //substituted vars have entry 0 in rw1
2059         }
2060      }
2061      Le=elimpart(endid);          //repeated application of elimpart
2062      q = q + size(Le[2]);
2063   }
2064   endphi = phi(endphi);
2065//---------- return -----------------------------------------------------------
2066// first the trivial case, where all variable have been eliminated
2067   if( nvars(newRing) == q )
2068   {
2069     ring lastRing = char(basering),T(1),dp;
2070     ideal endid = T(1);
2071     ideal endphi = T(1);
2072     for(ii=2; ii<=q; ii++ )
2073     {
2074        endphi[ii] = 0;
2075     }
2076     export(endid,endphi);
2077     list L = lastRing;
2078     setring newRing;
2079     return(L);
2080   }
2081
2082// in the homogeneous case put weights for the remaining vars correctly, i.e.
2083// delete from rw those weights for which the corresponding entry of rw1 is 0
2084
2085   if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 )
2086   {
2087      jj=1;
2088      for( ii=2; ii<size(rw1); ii++)
2089      {
2090         jj++;
2091         if( rw1[ii]==0 )
2092         {
2093            rw=rw[1..jj-1],rw[jj+1..size(rw)];
2094            jj=jj-1;
2095         }
2096      }
2097      if( rw1[1]==0 ) { rw=rw[2..size(rw)]; }
2098      if( rw1[size(rw1)]==0 ){ rw=rw[1..size(rw)-1]; }
2099
2100      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),(a(rw),dp);
2101   }
2102   else
2103   {
2104      ring lastRing = char(basering),(T(1..nvars(newRing)-q)),dp;
2105   }
2106   ideal lastmap;
2107   jj = 1;
2108
2109   for(ii=1; ii<=size(rw1); ii++ )
2110   {
2111      if ( rw1[ii]==1 ) { lastmap[ii] = T(jj); jj=jj+1; }
2112      if ( rw1[ii]==0 ) { lastmap[ii] = 0; }
2113   }
2114   map phi1 = newRing,lastmap;
2115   ideal endid  = phi1(endid);      //### bottelneck
2116   ideal endphi = phi1(endphi);
2117
2118/*
2119Versuch: subst statt phi
2120   for(ii=1; ii<=size(rw1); ii++ )
2121   {
2122      if ( rw1[ii]==1 ) { endid = subst(endid,var(ii),T(jj)); }
2123      if ( rw1[ii]==0 ) { endid = subst(endid,var(ii),0); }
2124   }
2125*/
2126   export(endid);
2127   export(endphi);
2128   list L = lastRing;
2129   setring newRing;
2130   return(L);
2131}
2132///////////////////////////////////////////////////////////////////////////////
2133static proc deltaP(ideal I)
2134{
2135   ASSUME(1, not isQuotientRing(basering) );
2136   def R=basering;
2137   int c,d,i;
2138   int n=nvars(R);
2139   list nor;
2140   if(size(I)>1){ERROR("no hypersurface");}
2141   ideal J=std(slocus(I));
2142   if(dim(J)<=0){return(0);}
2143   poly h;
2144   d=1;
2145   while((d)&&(i<n))
2146   {
2147      i++;
2148      h=var(i);
2149      d=dim(std(J+ideal(h)));
2150   }
2151   i=0;
2152   while(d)
2153   {
2154      i++;
2155      if(i>10){ERROR("delta not found, please inform the authors")};
2156      h=randomLast(100)[n];
2157      d=dim(std(J+ideal(h)));
2158   }
2159   I=I,h-1;
2160   if(char(R)<=19)
2161   {
2162      nor=normalP(I);
2163   }
2164   else
2165   {
2166      nor=normal(I);
2167   }
2168   return(nor[2][2]);
2169}
2170
2171proc genus(ideal I,list #)
2172"USAGE:   genus(I) or genus(I,<option>); I a 1-dimensional ideal over a perfect field
2173RETURN:  an integer, the geometric genus p_g = p_a - delta of the projective
2174         curve defined by i, where p_a is the arithmetic genus.
2175NOTE:    genus always treats projective curves and takes projective closure if input is affine 1-dim variety.
2176         delta is the sum of all local delta-invariants of the singularities,
2177         i.e. dim(R'/R), R' the normalization of the local ring R of the
2178         singularity. @*
2179         genus(I,"nor") uses the normalization to compute delta. Usually genus(I,"nor")
2180         is slower than genus(I) but sometimes not. @*
2181         genus(I,"pri") starts with a primary decompsition.
2182EXAMPLE: example genus; shows an example
2183"
2184{
2185   ASSUME(0, not isQuotientRing(basering) );
2186
2187   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
2188
2189   int ono,rpr,ll;
2190   if(size(#)>0)
2191   {
2192     if(typeof(#[1])=="string")
2193     {
2194        if(#[1]=="nor"){ono=1;}
2195        if(#[1]=="pri"){rpr=1;}
2196     }
2197     else { ERROR("invalid option for genus");}
2198   }
2199   def R0=basering;
2200   if((char(basering)>0)||(ono))
2201   {
2202     def R1=changeord(list(list("dp",1:nvars(basering))));
2203     setring R1;
2204     ideal I=imap(R0,I);
2205     I=radical(I);
2206     I=std(I);
2207     if(dim(I)!=1)
2208     {
2209       if(((homog(I))&&(dim(I)!=2))||(!homog(I)))
2210       {
2211         ERROR("This is not a curve");
2212       }
2213     }
2214     if(homog(I)&&(dim(I)==2))
2215     {
2216       def S=R0;
2217       setring S;
2218       ideal J=I;
2219     }
2220     else
2221     {
2222       def S=changevar(varstr(R0)+",@t");
2223       setring S;
2224       ideal J=imap(R1,I);
2225       J=homog(J,@t);
2226       J=std(J);
2227     }
2228
2229     list LL=normal(J,"prim");
2230     int pa,i;
2231     for(i=1;i<=size(LL[1]);i++)
2232     {
2233        def T=LL[1][i];
2234        setring T;
2235        pa=pa-hilbPoly(std(norid))[1];
2236        setring S;
2237        kill T;
2238     }
2239     pa=pa+1;
2240     setring R0;
2241     return(pa);
2242   }
2243   if(rpr)
2244   {
2245      list LZ=minAssGTZ(I);
2246      if(size(LZ)>1)
2247      {
2248         int p_g;
2249         for(ll=1;ll<=size(LZ);ll++)
2250         {
2251            p_g=p_g+genus(LZ[ll],"pri")-1;
2252         }
2253         return(p_g+1);
2254      }
2255      else
2256      {
2257         I=LZ[1];
2258      }
2259   }
2260   else
2261   {
2262     I=radical(I);
2263   }
2264   I=std(I);
2265   if(dim(I)!=1)
2266   {
2267      if(((homog(I))&&(dim(I)!=2))||(!homog(I)))
2268      {
2269        // ERROR("This is not a curve");
2270        if(w==1){"** WARNING: Input does not define a curve **"; "";}
2271      }
2272   }
2273   list L=elimpart(I);
2274   if(size(L[2])!=0)
2275   {
2276      map psi=R0,L[5];
2277      I=std(psi(I));
2278   }
2279   if(size(I)==0)
2280   {
2281      return(0);
2282   }
2283   list N=findvars(I,0);
2284   if(size(N[1])==1)
2285   {
2286
2287      poly p=I[1];
2288     // if(deg(squarefree(p))<deg(p)){ERROR("Curve is not reduced");}
2289      return(-deg(p)+1);
2290   }
2291   if(size(N[1]) < nvars(R0))
2292   {
2293      string newvar=string(N[1]);
2294      execute("ring R=("+charstr(R0)+"),("+newvar+"),dp;");
2295      ideal I =imap(R0,I);
2296      attrib(I,"isSB",1);
2297   }
2298   else
2299   {
2300     def R=changeord(list(list("dp",1:nvars(basering))));
2301     setring R;
2302     ideal I=imap(R0,I);
2303     I=std(I);
2304   }
2305   if(dim(I)==2)
2306   {
2307      def newR=basering;
2308   }
2309   else
2310   {
2311      if(dim(I)==0)
2312      {
2313         execute("ring Rhelp=("+charstr(R0)+"),(@s,@t),dp;");
2314      }
2315      else
2316      {
2317         execute("ring Rhelp=("+charstr(R0)+"),(@s),dp;");
2318      }
2319      def newR=R+Rhelp;
2320      setring newR;
2321      ideal I=imap(R,I);
2322      I=homog(I,@s);
2323      attrib(I,"isSB",1);
2324   }
2325
2326   if((nvars(basering)<=3)&&(size(I)>1))
2327   {
2328       ERROR("This is not equidimensional");
2329   }
2330
2331   intvec hp=hilbPoly(I);
2332   int p_a=1-hp[1];
2333   int d=hp[2];
2334
2335   if(w>=1)
2336   {
2337      "";"The ideal of the projective curve:";"";I;"";
2338      "The coefficients of the Hilbert polynomial";hp;
2339      "arithmetic genus:";p_a;
2340      "degree:";d;"";
2341   }
2342
2343   intvec v = hilb(I,1);
2344   int i,o;
2345   if(nvars(basering)>3)
2346   {
2347      map phi=newR,maxideal(1);
2348      int de;
2349      ideal K,L1;
2350      matrix M;
2351      poly m=var(4);
2352      poly he;
2353      for(i=5;i<=nvars(basering);i++){m=m*var(i);}
2354      K=eliminate(I,m,v);
2355      if(size(K)==1){de=deg(K[1]);}
2356      m=var(1);
2357      for(i=2;i<=nvars(basering)-3;i++){m=m*var(i);}
2358      i=0;
2359      while(d!=de)
2360      {
2361         o=1;
2362         i++;
2363         K=phi(I);
2364         K=eliminate(K,m,v);
2365         if(size(K)==1){de=deg(K[1]);}
2366         if((i==5)&&(d!=de))
2367         {
2368            K=reduce(equidimMax(I),I);
2369            if(size(K)!=0){ERROR("This is not equidimensional");}
2370         }
2371         if(i==10)
2372         {
2373            J;K;
2374            ERROR("genus: did not find a good projection for to
2375                           the plain");
2376         }
2377         if(i<5)
2378         {
2379            M=sparsetriag(nvars(newR),nvars(newR),80-5*i,i);
2380         }
2381         else
2382         {
2383            if(i<8)
2384            {
2385               M=transpose(sparsetriag(nvars(newR),nvars(newR),80-5*i,i));
2386            }
2387            else
2388            {
2389               he=0;
2390               while(he==0)
2391               {
2392                  M=randommat(nvars(newR),nvars(newR),ideal(1),20);
2393                  he=det(M);
2394               }
2395            }
2396         }
2397         L1=M*transpose(maxideal(1));
2398         phi=newR,L1;
2399      }
2400      I=K;
2401   }
2402   poly p=I[1];
2403
2404   execute("ring S=("+charstr(R)+"),(x,y,t),dp;");
2405   ideal L=maxideal(1);
2406   execute("ring C=("+charstr(R)+"),(x,y),ds;");
2407   ideal I;
2408   execute("ring A=("+charstr(R)+"),(x,t),dp;");
2409   map phi=S,1,x,t;
2410   map psi=S,x,1,t;
2411   poly g,h;
2412   ideal I,I1;
2413   execute("ring B=("+charstr(R)+"),(x,t),ds;");
2414
2415   setring S;
2416   if(o)
2417   {
2418     for(i=1;i<=nvars(newR)-3;i++){L[i]=0;}
2419     L=L,maxideal(1);
2420   }
2421   map sigma=newR,L;
2422   poly F=sigma(p);
2423   if(w>=1){"the projected curve:";"";F;"";}
2424
2425   kill newR;
2426
2427   int genus=(d-1)*(d-2) div 2;
2428   if(w>=1){"the arithmetic genus of the plane curve:";genus;pause();}
2429
2430   int delt,deltaloc,deltainf,tau,tauinf,cusps,iloc,iglob,l,nsing,
2431       tauloc,tausing,k,rat,nbranchinf,nbranch,nodes,cuspsinf,nodesinf;
2432   list inv;
2433
2434   if(w>=1)
2435     {"";"analyse the singularities at oo";"";"singular locus at (1,x,0):";"";}
2436   setring A;
2437   g=phi(F);
2438   h=psi(F);
2439   I=g,jacob(g),var(2);
2440   I=std(I);
2441   if(deg(I[1])>0)
2442   {
2443      list qr=minAssGTZ(I);
2444      if(w>=1){qr;"";}
2445
2446      for(k=1;k<=size(qr);k++)
2447      {
2448         if(w>=1){ nsing=nsing+vdim(std(qr[k]));}
2449         inv=deltaLoc(g,qr[k]);
2450         deltainf=deltainf+inv[1];
2451         tauinf=tauinf+inv[2];
2452         l=vdim(std(qr[k]));
2453         if(inv[2]==l){nodesinf=nodesinf+l;}
2454         if(inv[2]==2*l){cuspsinf=cuspsinf+l;}
2455         nbranchinf=nbranchinf+inv[3];
2456      }
2457   }
2458   else
2459   {
2460     if(w>=1){"            the curve is smooth at (1,x,0)";"";}
2461   }
2462   if(w>=1){"singular locus at (0,1,0):";"";}
2463   inv=deltaLoc(h,maxideal(1));
2464   if((w>=1)&&(inv[2]!=0)){ nsing++;}
2465   deltainf=deltainf+inv[1];
2466   tauinf=tauinf+inv[2];
2467   if(inv[2]==1){nodesinf++;}
2468   if(inv[2]==2){cuspsinf++;}
2469
2470   if((w>=1)&&(inv[2]==0)){"            the curve is smooth at (0,1,0)";"";}
2471   if(inv[2]>0){nbranchinf=nbranchinf+inv[3];}
2472
2473   if(w>=1)
2474   {
2475      if(tauinf==0)
2476      {
2477        "            the curve is smooth at oo";"";
2478      }
2479      else
2480      {
2481         "number of singularities at oo:";nsing;
2482         "nodes at oo:";nodesinf;
2483         "cusps at oo:";cuspsinf;
2484         "branches at oo:";nbranchinf;
2485         "Tjurina number at oo:";tauinf;
2486         "delta at oo:";deltainf;
2487         "Milnor number at oo:";2*deltainf-nbranchinf+nsing;
2488         pause();
2489      }
2490      "singularities at (x,y,1):";"";
2491   }
2492   execute("ring newR=("+charstr(R)+"),(x,y),dp;");
2493   //the singularities at the affine part
2494   map sigma=S,var(1),var(2),1;
2495   ideal I=sigma(F);
2496
2497   ideal I1=jacob(I);
2498   matrix Hess[2][2]=jacob(I1);
2499   ideal ID=I+I1+ideal(det(Hess));//singular locus of I+I1
2500
2501   ideal radID=std(radical(ID));//the non-nodal locus
2502   if(w>=1){"the non-nodal locus:";"";radID;pause();"";}
2503   if(deg(radID[1])==0)
2504   {
2505     ideal IDsing=1;
2506   }
2507   else
2508   {
2509     ideal IDsing=minor(jacob(ID),2)+radID;//singular locus of ID
2510   }
2511
2512   iglob=vdim(std(IDsing));
2513
2514   if(iglob!=0)//computation of the radical of IDsing
2515   {
2516      ideal radIDsing=reduce(IDsing,radID);
2517      if(size(radIDsing)==0)
2518      {
2519         radIDsing=radID;
2520         attrib(radIDsing,"isSB",1);
2521      }
2522      else
2523      {
2524         radIDsing=std(radical(IDsing));
2525      }
2526      iglob=vdim(radIDsing);
2527      if((w>=1)&&(iglob))
2528          {"the non-nodal-cuspidal locus:";radIDsing;pause();"";}
2529   }
2530   cusps=vdim(radID)-iglob;
2531   nsing=nsing+cusps;
2532
2533   if(iglob==0)
2534   {
2535      if(w>=1){"             there are only cusps and nodes";"";}
2536      tau=vdim(std(I+jacob(I)));
2537      tauinf=tauinf+tau;
2538      nodes=tau-2*cusps;
2539      delt=nodes+cusps;
2540      nbranch=2*tau-3*cusps;
2541      nsing=nsing+nodes;
2542   }
2543   else
2544   {
2545       if(w>=1){"the non-nodal-cuspidal singularities";"";}
2546       setring C;
2547       ideal I1=imap(newR,radIDsing);
2548       iloc=vdim(std(I1));
2549       if(iglob==iloc)
2550       {
2551          if(w>=1){"only cusps and nodes outside (0,0,1)";}
2552          setring newR;
2553          tau=vdim(std(I+jacob(I)));
2554          tauinf=tauinf+tau;
2555          inv=deltaLoc(I[1],maxideal(1));
2556          delt=inv[1];
2557          tauloc=inv[2];
2558          nodes=tau-tauloc-2*cusps;
2559          nsing=nsing+nodes;
2560          if (inv[2]!=0) { nsing++; }
2561          nbranch=inv[3]+ 2*nodes+cusps;
2562          delt=delt+nodes+cusps;
2563          if((w>=1)&&(inv[2]==0)){"smooth at (0,0,1)";}
2564        }
2565        else
2566        {
2567           setring newR;
2568           list pr=minAssGTZ(radIDsing);
2569           if(w>=1){pr;}
2570
2571           for(k=1;k<=size(pr);k++)
2572           {
2573              if(w>=1){nsing=nsing+vdim(std(pr[k]));}
2574              inv=deltaLoc(I[1],pr[k]);
2575              delt=delt+inv[1];
2576              tausing=tausing+inv[2];
2577              nbranch=nbranch+inv[3];
2578           }
2579           tau=vdim(std(I+jacob(I)));
2580           tauinf=tauinf+tau;
2581           nodes=tau-tausing-2*cusps;
2582           nsing=nsing+nodes;
2583           delt=delt+nodes+cusps;
2584           nbranch=nbranch+2*nodes+cusps;
2585        }
2586   }
2587   genus=genus-delt-deltainf;
2588   if(w>=1)
2589   {
2590      "The projected plane curve has locally:";"";
2591      "singularities:";nsing;
2592      "branches:";nbranch+nbranchinf;
2593      "nodes:"; nodes+nodesinf;
2594      "cusps:";cusps+cuspsinf;
2595      "Tjurina number:";tauinf;
2596      "Milnor number:";2*(delt+deltainf)-nbranch-nbranchinf+nsing;
2597      "delta of the projected curve:";delt+deltainf;
2598      "delta of the curve:";p_a-genus;
2599      "genus:";genus;
2600      "====================================================";
2601      "";
2602   }
2603   setring R0;
2604   return(genus);
2605}
2606example
2607{ "EXAMPLE:"; echo = 2;
2608   ring r=0,(x,y),dp;
2609   ideal i=y^9 - x^2*(x - 1)^9;
2610   genus(i);
2611   ring r7=7,(x,y),dp;
2612   ideal i=y^9 - x^2*(x - 1)^9;
2613   genus(i);
2614}
2615
2616///////////////////////////////////////////////////////////////////////////////
2617proc deltaLoc(poly f,ideal singL)
2618"USAGE:  deltaLoc(f,J);  f poly, J ideal
2619ASSUME: f is reduced bivariate polynomial; basering has exactly two variables;
2620        J is irreducible prime component of the singular locus of f (e.g., one
2621        entry of the output of @code{minAssGTZ(I);}, I = <f,jacob(f)>).
2622RETURN:  list L:
2623@texinfo
2624@table @asis
2625@item @code{L[1]}; int:
2626         the sum of (local) delta invariants of f at the (conjugated) singular
2627         points given by J.
2628@item @code{L[2]}; int:
2629         the sum of (local) Tjurina numbers of f at the (conjugated) singular
2630         points given by J.
2631@item @code{L[3]}; int:
2632         the sum of (local) number of branches of f at the (conjugated)
2633         singular points given by J.
2634@end table
2635@end texinfo
2636NOTE:    procedure makes use of @code{execute}; increasing printlevel displays
2637         more comments (default: printlevel=0).
2638SEE ALSO: delta, tjurina
2639KEYWORDS: delta invariant; Tjurina number
2640EXAMPLE: example deltaLoc;  shows an example
2641"
2642{
2643   ASSUME(0, not isQuotientRing(basering) );
2644
2645   intvec save_opt=option(get);
2646   option(redSB);
2647   def R=basering;
2648   execute("ring S=("+charstr(R)+"),(x,y),lp;");
2649   map phi=R,x,y;
2650   ideal singL=phi(singL);
2651   singL=simplify(std(singL),1);
2652   attrib(singL,"isSB",1);
2653   int d=vdim(singL);
2654   poly f=phi(f);
2655   int i;
2656   int w = printlevel-voice+2;  // w=printlevel (default: w=0)
2657   if(d==1)
2658   {
2659      map alpha=S,var(1)-singL[2][2],var(2)-singL[1][2];
2660      f=alpha(f);
2661      execute("ring C=("+charstr(S)+"),("+varstr(S)+"),ds;");
2662      poly f=imap(S,f);
2663      ideal singL=imap(S,singL);
2664      if((w>=1)&&(ord(f)>=2))
2665      {
2666        "local analysis of the singularities";"";
2667        basering;
2668        singL;
2669        f;
2670        pause();
2671      }
2672   }
2673   else
2674   {
2675      poly p;
2676      poly c;
2677      map psi;
2678      number co;
2679
2680      while((deg(lead(singL[1]))>1)&&(deg(lead(singL[2]))>1))
2681      {
2682         psi=S,x,y+random(-100,100)*x;
2683         singL=psi(singL);
2684         singL=std(singL);
2685          f=psi(f);
2686      }
2687
2688      if(deg(lead(singL[2]))==1)
2689      {
2690         p=singL[1];
2691         c=singL[2]-lead(singL[2]);
2692         co=leadcoef(singL[2]);
2693      }
2694      if(deg(lead(singL[1]))==1)
2695      {
2696         psi=S,y,x;
2697         f=psi(f);
2698         singL=psi(singL);
2699         p=singL[2];
2700         c=singL[1]-lead(singL[1]);
2701         co=leadcoef(singL[1]);
2702      }
2703
2704      execute("ring B=("+charstr(S)+"),a,dp;");
2705      map beta=S,a,a;
2706      poly p=beta(p);
2707
2708      execute("ring C=("+charstr(S)+",a),("+varstr(S)+"),ds;");
2709      number p=number(imap(B,p));
2710
2711      minpoly=p;
2712      map iota=S,a,a;
2713      number c=number(iota(c));
2714      number co=iota(co);
2715
2716      map alpha=S,x-c/co,y+a;
2717      poly f=alpha(f);
2718      f=cleardenom(f);
2719      if((w>=1)&&(ord(f)>=2))
2720      {
2721        "local analysis of the singularities";"";
2722        basering;
2723        alpha;
2724        f;
2725        pause();
2726        "";
2727      }
2728   }
2729   option(noredSB);
2730   ideal fstd=std(ideal(f)+jacob(f));
2731   poly hc=highcorner(fstd);
2732   int tau=vdim(fstd);
2733   int o=ord(f);
2734   int delt,nb;
2735
2736   if(tau==0)                 //smooth case
2737   {
2738      setring R;
2739      option(set,save_opt);
2740      return(list(0,0,1));
2741   }
2742   if((char(basering)>=181)||(char(basering)==0))
2743   {
2744      if(o==2)                //A_k-singularity
2745      {
2746        if(w>=1){"A_k-singularity";"";}
2747         setring R;
2748         delt=(tau+1) div 2;
2749         option(set,save_opt);
2750         return(list(d*delt,d*tau,d*(2*delt-tau+1)));
2751      }
2752      if((lead(f)==var(1)*var(2)^2)||(lead(f)==var(1)^2*var(2)))
2753      {
2754        if(w>=1){"D_k- singularity";"";}
2755
2756         setring R;
2757         delt=(tau+2) div 2;
2758         option(set,save_opt);
2759         return(list(d*delt,d*tau,d*(2*delt-tau+1)));
2760      }
2761
2762      int mu=vdim(std(jacob(f)));
2763
2764      poly g=f+var(1)^mu+var(2)^mu;  //to obtain a convenient Newton-polygon
2765
2766      list NP=newtonpoly(g);
2767      if(w>=1){"Newton-Polygon:";NP;"";}
2768      int s=size(NP);
2769
2770      if(is_NND(f,mu,NP))
2771      { // the Newton-polygon is non-degenerate
2772        // compute nb, the number of branches
2773        for(i=1;i<=s-1;i++)
2774        {
2775          nb=nb+gcd(NP[i][2]-NP[i+1][2],NP[i][1]-NP[i+1][1]);
2776        }
2777        if(w>=1){"Newton-Polygon is non-degenerated";"";}
2778        setring R;
2779        option(set,save_opt);
2780        return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
2781      }
2782
2783      if(w>=1){"Newton-Polygon is degenerated";"";}
2784/* need to re-consider the degree bound (de):
2785      // the following can certainly be made more efficient when replacing
2786      // 'hnexpansion' (used only for computing number of branches) by
2787      // successive blowing-up + test if Newton polygon degenerate:
2788      if(s>2)    //  splitting of f
2789      {
2790         if(w>=1){"Newton polygon can be used for splitting";"";}
2791         intvec v=NP[1][2]-NP[2][2],NP[2][1];
2792         int de=w_deg(g,v);
2793         //int st=w_deg(hc,v)+v[1]+v[2];
2794         int st=w_deg(var(1)^NP[size(NP)][1],v)+1;
2795         poly f1=var(2)^NP[2][2];
2796         poly f2=jet(g,de,v)/var(2)^NP[2][2];
2797         poly h=g-f1*f2;
2798         de=w_deg(h,v);
2799         poly k;
2800         ideal wi=var(2)^NP[2][2],f2;
2801         matrix li;
2802         while(de<st)
2803         {
2804           k=jet(h,de,v);
2805           li=lift(wi,k);
2806           f1=f1+li[2,1];
2807           f2=f2+li[1,1];
2808           h=g-f1*f2;
2809           de=w_deg(h,v);
2810         }
2811         nb=deltaLoc(f1,maxideal(1))[3]+deltaLoc(f2,maxideal(1))[3];
2812
2813         setring R;
2814         option(set,save_opt);
2815         return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
2816      }
2817*/
2818      f=jet(f,deg(hc)+2);
2819      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
2820      ideal fac=factorize(f,1);
2821      if(size(fac)>1)
2822      {
2823         nb=0;
2824         for(i=1;i<=size(fac);i++)
2825         {
2826            nb=nb+deltaLoc(fac[i],maxideal(1))[3];
2827         }
2828         setring R;
2829         option(set,save_opt);
2830         return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
2831      }
2832      list HNEXP=hnexpansion(f);
2833      if (typeof(HNEXP[1])=="ring")
2834      {
2835        def altring = basering;
2836        def HNEring = HNEXP[1]; setring HNEring;
2837        nb=size(hne);
2838        setring R;
2839        kill HNEring;
2840      }
2841      else
2842      {
2843        nb=size(HNEXP);
2844      }
2845      setring R;
2846      option(set,save_opt);
2847      return(list(d*(mu+nb-1) div 2,d*tau,d*nb));
2848   }
2849   else             //the case of small characteristic
2850   {
2851      f=jet(f,deg(hc)+2);
2852      if(w>=1){"now we have to use Hamburger-Noether (Puiseux) expansion";}
2853      delt=delta(f);
2854      setring R;
2855      option(set,save_opt);
2856      return(list(d*delt,d*tau,d));
2857   }
2858   option(set,save_opt);
2859}
2860example
2861{ "EXAMPLE:"; echo = 2;
2862  ring r=0,(x,y),dp;
2863  poly f=(x2+y^2-1)^3 +27x2y2;
2864  ideal I=f,jacob(f);
2865  I=std(I);
2866  list qr=minAssGTZ(I);
2867  size(qr);
2868  // each component of the singular locus either describes a cusp or a pair
2869  // of conjugated nodes:
2870  deltaLoc(f,qr[1]);
2871  deltaLoc(f,qr[2]);
2872  deltaLoc(f,qr[3]);
2873  deltaLoc(f,qr[4]);
2874  deltaLoc(f,qr[5]);
2875  deltaLoc(f,qr[6]);
2876}
2877///////////////////////////////////////////////////////////////////////////////
2878// compute the weighted degree of p;
2879// this code is an exact copy of the proc in paraplanecurves.lib
2880// (since we do not want to make it non-static)
2881static proc w_deg(poly p, intvec v)
2882{
2883   if(p==0){return(-1);}
2884   int d=0;
2885   while(jet(p,d,v)==0){d++;}
2886   d=(transpose(leadexp(jet(p,d,v)))*v)[1];
2887   return(d);
2888}
2889
2890//proc hilbPoly(ideal J)
2891//{
2892//   poly hp;
2893//   int i;
2894//   if(!attrib(J,"isSB")){J=std(J);}
2895//   intvec v = hilb(J,2);
2896//   for(i=1; i<=size(v); i++){ hp=hp+v[i]*(var(1)-i+2);}
2897//   return(hp);
2898//}
2899
2900
2901//////////////////////////////////////////////////////////////////////////////
2902
2903proc primeClosure (list L, list #)
2904"USAGE:    primeClosure(L [,c]); L a list of a ring containing a prime ideal
2905          ker, c an optional integer
2906RETURN:   a list L (of size n+1) consisting of rings L[1],...,L[n] such that
2907          - L[1] is a copy of (not a reference to!) the input ring L[1]
2908          - all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi
2909            such that
2910                    L[1]/ker --> ... --> L[n]/ker
2911            are injections given by the corresponding ideals phi, and L[n]/ker
2912            is the integral closure of L[1]/ker in its quotient field.
2913          - all rings L[i] contain a polynomial nzd such that elements of
2914            L[i]/ker are quotients of elements of L[i-1]/ker with denominator
2915            nzd via the injection phi.
2916            L[n+1] is the delta invariant
2917NOTE:     - L is constructed by recursive calls of primeClosure itself.
2918          - c determines the choice of nzd:
2919               - c not given or equal to 0: first generator of the ideal SL,
2920                 the singular locus of Spec(L[i]/ker)
2921               - c<>0: the generator of SL with least number of monomials.
2922EXAMPLE:  example primeClosure; shows an example
2923"
2924{
2925  //---- Start with a consistency check:
2926
2927  if (!(typeof(L[1])=="ring"))
2928  {
2929      "// Parameter must be a ring or a list containing a ring!";
2930      return(-1);
2931  }
2932
2933  int dblvl = printlevel-voice+2;
2934  list gnirlist = ringlist(basering);
2935
2936  //---- Some auxiliary variables:
2937  int delt;                      //finally the delta invariant
2938  if ( size(L) == 1 )
2939  {
2940      L[2] = delt;              //set delta to 0
2941  }
2942  int n = size(L)-1;            //L without delta invariant
2943
2944  //---- How to choose the non-zerodivisor later on?
2945
2946  int nzdoption=0;
2947  if (size(#)>0)
2948  {
2949      nzdoption=#[1];
2950  }
2951
2952// R0 below is the ring to work with, if we are in step one, make a copy of the
2953// input ring, so that all objects are created in the copy, not in the original
2954// ring (therefore a copy, not a reference is defined).
2955
2956  if (n==1)
2957  {
2958      def R = L[1];
2959      list Rlist = ringlist(R);
2960      def BAS = basering;
2961      setring R;
2962      if (!(typeof(ker)=="ideal"))
2963      {
2964          "// No ideal ker in the input ring!";
2965          return (-1);
2966      }
2967      ker=simplify(interred(ker),15);
2968      //execute ("ring R0="+charstr(R)+",("+varstr(R)+"),("+ordstr(R)+");");
2969      // Rlist may be not defined in this new ring, so we define it again.
2970      list Rlist2 = ringlist(R);
2971      def R0 = ring(Rlist2);
2972      setring R0;
2973      ideal ker=fetch(R,ker);
2974      // check whether we compute the normalization of the blow up of
2975      // an isolated singularity at the origin (checked in normalI)
2976
2977      if (typeof(attrib(L[1],"iso_sing_Rees"))=="int")
2978      {
2979        attrib(R0,"iso_sing_Rees",attrib(L[1],"iso_sing_Rees"));
2980      }
2981      L[1]=R0;
2982  }
2983  else
2984  {
2985      def R0 = L[n];
2986      setring R0;
2987  }
2988
2989// In order to apply HomJJ from normal.lib, we need the radical of the singular
2990// locus of ker, J:=rad(ker):
2991
2992   list SM=mstd(ker);
2993
2994// In the first iteration, we have to compute the singular locus "from
2995// scratch".
2996// In further iterations, we can fetch it from the previous one but
2997// have to compute its radical
2998// the next rings R1 contain already the (fetched) ideal
2999
3000  if (n==1)                              //we are in R0=L[1]
3001  {
3002      if (typeof(attrib(R0,"iso_sing_Rees"))=="int")
3003      {
3004        ideal J;
3005        for (int s=1;s<=attrib(R0,"iso_sing_Rees");s++)
3006        {
3007          J=J,var(s);
3008        }
3009        J = J,SM[2];
3010        list JM = mstd(J);
3011      }
3012      else
3013      {
3014        if ( dblvl >= 1 )
3015        {"";
3016           "// compute the singular locus";
3017        }
3018        //### Berechnung des singulaeren Orts geaendert (ist so schneller)
3019        ideal J = minor(jacob(SM[2]),nvars(basering)-dim(SM[1]),SM[1]);
3020        J = J,SM[2];
3021        list JM = mstd(J);
3022      }
3023
3024      if ( dblvl >= 1 )
3025      {"";
3026         "// dimension of singular locus is", dim(JM[1]);
3027         if (  dblvl >= 2 )
3028         {"";
3029            "// the singular locus is:"; JM[2];
3030         }
3031      }
3032
3033      if ( dblvl >= 1 )
3034      {"";
3035         "// compute radical of singular locus";
3036      }
3037
3038      J = simplify(radical(JM[2]),2);
3039      if ( dblvl >= 1 )
3040      {"";
3041         "// radical of singular locus is:"; J;
3042         pause();
3043      }
3044  }
3045  else
3046  {
3047      if ( dblvl >= 1 )
3048      {"";
3049         "// compute radical of test ideal in ideal of singular locus";
3050      }
3051      J = simplify(radical(J),2);
3052      if ( dblvl >= 1 )
3053      {"";
3054         "// radical of test ideal is:"; J;
3055         pause();
3056      }
3057  }
3058
3059  // having computed the radical J of/in the ideal of the singular locus,
3060  // we now need to pick an element nzd of J;
3061  // NOTE: nzd must be a non-zero divisor mod ker, i.e. not contained in ker
3062
3063  poly nzd = J[1];
3064  poly nzd1 = NF(nzd,SM[1]);
3065  if (nzd1 != 0)
3066  {
3067     if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) )
3068     {
3069        nzd = nzd1;
3070     }
3071  }
3072
3073  if (nzdoption || nzd1==0)
3074  {
3075    for (int ii=2;ii<=ncols(J);ii++)
3076    {
3077      nzd1 = NF(J[ii],SM[1]);
3078      if ( nzd1 != 0 )
3079      {
3080        if ( (deg(nzd)>=deg(J[ii])) && (size(nzd)>size(J[ii])) )
3081        {
3082          nzd=J[ii];
3083        }
3084        if ( deg(nzd)>=deg(nzd1) && size(nzd)>size(nzd1) )
3085        {
3086          nzd = nzd1;
3087        }
3088      }
3089    }
3090  }
3091
3092  export nzd;
3093  // In this case we do not eliminate variables, so that the maps
3094  // are well defined.
3095  list RR = SM[1],SM[2],J,nzd,1;
3096
3097  if ( dblvl >= 1 )
3098  {"";
3099     "// compute the first ring extension:";
3100     "RR: ";
3101     RR;
3102  }
3103
3104  list RS = HomJJ(RR);
3105  //NOTE: HomJJ creates new ring with variables X(i) and T(j)
3106//-------------------------------------------------------------------------
3107// If we've reached the integral closure (as determined by the result of
3108// HomJJ), then we are done, otherwise we have to prepare the next iteration.
3109
3110  if (RS[2]==1)     // we've reached the integral closure, we are still in R0
3111    {
3112      kill J;
3113      if ( n== 1)
3114      {
3115        def R1 = RS[1];
3116        setring R1;
3117        ideal norid, normap = endid, endphi;
3118        kill endid,  endphi;
3119
3120        //"//### case: primeClosure, final";
3121        //"size norid:", size(norid), size(string(norid));
3122        //"interred:";
3123        //norid = interred(norid);
3124        //"size norid:", size(norid), size(string(norid));
3125
3126        export (norid, normap);
3127        L[1] = R1;
3128      }
3129      return(L);
3130    }
3131  else                        // prepare the next iteration
3132    {
3133      if (n==1)               // In the first iteration: keep only the data
3134      {                       // needed later on.
3135         kill RR,SM;
3136         export(ker);
3137      }
3138      if ( dblvl >= 1 )
3139      {"";
3140         "// computing the next ring extension, we are in loop"; n+1;
3141      }
3142
3143      def R1 = RS[1];         // The data of the next ring R1:
3144      delt = RS[3];           // the delta invariant of the ring extension
3145      setring R1;             // keep only what is necessary and kill
3146      ideal ker=endid;        // everything else.
3147      export(ker);
3148      ideal norid=endid;
3149
3150      //"//### case: primeClosure, loop", n+1;
3151      //"size norid:", size(norid), size(string(norid));
3152      //"interred:";
3153      //norid = interred(norid);        //????
3154      //"size norid:", size(norid), size(string(norid));
3155
3156      export(norid);
3157      kill endid;
3158
3159      map phi = R0,endphi;                        // fetch the singular locus
3160      ideal J = mstd(simplify(phi(J)+ker,4))[2];  // ideal J in R1
3161      export(J);
3162      if(n>1)
3163      {
3164         ideal normap=phi(normap);
3165      }
3166      else
3167      {
3168         ideal normap=endphi;
3169      }
3170      export(normap);
3171      kill phi;              // we save phi as ideal, not as map, so that
3172      ideal phi=endphi;      // we have more flexibility in the ring names
3173      kill endphi;           // later on.
3174      export(phi);
3175      L=insert(L,R1,n);       // Add the new ring R1 and go on with the
3176                              // next iteration
3177      if ( L[size(L)] >= 0 && delt >= 0 )
3178      {
3179         delt = L[size(L)] + delt;
3180      }
3181      else
3182      {
3183         delt = -1;
3184      }
3185      L[size(L)] = delt;
3186
3187      if (size(#)>0)
3188      {
3189          return (primeClosure(L,#));
3190      }
3191      else
3192      {
3193          return(primeClosure(L));         // next iteration.
3194      }
3195    }
3196}
3197example
3198{
3199  "EXAMPLE:"; echo=2;
3200  ring R=0,(x,y),dp;
3201  ideal I=x4,y4;
3202  def K=ReesAlgebra(I)[1];        // K contains ker such that K/ker=R[It]
3203  list L=primeClosure(K);
3204  def R(1)=L[1];                  // L[4] contains ker, L[4]/ker is the
3205  def R(4)=L[4];                  // integral closure of L[1]/ker
3206  setring R(1);
3207  R(1);
3208  ker;
3209  setring R(4);
3210  R(4);
3211  ker;
3212}
3213
3214///////////////////////////////////////////////////////////////////////////////
3215
3216proc closureFrac(list L)
3217"USAGE:    closureFrac (L); L a list of size n+1 as in the result of
3218          primeClosure, L[n] contains an additional polynomial f
3219CREATE:   a list fraction of two elements of L[1], such that
3220          f=fraction[1]/fraction[2] via the injections phi L[i]-->L[i+1].
3221EXAMPLE:  example closureFrac; shows an example
3222"
3223{
3224// Define some auxiliary variables:
3225
3226  int n=size(L)-1;
3227  int i,j,k,l,n2,n3;
3228  intvec V;
3229  string mapstr;
3230  for (i=1; i<=n; i++)
3231  {
3232    ASSUME(0, not isQuotientRing( L[i] ) );
3233    def R(i) = L[i];
3234  }
3235
3236// The quotient representing f is computed as in 'closureGenerators' with
3237// the differences that
3238//   - the loop is done twice: for the numerator and for the denominator;
3239//   - the result is stored in the list fraction and
3240//   - we have to make sure that no more objects of the rings R(i) survive.
3241
3242  for (j=1; j<=2; j++)
3243    {
3244      setring R(n);
3245      if (j==1)
3246      {
3247         poly p=f;
3248      }
3249      else
3250      {
3251         p=1;
3252      }
3253
3254      for (k=n; k>1; k--)
3255      {
3256          if (j==1)
3257          {
3258             map phimap=R(k-1),phi;
3259          }
3260
3261          p=p*phimap(nzd);
3262
3263          if (j==2)
3264          {
3265            kill phimap;
3266          }
3267
3268          if (j==1)
3269          {
3270             //### noch abfragen ob Z(i) definiert ist
3271             list gnirlist = ringlist(R(k));
3272             n2 = size(gnirlist[2]);
3273             n3 = size(gnirlist[3]);
3274             for( i=1; i<=ncols(phi); i++)
3275             {
3276               gnirlist[2][n2+i] = "Z("+string(i)+")";
3277             }
3278             V=0;
3279             V[ncols(phi)]=0; V=V+1;
3280             gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1);
3281             def S(k) = ring(gnirlist);
3282             setring S(k);
3283
3284             //execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
3285             //          Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
3286             //          +"),dp("+string(ncols(phi))+"));");
3287
3288              ideal phi = imap(R(k),phi);
3289              ideal J = imap (R(k),ker);
3290              for (l=1;l<=ncols(phi);l++)
3291              {
3292                  J=J+(Z(l)-phi[l]);
3293              }
3294              J=groebner(J);
3295              poly h=NF(imap(R(k),p),J);
3296          }
3297          else
3298          {
3299              setring S(k);
3300              h=NF(imap(R(k),p),J);
3301              setring R(k);
3302              kill p;
3303          }
3304
3305          setring R(k-1);
3306
3307          if (j==1)
3308          {
3309              ideal maxi;
3310              maxi[nvars(R(k))] = 0;
3311              maxi = maxi,maxideal(1);
3312              map backmap = S(k),maxi;
3313
3314              //mapstr=" map backmap = S(k),";
3315              //for (l=1;l<=nvars(R(k));l++)
3316              //{
3317              //  mapstr=mapstr+"0,";
3318              //}
3319              //execute (mapstr+"maxideal(1);");
3320              poly p;
3321          }
3322          p=NF(backmap(h),std(ker));
3323          if (j==2)
3324          {
3325            kill backmap;
3326          }
3327        }
3328
3329      if (j==1)
3330        {
3331          if (defined(fraction))
3332            {
3333              kill fraction;
3334              list fraction=p;
3335            }
3336          else
3337            {
3338              list fraction=p;
3339            }
3340        }
3341      else
3342        {
3343          fraction=insert(fraction,p,1);
3344        }
3345    }
3346  export(fraction);
3347  return ();
3348}
3349example
3350{
3351  "EXAMPLE:"; echo=2;
3352  ring R=0,(x,y),dp;
3353  ideal ker=x2+y2;
3354  export ker;
3355  list L=primeClosure(R);          // We normalize R/ker
3356  for (int i=1;i<=size(L);i++) { def R(i)=L[i]; }
3357  setring R(2);
3358  kill R;
3359  phi;                             // The map R(1)-->R(2)
3360  poly f=T(2);                     // We will get a representation of f
3361  export f;
3362  L[2]=R(2);
3363  closureFrac(L);
3364  setring R(1);
3365  kill R(2);
3366  fraction;                        // f=fraction[1]/fraction[2] via phi
3367  kill R(1);
3368}
3369
3370///////////////////////////////////////////////////////////////////////////////
3371// closureGenerators is called inside proc normal (option "withGens" )
3372//
3373
3374// INPUT is the output of proc primeClosure (except for the last element, the
3375// delta invariant) : hence input is a list L consisting of rings
3376// L[1],...,L[n] (denoted R(1)...R(n) below) such that
3377// - L[1] is a copy of (not a reference to!) the input ring L[1]
3378// - all rings L[i] contain ideals ker, L[2],...,L[n] contain ideals phi
3379// such that
3380//                L[1]/ker --> ... --> L[n]/ker
3381// are injections given by the corresponding ideals phi, and L[n]/ker
3382// is the integral closure of L[1]/ker in its quotient field.
3383// - all rings L[i] contain a polynomial nzd such that elements of
3384// L[i]/ker are quotients of elements of L[i-1]/ker with denominator
3385// nzd via the injection phi.
3386
3387// COMPUTE: In the list L of rings R(1),...,R(n), compute representations of
3388// the ring variables of the last ring R(n) as fractions of elements of R(1):
3389// The proc returns an ideal preim s.t. preim[i]/preim[size(preim)] expresses
3390// the ith variable of R(n) as fraction of elements of the basering R(1)
3391// preim[size(preim)] is a non-zero divisor of basering/i.
3392
3393proc closureGenerators(list L);
3394{
3395  def Rees=basering;         // when called inside normalI (in reesclos.lib)
3396                             // the Rees Algebra is the current basering
3397
3398  // ------- First of all we need some variable declarations -----------
3399  int n = size(L);                // the number of rings R(1)-->...-->R(n)
3400  int length = nvars(L[n]);       // the number of variables of the last ring
3401  int j,k,l,n2,n3;
3402  intvec V;
3403  string mapstr;
3404  list preimages;
3405  //Note: the empty list belongs to no ring, hence preimages can be used
3406  //later in R(1)
3407  //this is not possible for ideals (belong always to some ring)
3408
3409  for (int i=1; i<=n; i++)
3410  {
3411     ASSUME(0, not isQuotientRing(L[i]) );
3412     def R(i) = L[i];          //give the rings from L a name
3413  }
3414
3415  // For each variable (counter j) and for each intermediate ring (counter k):
3416  // Find a preimage of var_j*phi(nzd(k-1)) in R(k-1).
3417  // Finally, do the same for nzd.
3418
3419  for (j=1; j <= length+1; j++ )
3420  {
3421      setring R(n);
3422
3423      if (j==1)
3424      {
3425        poly p;
3426      }
3427      if (j <= length )
3428      {
3429        p=var(j);
3430      }
3431      else
3432      {
3433        p=1;
3434      }
3435      //i.e. p=j-th var of R(n) for j<=length and p=1 for j=length+1
3436
3437      for (k=n; k>1; k--)
3438      {
3439
3440        if (j==1)
3441        {
3442          map phimap=R(k-1),phi;   //phimap:R(k-1)-->R(n), k=2..n, is the map
3443                                   //belonging to phi in R(n)
3444        }
3445
3446        p = p*phimap(nzd);
3447
3448          // Compute the preimage of [p mod ker(k)] under phi in R(k-1):
3449          // As p is an element of Image(phi), there is a polynomial h such
3450          // that h is mapped to [p mod ker(k)], and h can be computed as the
3451          // normal form of p w.r.t. a Groebner basis of
3452          // J(k) := <ker(k),Z(l)-phi(k)(l)> in R(k)[Z]=:S(k)
3453
3454        if (j==1)   // In the first iteration: Create S(k), fetch phi and
3455                    // ker(k) and construct the ideal J(k).
3456        {
3457         //### noch abfragen ob Z(i) definiert ist
3458         list gnirlist = ringlist(R(k));
3459         n2 = size(gnirlist[2]);
3460         n3 = size(gnirlist[3]);
3461         for( i=1; i<=ncols(phi); i++)
3462         {
3463            gnirlist[2][n2+i] = "Z("+string(i)+")";
3464         }
3465         V=0;
3466         V[ncols(phi)]=0;
3467         V=V+1;
3468         gnirlist[3] = insert(gnirlist[3],list("dp",V),n3-1);
3469         def S(k) = ring(gnirlist);
3470         setring S(k);
3471
3472        // execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",
3473        //           Z(1.."+string(ncols(phi))+")),(dp("+string(nvars(R(k)))
3474        //           +"),dp("+string(ncols(phi))+"));");
3475
3476          ideal phi = imap(R(k),phi);
3477          ideal J = imap (R(k),ker);
3478          for ( l=1; l<=ncols(phi); l++ )
3479          {
3480             J=J+(Z(l)-phi[l]);
3481          }
3482          J = groebner(J);
3483          poly h = NF(imap(R(k),p),J);
3484        }
3485        else
3486        {
3487           setring S(k);
3488           h = NF(imap(R(k),p),J);
3489        }
3490
3491        setring R(k-1);
3492
3493        if (j==1)  // In the first iteration: Compute backmap:S(k)-->R(k-1)
3494        {
3495           ideal maxi;
3496           maxi[nvars(R(k))] = 0;
3497           maxi = maxi,maxideal(1);
3498           map backmap = S(k),maxi;
3499
3500           //mapstr=" map backmap = S(k),";
3501           //for (l=1;l<=nvars(R(k));l++)
3502           //{
3503           //  mapstr=mapstr+"0,";
3504           //}
3505           //execute (mapstr+"maxideal(1);");
3506
3507           poly p;
3508        }
3509        p = NF(backmap(h),std(ker));
3510     }
3511     // Whe are down to R(1), store here the result in the list preimages
3512     preimages = insert(preimages,p,j-1);
3513  }
3514  ideal preim;                  //make the list preimages to an ideal preim
3515  for ( i=1; i<=size(preimages); i++ )
3516  {
3517     preim[i] = preimages[i];
3518  }
3519  // R(1) was a copy of Rees, so we have to get back to the basering Rees from
3520  // the beginning and fetch the result (the ideal preim) to this ring.
3521  setring Rees;
3522  return (fetch(R(1),preim));
3523}
3524
3525///////////////////////////////////////////////////////////////////////////////
3526//                From here: procedures for char p with Frobenius
3527///////////////////////////////////////////////////////////////////////////////
3528
3529proc normalP(ideal id,list #)
3530"USAGE:  normalP(id [,choose]); id = radical ideal, choose = optional list of
3531         strings.
3532         Optional parameters in list choose (can be entered in any order):@*
3533         \"withRing\", \"isPrim\", \"noFac\", \"noRed\", where@*
3534         - \"noFac\" -> factorization is avoided during the computation
3535         of the minimal associated primes.@*
3536         - \"isPrim\" -> assumes that the ideal is prime. If the assumption
3537         does not hold, output might be wrong.@*
3538         - \"withRing\" -> the ring structure of the normalization is
3539         computed. The number of variables in the new ring is reduced as much
3540         as possible.@*
3541         - \"noRed\" -> when computing the ring structure, no reduction on the
3542         number of variables is done, it creates one new variable for every
3543         new module generator of the integral closure in the quotient field.@*
3544ASSUME:  The characteristic of the ground field must be positive. If the
3545         option \"isPrim\" is not set, the minimal associated primes of id
3546         are computed first and hence normalP computes the normalization of
3547         the radical of id. If option \"isPrim\" is set, the ideal must be
3548         a prime ideal otherwise the result may be wrong.
3549RETURN:  a list, say 'nor' of size 2 (resp. 3 if \"withRing\" is set).@*
3550         ** If option \"withRing\" is not set: @*
3551         Only the module structure is computed: @*
3552         * nor[1] is a list of ideals Ii, i=1..r, in the basering R where r
3553         is the number of minimal associated prime ideals P_i of the input
3554         ideal id, describing the module structure:@*
3555         If Ii is given by polynomials g_1,...,g_k in R, then c:=g_k is
3556         non-zero in the ring R/P_i and g_1/c,...,g_k/c generate the integral
3557         closure of R/P_i as R-module in the quotient field of R/P_i.@*
3558         * nor[2] shows the delta invariants: it is a list of an intvec
3559         of size r, the delta invariants of the r components, and an integer,
3560         the total delta invariant of R/id
3561         (-1 means infinite, and 0 that R/P_i resp. R/id is normal). @*
3562         ** If option \"withRing\" is set: @*
3563         The ring structure is also computed, and in this case:@*
3564         * nor[1] is a list of r rings.@*
3565         Each ring Ri = nor[1][i], i=1..r, contains two ideals with given
3566         names @code{norid} and @code{normap} such that @*
3567         - Ri/norid is the normalization of R/P_i, i.e. isomorphic as
3568           K-algebra (K the ground field) to the integral closure of R/P_i in
3569           the field of fractions of R/P_i; @*
3570         - the direct sum of the rings Ri/norid is the normalization
3571           of R/id; @*
3572         - @code{normap} gives the normalization map from R to Ri/norid.@*
3573         * nor[2] gives the module generators of the normalization of R/P_i,
3574         it is the same as nor[1] if \"withRing\" is not set.@*
3575         * nor[3] shows the delta invariants, it is the same as nor[2] if
3576         \"withRing\" is not set.
3577THEORY:  normalP uses the Leonard-Pellikaan-Singh-Swanson algorithm (using the
3578         Frobenius) cf. [A. K. Singh, I. Swanson: An algorithm for computing
3579         the integral closure, arXiv:0901.0871].
3580         The delta invariant of a reduced ring A is dim_K(normalization(A)/A).
3581         For A=K[x1,...,xn]/id we call this number also the delta invariant of
3582         id. The procedure returns the delta invariants of the components P_i
3583         and of id.
3584NOTE:    To use the i-th ring type: @code{def R=nor[1][i]; setring R;}.
3585@*       Increasing/decreasing printlevel displays more/less comments
3586         (default: printlevel = 0).
3587@*       Not implemented for local or mixed orderings or quotient rings.
3588         For local or mixed orderings use proc 'normal'.
3589@*       If the input ideal id is weighted homogeneous a weighted ordering may
3590         be used (qhweight(id); computes weights).
3591@*       Works only in characteristic p > 0; use proc normal in char 0.
3592KEYWORDS: normalization; integral closure; delta invariant.
3593SEE ALSO: normal, normalC
3594EXAMPLE: example normalP; shows an example
3595"
3596{
3597   ASSUME(0, not isQuotientRing(basering) );
3598
3599   int i,j,y, sr, del, co;
3600   intvec deli;
3601   list resu, Resu, prim, Gens, mstdid;
3602   ideal gens;
3603
3604   // Default options
3605   int wring = 0;           // The ring structure is not computed.
3606   int noRed = 0;           // No reduction is done in the ring structure
3607   int isPrim = 0;          // Ideal is not assumed to be prime
3608   int noFac = 0;           // Use facstd when computing the decomposition
3609
3610
3611   y = printlevel-voice+2;
3612
3613   if ( attrib(basering,"global") != 1)
3614   {
3615     "";
3616     "// Not implemented for this ordering,";
3617     "// please change to global ordering!";
3618     return(resu);
3619   }
3620   if ( char(basering) <= 0)
3621   {
3622     "";
3623     "// Algorithm works only in positive characteristic,";
3624     "// use procedure 'normal' if the characteristic is 0";
3625     return(resu);
3626   }
3627
3628//--------------------------- define the method ---------------------------
3629   string method;                //make all options one string in order to use
3630                                 //all combinations of options simultaneously
3631   for ( i=1; i<= size(#); i++ )
3632   {
3633     if ( typeof(#[i]) == "string" )
3634     {
3635       method = method + #[i];
3636     }
3637   }
3638
3639   if ( find(method,"withring") or find(method,"withRing") )
3640   {
3641     wring=1;
3642   }
3643   if ( find(method,"noRed") or find(method,"nored") )
3644   {
3645     noRed=1;
3646   }
3647   if ( find(method,"isPrim") or find(method,"isprim") )
3648   {
3649     isPrim=1;
3650   }
3651   if ( find(method,"noFac") or find(method,"nofac"))
3652   {
3653     noFac=1;
3654   }
3655
3656   kill #;
3657   list #;
3658//--------------------------- start computation ---------------------------
3659   ideal II,K1,K2;
3660
3661   //----------- check first (or ignore) if input id is prime -------------
3662
3663   if ( isPrim )
3664   {
3665      prim[1] = id;
3666      if( y >= 0 )
3667      { "";
3668    "// ** WARNING: result is correct if ideal is prime (not checked) **";
3669    "// disable option \"isPrim\" to decompose ideal into prime components";"";
3670      }
3671   }
3672   else
3673   {
3674      if(y>=1)
3675      {  "// compute minimal associated primes"; }
3676
3677      if( noFac )
3678      { prim = minAssGTZ(id,1); }
3679      else
3680      { prim = minAssGTZ(id); }
3681
3682      if(y>=1)
3683      {
3684         prim;"";
3685         "// number of irreducible components is", size(prim);
3686      }
3687   }
3688
3689   //----------- compute integral closure for every component -------------
3690
3691      for(i=1; i<=size(prim); i++)
3692      {
3693         if(y>=1)
3694         {
3695            ""; pause(); "";
3696            "// start computation of component",i;
3697            "   --------------------------------";
3698         }
3699         if(y>=1)
3700         {  "// compute SB of ideal";
3701         }
3702         mstdid = mstd(prim[i]);
3703         if(y>=1)
3704         {  "// dimension of component is", dim(mstdid[1]);"";}
3705
3706      //------- 1-st main subprocedure: compute module generators ----------
3707         printlevel = printlevel+1;
3708         II = normalityTest(mstdid);
3709
3710      //------ compute also the ringstructure if "withRing" is given -------
3711         if ( wring )
3712         {
3713         //------ 2-nd main subprocedure: compute ring structure -----------
3714           if(noRed == 0){
3715             resu = list(computeRing(II,prim[i])) + resu;
3716           }
3717           else
3718           {
3719             resu = list(computeRing(II,prim[i], "noRed")) + resu;
3720           }
3721         }
3722         printlevel = printlevel-1;
3723
3724      //----- rearrange module generators s.t. denominator comes last ------
3725         gens=0;
3726         for( j=2; j<=size(II); j++ )
3727         {
3728            gens[j-1]=II[j];
3729         }
3730         gens[size(gens)+1]=II[1];
3731         Gens = list(gens) + Gens;
3732      //------------------------------ compute delta -----------------------
3733         K1 = mstdid[1]+II;
3734         K1 = std(K1);
3735         K2 = mstdid[1]+II[1];
3736         K2 = std(K2);
3737         // K1 = std(mstdid[1],II);      //### besser
3738         // K2 = std(mstdid[1],II[1]);   //### besser: Hannes, fixen!
3739         co = codim(K1,K2);
3740         deli = co,deli;
3741         if ( co >= 0 && del >= 0 )
3742         {
3743            del = del + co;
3744         }
3745         else
3746         { del = -1; }
3747      }
3748
3749      if ( del >= 0 )
3750      {
3751         int mul = iMult(prim);
3752         del = del + mul;
3753      }
3754      else
3755      { del = -1; }
3756
3757      deli = deli[1..size(deli)-1];
3758      if ( wring )
3759      { Resu = resu,Gens,list(deli,del); }
3760      else
3761      { Resu = Gens,list(deli,del); }
3762
3763   sr = size(prim);
3764
3765//-------------------- Finally print comments and return --------------------
3766   if(y >= 0)
3767   {"";
3768     if ( wring )
3769     {
3770"// 'normalP' created a list, say nor, of three lists:
3771// To see the result, type
3772     nor;
3773
3774// * nor[1] is a list of",sr,"ring(s):
3775// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
3776     def R1 = nor[1][1]; setring R1;  norid; normap;
3777// for the other rings type first setring R; (if R is the name of your
3778// original basering) and then continue as for R1;
3779// Ri/norid is the affine algebra of the normalization of the i-th
3780// component R/P_i (where P_i is a min. associated prime of the input ideal)
3781// and normap the normalization map from R to Ri/norid;
3782
3783// * nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of
3784// elements g1..gk of r such that the gj/gk generate the integral
3785// closure of R/P_i as R-module in the quotient field of R/P_i.
3786
3787// * nor[3] shows the delta-invariant of each component and of the input
3788// ideal (-1 means infinite, and 0 that r/P_i is normal).";
3789     }
3790     else
3791     {
3792"// 'normalP' computed a list, say nor, of two lists:
3793// To see the result, type
3794     nor;
3795
3796// * nor[1] is a list of",sr,"ideal(s), where each ideal nor[1][i] consists
3797// of elements g1..gk of the basering R such that gj/gk generate the integral
3798// closure of R/P_i (where P_i is a min. associated prime of the input ideal)
3799// as R-module in the quotient field of R/P_i;
3800
3801// * nor[2] shows the delta-invariant of each component and of the input ideal
3802// (-1 means infinite, and 0 that R/P_i is normal).";
3803     }
3804   }
3805
3806   return(Resu);
3807}
3808example
3809{ "EXAMPLE:"; echo = 2;
3810   ring r  = 11,(x,y,z),wp(2,1,2);
3811   ideal i = x*(z3 - xy4 + x2);
3812   list nor= normalP(i); nor;
3813   //the result says that both components of i are normal, but i itself
3814   //has infinite delta
3815   pause("hit return to continue");
3816
3817   ring s = 2,(x,y),dp;
3818   ideal i = y*((x-y^2)^2 - x^3);
3819   list nor = normalP(i,"withRing"); nor;
3820
3821   def R2  = nor[1][2]; setring R2;
3822   norid; normap;
3823}
3824
3825///////////////////////////////////////////////////////////////////////////////
3826// Assume: mstdid is the result of mstd(prim[i]), prim[i] a prime component of
3827// the input ideal id of normalP.
3828// Output is an ideal U s.t. U[i]/U[1] are module generators.
3829
3830static proc normalityTest(list mstdid)
3831{
3832   ASSUME(1, not isQuotientRing(basering) );
3833
3834   int y = printlevel-voice+2;
3835   intvec op = option(get);
3836   option(redSB);
3837   def R = basering;
3838   int n, p = nvars(R), char(R);
3839   int ii;
3840
3841   ideal J = mstdid[1];         //J is the SB of I
3842   ideal I = mstdid[2];
3843   int h = n-dim(J);            //codimension of V(I), I is a prime ideal
3844
3845   //-------------------------- compute singular locus ----------------------
3846   qring Q = J;                 //pass to quotient ring
3847   ideal I = imap(R,I);
3848   ideal J = imap(R,J);
3849   attrib(J,"isSB",1);
3850   if ( y >= 1)
3851   { "start normality test";  "compute singular locus";}
3852
3853   ideal M = minor(jacob(I),h,J); //use the command minor modulo J (hence J=0)
3854   M = std(M);                    //this makes M much smaller
3855   //keep only minors which are not 0 mod I (!) this is important since we
3856   //need a nzd mod I
3857
3858   //---------------- choose nzd from ideal of singular locus --------------
3859   ideal D = M[1];
3860   for( ii=2; ii<=size(M); ii++ )            //look for the shortest one
3861   {
3862      if( size(M[ii]) < size(D[1]) )
3863      {
3864          D = M[ii];
3865      }
3866   }
3867
3868   //--------------- start p-th power algorithm and return ----------------
3869   ideal F = var(1)^p;
3870   for(ii=2; ii<=n; ii++)
3871   {
3872      F=F,var(ii)^p;
3873   }
3874
3875   ideal Dp=D^(p-1);
3876   ideal U=1;
3877   ideal K,L;
3878   map phi=Q,F;
3879   if ( y >= 1)
3880   {  "compute module generators of integral closure";
3881      "denominator D is:";  D;
3882      pause();
3883   }
3884
3885   ii=0;
3886   list LK;
3887   while(1)
3888   {
3889      ii=ii+1;
3890      if ( y >= 1)
3891      { "iteration", ii; }
3892      L = U*Dp + I;
3893      //### L=interred(L) oder mstd(L)[2]?
3894      //Wird dadurch kleiner aber string(L) wird groesser
3895      K = preimage(Q,phi,L);    //### Improvement by block ordering?
3896      option(returnSB);
3897      K = intersect(U,K);          //K is the new U, it is a SB
3898      LK = mstd(K);
3899      K = LK[2];
3900
3901   //---------------------------- simplify output --------------------------
3902      if(size(reduce(U,LK[1]))==0)  //previous U coincides with new U
3903      {                             //i.e. we reached the integral closure
3904         U=simplify(reduce(U,groebner(D)),2);
3905         U = D,U;
3906         poly gg = gcd(U[1],U[size(U)]);
3907         for(ii=2; ii<=size(U)-1 ;ii++)
3908         {
3909            gg = gcd(gg,U[ii]);
3910         }
3911         for(ii=1; ii<=size(U); ii++)
3912         {
3913            U[ii]=U[ii]/gg;
3914         }
3915         U = simplify(U,6);
3916         //if ( y >= 1)
3917         //{ "module generators are U[i]/U[1], with U:"; U;
3918         //  ""; pause(); }
3919         setring R;
3920         option(set,op);
3921         ideal U = imap(Q,U);
3922         return(U);
3923      }
3924      U=K;
3925   }
3926}
3927
3928///////////////////////////////////////////////////////////////////////////////
3929
3930static proc substpartSpecial(ideal endid, ideal endphi)
3931{
3932   ASSUME(1, not isQuotientRing(basering) );
3933
3934   //Note: newRing is of the form (R the original basering):
3935   //char(R),(T(1..N),X(1..nvars(R))),(dp(N),...);
3936
3937   int ii,jj,kk;
3938   def BAS = basering;
3939   int n = nvars(basering);
3940
3941   list Le = elimpart(endid);
3942   int q = size(Le[2]);                   //q variables have been substituted
3943//Le;"";
3944   if ( q == 0 )
3945   {
3946      ideal normap = endphi;
3947      ideal norid = endid;
3948      export(norid);
3949      export(normap);
3950      list L = BAS;
3951      return(L);
3952   }
3953
3954      list gnirlist = ringlist(basering);
3955      endid = Le[1];
3956//endphi;"";
3957      for( ii=1; ii<=n; ii++)
3958      {
3959         if( Le[4][ii] == 0 )            //ii=index of substituted var
3960         {
3961            endphi = subst(endphi,var(ii),Le[5][ii]);
3962         }
3963      }
3964//endphi;"";
3965      list g2 = gnirlist[2];             //the varlist
3966      list g3 = gnirlist[3];             //contains blocks of orderings
3967      int n3 = size(g3);
3968
3969   //----------------- first identify module ordering ------------------
3970      if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
3971      {
3972         list gm = g3[n3];              //last blockis module ordering
3973         g3 = delete(g3,n3);
3974         int m = 0;
3975      }
3976      else
3977      {
3978         list gm = g3[1];              //first block is module ordering
3979         g3 = delete(g3,1);
3980         int m = 1;
3981      }
3982   //---- delete variables which were substituted and weights  --------
3983      intvec V;
3984      int n(0);
3985      list newg2;
3986      list newg3;
3987      for ( ii=1; ii<=n3-1; ii++ )
3988      {
3989        // If order is a matrix ordering, it is replaced by dp ordering.
3990        // TODO: replace it only when some of the original
3991        //       variables are eliminated.
3992        if(g3[ii][1] == "M"){
3993          g3[ii][1] = "dp";
3994          g3[ii][2] = (1..sqroot(size(g3[ii][2])))*0+1;
3995        }
3996        V = V,g3[ii][2];           //copy weights for ordering in each block
3997        if ( ii==1 )               //into one intvector
3998        {
3999           V = V[2..size(V)];
4000        }
4001        // int n(ii) = size(g3[ii][2]);
4002        int n(ii) = size(V);
4003        intvec V(ii);
4004
4005        for ( jj = n(ii-1)+1; jj<=n(ii); jj++)
4006        {
4007          if( Le[4][jj] !=0 )     //jj=index of var which was not substituted
4008          {
4009            kk=kk+1;
4010            newg2[kk] = g2[jj];   //not substituted var from varlist
4011            V(ii)=V(ii),V[jj];    //weight of not substituted variable
4012          }
4013        }
4014        if ( size(V(ii)) >= 2 )
4015        {
4016           V(ii) = V(ii)[2..size(V(ii))];
4017           list g3(ii)=g3[ii][1],V(ii);
4018           newg3 = insert(newg3,g3(ii),size(newg3));
4019//"newg3"; newg3;
4020        }
4021      }
4022//"newg3"; newg3;
4023      //newg3 = delete(newg3,1);    //delete empty list
4024
4025/*
4026//### neue Ordnung, 1 Block fuer alle vars, aber Gewichte erhalten;
4027//vorerst nicht realisiert, da bei leonhard1 alte Version (neue Variable T(i)
4028//ein neuer Block) ein kuerzeres Ergebnis liefert
4029      kill g3;
4030      list g3;
4031      V=0;
4032      for ( ii= 1; ii<=n3-1; ii++ )
4033      {
4034        V=V,V(ii);
4035      }
4036      V = V[2..size(V)];
4037
4038      if ( V==1 )
4039      {
4040         g3[1] = list("dp",V);
4041      }
4042      else
4043      {
4044         g3[1] = lis("wp",V);
4045      }
4046      newg3 = g3;
4047
4048//"newg3";newg3;"";
4049//### Ende neue Ordnung
4050*/
4051
4052      if ( m == 0 )
4053      {
4054         newg3 = insert(newg3,gm,size(newg3));
4055      }
4056      else
4057      {
4058         newg3 = insert(newg3,gm);
4059      }
4060      gnirlist[2] = newg2;
4061      gnirlist[3] = newg3;
4062
4063//gnirlist;
4064      def newBAS = ring(gnirlist);            //change of ring to less vars
4065      setring newBAS;
4066      ideal normap = imap(BAS,endphi);
4067      //normap = simplify(normap,2);
4068      ideal norid =  imap(BAS,endid);
4069      export(norid);
4070      export(normap);
4071      list L = newBAS;
4072      return(L);
4073
4074   //Hier scheint interred gut zu sein, da es Ergebnis relativ schnell
4075   //verkleinert. Hier wird z.B. bei leonard1 size(norid) verkleinert aber
4076   //size(string(norid)) stark vergroessert, aber es hat keine Auswirkungen
4077   //da keine map mehr folgt.
4078   //### Bei Leonard2 haengt interred (BUG)
4079   //mstd[2] verkleinert norid nocheinmal um die Haelfte, dauert aber 3.71 sec
4080   //### Ev. Hinweis auf mstd in der Hilfe?
4081
4082}
4083
4084///////////////////////////////////////////////////////////////////////////////
4085// Computes the ring structure of a ring given by module generators.
4086// Assume: J[i]/J[1] are the module generators in the quotient field
4087// with J[1] as universal denominator.
4088// If option "noRed" is not given, a reduction in the number of variables is
4089// attempted.
4090static proc computeRing(ideal J, ideal I, list #)
4091{
4092  ASSUME(1, not isQuotientRing(basering) );
4093
4094  int i, ii,jj;
4095  intvec V;                          // to be used for variable weights
4096  int y = printlevel-voice+2;
4097  def R = basering;
4098  poly c = J[1];                     // the denominator
4099  list gnirlist = ringlist(basering);
4100  string svars = varstr(basering);
4101  int nva = nvars(basering);
4102  string svar;
4103  ideal maxid = maxideal(1);
4104
4105  int noRed = 0;     // By default, we try to reduce the number of generators.
4106  if(size(#) > 0){
4107    if ( typeof(#[1]) == "string" )
4108    {
4109      if (#[1] == "noRed"){noRed = 1;}
4110    }
4111  }
4112
4113  if ( y >= 1){"// computing the ring structure...";}
4114
4115  if(c==1)
4116  {
4117/*    if( defined(norid) )  { kill norid; }
4118      if( defined(normap) ) { kill normap; }
4119      ideal norid = I;
4120      ideal normap =  maxid;  */
4121
4122    def R1 = ring(gnirlist);
4123    setring R1;
4124    ideal norid = imap(R, I);
4125    ideal normap = imap(R, maxid);
4126    export norid;
4127    export normap;
4128
4129    if(noRed == 1){
4130      setring R;
4131      return(R1);
4132    }
4133    else
4134    {
4135      list L = substpartSpecial(norid,normap);
4136      def lastRing = L[1];
4137      setring R;
4138      return(lastRing);
4139    }
4140  }
4141
4142
4143  //-------------- Enlarge ring by creating new variables ------------------
4144  //check first whether variables T(i) and then whether Z(i),...,A(i) exist
4145  //old variable names are not touched
4146
4147  if ( find(svars,"T(") == 0 )
4148  {
4149    svar = "T";
4150  }
4151  else
4152  {
4153    for (ii=90; ii>=65; ii--)
4154    {
4155      if ( find(svars,ASCII(ii)+"(") == 0 )
4156      {
4157        svar = ASCII(ii);  break;
4158      }
4159    }
4160  }
4161
4162  int q = size(J)-1;
4163  if ( size(svar) != 0 )
4164  {
4165    for ( ii=q; ii>=1; ii-- )
4166    {
4167      gnirlist[2] = insert(gnirlist[2],svar+"("+string(ii)+")");
4168    }
4169  }
4170  else
4171  {
4172    for ( ii=q; ii>=1; ii-- )
4173    {
4174      gnirlist[2] = insert(gnirlist[2],"T("+string(100*nva+ii)+")");
4175    }
4176  }
4177
4178  V[q]=0;                        //create intvec of variable weights
4179  V=V+1;
4180  gnirlist[3] = insert(gnirlist[3],list("dp",V));
4181
4182  //this is a block ordering with one dp-block (1st block) for new vars
4183  //the remaining weights and blocks for old vars are kept
4184  //### perhaps better to make only one block, keeping weights ?
4185  //this might effect syz below
4186  //alt: ring newR = char(R),(X(1..nvars(R)),T(1..q)),dp;
4187  //Reihenfolge geaendert:neue Variablen kommen zuerst, Namen ev. nicht T(i)
4188
4189  def newR = ring(gnirlist);
4190  setring newR;                //new extended ring
4191  ideal I = imap(R,I);
4192
4193  //------------- Compute linear and quadratic relations ---------------
4194  if(y>=1)
4195  {
4196     "// compute linear relations:";
4197  }
4198  qring newQ = std(I);
4199
4200  ideal f = imap(R,J);
4201  module syzf = syz(f);
4202  ideal pf = f[1]*f;
4203  //f[1] is the denominator D from normalityTest, a non zero divisor of R/I
4204
4205  ideal newT = maxideal(1);
4206  newT = 1,newT[1..q];
4207  //matrix T = matrix(ideal(1,T(1..q)),1,q+1);   //alt
4208  matrix T = matrix(newT,1,q+1);
4209  ideal Lin = ideal(T*syzf);
4210  //Lin=interred(Lin);
4211  //### interred reduziert ev size aber size(string(LIN)) wird groesser
4212
4213  if(y>=1)
4214  {
4215    if(y>=3)
4216    {
4217      "//   the linear relations:";  Lin; pause();"";
4218    }
4219      "// the ring structure of the normalization as affine algebra";
4220      "//   number of linear relations:", size(Lin);
4221  }
4222
4223  if(y>=1)
4224  {
4225    "// compute quadratic relations:";
4226  }
4227  matrix A;
4228  ideal Quad;
4229  poly ff;
4230  newT = newT[2..size(newT)];
4231  matrix u;  // The units for non-global orderings.
4232
4233  // Quadratic relations
4234  for (ii=2; ii<=q+1; ii++ )
4235  {
4236    for ( jj=2; jj<=ii; jj++ )
4237    {
4238      ff = NF(f[ii]*f[jj],std(0));     // this makes lift much faster
4239      // For non-global orderings, we have to take care of the units.
4240      if(attrib(basering,"global") != 1)
4241      {
4242        A = lift(pf, ff, u);
4243        Quad = Quad,ideal(newT[jj-1]*newT[ii-1] * u[1, 1]- T*A);
4244      }
4245      else
4246      {
4247        A = lift(pf,ff);              // ff lin. comb. of elts of pf mod I
4248        Quad = Quad,ideal(newT[jj-1]*newT[ii-1] - T*A);
4249      }
4250      //A = lift(pf, f[ii]*f[jj]);
4251      //Quad = Quad, ideal(T(jj-1)*T(ii-1) - T*A);
4252    }
4253  }
4254  Quad = Quad[2..ncols(Quad)];
4255
4256  if(y>=1)
4257  {
4258    if(y>=3)
4259    {
4260      "//   the quadratic relations"; Quad; pause();"";
4261    }
4262      "//   number of quadratic relations:", size(Quad);
4263  }
4264  ideal Q1 = Lin,Quad;     //elements of Q1 are in NF w.r.t. I
4265
4266  //Q1 = mstd(Q1)[2];
4267  //### weglassen, ist sehr zeitaufwendig.
4268  //Ebenso interred, z.B. bei Leonard1 (1. Komponente von Leonard):
4269  //"size Q1:", size(Q1), size(string(Q1));   //75 60083
4270  //Q1 = interred(Q1);
4271  //"size Q1:", size(Q1), size(string(Q1));   //73 231956 (!)
4272  //### Speicherueberlauf bei substpartSpecial bei 'ideal norid  = phi1(endid)'
4273  //Beispiel fuer Hans um map zu testen!
4274
4275  setring newR;
4276  ideal endid  = imap(newQ,Q1),I;
4277  ideal endphi = imap(R,maxid);
4278
4279  if(noRed == 0){
4280    list L=substpartSpecial(endid,endphi);
4281    def lastRing=L[1];
4282    if(y>=1)
4283    {
4284      "//   number of substituted variables:", nvars(newR)-nvars(lastRing);
4285      pause();"";
4286    }
4287    setring R;
4288    return(lastRing);
4289  }
4290  else
4291  {
4292    ideal norid = endid;
4293    ideal normap = endphi;
4294    export(norid);
4295    export(normap);
4296    setring R;
4297    return(newR);
4298  }
4299}
4300
4301//                Up to here: procedures for char p with Frobenius
4302///////////////////////////////////////////////////////////////////////////////
4303
4304///////////////////////////////////////////////////////////////////////////////
4305//                From here: subprocedures for normal
4306
4307// inputJ is used in parametrization of rational curves algorithms, to specify
4308// a different test ideal.
4309
4310static proc normalM(ideal I, int decomp, int withDelta, int denomOption, ideal inputJ, ideal inputC){
4311// Computes the normalization of a ring R / I using the module structure as far
4312// as possible.
4313// The ring R is the basering.
4314// Input: ideal I
4315// Output: a list of 3 elements (resp 4 if withDelta = 1), say nor.
4316// - nor[1] = U, an ideal of R.
4317// - nor[2] = c, an element of R.
4318// U and c satisfy that 1/c * U is the normalization of R / I in the
4319// quotient field Q(R/I).
4320// - nor[3] = ring say T, containing two ideals norid and normap such that
4321// normap gives the normalization map from R / I to T / norid.
4322// - nor[4] = the delta invariant, if withDelta = 1.
4323
4324// Details:
4325// --------
4326// Computes the ideal of the minors in the first step and then reuses it in all
4327// steps.
4328// In step s, the denominator is D^s, where D is a nzd of the original quotient
4329// ring, contained in the radical of the singular locus.
4330// This denominator is used except when the degree of D^i is greater than the
4331// degree of a universal denominator.
4332// The nzd is taken as the smallest degree polynomial in the radical of the
4333// singular locus.
4334
4335// It assumes that the ideal I is equidimensional radical. This is not checked
4336// in the procedure!
4337// If decomp = 0 or decomp = 3 it assumes further that I is prime. Therefore
4338// any non-zero element in the jacobian ideal is assumed to be a
4339// non-zerodivisor.
4340
4341// It works over the given basering.
4342// If it has a non-global ordering, it changes it to dp global only for
4343// computing radical.
4344
4345// The delta invariant is only computed if withDelta = 1, and decomp = 0 or
4346// decomp = 3 (meaning that the ideal is prime).
4347
4348// denomOption = 0      -> Uses the smallest degree polynomial
4349// denomOption = i > 0  -> Uses a polynomial in the i-th variable
4350
4351  ASSUME(1, not isQuotientRing(basering) );
4352
4353  intvec save_opt=option(get);
4354  option(redSB);
4355  option(returnSB);
4356  int step = 0;                       // Number of steps. (for debugging)
4357  int dbg = printlevel - voice + 2;   // dbg = printlevel (default: dbg = 0)
4358  int i;                              // counter
4359  int isGlobal = attrib(basering,"global");
4360
4361  poly c;                     // The common denominator.
4362
4363  def R = basering;
4364
4365//------------------------ Groebner bases and dimension of I-----------------
4366  if(isGlobal == 1)
4367  {
4368    list IM = mstd(I);
4369    I = IM[1];
4370    ideal IMin = IM[2];   // A minimal set of generators in the groebner basis.
4371  }
4372  else
4373  {
4374    // The minimal set of generators is not computed by mstd for
4375    // non-global orderings.
4376    I = groebner(I);
4377    ideal IMin = I;
4378  }
4379  int d = dim(I);
4380
4381  // ---------------- computation of the singular locus ---------------------
4382  // We compute the radical of the ideal of minors modulo the original ideal.
4383  // This is done only in the first step.
4384  qring Q = I;   // We work in the quotient by the groebner base of the ideal I
4385  option(redSB);
4386  option(returnSB);
4387
4388  // If a conductor ideal was given as input, we use it instead of the
4389  // singular locus. If a test ideal was given as input, we do not compute the
4390  // singular locus.
4391  ideal inputC = fetch(R, inputC);
4392  ideal inputJ = fetch(R, inputJ);
4393  if((inputC == 0) && (inputJ == 0))
4394  {
4395    // We compute the radical of the ideal of minors modulo the original ideal.
4396    // This is done only in the first step.
4397    ideal I = fetch(R, I);
4398    attrib(I, "isSB", 1);
4399    ideal IMin = fetch(R, IMin);
4400
4401    dbprint(dbg, "Computing the jacobian ideal...");
4402
4403    // If a given conductor ideal is given, we use it.
4404    // If a given test ideal is given, we don't need to compute the jacobian
4405
4406    // reduction mod I in 'minor' is not working for local orderings!
4407    if(attrib(basering,"global"))
4408    {
4409      ideal J = minor(jacob(IMin), nvars(basering) - d, I);
4410    }
4411    else
4412    {
4413      ideal J = minor(jacob(IMin), nvars(basering) - d);
4414      J = reduce(J, groebner(I));
4415    }
4416    J = groebner(J);
4417  }
4418  else
4419  {
4420    ideal J = fetch(R, inputC);
4421    J = groebner(J);
4422  }
4423
4424  //------------------ We check if the singular locus is empty -------------
4425  if(J[1] == 1)
4426  {
4427    // The original ring R/I was normal. Nothing to do.
4428    // We define anyway a new ring, equal to R, to be able to return it.
4429    setring R;
4430    list lR = ringlist(R);
4431    def ROut = ring(lR);
4432    setring ROut;
4433    ideal norid = fetch(R, I);
4434    ideal normap = maxideal(1);
4435    export norid;
4436    export normap;
4437    setring R;
4438    if(withDelta)
4439    {
4440      list output = ideal(1), poly(1), ROut, 0;
4441    }
4442    else
4443    {
4444      list output = ideal(1), poly(1), ROut;
4445    }
4446    option(set,save_opt);
4447    return(list(output));
4448  }
4449
4450
4451  // -------------------- election of the universal denominator----------------
4452  // We first check if a conductor ideal was computed. If not, we don't
4453  // compute a universal denominator.
4454  ideal Id1;
4455  if(J != 0)
4456  {
4457    if(denomOption == 0)
4458    {
4459      poly condu = getSmallest(J);   // Choses the polynomial of smallest degree
4460                                     // of J as universal denominator.
4461    }
4462    else
4463    {
4464      poly condu = getOneVar(J, denomOption);
4465    }
4466    if(dbg >= 1)
4467    {
4468      "";
4469      "The universal denominator is ", condu;
4470    }
4471
4472    // ----- splitting the ideal by the universal denominator (if possible) -----
4473    // If the ideal is equidimensional, but not necessarily prime, we check if
4474    // the universal denominator is a non-zerodivisor of R/I.
4475    // If not, we split I.
4476    if((decomp == 1) or (decomp == 2))
4477    {
4478      Id1 = quotient(0, condu);
4479      if(size(Id1) > 0)
4480      {
4481        // We have to split.
4482        if(dbg >= 1)
4483        {
4484          "A zerodivisor was found. We split the ideal. The zerodivisor is ", condu;
4485        }
4486        setring R;
4487        ideal Id1 = fetch(Q, Id1), I;
4488        Id1 = groebner(Id1);
4489        ideal Id2 = quotient(I, Id1);
4490        // I = I1 \cap I2
4491        printlevel = printlevel + 1;
4492        ideal JDefault = 0; // Now it uses the default J;
4493        list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault);
4494        list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault);
4495        printlevel = printlevel - 1;
4496        option(set,save_opt);
4497        list res = nor1 + nor2;
4498        return(res);
4499      }
4500    }
4501  }
4502  else
4503  {
4504    poly condu = 0;
4505  }
4506
4507  // --------------- computation of the first test ideal ---------------------
4508  // To compute the radical we go back to the original ring.
4509  // If we are using a non-global ordering, we must change to the global
4510  // ordering.
4511  setring R;
4512  // If a test ideal is given at the input, we use it.
4513  if(inputJ == 0)
4514  {
4515    if(isGlobal == 1)
4516    {
4517      ideal J = fetch(Q, J);
4518      J = J, I;
4519      if(dbg >= 1)
4520      {
4521        "The original singular locus is";
4522        groebner(J);
4523        if(dbg >= 2){pause();}
4524        "";
4525      }
4526      // We check if the only singular point is the origin.
4527      // If so, the radical is the maximal ideal at the origin.
4528      J = groebner(J);
4529      if(locAtZero(J))
4530      {
4531        J = maxideal(1);
4532      }
4533      else
4534      {
4535        J = radical(J);
4536      }
4537    }
4538    else
4539    {
4540      // We change to global dp ordering.
4541      list rl = ringlist(R);
4542      list origOrd = rl[3];
4543      list newOrd = list("dp", intvec(1:nvars(R))), list("C", 0);
4544      rl[3] = newOrd;
4545      def globR = ring(rl);
4546      setring globR;
4547      ideal J = fetch(Q, J);
4548      ideal I = fetch(R, I);
4549      J = J, I;
4550      if(dbg >= 1)
4551      {
4552        "The original singular locus is";
4553        groebner(J);
4554        if(dbg>=2){pause();}
4555        "";
4556      }
4557      J = radical(J);
4558      setring R;
4559      ideal J = fetch(globR, J);
4560    }
4561  }
4562  else
4563  {
4564    ideal J = inputJ;
4565  }
4566
4567  if(dbg >= 1)
4568  {
4569    "The radical of the original singular locus is";
4570    J;
4571    if(dbg>=2){pause();}
4572  }
4573
4574  // ---------------- election of the non zero divisor ---------------------
4575  setring Q;
4576  J = fetch(R, J);
4577  J = interred(J);
4578  if(denomOption == 0)
4579  {
4580    poly D = getSmallest(J);    // Chooses the polynomial of smallest degree as
4581                                // non-zerodivisor.
4582  }
4583  else
4584  {
4585    poly D = getOneVar(J, denomOption);
4586  }
4587  if(dbg >= 1)
4588  {
4589    "The non zero divisor is ", D;
4590    "";
4591  }
4592
4593  // ------- splitting the ideal by the non-zerodivisor (if possible) --------
4594  // If the ideal is equidimensional, but not necessarily prime, we check if D
4595  // is actually a non-zerodivisor of R/I.
4596  // If not, we split I.
4597  if((decomp == 1) or (decomp == 2))
4598  {
4599    // We check if D is actually a non-zerodivisor of R/I.
4600    // If not, we split I.
4601    Id1 = quotient(0, D);
4602    if(size(Id1) > 0)
4603    {
4604      // We have to split.
4605      if(dbg >= 1)
4606      {
4607        "A zerodivisor was found. We split the ideal. The zerodivisor is ", D;
4608      }
4609      setring R;
4610      ideal Id1 = fetch(Q, Id1), I;
4611      Id1 = groebner(Id1);
4612      ideal Id2 = quotient(I, Id1);
4613      // I = Id1 \cap Id2
4614      printlevel = printlevel + 1;
4615
4616      ideal JDefault = 0;  // Now it uses the default J;
4617      list nor1 = normalM(Id1, decomp, withDelta, denomOption, JDefault, JDefault);
4618      list nor2 = normalM(Id2, decomp, withDelta, denomOption, JDefault, JDefault);
4619      printlevel = printlevel - 1;
4620      option(set,save_opt);
4621      list res = nor1 + nor2;
4622      return(res);
4623    }
4624  }
4625
4626  // --------------------- normalization ------------------------------------
4627  // We call normalMEqui to compute the normalization.
4628  setring R;
4629  poly D = fetch(Q, D);
4630  poly condu = fetch(Q, condu);
4631  J = fetch(Q, J);
4632  printlevel = printlevel + 1;
4633  list result = normalMEqui(I, J, condu, D, withDelta, denomOption);
4634  printlevel = printlevel - 1;
4635  option(set,save_opt);
4636  return(list(result));
4637}
4638
4639///////////////////////////////////////////////////////////////////////////////
4640
4641static proc normalMEqui(ideal I, ideal origJ, poly condu, poly D, int withDelta)
4642// Here is where the normalization is actually computed.
4643
4644// Computes the normalization of R/I. (basering is R)
4645// I is assumed to be radical and equidimensional.
4646// origJ is the first test ideal.
4647// D is a non-zerodivisor of R/I.
4648// condu is a non-zerodivisor in the conductor or 0 if it was not computed.
4649// If withDelta = 1, computes the delta invariant.
4650{
4651  ASSUME(1, not isQuotientRing(basering) );
4652  int step = 0;                       // Number of steps. (for debugging)
4653  int dbg = printlevel - voice + 2;   // dbg = printlevel (default: dbg = 0)
4654  int i;                              // counter
4655  int isNormal = 0;                   // check for exiting the loop
4656  int isGlobal = attrib(basering,"global");
4657  int delt;
4658
4659  def R = basering;
4660  poly c = D;
4661  ideal U;
4662  ideal cJ;
4663  list testOut;                 // Output of proc testIdeal
4664                                // (the test ideal and the ring structure)
4665
4666  qring Q = groebner(I);
4667  intvec save_opt=option(get);
4668  option(redSB);
4669  option(returnSB);
4670  ideal J = imap(R, origJ);
4671  poly c = imap(R, c);
4672  poly D = imap(R, D);
4673  poly condu = imap(R, condu);
4674  ideal cJ;
4675  ideal cJMod;
4676
4677  dbprint(dbg, "Preliminar step begins.");
4678
4679  // --------------------- computation of A1 -------------------------------
4680  dbprint(dbg, "Computing the quotient (DJ : J)...");
4681  ideal U = groebner(quotient(D*J, J));
4682  ideal oldU = 1;
4683
4684  if(dbg >= 2) { "The quotient is"; U; }
4685
4686  // ----------------- Grauer-Remmert criterion check -----------------------
4687  // We check if the equality in Grauert - Remmert criterion is satisfied.
4688  isNormal = checkInclusions(D*oldU, U);
4689  if(isNormal == 0)
4690  {
4691    if(dbg >= 1)
4692    {
4693      "In this step, we have the ring 1/c * U, with c =", c;
4694      "and U = "; U;
4695    }
4696  }
4697  else
4698  {
4699    // The original ring R/I was normal. Nothing to do.
4700    // We define anyway a new ring, equal to R, to be able to return it.
4701    setring R;
4702    list lR = ringlist(R);
4703    def ROut = ring(lR);
4704    setring ROut;
4705    ideal norid = fetch(R, I);
4706    ideal normap = maxideal(1);
4707    export norid;
4708    export normap;
4709    setring R;
4710    if(withDelta)
4711    {
4712      list output = ideal(1), poly(1), ROut, 0;
4713    }
4714    else
4715    {
4716      list output = ideal(1), poly(1), ROut;
4717    }
4718    option(set,save_opt);
4719    return(output);
4720  }
4721
4722  // ----- computation of the chain of ideals A1 c A2 c ... c An ------------
4723  while(isNormal == 0)
4724  {
4725    step++;
4726    if(dbg >= 1) { ""; "Step ", step, " begins."; }
4727    dbprint(dbg, "Computing the test ideal...");
4728
4729    // --------------- computation of the test ideal ------------------------
4730    // Computes a test ideal for the new ring.
4731    // The test ideal will be the radical in the new ring of the original
4732    // test ideal.
4733    setring R;
4734    U = imap(Q, U);
4735    c = imap(Q, c);
4736    testOut = testIdeal(I, U, origJ, c, D);
4737    cJ = testOut[1];
4738
4739    setring Q;
4740    cJ = imap(R, cJ);
4741    cJ = groebner(cJ);
4742
4743    // cJ / c is now the ideal mapped back.
4744    // We have the generators as an ideal in the new ring,
4745    // but we want the generators as an ideal in the original ring.
4746    cJMod = getGenerators(cJ, U, c);
4747
4748    if(dbg >= 2) { "The test ideal in this step is "; cJMod; }
4749
4750    cJ = cJMod;
4751
4752    // ------------- computation of the quotient (DJ : J)--------------------
4753    oldU = U;
4754    dbprint(dbg, "Computing the quotient (c*D*cJ : cJ)...");
4755    U = quotient(c*D*cJ, cJ);
4756    if(dbg >= 2){"The quotient is "; U;}
4757
4758    // ------------- Grauert - Remmert criterion check ----------------------
4759    // We check if the equality in Grauert - Remmert criterion is satisfied.
4760    isNormal = checkInclusions(D*oldU, U);
4761
4762    if(isNormal == 1)
4763    {
4764      // We go one step back. In the last step we didnt get antyhing new,
4765      // we just verified that the ring was already normal.
4766      dbprint(dbg, "The ring in the previous step was already normal.");
4767      dbprint(dbg, "");
4768      U = oldU;
4769    }
4770    else
4771    {
4772      // ------------- preparation for next iteration ----------------------
4773      // We have to go on.
4774      // The new denominator is chosen.
4775      c = D * c;
4776
4777      // If we have a universal denominator of smaller degree than c,
4778      // we replace c by it.
4779      if(condu != 0)
4780      {
4781        if(deg(c) > deg(condu))
4782        {
4783          U = changeDenominatorQ(U, c, condu);
4784          c = condu;
4785        }
4786      }
4787      if(dbg >= 1)
4788      {
4789        "In this step, we have the ring 1/c * U, with c =", c;
4790        "and U = ";
4791        U;
4792        if(dbg>=2){pause();}
4793      }
4794    }
4795  }
4796
4797  // ------------------------- delta computation ----------------------------
4798  if(withDelta)
4799  {
4800    ideal UD = groebner(U);
4801    delt = vdim(std(modulo(UD, c)));
4802  }
4803
4804  // -------------------------- prepare output -----------------------------
4805  setring R;
4806  U = fetch(Q, U);
4807  c = fetch(Q, c);
4808
4809  // Ring structure of the new ring
4810  def ere = testOut[2];
4811  if(withDelta)
4812  {
4813    list output = U, c, ere, delt;
4814  }
4815  else
4816  {
4817    list output = U, c, ere;
4818  }
4819  option(set,save_opt);
4820  return(output);
4821}
4822
4823///////////////////////////////////////////////////////////////////////////////
4824
4825static proc lineUpLast(ideal U, poly c)
4826// Sets c as the last generator of U.
4827{
4828  int i;
4829  ideal newU;
4830  for (i = 1; i <= ncols(U); i++)
4831  {
4832    if(U[i] != c)
4833    {
4834      if(size(newU) == 0)
4835      { newU = U[i]; }
4836      else
4837      { newU = newU, U[i]; }
4838    }
4839  }
4840  if(size(newU) == 0)
4841  { newU = c; }
4842  else
4843  { newU = newU, c; }
4844  return(newU);
4845}
4846
4847///////////////////////////////////////////////////////////////////////////////
4848
4849static proc lineUp(ideal U, poly c)
4850// Sets c as the first generator of U.
4851{
4852  int i;
4853  ideal newU = c;
4854  for (i = 1; i <= ncols(U); i++)
4855  {
4856    if(U[i] != c)
4857    {
4858      newU = newU, U[i];
4859    }
4860  }
4861  return(newU);
4862}
4863
4864///////////////////////////////////////////////////////////////////////////////
4865
4866//WARNING - elim is not working here!! Check!!
4867//It is now replaced by computing an eliminating groebner basis.
4868proc getOneVar(ideal J, int vari)
4869"USAGE:   getOneVar(J, vari); J is a 0-dimensional ideal, vari is an integer.
4870RETURN:  a polynomial of J in the variable indicated by vari of smallest
4871         degree.@*
4872NOTE:    Works only over rings of two variables.@*
4873         It is intended mainly as an auxiliary procedure for computing
4874         integral bases. @*
4875EXAMPLE: example getOneVar; shows an example
4876"
4877{
4878  ASSUME(0, nvars(basering)==2 );
4879  ASSUME(0, (vari==2) || (vari==1) );
4880
4881  def R = basering;
4882  list RL = ringlist(R);
4883  // We keep everything from R but we change the ordering to lp, and we
4884  // order the variables as needed.
4885  RL[3] = list(list("lp", 1:2), list("C", 0:1));
4886  RL[2] = list(var(3-vari), var(vari));
4887  RL[4]=ideal(0); // does not work with qrings: Ex.7 of paraplanecurves
4888  def RR = ring(RL);
4889  setring RR;
4890  ideal J = imap(R, J);
4891  J = groebner(J);
4892  poly g = J[1];
4893  setring R;
4894  poly g = imap(RR, g);
4895  return(g);
4896}
4897example
4898{ "EXAMPLE:";
4899  printlevel = printlevel+1;
4900  echo = 2;
4901  ring s = 0,(x,y),dp;
4902  ideal J = x3-y, y3;
4903  getOneVar(J, 1);
4904
4905  echo = 0;
4906  printlevel = printlevel-1;
4907}
4908///////////////////////////////////////////////////////////////////////////////
4909
4910proc getSmallest(ideal J)
4911"USAGE:   getSmallest(J); J is an ideal.
4912RETURN:  the generator of J of smallest degree. If there are more than one, it
4913         chooses the one with smallest number of monomials.@*
4914NOTE:    It looks only at the generator of J, not at all the polynomials in
4915         the ideal.@*
4916         It is intended maninly to compute a good universal denominator in the
4917         normalization algorithms.@*
4918EXAMPLE: example getSmallest; shows an example
4919"
4920{
4921
4922// Computes the polynomial of smallest degree of J.
4923//
4924  int i;
4925  poly p = J[1];
4926  int d = deg(p);
4927  int di;
4928  for(i = 2; i <= ncols(J); i++)
4929  {
4930    if(J[i] != 0)
4931    {
4932      di = deg(J[i]);
4933      if(di < d)
4934      {
4935        p = J[i];
4936        d = di;
4937      }
4938      else
4939      {
4940        if(di == d)
4941        {
4942          if(size(J[i]) < size(p))
4943          {
4944            p = J[i];
4945          }
4946        }
4947      }
4948    }
4949  }
4950  return(p);
4951}
4952example
4953{ "EXAMPLE:";
4954  printlevel = printlevel+1;
4955  echo = 2;
4956  ring s = 0,(x,y),dp;
4957  ideal J = x3-y, y5, x2-y2+1;
4958  getSmallest(J);
4959
4960  echo = 0;
4961  printlevel = printlevel-1;
4962}
4963
4964///////////////////////////////////////////////////////////////////////////////
4965
4966static proc getGenerators(ideal J, ideal U, poly c)
4967{
4968
4969// Computes the generators of J as an ideal in the original ring,
4970// where J is given by generators in the new ring.
4971
4972// The new ring is given by 1/c * U in the total ring of fractions.
4973
4974  int i, j;                             // counters;
4975  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)
4976  poly p;                               // The lifted polynomial
4977  ideal JGr = groebner(J);              // Groebner base of J
4978
4979  if(dbg>1){"Checking for new generators...";}
4980  for(i = 1; i <= ncols(J); i++)
4981  {
4982    for(j = 1; j <= ncols(U); j++)
4983    {
4984      p = lift(c, J[i]*U[j])[1,1];
4985      p = reduce(p, JGr);
4986      if(p != 0)
4987      {
4988        if(dbg>1)
4989        {
4990          "New polynoial added:", p;
4991          if(dbg>4) {pause();}
4992        }
4993        JGr = JGr, p;
4994        JGr = groebner(JGr);
4995        J = J, p;
4996      }
4997    }
4998  }
4999  return(J);
5000}
5001
5002///////////////////////////////////////////////////////////////////////////////
5003
5004static proc testIdeal(ideal I, ideal U, ideal origJ, poly c, poly D)
5005{
5006
5007  ASSUME(1, not isQuotientRing(basering) );
5008
5009// Internal procedure, used in normalM.
5010// Computes the test ideal in the new ring.
5011// It takes the original test ideal and computes the radical of it in the
5012// new ring.
5013
5014// The new ring is 1/c * U.
5015// The original test ideal is origJ.
5016// The original ring is R / I, where R is the basering.
5017  intvec save_opt=option(get);
5018  int i;                                // counter
5019  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)
5020  def R = basering;                      // We dont work in the quo
5021  ideal J = origJ;
5022
5023  // ---------- computation of the ring structure of 1/c * U ----------------
5024  U = lineUp(U, c);
5025
5026  if(dbg > 1){"Computing the new ring structure...";}
5027  list ele = computeRing(U, I, "noRed");
5028
5029  def origEre = ele[1];
5030  setring origEre;
5031  if(dbg > 1){"The relations are"; norid;}
5032
5033  // ---------------- setting the ring to work in  --------------------------
5034  int isGlobal = attrib(origEre,"global");      // Checks if the original ring has
5035                                         // global ordering.
5036  if(isGlobal != 1)
5037  {
5038    list rl = ringlist(origEre);
5039    list origOrd = rl[3];
5040    list newOrd = list("dp", intvec(1:nvars(origEre))), list("C", 0);
5041    rl[3] = newOrd;
5042    def ere = ring(rl);     // globR is the original ring but
5043                            // with a global ordering.
5044    setring ere;
5045    ideal norid = imap(origEre, norid);
5046  }
5047  else
5048  {
5049    def ere = origEre;
5050  }
5051
5052  ideal I = imap(R, I);
5053  ideal J = imap(R, J);
5054  J = J, norid, I;
5055
5056
5057  // ----- computation of the test ideal using the ring structure of Ai -----
5058
5059  option(redSB);
5060  option(returnSB);
5061
5062  if(dbg > 1){"Computing the radical of J...";}
5063  J = radical(J);
5064  if(dbg > 1){"Computing the interreduction of the radical...";}
5065  J = groebner(J);
5066  //J = interred(J);
5067  if(dbg > 1)
5068  {
5069    "The radical in the generated ring is";
5070    J;
5071    if(dbg>4){pause();}
5072  }
5073
5074  setring ere;
5075
5076  // -------------- map from Ai to the total ring of fractions ---------------
5077  // Now we must map back this ideal J to U_i / c in the total ring of
5078  // fractions.
5079  // The map sends T_j -> u_j / c.
5080  // The map is built by the following steps:
5081  // 1) We compute the degree of the generators of J with respect to the
5082  //    new variables T_j.
5083  // 2) For each generator, we multiply each term by a power of c, as if
5084  //    taking c^n as a common denominator (considering the new variables as
5085  //    a polynomial in the old variables divided by c).
5086  // 3) We replace the new variables T_j by the corresponding numerator u_j.
5087  // 4) We lift the resulting polynomial to change the denominator
5088  //    from c^n to c.
5089  int nNewVars = nvars(ere) - nvars(R);      // Number of new variables
5090  poly c = imap(R, c);
5091  intvec @v = 1..nNewVars;    // Vector of the new variables.
5092                              // They must be the first ones.
5093  if(dbg > 1){"The indices of the new variables are", @v;}
5094
5095  // ---------------------- step 1 of the mapping ---------------------------
5096  intvec degs;
5097  for(i = 1; i<=ncols(J); i++)
5098  {
5099    degs[i] = degSubring(J[i], @v);
5100  }
5101  if(dbg > 1)
5102  {
5103    "The degrees with respect to the new variables are";
5104    degs;
5105  }
5106
5107  // ---------------------- step 2 of the mapping ---------------------------
5108  ideal mapJ = mapBackIdeal(J, c, @v);
5109
5110  setring R;
5111
5112  // ---------------------- step 3 of the mapping ---------------------------
5113  ideal z;                    // The variables of the original ring in order.
5114  for(i = 1; i<=nvars(R); i++)
5115  {
5116    z[i] = var(i);
5117  }
5118
5119  map f = ere, U[2..ncols(U)], z[1..ncols(z)]; // The map to the original ring.
5120  if(dbg > 1)
5121  {
5122    "The map is ";
5123    f;
5124    if(dbg>4){pause();}
5125  }
5126
5127  if(dbg > 1){ "Computing the map..."; }
5128
5129  J = f(mapJ);
5130  if(dbg > 1)
5131  {
5132    "The ideal J mapped back (before lifting) is";
5133    J;
5134    if(dbg>4){pause();}
5135  }
5136
5137  // ---------------------- step 4 of the mapping ---------------------------
5138  qring Q = groebner(I);
5139  ideal J = imap(R, J);
5140  poly c = imap(R, c);
5141  for(i = 1; i<=ncols(J); i++)
5142  {
5143    if(degs[i]>1)
5144    {
5145      J[i] = lift(c^(degs[i]-1), J[i])[1,1];
5146    }
5147    else
5148    {
5149      if(degs[i]==0) { J[i] = c*J[i]; }
5150    }
5151  }
5152
5153  if(dbg > 1)
5154  {
5155    "The ideal J lifted is";
5156    J;
5157    if(dbg>4){pause();}
5158  }
5159
5160  // --------------------------- prepare output ----------------------------
5161  J = groebner(J);
5162
5163  setring R;
5164  J = imap(Q, J);
5165
5166  option(set,save_opt);
5167  return(list(J, ele[1]));
5168}
5169
5170///////////////////////////////////////////////////////////////////////////////
5171
5172proc changeDenominator(ideal U1, poly c1, poly c2, ideal I)
5173"USAGE:   changeDenominator(U1, c1, c2, I); U1 and I ideals, c1 and c2
5174         polynomials.@*
5175RETURN:  an ideal U2 such that the A-modules 1/c1 * U1 and 1/c2 * U2 are equal,
5176         where A = R/I and R is the basering.@*
5177NOTE:    It assumes that such U2 exists. It is intended maninly as an auxiliary
5178         procedure in the normalization algorithms.@*
5179EXAMPLE: example changeDenominator; shows an example
5180"
5181{
5182  ASSUME(0, not isQuotientRing(basering) );
5183// Let A = R / I. Given an A-module in the form 1/c1 * U1 (U1 ideal of A), it
5184// computes a new ideal U2 such that the the A-module is 1/c2 * U2.
5185// The base ring is R, but the computations are to be done in R / I.
5186  int a;      // counter
5187  def R = basering;
5188  qring Q = I;
5189  ideal U1 = fetch(R, U1);
5190  poly c1 = fetch(R, c1);
5191  poly c2 = fetch(R, c2);
5192  ideal U2 = changeDenominatorQ(U1, c1, c2);
5193  setring R;
5194  ideal U2 = fetch(Q, U2);
5195  return(U2);
5196}
5197example
5198{
5199  "EXAMPLE:";
5200  echo = 2;
5201  ring s = 0,(x,y),dp;
5202  ideal I = y5-y4x+4y2x2-x4;
5203  ideal U1 = normal(I)[2][1];
5204  poly c1 = U1[4];
5205  U1;c1;
5206  // 1/c1 * U1 is the normalization of I.
5207  ideal U2 = changeDenominator(U1, c1, x3, I);
5208  U2;
5209  // 1/x3 * U2 is also the normalization of I, but with a different denominator.
5210  echo = 0;
5211}
5212
5213///////////////////////////////////////////////////////////////////////////////
5214
5215static proc changeDenominatorQ(ideal U1, poly c1, poly c2)
5216{
5217// Given a ring in the form 1/c1 * U, it computes a new U2 st the ring
5218// is 1/c2 * U2.
5219// The base ring is already a quotient ring R / I.
5220  int a;      // counter
5221  ideal U2;
5222  poly p;
5223  for(a = 1; a <= ncols(U1); a++)
5224  {
5225    p = lift(c1, c2*U1[a])[1,1];
5226    U2[a] = p;
5227  }
5228  return(U2);
5229}
5230
5231///////////////////////////////////////////////////////////////////////////////
5232
5233static proc checkInclusions(ideal U1, ideal U2)
5234{
5235// Checks if the identity A = Hom(J, J) of Grauert-Remmert criterion is
5236// satisfied.
5237  int dbg = printlevel - voice + 2;     // dbg = printlevel (default: dbg = 0)
5238  list reduction1;
5239  list reduction2;
5240
5241  // ---------------------- inclusion Hom(J, J) c A -------------------------
5242  if(dbg > 1){"Checking the inclusion Hom(J, J) c A:";}
5243  // This interred is used only because a bug in groebner!
5244  U1 = groebner(U1);
5245  reduction1 = reduce(U2, U1);
5246  if(dbg > 1){reduction1[1];}
5247
5248  // ---------------------- inclusion A c Hom(J, J) -------------------------
5249  // The following check should always be satisfied.
5250  // This is only used for debugging.
5251  if(dbg > 1)
5252  {
5253    "and the inclusion A c Hom(J, J): (this should always be satisfied)";
5254    // This interred is used only because a bug in groebner!
5255    U2 = groebner(U2);
5256    reduction2 = reduce(U1, groebner(U2));
5257    reduction2[1];
5258    if(size(reduction2[1]) > 0)
5259    {
5260      "Something went wrong... (this inclusion should always be satisfied)";
5261      ~;
5262    }
5263    else
5264    {
5265      if(dbg>4){pause();}
5266    }
5267  }
5268
5269  if(size(reduction1[1]) == 0)
5270  {
5271    // We are done! The ring computed in the last step was normal.
5272    return(1);
5273  }
5274  else
5275  {
5276    return(0);
5277  }
5278}
5279
5280///////////////////////////////////////////////////////////////////////////////
5281
5282static proc degSubring(poly p, intvec @v)
5283{
5284  ASSUME(1, not isQuotientRing(basering) );
5285// Computes the degree of a polynomial taking only some variables as variables
5286// and the others as parameters.
5287
5288// The degree is taken w.r.t. the variables indicated in v.
5289  int i;      // Counter
5290  int d = 0;  // The degree
5291  int e;      // Degree (auxiliar variable)
5292  for(i = 1; i <= size(p); i++)
5293  {
5294    e = sum(leadexp(p[i]), @v);
5295    if(e > d){d = e;}
5296  }
5297  return(d);
5298}
5299
5300///////////////////////////////////////////////////////////////////////////////
5301
5302static proc mapBackIdeal(ideal I, poly c, intvec @v)
5303{
5304   ASSUME(1, not isQuotientRing(basering) );
5305
5306// Modifies all polynomials in I so that a map x(i) -> y(i)/c can be
5307// carried out.
5308
5309// v indicates wicih variables x(i) of the ring will be mapped to y(i)/c.
5310
5311  int i;  // counter
5312  for(i = 1; i <= ncols(I); i++)
5313  {
5314    I[i] = mapBackPoly(I[i], c, @v);
5315  }
5316  return(I);
5317}
5318
5319///////////////////////////////////////////////////////////////////////////////
5320
5321static proc mapBackPoly(poly p, poly c, intvec @v)
5322{
5323  ASSUME(1, not isQuotientRing(basering) );
5324
5325// Multiplies each monomial of p by a power of c so that a map x(i) -> y(i)/c
5326// can be carried out.
5327
5328// v indicates wicih variables x(i) of the ring will be mapped to y(i)/c.
5329  int i;  // counter
5330  int e;  // exponent
5331  int d = degSubring(p, @v);
5332  poly g = 0;
5333  int size_p=size(p);
5334  for(i = 1; i <= size_p; i++)
5335  {
5336    e = sum(leadexp(p[i]), @v);
5337    g = g + p[i] * c^(d-e);
5338  }
5339  return(g);
5340}
5341
5342//                    End procedures for normal
5343///////////////////////////////////////////////////////////////////////////////
5344
5345
5346///////////////////////////////////////////////////////////////////////////////
5347//                  Begin procedures for normalC
5348
5349// We first define resp. copy some attributes to be used in proc normal and
5350// static proc normalizationPrimes, and ..., to speed up computation in
5351// special cases
5352//NOTE:  We use the following attributes:
5353// 1     attrib(id,"isCohenMacaulay");         //--- Cohen Macaulay
5354// 2     attrib(id,"isCompleteIntersection");  //--- complete intersection
5355// 3     attrib(id,"isHypersurface");          //--- hypersurface
5356// 4     attrib(id,"isEquidimensional");       //--- equidimensional ideal
5357// 5     attrib(id,"isPrim");                  //--- prime ideal
5358// 6     attrib(id,"isRegInCodim2");           //--- regular in codimension 2
5359// 7     attrib(id,"isIsolatedSingularity");   //--- isolated singularities
5360// 8     attrib(id,"onlySingularAtZero");      //--- only singular at 0
5361// 9     attrib(id,"isRadical");               //--- radical ideal
5362//Recall: (attrib(id,"xy"),1) sets attrib xy to TRUE and
5363//        (attrib(id,"xy"),0) to FALSE
5364
5365static proc getAttrib (ideal id)
5366"USAGE:   getAttrib(id);  id=ideal
5367COMPUTE: check attributes for id. If the attributes above are defined,
5368         take its value, otherwise define it and set it to 0
5369RETURN:  intvec of size 9, with entries 0 or 1,  values of attributes defined
5370         above (in this order)
5371EXAMPLE: no example
5372"
5373{
5374  int isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad;
5375
5376  if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
5377  {
5378    if( attrib(id,"isCohenMacaulay")==1 )
5379    { isCoM=1; isEq=1; }
5380  }
5381
5382  if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
5383  {
5384    if(attrib(id,"isCompleteIntersection")==1)
5385    { isCoI=1; isCoM=1; isEq=1; }
5386  }
5387
5388  if( typeof(attrib(id,"isHypersurface"))=="int" )
5389  {
5390    if(attrib(id,"isHypersurface")==1)
5391    { isHy=1; isCoI=1; isCoM=1; isEq=1; }
5392  }
5393
5394  if( typeof(attrib(id,"isEquidimensional"))=="int" )
5395  {
5396    if(attrib(id,"isEquidimensional")==1)
5397    { isEq=1; }
5398  }
5399
5400  if( typeof(attrib(id,"isPrim"))=="int" )
5401  {
5402    if(attrib(id,"isPrim")==1)
5403    { isPr=1; }
5404  }
5405
5406  if( typeof(attrib(id,"isRegInCodim2"))=="int" )
5407  {
5408    if(attrib(id,"isRegInCodim2")==1)
5409    { isReg=1; }
5410  }
5411
5412  if( typeof(attrib(id,"isIsolatedSingularity"))=="int" )
5413  {
5414    if(attrib(id,"isIsolatedSingularity")==1)
5415    { isIso=1; }
5416  }
5417
5418  if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
5419  {
5420    if(attrib(id,"onlySingularAtZero")==1)
5421    { oSAZ=1; }
5422  }
5423
5424  if( typeof(attrib(id,"isRad"))=="int" )
5425  {
5426    if(attrib(id,"isRad")==1)
5427    { isRad=1; }
5428  }
5429
5430  intvec atr = isCoM,isCoI,isHy,isEq,isPr,isReg,isIso,oSAZ,isRad;
5431  return(atr);
5432}
5433
5434///////////////////////////////////////////////////////////////////////////////
5435
5436static proc setAttrib (ideal id, intvec atr)
5437"USAGE:   setAttrib(id,atr);  id ideal, atr intvec
5438COMPUTE: set attributes to id specified by atr
5439RETURN:  id, with assigned attributes from atr
5440EXAMPLE: no example
5441"
5442{
5443  attrib(id,"isCohenMacaulay",atr[1]);         //--- Cohen Macaulay
5444  attrib(id,"isCompleteIntersection",atr[2]);  //--- complete intersection
5445  attrib(id,"isHypersurface",atr[3]);          //--- hypersurface
5446  attrib(id,"isEquidimensional",atr[4]);       //--- equidimensional ideal
5447  attrib(id,"isPrim",atr[5]);                  //--- prime ideal
5448  attrib(id,"isRegInCodim2",atr[6]);           //--- regular in codimension 2
5449  attrib(id,"isIsolatedSingularity",atr[7]);   //--- isolated singularities
5450  attrib(id,"onlySingularAtZero",atr[8]);      //--- only singular at 0
5451  attrib(id,"isRadical",atr[9]);               //--- radical ideal
5452
5453  return(id);
5454}
5455
5456///////////////////////////////////////////////////////////////////////////////
5457// copyAttribs is not used anywhere so far
5458
5459static proc copyAttribs (ideal id1, ideal id)
5460"USAGE:   copyAttribs(id1,id);  id1, id ideals
5461COMPUTE: copy attributes from id1 to id
5462RETURN:  id, with assigned attributes from id1
5463EXAMPLE: no example
5464"
5465{
5466  if( typeof(attrib(id1,"isCohenMacaulay"))=="int" )
5467  {
5468    if( attrib(id1,"isCohenMacaulay")==1 )
5469    {
5470      attrib(id,"isEquidimensional",1);
5471    }
5472  }
5473  else
5474  {
5475    attrib(id,"isCohenMacaulay",0);
5476  }
5477
5478  if( typeof(attrib(id1,"isCompleteIntersection"))=="int" )
5479  {
5480    if(attrib(id1,"isCompleteIntersection")==1)
5481    {
5482      attrib(id,"isCohenMacaulay",1);
5483      attrib(id,"isEquidimensional",1);
5484    }
5485  }
5486  else
5487  {
5488    attrib(id,"isCompleteIntersection",0);
5489  }
5490
5491  if( typeof(attrib(id1,"isHypersurface"))=="int" )
5492  {
5493    if(attrib(id1,"isHypersurface")==1)
5494    {
5495      attrib(id,"isCompleteIntersection",1);
5496      attrib(id,"isCohenMacaulay",1);
5497      attrib(id,"isEquidimensional",1);
5498    }
5499  }
5500  else
5501  {
5502    attrib(id,"isHypersurface",0);
5503  }
5504
5505  if( (typeof(attrib(id1,"isEquidimensional"))=="int") )
5506  {
5507    if(attrib(id1,"isEquidimensional")==1)
5508    {
5509      attrib(id,"isEquidimensional",1);
5510    }
5511  }
5512  else
5513  {
5514    attrib(id,"isEquidimensional",0);
5515  }
5516
5517  if( typeof(attrib(id1,"isPrim"))=="int" )
5518  {
5519    if(attrib(id1,"isPrim")==1)
5520    {
5521      attrib(id,"isEquidimensional",1);
5522    }
5523  }
5524  else
5525  {
5526    attrib(id,"isPrim",0);
5527  }
5528
5529  if( (typeof(attrib(id1,"isRegInCodim2"))=="int") )
5530  {
5531    if(attrib(id1,"isRegInCodim2")==1)
5532    {
5533      attrib(id,"isRegInCodim2",1);
5534    }
5535  }
5536  else
5537  {
5538    attrib(id,"isRegInCodim2",0);
5539  }
5540
5541  if( (typeof(attrib(id1,"isIsolatedSingularity"))=="int") )
5542  {
5543    if(attrib(id1,"isIsolatedSingularity")==1)
5544    {
5545      attrib(id,"isIsolatedSingularity",1);
5546    }
5547  }
5548  else
5549  {
5550    attrib(id,"isIsolatedSingularity",0);
5551  }
5552
5553  if( typeof(attrib(id1,"onlySingularAtZero"))=="int" )
5554  {
5555    if(attrib(id1,"onlySingularAtZero")==1)
5556    {
5557      attrib(id,"isIsolatedSingularity",1);
5558    }
5559  }
5560  else
5561  {
5562    attrib(id,"onlySingularAtZero",0);
5563  }
5564
5565  if( typeof(attrib(id1,"isRad"))=="int" )
5566  {
5567    if(attrib(id1,"isRad")==1)
5568    {
5569      attrib(id,"isRad",1);
5570    }
5571  }
5572  else
5573  {
5574    attrib(id,"isRad",0);
5575  }
5576  return(id);
5577}
5578///////////////////////////////////////////////////////////////////////////////
5579
5580proc normalC(ideal id, list #)
5581"USAGE:  normalC(id [,choose]);  id = radical ideal, choose = optional list
5582         of string.
5583         Optional parameters in list choose (can be entered in any order):@*
5584         Decomposition:@*
5585         - \"equidim\" -> computes first an equidimensional decomposition,
5586         and then the normalization of each component (default).@*
5587         - \"prim\" -> computes first the minimal associated primes, and then
5588         the normalization of each prime. @*
5589         - \"noDeco\" -> no preliminary decomposition is done. If the ideal is
5590         not equidimensional radical, output might be wrong.@*
5591         - \"isPrim\" -> assumes that the ideal is prime. If the assumption does
5592         not hold, output might be wrong.@*
5593         - \"noFac\" -> factorization is avoided in the computation of the
5594         minimal associated primes;
5595         Other:@*
5596         - \"withGens\" -> the minimal associated primes P_i of id are
5597         computed and for each P_i, algebra generators of the integral closure
5598         of basering/P_i are computed as elements of its quotient field;@*
5599         If choose is not given or empty, the default options are used.@*
5600ASSUME:  The ideal must be radical, for non-radical ideals the output may
5601         be wrong (id=radical(id); makes id radical). However, if option
5602         \"prim\" is set the minimal associated primes are computed first
5603         and hence normalC computes the normalization of the radical of id.
5604         \"isPrim\" should only be used if id is known to be irreducible.
5605RETURN:  a list, say nor, of size 2 (resp. 3 if option \"withGens\" is set).@*
5606         * nor[1] is always a of r rings, where r is the number of associated
5607         primes with option \"prim\" (resp. >= no of equidimenensional
5608         components with option  \"equidim\").@*
5609         Each ring Ri=nor[1][i], i=1..r, contains two ideals with given
5610         names @code{norid} and @code{normap} such that @*
5611         - Ri/norid is the normalization of the i-th component, i.e. the
5612          integral closure in its field of fractions as affine ring, i.e. Ri is
5613          given in the form K[X(1..p),T(1..q)], where K is the ground field;
5614         - normap gives the normalization map from basering/id to
5615           Ri/norid for each i (the j-th element of normap is mapped to the
5616           j-th variable of R).@*
5617         - the direct sum of the rings Ri/norid is the normalization
5618           of basering/id; @*
5619         ** If option \"withGens\" is not set: @*
5620         * nor[2] shows the delta invariants: nor[2] is a list of an intvec of
5621         size r, the delta invariants of the r components, and an integer, the
5622         delta invariant of basering/id. (-1 means infinite, 0 that basering/P_i
5623         resp. basering/input is normal, -2 means that delta resp. delta of one
5624         of the components is not computed (which may happen if \"equidim\" is
5625         given). @*
5626         ** If option \"withGens\" is set:
5627         * nor[2] is a list of ideals Ii=nor[2][i], i=1..r, in the basering,
5628         generating the integral closure of basering/P_i in its quotient field
5629         as K-algebra (K the ground field):@*
5630         If Ii is given by polynomials g_1,...,g_k, then c:=g_k is a non-zero
5631         divisor and the j-th variables of the ring Ri satisfies var(j)=g_j/c,
5632         j=1..k-1, as element in the quotient field of basering/P_i. The
5633         g_j/g_k+1 are K-algebra generators  of the integral closure of
5634         basering/P_i.@*
5635         * nor[3] shows the delta invariant as above.
5636THEORY:  We use the Grauert-Remmert-de Jong algorithm [c.f. G.-M. Greuel,
5637         G. Pfister: A SINGULAR Introduction to Commutative Algebra, 2nd Edition.
5638         Springer Verlag (2007)].
5639         The procedure computes the algebra structure and the delta invariant of
5640         the normalization of R/id:@*
5641         The normalization is an affine algebra over the ground field K
5642         and nor[1] presents it as such: Ri = K[X(1..p),T(1..q)] and Ri/norid
5643         is the integral closure of R/P_i; if option \"withGens\" is set the
5644         X(j) and T(j) are expressed as quotients in the total ring of
5645         fractions. Note that the X(j) and T(j) generate the integral closure
5646         as K-algebra, but not necessarily as R-module (since relations of the
5647         form X(1)=T(1)*T(2) may have been eliminated). Geometrically the
5648         algebra structure is relevant since the variety of the ideal norid in
5649         Ri is the normalization of the variety of the ideal P_i in R.@*
5650         The delta invariant of a reduced ring A is dim_K(normalization(A)/A).
5651         For A=K[x1,...,xn]/id we call this number also the delta invariant of
5652         id. nor[3] returns the delta invariants of the components P_i and of
5653         id.
5654NOTE:    To use the i-th ring type: @code{def R=nor[1][i]; setring R;}.
5655@*       Increasing/decreasing printlevel displays more/less comments
5656         (default: printlevel=0).
5657@*       Not implemented for local or mixed orderings or quotient rings.
5658         For local or mixed orderings use proc 'normal'.
5659@*       If the input ideal id is weighted homogeneous a weighted ordering may
5660         be used (qhweight(id); computes weights).
5661KEYWORDS: normalization; integral closure; delta invariant.
5662SEE ALSO: normal, normalP.
5663EXAMPLE: example normalC; shows an example
5664"
5665{
5666   ASSUME(0, not isQuotientRing(basering) );
5667
5668   int i,j;
5669   int withGens, withEqui, withPrim, isPrim, noFac;
5670   int dbg = printlevel-voice+2;
5671   int nvar = nvars(basering);
5672   int chara  = char(basering);
5673   list result, prim, keepresult;
5674
5675  int decomp;   // Preliminar decomposition:
5676                // 0 -> no decomposition (id is assumed to be prime)
5677                // 1 -> no decomposition
5678                //      (id is assumed to be equidimensional radical)
5679                // 2 -> equidimensional decomposition
5680                // 3 -> minimal associated primes
5681
5682   // Default methods:
5683   noFac = 0;         // Use facstd when computing minimal associated primes
5684   decomp = 2;        // Equidimensional decomposition for nvar > 2
5685   if (nvar <= 2)
5686   { decomp = 3; }    // Compute minimal associated primes if nvar <= 2
5687
5688   if ( attrib(basering,"global") != 1 )
5689   {
5690     "";
5691     "// Not implemented for this ordering,";
5692     "// please change to global ordering or use proc normal";
5693     return(result);
5694   }
5695
5696//--------------------------- define the method ---------------------------
5697   string method;                //make all options one string in order to use
5698                                 //all combinations of options simultaneously
5699   for ( i=1; i <= size(#); i++ )
5700   {
5701     if ( typeof(#[i]) == "string" )
5702     {
5703       method = method + #[i];
5704     }
5705   }
5706
5707   //--------------------------- choosen methods -----------------------
5708   // "withGens": computes algebra generators for each irreducible component
5709   // ### the extra code for withGens should be incorporated in the general case
5710
5711   if ( find(method,"withgens") or find(method,"withGens"))
5712   {
5713     withGens = 1;
5714   }
5715
5716   // the general case: either equidim or minAssGTZ or no decomposition
5717
5718   if ( find(method,"isprim") or find(method,"isPrim") )
5719   {decomp = 0; isPrim=1;}
5720
5721   if ( find(method,"nodeco") or find(method,"noDeco") )
5722   {decomp = 1;}
5723
5724   if ( find(method,"equidim") )
5725   { decomp = 2; }
5726
5727   if ( find(method,"prim") )
5728   { decomp = 3; }
5729
5730   if ( find(method,"nofac") or find(method,"noFac") )
5731   { noFac = 1; }
5732
5733   kill #;
5734   list #;
5735
5736//------- Special algorithm with computation of the generators, RETURN -------
5737   //--------------------- method "withGens" ----------------------------------
5738   //the integral closure is computed in proc primeClosure. In the general case
5739   //it is computed in normalizationPrimes. The main difference is that in
5740   //primeClosure the singular locus is only computed in the first iteration,
5741   //that no attributes are used, and that the generators are computed.
5742   //In primeClosure the (algebra) generators for each irreducible component
5743   //are computed in the static proc closureGenerators
5744
5745   if( withGens )
5746   {
5747      if( dbg >= 1 )
5748      {  "";
5749         "// We use method 'withGens'";
5750      }
5751      if ( decomp == 0 or decomp == 1 )
5752      {
5753         prim[1] = id;
5754         if( dbg >= 0 )
5755         {
5756           "";
5757           "// ** WARNING: result is correct if ideal is prime (not checked) **";
5758           "// if procedure is called with string \"prim\", primality is checked";
5759         }
5760      }
5761      else
5762      {
5763         if(dbg >= 1)
5764         {  "// Computing minimal associated primes..."; }
5765
5766         if( noFac )
5767         { prim = minAssGTZ(id,1); }
5768         else
5769         { prim = minAssGTZ(id); }
5770
5771         if(dbg >= 2)
5772         {  prim;""; }
5773         if(dbg >= 1)
5774         {
5775            "// number of irreducible components is", size(prim);
5776         }
5777      }
5778   //----------- compute integral closure for every component -------------
5779      int del;
5780      intvec deli;
5781      list Gens,l,resu,Resu;
5782      ideal gens;
5783      def R = basering;
5784      poly gg;
5785
5786      for(i=1; i<=size(prim); i++)
5787      {
5788         if(dbg>=1)
5789         {
5790            ""; pause(); "";
5791            "// Computing normalization of component",i;
5792            "   ---------------------------------------";
5793         }
5794
5795         if( defined(ker) ) { kill ker; }
5796         ideal ker = prim[i];
5797         export(ker);
5798         l = R;
5799         l = primeClosure(l,1);              //here the work is done
5800         // primeClosure is called with list l consisting of the basering
5801         //### ausprobieren ob primeClosure(l,1) schneller als primeClosure(l)
5802         // 1 bedeutet: kuerzester nzd
5803         // l[size(l)] is the delta invariant
5804
5805         if ( l[size(l)] >= 0 && del >= 0 )
5806         {
5807            del = del + l[size(l)];
5808         }
5809         else
5810         { del = -1; }
5811         deli = l[size(l)],deli;
5812
5813         l = l[1..size(l)-1];
5814         resu = list(l[size(l)]) + resu;
5815         gens = closureGenerators(l);         //computes algebra(!) generators
5816
5817         //NOTE: gens[i]/gens[size(gens)] expresses the ith variable of resu[1]
5818         //(the normalization) as fraction of elements of the basering;
5819         //the variables of resu[1] are algebra generators.
5820         //gens[size(gens)] is a non-zero divisor of basering/i
5821
5822         //divide by the greatest common divisor:
5823         gg = gcd( gens[1],gens[size(gens)] );
5824         for(j=2; j<=size(gens)-1; j++)
5825         {
5826            gg=gcd(gg,gens[j]);
5827         }
5828         for(j=1; j<=size(gens); j++)
5829         {
5830            gens[j]=gens[j]/gg;
5831         }
5832         Gens = list(gens) + Gens;
5833
5834/*       ### Da die gens Algebra-Erzeuger sind, ist reduce nach Bestimmung
5835         der Algebra-Variablen T(i) nicht zulaessig!
5836         for(i=1;i<=size(gens)-1;i++)
5837         {
5838            gens[i]= reduce(gens[i],std(gens[size(gens)]));
5839         }
5840         for(i=size(gens)-1; i>=1; i--)
5841         {
5842            if(gens[i]==0)
5843            { gens = delete(gens,i); }
5844         }
5845*/
5846         if( defined(ker) ) { kill ker; }
5847      }
5848
5849      if ( del >= 0 )
5850      {
5851         int mul = iMult(prim);
5852         del = del + mul;
5853      }
5854      else
5855      { del = -1; }
5856      deli = deli[1..size(deli)-1];
5857      Resu = resu,Gens,list(deli,del);
5858      int sr = size(resu);
5859
5860      if ( dbg >= 0 )
5861      {"";
5862"// 'normalC' created a list, say nor, of three lists:
5863// To see the list type
5864      nor;
5865
5866// * nor[1] is a list of",sr,"ring(s)
5867// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
5868     def R1 = nor[1][1]; setring R1;  norid; normap;
5869// For the other rings type first (if R is the name of your original basering)
5870     setring R;
5871// and then continue as for R1.
5872// Ri/norid is the affine algebra of the normalization of the i-th
5873// component R/P_i (where P_i is an associated prime of the input ideal id)
5874// and normap the normalization map from R to Ri/norid.
5875
5876// * nor[2] is a list of",sr,"ideal(s), each ideal nor[2][i] consists of
5877// elements g1..gk of R such that the gj/gk generate the integral
5878// closure of R/P_i as sub-algebra in the quotient field of R/P_i, with
5879// gj/gk being mapped by normap to the j-th variable of Ri;
5880
5881// * nor[3] shows the delta-invariant of each component and of id
5882// (-1 means infinite, and 0 that R/P_i resp. R/id is normal).";
5883      }
5884      return(Resu);
5885   }
5886   //----------------- end method "withGens" --------------------------------
5887
5888//-------- The general case without computation of the generators -----------
5889// (attrib(id,"xy"),1) sets attrib xy to TRUE and (attrib(id,"xy"),0) to FALSE
5890// We use the following attributes:
5891//   attrib(id,"isCohenMacaulay");         //--- Cohen Macaulay
5892//   attrib(id,"isCompleteIntersection");  //--- complete intersection
5893//   attrib(id,"isHypersurface");          //--- hypersurface
5894//   attrib(id,"isEquidimensional",-1);    //--- equidimensional ideal
5895//   attrib(id,"isPrim");                  //--- prime ideal
5896//   attrib(id,"isRegInCodim2");           //--- regular in codimension 2
5897//   attrib(id,"isIsolatedSingularity";    //--- isolated singularities
5898//   attrib(id,"onlySingularAtZero");      //--- only singular at 0
5899
5900 //------------------- first set the attributes ----------------------
5901   if( typeof(attrib(id,"isCohenMacaulay"))=="int" )
5902   {
5903      if( attrib(id,"isCohenMacaulay")==1 )
5904      {
5905         attrib(id,"isEquidimensional",1);
5906      }
5907   }
5908   else
5909   {
5910      attrib(id,"isCohenMacaulay",0);
5911   }
5912
5913   if( typeof(attrib(id,"isCompleteIntersection"))=="int" )
5914   {
5915      if(attrib(id,"isCompleteIntersection")==1)
5916      {
5917         attrib(id,"isCohenMacaulay",1);
5918         attrib(id,"isEquidimensional",1);
5919      }
5920   }
5921   else
5922   {
5923      attrib(id,"isCompleteIntersection",0);
5924   }
5925
5926   if( typeof(attrib(id,"isHypersurface"))=="int" )
5927   {
5928      if(attrib(id,"isHypersurface")==1)
5929      {
5930         attrib(id,"isCompleteIntersection",1);
5931         attrib(id,"isCohenMacaulay",1);
5932         attrib(id,"isEquidimensional",1);
5933      }
5934   }
5935   else
5936   {
5937      attrib(id,"isHypersurface",0);
5938   }
5939
5940   if( ! (typeof(attrib(id,"isEquidimensional"))=="int") )
5941   {
5942         attrib(id,"isEquidimensional",0);
5943   }
5944
5945   if( typeof(attrib(id,"isPrim"))=="int" )
5946   {
5947      if(attrib(id,"isPrim")==1)
5948      {
5949         attrib(id,"isEquidimensional",1);
5950      }
5951   }
5952   else
5953   {
5954      attrib(id,"isPrim",0);
5955   }
5956
5957   if( ! (typeof(attrib(id,"isRegInCodim2"))=="int") )
5958   {
5959         attrib(id,"isRegInCodim2",0);
5960   }
5961
5962   if( ! (typeof(attrib(id,"isIsolatedSingularity"))=="int") )
5963   {
5964         attrib(id,"isIsolatedSingularity",0);
5965   }
5966
5967   if( typeof(attrib(id,"onlySingularAtZero"))=="int" )
5968   {
5969      if(attrib(id,"onlySingularAtZero")==1)
5970      {
5971         attrib(id,"isIsolatedSingularity",1);
5972      }
5973   }
5974   else
5975   {
5976      attrib(id,"onlySingularAtZero",0);
5977   }
5978
5979   //-------------- compute equidimensional decomposition --------------------
5980   //If the method "equidim" is given, compute the equidim decomposition
5981   //and goto the next step (no normalization
5982   //ACHTUNG: equidim berechnet bei nicht reduzierten id die eingebetteten
5983   //Komponenten als niederdim Komponenten, waehrend diese bei primdecGTZ
5984   //nicht auftauchen: ideal(x,y)*xy
5985   //this is default for nvars > 2
5986
5987   if( decomp == 2 )
5988   {
5989      withPrim = 0;                 //this is used to check later that prim
5990                                    //contains equidim but not prime components
5991      if( dbg >= 1 )
5992      {
5993         "// We use method 'equidim'";
5994      }
5995      if( typeof(attrib(id,"isEquidimensional"))=="int" )
5996      {
5997         if(attrib(id,"isEquidimensional")==1)
5998         {
5999            prim[1] = id;
6000         }
6001         else
6002         {
6003            prim = equidim(id);
6004         }
6005      }
6006      else
6007      {
6008         prim = equidim(id);
6009      }
6010      if(dbg>=1)
6011      {  "";
6012         "// number of equidimensional components:", size(prim);
6013      }
6014      if ( !noFac )
6015      {
6016        intvec opt = option(get);
6017        option(redSB);
6018        for(j=1; j<=size(prim); j++)
6019        {
6020           keepresult = keepresult+facstd(prim[j]);
6021        }
6022        prim = keepresult;
6023        if ( size(prim) == 0 )
6024        {
6025          prim=ideal(0);     //Bug in facstd, liefert leere Liste bei 0-Ideal
6026        }
6027
6028        if(dbg>=1)
6029        {  "";
6030         "// number of components after application of facstd:", size(prim);
6031        }
6032        option(set,opt);
6033      }
6034   }
6035
6036   //------------------- compute associated primes -------------------------
6037   //the case where withEqui = 0, here the min. ass. primes are computed
6038   //start with the computation of the minimal associated primes:
6039
6040   else
6041   {
6042    if( isPrim )
6043    {
6044      if( dbg >= 0 )
6045      {
6046         "// ** WARNING: result is correct if ideal is prime";
6047         "// or equidimensional (not checked) **";
6048         "// disable option \"isPrim\" to decompose ideal into prime";
6049         "// or equidimensional components";"";
6050      }
6051      if( dbg >= 1 )
6052      {
6053        "// We use method 'isPrim'";"";
6054      }
6055      prim[1]=id;
6056    }
6057    else
6058    {
6059      withPrim = 1;                 //this is used to check later that prim
6060                                    //contains prime but not equidim components
6061      if( dbg >= 1 )
6062      {
6063         "// We use method 'prim'";
6064      }
6065
6066      if( typeof(attrib(id,"isPrim"))=="int" )
6067      {
6068         if(attrib(id,"isPrim")==1)
6069         {
6070            prim[1]=id;
6071         }
6072         else
6073         {
6074            if( noFac )
6075            { prim=minAssGTZ(id,1); }     //does not use factorizing groebner
6076            else
6077            { prim=minAssGTZ(id); }       //uses factorizing groebner
6078         }
6079      }
6080      else
6081      {
6082            if( noFac )
6083            { prim=minAssGTZ(id,1); }
6084            else
6085            { prim=minAssGTZ(id); }
6086      }
6087      if(dbg>=1)
6088      {  "";
6089         "// number of irreducible components:", size(prim);
6090      }
6091    }
6092   }
6093
6094   //----- for each component (equidim or irred) compute normalization -----
6095   int sr, skr, del;
6096   intvec deli;
6097   int sp = size(prim);     //size of list prim (# irred or equidim comp)
6098
6099   for(i=1; i<=sp; i++)
6100   {
6101      if(dbg>=1)
6102      {  "";
6103         "// computing the normalization of component",i;
6104         "   ----------------------------------------";
6105      }
6106      //-------------- first set attributes for components ------------------
6107      attrib(prim[i],"isEquidimensional",1);
6108      if( withPrim )
6109      {
6110         attrib(prim[i],"isPrim",1);
6111      }
6112      else
6113      { attrib(prim[i],"isPrim",0); }
6114
6115      if(attrib(id,"onlySingularAtZero")==1)
6116      { attrib(prim[i],"onlySingularAtZero",1); }
6117      else
6118      { attrib(prim[i],"onlySingularAtZero",0); }
6119
6120      if(attrib(id,"isIsolatedSingularity")==1)
6121      { attrib(prim[i],"isIsolatedSingularity",1); }
6122      else
6123      { attrib(prim[i],"isIsolatedSingularity",0); }
6124
6125      if( attrib(id,"isHypersurface")==1 )
6126      {
6127         attrib(prim[i],"isHypersurface",1);
6128         attrib(prim[i],"isCompleteIntersection",1);
6129         attrib(prim[i],"isCohenMacaulay",1);
6130      }
6131      else
6132      { attrib(prim[i],"isHypersurface",0); }
6133
6134      if ( sp == 1)         //the case of one component: copy attribs from id
6135      {
6136        if(attrib(id,"isRegInCodim2")==1)
6137        {attrib(prim[i],"isRegInCodim2",1); }
6138        else
6139        {attrib(prim[i],"isRegInCodim2",0); }
6140
6141        if(attrib(id,"isCohenMacaulay")==1)
6142        {attrib(prim[i],"isCohenMacaulay",1); }
6143        else
6144        {attrib(prim[i],"isCohenMacaulay",0); }
6145
6146        if(attrib(id,"isCompleteIntersection")==1)
6147        {attrib(prim[i],"isCompleteIntersection",1); }
6148        else
6149        {attrib(prim[i],"isCompleteIntersection",0); }
6150      }
6151      else
6152      {
6153        attrib(prim[i],"isRegInCodim2",0);
6154        attrib(prim[i],"isCohenMacaulay",0);
6155        attrib(prim[i],"isCompleteIntersection",0);
6156      }
6157
6158      //------ Now compute the normalization of each component ---------
6159      //note: for equidimensional components the "splitting tools" can
6160      //create further decomposition
6161      //We now start normalizationPrimes with
6162      //ihp = partial normalisation map = identity map = maxideal(1)
6163      //del = partial delta invariant = 0
6164      //deli= intvec of partial delta invariants of components
6165      //in normalizationPrimes all the work is done:
6166
6167      keepresult = normalizationPrimes(prim[i],maxideal(1),0,0);
6168
6169      for(j=1; j<=size(keepresult)-1; j++)
6170      {
6171         result=insert(result,keepresult[j]);
6172      }
6173      skr = size(keepresult);
6174
6175      //compute delta:
6176      if( del >= 0 && keepresult[skr][1] >=0 )
6177      {
6178         del = del + keepresult[skr][1];
6179      }
6180      else
6181      {
6182         del = -1;
6183      }
6184      deli = keepresult[skr][2],deli;
6185
6186      if ( dbg>=1 )
6187      {
6188           "// delta of component",i; keepresult[skr][1];
6189      }
6190   }
6191   sr = size(result);
6192
6193   // -------------- Now compute intersection multiplicities -------------
6194   //intersection multiplicities of list prim, sp=size(prim).
6195      if ( dbg>=1 )
6196      {
6197        "// Sum of delta for all components"; del;
6198        if ( sp>1 )
6199        {
6200           "// Compute intersection multiplicities of the components";
6201        }
6202      }
6203
6204      if ( sp > 1 )
6205      {
6206        int mul = iMult(prim);
6207        if ( mul < 0 )
6208        {
6209           del = -1;
6210        }
6211        else
6212        {
6213           del = del + mul;
6214        }
6215      }
6216   deli = deli[1..size(deli)-1];
6217   result = result,list(deli,del);
6218
6219//--------------- Finally print comments and return ------------------
6220   if ( dbg >= 0)
6221   {"";
6222"// 'normalC' created a list, say nor, of two lists:
6223// To see the result, type
6224      nor;
6225
6226// * nor[1] is a list of",sr,"ring(s).
6227// To access the i-th ring nor[1][i] give it a name, say Ri, and type e.g.
6228      def R1 = nor[1][1];  setring R1;  norid;  normap;
6229// and similair for the other rings nor[1][i];
6230// Ri/norid is the affine algebra of the normalization of r/P_i  (where P_i
6231// is an associated prime or an equidimensional part of the input ideal id)
6232// and normap the normalization map from the basering to Ri/norid;
6233
6234// * nor[2] shows the delta-invariant of each component and of id
6235// (-1 means infinite, 0 that r/P_i resp. r/id is normal, and -2 that delta
6236// of a component was not computed).";
6237   }
6238   return(result);
6239}
6240
6241example
6242{ "EXAMPLE:";
6243   printlevel = printlevel+1;
6244   echo = 2;
6245   ring s = 0,(x,y),dp;
6246   ideal i = (x2-y3)*(x2+y2)*x;
6247
6248   list nor = normalC(i);
6249
6250   nor;
6251   // 2 branches have delta = 1, and 1 branch has delta = 0
6252   // the total delta invariant is 13
6253
6254   def R2 = nor[1][2];  setring R2;
6255   norid; normap;
6256
6257   echo = 0;
6258   printlevel = printlevel-1;
6259   pause("   hit return to continue"); echo=2;
6260
6261   ring r = 2,(x,y,z),dp;
6262   ideal i = z3-xy4;
6263   nor = normalC(i);  nor;
6264   // the delta invariant is infinite
6265   // xy2z/z2 and xy3/z2 generate the integral closure of r/i as r/i-module
6266   // in its quotient field Quot(r/i)
6267
6268   // the normalization as affine algebra over the ground field:
6269   def R = nor[1][1]; setring R;
6270   norid; normap;
6271
6272   echo = 0;
6273   pause("   hit return to continue");echo = 2;
6274
6275   setring r;
6276   nor = normalC(i, "withGens", "prim");    // a different algorithm
6277   nor;
6278}
6279
6280//////////////////////////////////////////////////////////////////////////////
6281//closureRingtower seems not to be used anywhere
6282static proc closureRingtower(list L)
6283"USAGE:    closureRingtower(list L); L a list of rings
6284CREATE:   rings R(1),...,R(n) such that R(i)=L[i] for all i
6285EXAMPLE:  example closureRingtower; shows an example
6286"
6287{
6288  int n=size(L);
6289  for (int i=1;i<=n;i++)
6290    {
6291      if (defined(R(i)))
6292      {
6293        string s="Fixed name R("+string(i)+") leads to conflict with existing "
6294              +"object having this name";
6295        ERROR(s);
6296      }
6297      def R(i)=L[i];
6298      export R(i);
6299    }
6300
6301  return();
6302}
6303example
6304{
6305  "EXAMPLE:"; echo=2;
6306  ring R=0,(x,y),dp;
6307  ideal I=x4,y4;
6308  list L=primeClosure(ReesAlgebra(I)[1]);
6309  L=delete(L,size(L));
6310  L;
6311  closureRingtower(L);
6312  R(1);
6313  R(4);
6314  kill R(1),R(2),R(3),R(4);
6315}
6316
6317//                Up to here: procedures for normalC
6318///////////////////////////////////////////////////////////////////////////////
6319
6320///////////////////////////////////////////////////////////////////////////////
6321//                From here: miscellaneous procedures
6322
6323// Used for timing and comparing the different normalization procedures.
6324// Option (can be entered in any order)
6325// "normal"   -> uses the new algortihm (normal)
6326// "normalP"  -> uses normalP
6327// "normalC"  -> uses normalC, without "withGens" option
6328// "primCl"   -> uses normalC, with option "withGens".
6329// "111"      -> checks the output of normalM using norTest.
6330// "p"        -> compares the output of norM with the output of normalP
6331//               ("normalP" option must also be set).
6332// "pc"       -> compares the output of norM with the output of normalC with
6333//               option "withGens"
6334//               ("primCl" option must also be set).
6335
6336proc timeNormal(ideal I, list #)
6337{
6338  ASSUME(0, not isQuotientRing(basering) );
6339
6340  def r = basering;
6341
6342  //--------------------------- define the method ---------------------------
6343  int isPrim, useRing;
6344  int decomp = -1;
6345  int norM, norC, norP, primCl;
6346  int checkP, check111, checkPC;
6347  int i;
6348  ideal U1, U2, W;
6349  poly c1, c2;
6350  int ch;
6351  string check;
6352  string method;                //make all options one string in order to use
6353                                //all combinations of options simultaneously
6354  for ( i=1; i <= size(#); i++ )
6355  {
6356    if ( typeof(#[i]) == "string" )
6357    {
6358      method = method + #[i];
6359    }
6360  }
6361  if ( find(method, "normal"))
6362  {norM = 1;}
6363  if ( find(method, "normalP") and (char(basering) > 0))
6364  {norP = 1;}
6365  if ( find(method, "normalC"))
6366  {norC = 1;}
6367  if ( find(method, "primCl"))
6368  {primCl = 1;}
6369  if ( find(method, "isprim") or find(method,"isPrim") )
6370  {decomp = 0;}
6371  if ( find(method, "p") )
6372  {checkP = 1;}
6373  if ( find(method, "pc") )
6374  {checkPC = 1;}
6375  if ( find(method, "111") )
6376  {check111 = 1;}
6377
6378  int tt;
6379  if(norM)
6380  {
6381    tt = timer;
6382    if(decomp == 0)
6383    {
6384      "Running normal(useRing, isPrim)...";
6385      list a1 = normal(I, "useRing", "isPrim");
6386      "Time normal(useRing, isPrim): ", timer - tt;
6387    }
6388    else
6389    {
6390      "Running normal(useRing)...";
6391      list a1 = normal(I, "useRing");
6392      "Time normal(useRing): ", timer - tt;
6393    }
6394    "";
6395  }
6396  if(norP)
6397  {
6398    tt = timer;
6399    if(decomp == 0)
6400    {
6401      "Running normalP(isPrim)...";
6402      list a2 = normalP(I, "isPrim");
6403      "Time normalP(isPrim): ", timer - tt;
6404    }
6405    else
6406    {
6407      "Running normalP()...";
6408      list a2 = normalP(I);
6409      "Time normalP(): ", timer - tt;
6410    }
6411    "";
6412  }
6413
6414  if(norC)
6415  {
6416    tt = timer;
6417    if(decomp == 0)
6418    {
6419      "Running normalC(isPrim)...";
6420      list a3 = normalC(I, "isPrim");
6421      "Time normalC(isPrim): ", timer - tt;
6422    }
6423    else
6424    {
6425      "Running normalC()...";
6426      list a3 = normalC(I);
6427      "Time normalC(): ", timer - tt;
6428    }
6429    "";
6430  }
6431
6432  if(primCl)
6433  {
6434    tt = timer;
6435    if(decomp == 0)
6436    {
6437      "Running normalC(withGens, isPrim)...";
6438      list a4 = normalC(I, "isPrim", "withGens");
6439      "Time normalC(withGens, isPrim): ", timer - tt;
6440    }
6441    else
6442    {
6443      "Running normalC(withGens)...";
6444      list a4 = normalC(I, "withGens");
6445      "Time normalC(withGens): ", timer - tt;
6446    }
6447    "";
6448  }
6449
6450  if(check111 and norM)
6451  {
6452    "Checking output with norTest...";
6453    "WARNING: this checking only works if the original ideal was prime.";
6454    norTest(I, a1);
6455    "";
6456  }
6457
6458  if(checkP and norP and norM)
6459  {
6460    "Comparing with normalP output...";
6461    if(size(a2) > 0)
6462    {
6463      "WARNING: this checking only works if the original ideal was prime.";
6464      U1 = a1[2][1];
6465      c1 = U1[size(U1)];
6466      U2 = a2[1][1];
6467      c2 = a2[1][1][size(a2[1][1])];
6468      W = changeDenominator(U1, c1, c2, groebner(I));
6469      qring q = groebner(I);
6470      ideal U2 = fetch(r, U2);
6471      ideal W = fetch(r, W);
6472      ch = 0;
6473      if(size(reduce(U2, groebner(W))) == 0)
6474      {
6475        "U2 c U1";
6476        ch = 1;
6477      }
6478      if(size(reduce(W, groebner(U2))) == 0)
6479      {
6480        "U1 c U2";
6481        ch = ch + 1;
6482      }
6483      if(ch == 2)
6484      {
6485        "Output of normalP is equal.";
6486      }
6487      else
6488      {
6489        "ERROR: Output of normalP is different.";
6490      }
6491      setring r;
6492      kill q;
6493    }
6494    else
6495    {
6496      "normalP returned no output. Comparison is not possible.";
6497    }
6498    "";
6499  }
6500
6501  if(checkPC and norM and primCl)
6502  {
6503    "Comparing with primeClosure output...";
6504    if(size(a4) > 0)
6505    {
6506      "WARNING: this checking only works if the original ideal was prime.";
6507      // primeClosure check
6508      U1 = a1[2][1];
6509      c1 = U1[size(U1)];
6510      U2 = a4[2][1];
6511      c2 = a4[2][1][size(a4[2][1])];
6512      W = changeDenominator(U1, c1, c2, groebner(I));
6513      qring q = groebner(I);
6514      ideal U2 = fetch(r, U2);
6515      ideal W = fetch(r, W);
6516      ch = 0;
6517      if(size(reduce(U2, groebner(W))) == 0)
6518      {
6519        "U2 c U1";
6520        ch = 1;
6521      }
6522      if(size(reduce(W, groebner(U2))) == 0)
6523      {
6524        "U1 c U2";
6525        ch = ch + 1;
6526      }
6527      if(ch == 2)
6528      {
6529        "Output of normalC(withGens) is equal.";
6530      }
6531      else
6532      {
6533        "ERROR: Output of normalC(withGens) is different.";
6534      }
6535      setring r;
6536      kill q;
6537    }
6538    else
6539    {
6540      "normalC(withGens) returned no output. Comparison is not possible.";
6541    }
6542    "";
6543  }
6544}
6545
6546///////////////////////////////////////////////////////////////////////////
6547static proc sqroot(int n);
6548{
6549  int s = 1;
6550  while(s*s < n) { s++; }
6551  return(s);
6552}
6553
6554///////////////////////////////////////////////////////////////////////////
6555proc norTest (ideal i, list nor, list #)
6556"USAGE:   norTest(i,nor,[n]); i=prime ideal, nor=list, n=optional integer
6557ASSUME:  nor is the output of normal(i) (any options) or
6558         normalP(i,"withRing") or normalC(i) (any options).
6559         In particular, the ring nor[1][1] contains the ideal norid
6560         and the map normap: basering/i --> nor[1][1]/norid.
6561RETURN:  an intvec v such that:
6562@format
6563         v[1] = 1 if the normap is injective and 0 otherwise
6564         v[2] = 1 if the normap is finite and 0 otherwise
6565         v[3] = 1 if nor[1][1]/norid is normal and 0 otherwise
6566@end format
6567         If n=1 (resp n=2) only v[1] (resp. v[2]) is computed and returned
6568THEORY:  The procedure can be used to test whether the computation of the
6569         normalization was correct: basering/i --> nor[1][1]/norid is the
6570         normalization of basering/i if and only if v=1,1,0.
6571NOTE:    For big examples it can be hard to fully test correctness; the
6572         partial test norTest(i,nor,2) is usually fast
6573EXAMPLE: example norTest; shows an example
6574"
6575{
6576   ASSUME(0, not isQuotientRing(basering) );
6577//### Sollte erweitert werden auf den reduziblen Fall: einen neuen affinen
6578// Ring nor[1][1]+...+nor[1][r] (direkte Summe) erzeugen, map dorthin
6579// definieren und dann testen.
6580
6581    int prl = printlevel - voice + 2;
6582    int a,b,d;
6583    int n,ii;
6584    if (size(#) > 0) {  n = #[1];  }
6585
6586    def BAS = basering;
6587
6588    //### make a copy of nor to have a cpoy of nor[1][1]  (not a reference to)
6589    // in order not to override norid and normap.
6590    // delete nor[2] (if it contains the module generators, which are not used)
6591    // s.t. newnor does not belong to a ring.
6592
6593    list newnor = nor;
6594    if ( size(newnor) == 3 )
6595    {
6596       newnor = delete(newnor,2);
6597    }
6598    def R = newnor[1][1];
6599    qring QAS = std(i);
6600
6601
6602    setring R;
6603    int nva = nvars(R);
6604    string svars = varstr(R);
6605    string svar;
6606
6607    norid = interred(norid);
6608
6609    //--------- create new ring with one dp block keeping weights ------------
6610    list LR = ringlist(R);
6611    list g3 = LR[3];
6612    int n3 = size(g3);
6613    list newg3;
6614    intvec V;
6615
6616    //--------- check first whether variables Z(i),...,A(i) exist -----------
6617    for (ii=90; ii>=65; ii--)
6618    {
6619       if ( find(svars,ASCII(ii)+"(") == 0 )
6620       {
6621          svar = ASCII(ii);  break;
6622       }
6623    }
6624    if ( size(svar) != 0 )
6625    {
6626        for ( ii = 1; ii <= nva; ii++ )
6627        {
6628            LR[2][ii] = svar+"("+string(ii)+")";
6629            V[ii] = 1;
6630        }
6631    }
6632    else
6633    {
6634        for ( ii = 1; ii <= nva; ii++ )
6635        {
6636           LR[2][ii] = "Z("+string(100*nva+ii)+")";
6637           V[ii] = 1;
6638        }
6639    }
6640
6641    if ( g3[n3][1]== "c" or g3[n3][1] == "C" )
6642    {
6643       list gm = g3[n3];       //last blockis module ordering
6644       newg3[1] = list("dp",V);
6645       newg3 = insert(newg3,gm,size(newg3));
6646    }
6647    else
6648    {
6649       list gm = g3[1];              //first block is module ordering
6650       newg3[1] = list("dp",V);
6651       newg3 = insert(newg3,gm);
6652    }
6653    LR[3] = newg3;
6654//LR;"";
6655    def newR = ring(LR);
6656
6657    setring newR;
6658    ideal norid = fetch(R,norid);
6659    ideal normap = fetch(R,normap);
6660    if( defined(lnorid) )  { kill lnorid; }     //um ** redefinig zu beheben
6661    if( defined(snorid) )  { kill snorid; }     //sollte nicht noetig sein
6662
6663    //----------- go to quotient ring for checking injectivity -------------
6664//"mstd";
6665    list lnorid = mstd(norid);
6666    ideal snorid = lnorid[1];
6667//"size mstdnorid:", size(snorid),size(lnorid[2]);
6668//"size string mstdnorid:", size(string(snorid)),size(string(lnorid[2]));
6669    qring QR = snorid;
6670    ideal qnormap = fetch(newR,normap);
6671    //ideal qnormap = imap(newR,normap);
6672    //ideal qnormap = imap(R,normap);
6673    map Qnormap = QAS,qnormap;    //r/id --> R/norid
6674
6675    //------------------------ check injectivity ---------------------------
6676//"injective:";
6677    a = is_injective(Qnormap,QAS);          //a. Test for injectivity of Qnormap
6678    dbprint ( prl, "injective: "+string(a) );
6679    if ( n==1 )
6680    {
6681     intvec result = intvec(a);
6682     setring BAS;
6683     return (result);
6684   }
6685   a;
6686
6687    //------------------------ check finiteness ---------------------------
6688    setring newR;
6689    b = mapIsFinite(normap,BAS,lnorid[2]);  //b. Test for finiteness of normap
6690    dbprint ( prl, "finite: "+string(b) );
6691    if ( n==2 )
6692    {
6693       intvec result = intvec(a,b);
6694       setring BAS;
6695       return (result);
6696    }
6697   b;
6698
6699    //------------------------ check normality ---------------------------
6700    list testnor = normal(lnorid[2],"isPrim","noFac", "withDelta");
6701    //### Problem: bei mehrfachem Aufruf von norTest gibt es
6702    // ** redefining norid & ** redefining normap
6703    //Dies produziert Fehler, da alte norid und normap ueberschrieben werden
6704    //norid und normap werden innnerhalb von proc computeRing ueberschrieben
6705    //Die Kopie newR scheint das Problem zu loesen
6706
6707
6708    d = testnor[3][2];             //d = delta
6709    kill testnor;                              //### sollte ueberfluessig sein
6710    int d1 = (d==0);                           //d1=1 if delta=0
6711    dbprint ( prl, "delta: "+string(d) );
6712    intvec result = intvec(a,b,d1);
6713    setring BAS;
6714    return(result);
6715}
6716example
6717{ "EXAMPLE:"; echo = 2;
6718   int prl = printlevel;
6719   printlevel = -1;
6720   ring r = 0,(x,y),dp;
6721   ideal i = (x-y^2)^2 - y*x^3;
6722   list nor = normal(i);
6723   norTest(i,nor);                //1,1,1 means that normal was correct
6724
6725   nor = normalC(i);
6726   norTest(i,nor);                //1,1,1 means that normal was correct
6727
6728   ring s = 2,(x,y),dp;
6729   ideal i = (x-y^2)^2 - y*x^3;
6730   nor = normalP(i,"withRing");
6731   norTest(i,nor);               //1,1,1 means that normalP was correct
6732   printlevel = prl;
6733}
6734
6735///////////////////////////////////////////////////////////////////////////
6736//
6737//                            EXAMPLES
6738//
6739///////////////////////////////////////////////////////////////////////////
6740/*
6741//commands for computing the normalization:
6742// options for normal:  "equidim", "prim"
6743//                      "noDeco", "isPrim", "noFac"
6744//                       (prim by default)
6745// options for normalP: "withRing", "isPrim" or "noFac"
6746// options for normalC: "equidim", "prim", "withGens",
6747//                      "noDeco", "isPrim", "noFac"
6748
6749//Commands for testing 'normal'
6750 list nor = normal(i); nor;
6751 list nor = normal(i,"isPrim");nor;
6752 list nor = normal(i,"equidim");nor;
6753 list nor = normal(i,"prim");nor;
6754 list nor = normal(i,"equidim","noFac");nor;
6755 list nor = normal(i,"prim","noFac");nor;
6756
6757//Commands for testing 'normalP' in positive char
6758 list nor = normalP(i);nor;              //withGens but no ringstructure
6759 list nor = normalP(i,"withRing"); nor;  //compute the ringstructure
6760 list nor = normalP(i,"isPrim"); nor;    //if i is known to be prime
6761
6762//Commands for testing 'normalC'
6763 list nor = normal(i); nor;
6764 list nor = normal(i,"withGens");nor;
6765 list nor = normal(i,"isPrim");nor;
6766 list nor = normal(i,"equidim");nor;
6767 list nor = normal(i,"prim");nor;
6768 list nor = normal(i,"equidim","noFac");nor;
6769 list nor = normal(i,"prim","noFac");nor;
6770
6771//Commands for testing correctness (i must be prime):
6772list nor = normalP(i,"withRing","isPrim");
6773list nor = normal(i,"isPrim");
6774norTest(i,nor);       //full test for not too big examples (1,1,1 => ok)
6775norTest(i,nor,2);     //partial test for big examples (1,1 => ok)
6776factorize(i[1]);      //checks for irreducibility
6777
6778/////////////////////////////////////////////////////////////////////////////
6779
6780//----------------------Examples for normal (new algorithm)------------------
6781// Timings with Computeserver Dual AMD Opteron 242 1.60GHz.
6782// Examples from "Normalization of Rings" paper.
6783
6784// Example 1
6785// char 0 : normal = 0 secs (7 steps) - normalC = 75 secs
6786// char 2 : normal = 0 secs (7 steps) - normalP = 0 secs - normalC = 0 secs
6787// char 5 : normal = 1 secs (7 steps) - normalP = 71 - normalC = 1 secs
6788// char 11 : normal = 2 secs (7 steps) - normalP = 12 secs - normalC doesn't finish
6789// char 32003 : normal = 1 secs (7 steps) - normalP doesn't finish - normalC = 1 sec
6790LIB"normal.lib";
6791ring r = 2, (x, y), dp;
6792ideal i = (x-y)*x*(y+x^2)^3-y^3*(x^3+x*y-y^2);
6793timeNormal(i, "normal", "normalC", "normalP", "isPrim", "p");
6794
6795// Example 2
6796// char 0  : normal = 1 sec (7 steps) - normalC doesn't finish
6797// char 3 : normal = 1 secs (8 steps) - normalP = 0 secs - normalC = 4 secs
6798// char 13 : normal = 1 sec (7 steps) - normalP doesn't finish - normalC = 13 secs
6799// char 32003 : normal = 1 secs (7 steps) - normalP doesn't finish - normalC = 10 sec
6800//Example is reducible in char 5 and 7
6801LIB"normal.lib";
6802ring r = 3, (x, y), dp;
6803ideal 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;
6804timeNormal(i, "normal", "normalC", "normalP", "p", "isPrim");
6805
6806// Example 3
6807// char 0 : normal = 3 secs (6 steps) - normalC doesn't finish
6808// char 2 : normal = 1 secs (13 steps) - normalP = 0 secs - normalC doesn't finish
6809// char 5 : normal = 0 secs (6 steps) - normalP = 8 secs - normalC doesn't finish
6810LIB"normal.lib";
6811ring r=5,(x, y),dp;
6812ideal i=y9+y8x+y8+y5+y4x+y3x2+y2x3+yx8+x9;
6813timeNormal(i, "normal", "normalC", "normalP", "isPrim");
6814
6815// Example 4
6816// char 0 : normal = 0 secs (1 step) - normalC = 0 secs
6817// char 5 : normal = 0 secs (1 step) - normalP = 3 secs - normalC = 0 secs
6818// char 11 : normal = 0 secs (1 step) - normalP doesn't finish - normalC = 0 secs
6819// char 32003 : normal = 0 secs (1 step) - normalP doesn't finish - normalC = 0 secs
6820LIB"normal.lib";
6821ring r=5,(x,y),dp;   // genus 0 4 nodes and 6 cusps im P2
6822ideal i=(x2+y^2-1)^3 +27x2y2;
6823timeNormal(i, "normal", "normalC", "normalP", "isPrim");
6824
6825// Example 5
6826// char 0 : normal = 0 secs (1 step) - normalC = 0 secs
6827// char 5 : normal = 1 secs (3 step) - normalP doesn't finish - normalC doesn't finish
6828// char 11 : normal = 0 secs (1 step) - normalP 0 secs - normalC = 0 secs
6829// char 32003 : normal = 0 secs (1 step) - normalP doesn't finish - normalC = 0 secs
6830LIB"normal.lib";
6831ring r=11,(x,y),dp;    //24 sing, delta 24
6832ideal i=-x10+x8y2-x6y4-x2y8+2y10-x8+2x6y2+x4y4-x2y6-y8+2x6-x4y2+x2y4+2x4+2x2y2-y4-x2+y2-1;
6833timeNormal(i, "normal", "normalC", "normalP", "isPrim", "p");
6834
6835// Example 6
6836// char 2 : normal = 5 secs (2 steps) - normalP = 25 secs - normalC = 166 secs
6837LIB"normal.lib";
6838ring r=2,(v,u,z,y,x),dp;
6839ideal i = z3+zyx+y3x2+y2x3, uyx+z2,uz+z+y2x+yx2, u2+u+zy+zx, v3+vux+vz2+vzyx+vzx+uz3+uz2y+z3+z2yx2;
6840timeNormal(i, "normal", "normalC", "normalP", "isPrim", "p");
6841
6842// Example 7
6843// char 0 : normal = 11 secs (6 steps) - normalC = 11 secs
6844// char 2 : normal = 11 secs (6 steps) - normalP = 0 secs - normalC = 11 secs
6845// char 5 : normal = 11 secs (6 steps) - normalP = 3 secs - normalC = 11 secs
6846// char 11 : normal = 11 secs (6 steps) - normalP = 43 secs - normalC = 11 secs
6847// char 32003 : normal = 11 secs (6 steps) - normalP doesn't finish - normalC = 11 secs
6848LIB"normal.lib";
6849ring r=11,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
6850ideal i= x2+zw, y3+xwt, xw3+z3t+ywt2, y2w4-xy2z2t-w3t3;
6851timeNormal(i, "normal", "normalC", "normalP", "isPrim");
6852
6853////////////////////////////////////////////////////////////////////////////////
6854
6855// Other examples with new algorithm
6856
6857// Example 1
6858// char 0 : normal = 1 secs (13 steps) - normalC doesn't finish
6859// char 2 : normal = 1 secs (13 steps) - normalP = 0 secs - normalC doesn't finish
6860// char 5 : normal = 1 secs (13 steps) - normalP = 29 secs - normalC doesn't finish
6861ring r=2,(x,y),dp;  //genus 35
6862ideal i=y30+y13x+x4y5+x3*(x+1)^2;
6863timeNormal(i, "normal", "normalC", "normalP");
6864
6865// Example 2
6866// char 0 : normal = 1 secs (13 steps) - normalC doesn't finish
6867// char 3 : normal = 2 secs (13 steps) - normalP = 0 secs - normalC doesn't finish
6868ring r=3,(x,y),dp;  //genus 19, delta 21
6869ideal i=y20+y13x+x4y5+x3*(x+1)^2;
6870timeNormal(i, "normal", "normalC", "normalP");
6871
6872// Example 3
6873// Very fast with all algorithms
6874ring r = 3, (x, y), dp;
6875ideal I = (x-y^2)^2-x*y^3;
6876timeNormal(I, "normal", "normalC", "normalP", "primCl", "111", "p", "pc");
6877
6878
6879
6880//----------------------Test Example for charp -------------------
6881//Zu tun:
6882//### nach minor nur std statt mstd verwenden
6883//***hat bei keinem Beisp etwas gebracht -> wieder zurueck
6884//### wenn interred ok, dann wieder einsetzen (am Schluss)
6885//### bottelnecks bei maps beheben
6886//### minor verbessern
6887//### preimage verbessern (Ist imm Kern map oder imap verwendet?)
6888//### Gleich in Ordnung dp wechseln, ringlist verwenden
6889//### interred ev nur zum Schluss
6890//    (z.B. wenn nacher std; wenn nacher minor: testen )
6891
6892//Zeiten mit normalV5.lib (mstd aktiv, interred inaktiv)
6893
6894//SWANSON EXAMPLES: (Macaulay2, icFracP=normalP, icFractions<->normal)
6895//---------------------------------------------------------------------
6896//1. Series Fp[x,y,u,v]/(x2v-y2u)
6897//-------------------------------
6898//characteristic p   2   3    5    7    11   13   17   37   97
6899//icFracP          0.04 0.03 0.04 0.04 0.04 0.05 0.05 0.13 0.59  Mac
6900//normalP           0   0    0    0     0    0    0    0   1    Sing
6901//icFractions      0.08 0.09 0.09 0.09 0.14 0.15 0.15 0.15 0.15  Mac
6902//normal             0   0    0    0     0    0    0    0    0   Sing
6903
69042. Series Fp[u, v, w, x, y, z]/u2x4+uvy4+v2z4
6905//--------------------------------------------
6906//characteristic p 2    3    5    7   11
6907//icFracP         0.07 0.22 9.67 143 12543
6908//normalP          0    0    5   42  1566
6909//icFractions     1.16   *    *   *    *       *: > 6h
6910//normal            0    0    0   0    0
6911
6912//3. Series Fp[u, v, w, x, y, z]/(u2xp+uvyp+v2zp)
6913//-----------------------------------------------
6914//characteristic p  2    3    5    7    11   13  17 19 23
6915//icFracP          0.06 0.07 0.09 0.27 1.81 4.89 26 56 225
6916//normalP          0     0    0    0    1    2  6  10  27
6917//icFractions      0.16 1.49 75.00 4009 *    *   *  *  *
6918//normal            0     0    2   836
6919//normal(neu)       0     0    1   2    10  155
6920//### p=7 normal braucht 807 sec in:
6921// ideal endid  = phi1(endid);      //### bottelneck'
6922
6923//1.
6924int p = 2;  ring r = p,(u,v,x,y,z),dp; ideal i = x2v-y2u;
6925//2.
6926int p = 7; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2x4+uvy4+v2z4;
6927//3.
6928int p=11; ring r=p,(u,v,w,x,y,z),dp; ideal i=u2*x^p+uv*y^p+v2*z^p;
6929
6930//IRREDUCIBLE EXAMPLES:
6931//---------------------
6932//timing for MacBookPro 2.2GHz Intel Core 2 Duo, 4GB Ram
6933//Sing. ix86Mac-darwin version 3-1-0 (3100-2008101314)  Oct 13 2008 14:46:59
6934//if no time is given: < 1  sec
6935
6936//Apply:
6937list nor = normal(i,"isPrim"); nor;
6938list nor = normalP(i,"withRing","isPrim"); nor;
6939def R=nor[1][1]; setring R; norid; normap;
6940setring r;
6941norTest(i,nor);
6942
6943int tt = timer;
6944list nor = normalP(i,"withRing","isPrim"); nor;
6945timer-tt;
6946int tt = timer;
6947list nor = normal(i,"isPrim");
6948timer-tt;
6949
6950ring r=19,(x,y,u,v),dp;    //delta -1
6951ideal i=x2v-y2u;
6952//norTest 2 sec
6953
6954ring r=2,(y,x2,x1),lp;     //delta -1
6955ideal i=y^4+y^2*x2*x1+x2^3*x1^2+x2^2*x1^3;
6956//### norid hat 1 Element nach interred
6957
6958ring r  = 11,(x,y,z),wp(2,1,2); //alles < 1 sec
6959ideal i=z3 - xy4 + x2;          //not reduced, delta =0 ok
6960ideal i=y4+x5+y2x;              //not reduced, delta -1
6961//interred verkleinert norid
6962
6963ring r=3,(u,v,x,y,z),dp;   //delta -1
6964ideal i=u2x3+uvy3+v2z3;
6965
6966ring r=3,(u,v,x,y,z),dp;   //delta -1
6967ideal i=u2x4+uvy4+v2z4;
6968//norTest(i,nor);  0 sec, norTest(i,nor) haengt!
6969
6970ring r=5,(u,v,x,y,z),dp;   //delta -1
6971ideal i=u2x6+uvy6+v2z6;
6972//normalP 5sec, normalC 1sec
6973//V5: norTest(i,nor); 45 sec bei normalP, V6 12 sec
6974//28 sec bei normal
6975
6976ring r=5,(u,v,x,y,z),dp;   //delta -1
6977ideal i=u2x5+uvy5+v2z5;
6978//normalP 1sec, normalC 1 sec,
6979//norTest lange: minor(jacob(I),h,J) 193 (308)sec, haengt dann bei M = std(M);
6980//norTest(i,nor,2); verwenden!
6981//Sing 3.0-4 orig  >9h! haengt bei Q = mstd(Q)[2];
6982
6983ring r=2,(y,x),wp(12,5);  //delta 3
6984ideal i=y5+y2x4+y2x+yx2+x12;
6985//normalP 0 sec (Test 0 sec), normalC 2 sec (Test 2 sec)
6986//normalC withGens (ohne interred) 0sec
6987
6988ring r=2,(y,x),dp;       //delta= 22
6989ideal i=y9+y8x+y8+y5+y4x+y3x2+y2x3+yx8+x9;
6990//normalP 1sec, interred verkleinert norid betraechtlich
6991//normalC haengt bei minor, ideal im loop wird zu gross ###
6992//interred bei normalC vergroeesert string um Faktor 4000!
6993//withGens haengt bei interred in loop 4 (> 10 h) oder
6994//(nach Ausschalten von interred) bei
6995//int delt=vdim(std(modulo(f,ideal(p)))); (>?h)
6996
6997//Leonard1: (1. Komponente von Leonard),  delta -1
6998ring r=2,(v,u,z,y,x),dp;
6999ideal i = z3+zyx+y3x2+y2x3, uyx+z2,uz+z+y2x+yx2, u2+u+zy+zx,
7000          v3+vux+vz2+vzyx+vzx+uz3+uz2y+z3+z2yx2;
7001//normalP 5 sec (withRing 9 sec), norTest(i,nor,2); 45 sec
7002//normalC 102sec, 99sec
7003//### Zeit wird bei ideal Ann = quotient(SM[2],SL[1]); und bei
7004// f  = quotient(p*J,J); verbraucht
7005//withGens (ohne interred) 131sec, norTest(i,nor,2); 2min25sec
7006//norTest(i,nor,2);  45 sec
7007
7008 ring r=2,(y,x),wp(25,21); //Leonard2, delta 232
7009 ring r=2,(y,x),dp;
7010 ideal i=
7011 y^21+y^20*x +y^18*(x^3+x+1) +y^17*(x^3+1) +y^16*(x^4+x)
7012 +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)
7013 +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)
7014 +y^10*(x^12+x^9+x^8+x^7+x^5+x^3+x+1)
7015 +y^9*(x^14+x^13+x^10+x^9+x^8+x^7+x^6+x^3+x^2+1)
7016 +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)
7017 +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)
7018 +y^4*(x^19+x^16+x^15+x^12+x^6+x^5+x^3+1)
7019 +y^3*(x^18+x^15+x^12+x^10+x^9+x^7+x^4+x)
7020 +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)
7021 +y*(x^23+x^22+x^20+x^17+x^15+x^14+x^12+x^9)
7022 +(x^25+x^23+x^19+x^17+x^15+x^13+x^11+x^5);
7023//normalP: dp 2sec withRing 8sec,
7024//wp 4sec, withRing:51sec Zeit in lin = subst(lin, var(ii), vip); in elimpart ),
7025//norTest(i,nor,2): haengt bei mstd(norid);
7026//### normalC: (m. interred): haengt bei endid = interred(endid);
7027//GEFIXTES INTERRED ABWARTEN. Dann interred aktivieren
7028//interred(norid) haengt u. mst(norid) zu lange
7029//(o.interred): haengt bei  haengt bei list SM = mstd(i);
7030//ideal in der Mitte zu gross
7031//i = Ideal (size 118, 13 var) fuer die neue Normalisierung
7032//normal(neu) haengt bei return(std(i)) (offensichtlich in eineranderen lib)
7033
7034REDUCIBLE EXAMPLES:
7035------------------
7036//Apply:
7037int tt = timer;
7038list nor=normalP(i,"isPrim","withRing");
7039timer-tt;
7040
7041list nor = normal(i); nor;
7042list nor = normalC(i); nor;
7043list nor = normalC(i, "withGens"); nor;
7044list nor = normalP(i,"withRing"); nor;
7045list nor = normalP(i); nor;
7046def R=nor[1][1]; setring R; norid; normap;
7047
7048//Leonhard 4 Komponenten, dim=2, delta: 0,0,0,-1
7049ring r=2,(v,u,z,y,x),dp;      //lp zu lange
7050ideal i=z3+zyx+y3x2+y2x3, uyx+z2, v3+vuyx+vux+vzyx+vzx+uy3x2+uy2x+zy3x+zy2x2;
7051//normalP: 19 sec, withRing: 22 sec
7052//normalC ohne (mit) interred: 112 (113)sec, equidim: 99sec
7053//normalC 1. mal 111 sec, (2.mal) 450sec!! 3.mal 172 sec
7054//(unterschiedlich lange primdec, mit Auswirkungen)
7055//char 19: normalC: 15sec , withGens: 14sec (o.interr.)
7056
7057//----------------------Test Example for special cases -------------------
7058int tt = timer;
7059list nor=normalP(i,"withRing");nor;
7060//list nor=normalP(i,"withRing", "isPrim");nor;
7061timer-tt;
7062def R1 = nor[1][1]; setring R1;  norid; normap; interred(norid);
7063setring r;
7064
7065int tt = timer;
7066list nor=normal(i,"isPrim");nor;
7067timer-tt;
7068
7069ring r = 29,(x,y,z),dp;
7070ideal i = x2y2,x2z2;       //Nicht equidimensional, equidim reduziert nicht, ok
7071ideal i  = xyz*(z3-xy4);   //### interred(norid) verkuerzt
7072//je 0 sec
7073
7074ideal j = x,y;
7075ideal i = j*xy;
7076equidim(i);
7077//hat eingebettete Komponente, equidim rechnet wie in Beschreibung (ok)
7078
7079ring r  = 19,(x,y),dp;
7080   ideal i = x3-y4;                   //delta = 3
7081   ideal i = y*x*(x3-y4);             //delta = 11; 0,0,3
7082   ideal i = (x2-y3)*(x3-y4);         //delta = 13; 1,3
7083   ideal i = (x-y)*(x3+y2)*(x3-y4);   //delta = 23; 0,1,3
7084   ideal i = (x-1)*(x3+y2)*(x2-y3);   //delta = 16; 0,1,1
7085   ideal i = (x-y^2)^2 - y*x^3;       //delta = 3
7086   //singularities at not only at 0, hier rechnet equidim falsch
7087
7088// -------------------------- General Examples  ---------------------------//Huneke, irred., delta=2 (Version 3-0-4: < 1sec)
7089//Version 3-0-6 default: 1sec, mit gens 2sec, mit delta 5 sec
7090//(prim,noFac):ca 7 Min, prim:ca 10 min(wg facstd)
7091//
7092// "equidim" < 1sec irred. 5sec
7093// ring r=31991,(a,b,c,d,e),dp;
7094ring r=2,(a,b,c,d,e),dp;                    //delta=2
7095ideal i=
70965abcde-a5-b5-c5-d5-e5,
7097ab3c+bc3d+a3be+cd3e+ade3,
7098a2bc2+b2cd2+a2d2e+ab2e2+c2de2,
7099abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5,
7100ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5,
7101a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6,
7102a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4,
7103b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5;
7104//normalC: char 2, 31991: 0 sec (isPrim); char 2, equidim: 7 sec
7105//norTest(i,nor,2); 1sec
7106//normalP char 2: 1sec (isPrim)
7107//size(norid); size(string(norid));21 1219 interred(norid): 21 1245 (0 sec)
7108
7109int tt = timer;
7110list nor=normalC(i);nor;
7111timer-tt;
7112
7113list nor = normalP(i,"isPrim");
7114
7115//Vasconcelos irred., delta -1 (dauert laenger)
7116//auf macbook pro = 20 sec mit alter Version,
7117//Sing 3-0-6:
7118// Char 32003: "equidim" 30 sec, "noFac": 30sec
7119//gens: nach 9 min abgebr (haengt in Lin = ideal(T*syzf);) !!!! Hans zu tun
7120//Char 2: default (charp) 2 sec, normalC ca 30 sec
7121//ring r=32003,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
7122ring r=2,(x,y,z,w,t),dp;   //dim 2, dim s_locus 1
7123ideal i= x2+zw, y3+xwt, xw3+z3t+ywt2, y2w4-xy2z2t-w3t3;
7124//normalC: char 2: 22,  sec (mit und ohne isPrim)
7125//normalP char 2: 0sec (isPrim)      o. interred
7126//char 32003: ### haengt in ideal endid  = phi1(endid);
7127
7128//-------------------------------------------------------
7129//kleine Beispiele:
7130
7131//Theo1, irred, delta=-1
7132//normalC: 1sec, normalP: 3 sec
7133ring r=32003,(x,y,z),wp(2,3,6); //dim 2,dim slocus 1
7134ideal i=zy2-zx3-x6;
7135//normalC: char 2,19,32003: 0  sec (isPrim)
7136//normalP (isPrim) char 2,19: 0sec, char 29: 1sec
7137
7138//Theo1a, CohenMacaulay regular in codim 2, dim slocus=1, delta=0
7139//normalC: 0 sec, normalP: haegt in K=preimage(R,phi,L);
7140ring r=32003,(x,y,z,u),dp;
7141ideal i=zy2-zx3-x6+u2;
7142//normalC: char 2,32003: 0  sec (isPrim)
7143//normalP (isPrim) char 2: 0sec, char 19: haengt in K = preimage(Q,phi,L);
7144
7145//Theo2, irreduzibel, reduziert, < 1sec, delta -1
7146ring r=0,(x,y,z),wp(3,4,12);
7147ideal i=z*(y3-x4)+x8;
7148//normalC: char 2,32003,0: 0  sec (isPrim)
7149//normalP (isPrim) char 2: 0 1sec, char 19: 1sec char 29: 7 sec
7150
7151//Theo2a, reduiziert, 2-dim, dim_slocus=1, alte Version 3 sec,
7152//normalP ca 30 sec, normalC ca 4sec, delta -1
7153//ring r=32003,(T(1..4)),wp(3,4,12,17);
7154//ring r=11,(T(1..4)),dp;
7155ring r=11,(T(1..4)),wp(3,4,12,17);
7156ideal i=
7157T(1)^8-T(1)^4*T(3)+T(2)^3*T(3),
7158T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4),
7159T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4),
7160T(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;
7161//normalC: char 2,32003: 0  sec (isPrim)
7162//normalP (isPrim) char 2: 0sec, char 11 2se, char 19: 13sec
7163//norTest 48sec in char11
7164//### interred verkuerzt
7165//char 29: haengt in K = preimage(Q,phi,L);
7166
7167//Theo3, irred, 2-dim, 1-dim sing, < 1sec
7168ring r=11,(x,y,z),wp(3,5,15);
7169ideal i=z*(y3-x5)+x10;
7170//normalC: char 2,0: 0  sec (withRing)
7171//normalP (withRing) char 2,11: 0sec, char 19: 13sec norTest 12sec(char 11)
7172
7173//Theo4 reducible, delta (0,0,0) -1
7174ring r=29,(x,y,z),dp;
7175ideal i=(x-y)*(x-z)*(y-z);
7176//normalC: char 2,32003: 0  sec
7177//normalP char withRing 2, 29: 0sec, 6sec
7178
7179//Theo6
7180ring r=32003,(x,y,z),dp;
7181ideal i=x2y2+x2z2+y2z2;
7182//normalC: char 2,32003: 0  sec
7183//normalP char withRing 2, 29: 0sec, 4sec
7184
7185//Sturmfels, CM, 15 componenten, alle glatt
7186ring r=0,(b,s,t,u,v,w,x,y,z),dp;
7187ideal i= bv+su, bw+tu, sw+tv, by+sx, bz+tx, sz+ty,uy+vx,uz+wx,vz+wy,bvz;
7188//normalC car 11, 0: 1sec, normalP 0 sec
7189
7190//riemenschneider, , dim 3, 5 Komp. delta (0,0,0,0,0), -1
7191ring r=2,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1);
7192ideal i=xz,vx,ux,su,qu,txy,stx,qtx,uv2z-uwz,uv3-uvw,puv2-puw;
7193//alles 0 sec in char 2
7194
7195//4 Komponenten, alle glatt, 0sec
7196ring r=11,(x,y,z,t),dp;
7197ideal i=x2z+xzt,xyz,xy2-xyt,x2y+xyt;
7198
7199//dim 3, 2 Komponenten delta (-1,0), -1
7200ring r=2,(u,v,w,x,y,z),wp(1,1,1,3,2,1);
7201ideal i=wx,wy,wz,vx,vy,vz,ux,uy,uz,y3-x2;
7202//alles 0 sec in char 2
7203//---------------------------------------------------------
7204int tt = timer;
7205list nor=normalP(i,"normalC","withRing");nor;
7206timer-tt;
7207
7208//St_S/Y, 3 Komponenten, 2 glatt, 1 normal
7209//charp haengt (in char 20) in K=preimage(R,phi,L);
7210//ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
7211ring r=11,(b,s,t,u,v,w,x,y,z),dp;
7212ideal i=wy-vz,vx-uy,tv-sw,su-bv,tuy-bvz;
7213//normalC: char 2,32003: 0  sec
7214//normalP char withRing 2: 1sec, char 11: 40sec
7215
7216//Horrocks: cahr 0: 17 (8 in char 11) Komponenten alle normal, delta 1
7217//char 11: 8 Komponenten alle normal, delta -1
7218ring r=0,(a,b,c,d,e,f),dp;
7219//ring r=11,(a,b,c,d,e,f),dp; //Charp bis p = 200 ca 3sec
7220ideal i=
7221adef-16000be2f+16001cef2, ad2f+8002bdef+8001cdf2, abdf-16000b2ef+16001bcf2,
7222a2df+8002abef+8001acf2, ad2e-8000bde2-7999cdef, acde-16000bce2+16001c2ef,
7223a2de-8000abe2-7999acef, acd2+8002bcde+8001c2df, abd2-8000b2de-7999bcdf,
7224a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2,
7225abde-8000b2e2-acdf-16001bcef-8001c2f2, abcd-16000b2ce+16001bc2f,
7226a2cd+8002abce+8001ac2f, a2bd-8000ab2e-7999abcf, ab3f-3bdf3,
7227a2b2f-2adf3-16000bef3+16001cf4, a3bf+4aef3, ac3e-10668cde3,
7228a2c2e+10667ade3+16001be4+5334ce3f, a3ce+10669ae3f, bc3d+8001cd3e,
7229ac3d+8000bc3e+16001cd2e2+8001c4f, b2c2d+16001ad4+4000bd3e+12001cd3f,
7230b2c2e-10668bc3f-10667cd2ef, abc2e-cde2f, b3cd-8000bd3f, b3ce-10668b2c2f-10667bd2ef, abc2f-cdef2, a2bce-16000be3f+16001ce2f2,
7231ab3d-8000b4e-8001b3cf+16000bd2f2, ab2cf-bdef2,
7232a2bcf-16000be2f2+16001cef3, a4d-8000a3be+8001a3cf-2ae2f2;
7233//normalC: char 0: 1sec char 11: 0sec
7234//normalP: char 11: 0sec
7235
7236//2sec mit normalC, in char 2 ebenfalls (char 20 mit charp >1 min)
7237//4 Komp. in char 2, delta (0,0,0,0) -1, char 11:delta (-1,0,0,0) -1
7238ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
7239ideal i=
7240wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3,
7241vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2,
7242tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2,
7243sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2,
7244tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3;
7245//normalC: char 2,32003: 1 sec
7246//normalP char withRing 2: 1sec, char 11: 40sec
7247
7248//---------------------------------------------------------
7249//genus:
7250int tt = timer;
7251list nor=normal(i, "noFac");nor;
7252timer-tt;
7253
7254//Yoshihiko Sakai, irred, 0sec, delta = 8
7255ring r=0,(x,y),dp;                    //genus 0 4 nodes and 6 cusps im P2
7256//ring r=7,(x,y),dp;                  //charp haengt in K = preimage(Q,phi,L)
7257ideal i=(x2+y^2-1)^3 +27x2y2;
7258
7259ring r=0,(x,y),dp;   //genus 0
7260ideal i=(x-y^2)^2 - y*x^3;
7261
7262ring r=0,(x,y),dp;  //genus 4
7263ideal i=y3-x6+1;
7264
7265int m=9;           // q=9: genus 0
7266int p=2;
7267int q=9;//2,...,9
7268ring r=0,(x,y),dp;
7269ideal i=y^m - x^p*(x - 1)^q;
7270
7271ring r=0,(x,y),dp;  //genus 19
7272ideal 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;
7273
7274ring r=23,(x,y),dp;  //genus 34, delta 2
7275ideal i=y10+(-2494x2+474)*y8+(84366+2042158x4-660492)*y6
7276        +(128361096x4-47970216x2+6697080-761328152x6)*y4
7277        +(-12024807786x4-506101284x2+15052058268x6+202172841-3212x8)*y2
7278        +34263110700x4-228715574724x6+5431439286x2+201803238
7279        -9127158539954x10-3212722859346x8;
7280//normalC, normalP 0 sec
7281
7282//Rob Koelman
7283//ring r=0,(x,y,z),dp;      //dim sing = 1 (nach ca 15 min abgebrochen)
7284ring r=32003,(x,y,z),dp;
7285ideal i=
7286761328152*x^6*z^4-5431439286*x^2*y^8+2494*x^2*z^8+228715574724*x^6*y^4+
7287 9127158539954*x^10-15052058268*x^6*y^2*z^2+3212722859346*x^8*y^2-
7288 134266087241*x^8*z^2-202172841*y^8*z^2-34263110700*x^4*y^6-6697080*y^6*z^4-
7289 2042158*x^4*z^6-201803238*y^10+12024807786*x^4*y^4*z^2-128361096*x^4*y^2*z^4+
7290 506101284*x^2*z^2*y^6+47970216*x^2*z^4*y^4+660492*x^2*z^6*y^2-
7291 z^10-474*z^8*y^2-84366*z^6*y^4;
7292//normalC char 32003: 10 sec, char 0 :
7293
7294//ring r=0,(x,y),dp;//genus 10  with 26 cusps (nach ca 4 min abgebrochen)
7295ring r=32003,(x,y),dp;    //24 sing, delta 24
7296ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6
7297-5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4
7298+506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6
7299-2042158x4+660492x2y2-84366y4+2494x2-474y2-1;
7300//normalC 32003: 4 sec, char 0: abgebrochen bei pr = facstd(i); ###
7301
7302ring r=0,(x,y),dp;   //irred, genus 1  with 5 cusps, delta 5
7303ideal i=57y5+516x4y-320x4+66y4-340x2y3+73y3+128x2-84x2y2-96x2y;
7304//normalC 0 sec
7305
7306ring r=2,(x,y),dp;  //genus 4, 2 Zweige, delta (13,9) 89
7307ideal i=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
7308//normalC: char 2 : 1sec, char 0: lange
7309//normalP char 2 withRing: 0sec
7310
7311ring r=2,(y,z,w,u),dp; //2 Komp. genus -5
7312ideal i=y2+z2+w2+u2,w4-u4;
7313//normalC: char 2 : 0sec, char 0: 1sec
7314//normalP char 2 withRing: 0sec
7315
7316ring r=0,(y,z,w,u),dp; //irred. genus 9
7317ideal i=y2+z2+w2+u2,z4+w4+u4;
7318//char 0: 0sec
7319
7320ring r=0,(x,y,t),dp;  //irred, delta -1
7321ideal i= 25x8+200x7y+720x6y2+1520x5y3+2064x4y4+1856x3y5+1088x2y6+384xy7+64y8-12x6t2-72x5yt2-184x4y2t2-256x3y3t2-192x2y4t2-64xy5t2-2x4t4-8x3yt4+16xy3t4+16y4t4+4x2t6+8xyt6+8y2t6+t8;
7322//char 0: 0sec
7323
7324ring r=0,(x,y,z,w,u),dp;
7325ideal i=x2+y2+z2+w2+u2,x3+y3+z3,z4+w4+u4;
7326//char 0: 0sec
7327
7328//---------------------------------------------------------
7329//Probleme mit normalC in char 2 und char 0
7330
7331int tt = timer;
7332list nor=normalC(i,"withRing");nor;
7333timer-tt;
7334
7335//Mark van Hoeij
7336ring r=3,(x,y),dp;  //genus 19, delta 21
7337ideal i=y20+y13x+x4y5+x3*(x+1)^2;
7338//normalC: char 2 > 10 min   bei list SM = mstd(i);###
7339//normalP char 2 withRing: 0sec, char 11: haengt bei K = preimage(Q,phi,L);
7340
7341ring r=2,(x,y),dp;  //genus 35
7342ideal i=y30+y13x+x4y5+x3*(x+1)^2;
7343//char 0 abgebrochen bei list SM = mstd(i); ###
7344//char 2 nach ca 30 min
7345//normalC: char 2: abgebr. bei list SM = mstd(i);  //Now the work starts'
7346//normalC, withGens, char 2: abgebrochen bei Q=mstd(Q)[2];
7347//normalP char 2 withRing: 0sec
7348
7349ring r=0,(x,y),dp;   //irred, genus 55, delta 21
7350ideal i=y40+y13x+x4y5+x3*(x+1)^2;
7351//normalC: char 2 lange
7352//normalP char 2 withRing: 0sec
7353
7354ring r=29,(x,y,t),dp; //char 0: genus -5, 4 Komp, delta (-1,-1,0,0), -1
7355ideal 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;
7356//normalC: char 29 : 0sec, char 0: 0sec  //char 29 6 Komponenten
7357//normalP char 29 withRing: 1sec
7358
7359//-------------------------- problematic examples ------------------------
7360//ring r=0,(x,y,t),dp;
7361ring r=32003,(x,y,t),dp;
7362ideal i=
736332761x8+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;
7364//char 0: lange (es liegt an den grossen Zahlen), char 32003: 0 sec
7365
7366//dasselbe Beipiel in char 19: irred
7367ring r=0,(x,y,t),dp;
7368ideal i=
73695x8+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;
7370//normalP: char 2,3: 0sec, norTest 0,2 sec, char 11 haengt bei peimage
7371//normalC: char 3: 0 sec, char 0: 1sec
7372
7373//ring r=0,(x,y),dp;
7374ring r=32003,(x,y),dp;
7375ideal i=
7376x30y21+21x29y20+210x28y19+10x27y19+1330x27y18+190x26y18+5985x26y17
7377+1710x25y17+20349x25y16+45x24y17+9690x24y16+54264x24y15+765x23y16
7378+38760x23y15+116280x23y14+6120x22y15+116280x22y14+120x21y15
7379+203490x22y13+30600x21y14+271320x21y13+1799x20y14+293930x21y12+107100x20y13
7380+503880x20y12+12586x19y13+352716x20y11+278460x19y12+210x18y13+755820x19y11
7381+54509x18y12+352716x19y10+556920x18y11+2723x17y12+923780x18y10+163436x17y11
7382+293930x18y9+875160x17y10+16296x16y11+923780x17y9+359359x16y10+252x15y11
7383+203490x17y8+1093950x16y9+59598x15y10+755820x16y8+598598x15y9+2751x14y10
7384+116280x16y7+1093950x15y8+148610x14y9+503880x15y7+769197x14y8+13650x13y9
7385+54264x15y6+875160x14y7+266805x13y8+210x12y9+271320x14y6+768768x13y7
7386+40635x12y8+20349x14y5+556920x13y6+354816x12y7+1855x11y8+116280x13y5
7387+597597x12y6+80640x11y7+5985x13y4+278460x12y5+353892x11y6+7280x10y7+38760x12y4
7388+358358x11y5+112014x10y6+120x9y7+1330x12y3+107100x11y4+264726x10y5+16660x9y6
7389+9690x11y3+162799x10y4+111132x9y5+805x8y6+210x11y2+30600x10y3+146685x9y4
7390+24500x8y5+1710x10y2+54236x9y3+78750x8y4+2310x7y5+21x10y+6120x9y2+58520x8y3
7391+24010x7y4+45x6y5+190x9y+12509x8y2+39060x7y3+3675x6y4+x9+765x8y+15918x7y2
7392+15680x6y3+204x5y4+10x8+1786x7y+12915x6y2+3500x5y3+45x7+2646x6y+6580x5y2
7393+366x4y3+119x6+2562x5y+1995x4y2+10x3y3+203x5+1610x4y+324x3y2+231x4+630x3y
7394+23x2y2+175x3+141x2y+85x2+16xy+24x+y+4;
7395list nor = normal(i);
7396//normalC: char 0: ### haengt in SB of singular locus JM = mstd(J);
7397//normalC: char 32003,"noFac","equidim": 0sec, "noFac": 1sec
7398// ev neues interred
7399genus(i);         // haengt bei int mu=vdim(std(jacob(f)));
7400                  //### ist das noetig?
7401
7402//Singular rechnet genus richtig, auch im Fall, dass Kurve irreduzibel,
7403//aber nicht absolut irreduzibel ist:
7404ring r = 0,(x,y),dp;
7405ideal i = x2+y2;      //irreduzibel /Q aber reduzibel /C (x-iy)*(x+iy)
7406factorize( x2+y2);    //liefert irreduzibel
7407genus(i);             //sollte 0+0-2+1= -1 sein
7408genus(i,1);           //beides ist korrekt in Singular
7409
7410*/
Note: See TracBrowser for help on using the repository browser.