source: git/Singular/LIB/normal.lib @ e05cb2

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