source: git/Singular/LIB/normal.lib @ 6e80ab

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