source: git/Singular/LIB/normal.lib @ 81f4bd5

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