# source:git/Singular/LIB/reesclos.lib@fbb23b

spielwiese
Last change on this file since fbb23b was fbb23b, checked in by Hans Schoenemann <hannes@…>, 13 years ago
fix version string git-svn-id: file:///usr/local/Singular/svn/trunk@14005 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 18.1 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="\$Id\$";
3category="Commutative Algebra";
4
5info="
6LIBRARY:     reesclos.lib   PROCEDURES TO COMPUTE THE INT. CLOSURE OF AN IDEAL
7AUTHOR:      Tobias Hirsch, email: hirsch@math.tu-cottbus.de
8
9OVERVIEW:
10 A library to compute the integral closure of an ideal I in a polynomial ring
11 R=k[x(1),...,x(n)] using the Rees Algebra R[It] of I. It computes the integral
12 closure of R[It] (in the same manner as done in the library 'normal.lib'),
13 which is a graded subalgebra of R[t]. The degree-k-component is the integral
14 closure of the k-th power of I.
15
16PROCEDURES:
17 ReesAlgebra(I);        computes the Rees Algebra of an ideal I
18 normalI(I[,p[,r]]);    computes the integral closure of an ideal I using R[It]
19";
20
21LIB "normal.lib";       // for HomJJ
22LIB "standard.lib";     // for groebner
23
24///////////////////////////////////////////////////////////////////////////////
25
26proc ReesAlgebra (ideal I)
27"USAGE:    ReesAlgebra (I); I = ideal
28RETURN:   The Rees algebra R[It] as an affine ring, where I is an ideal in R.
29          The procedure returns a list containing two rings:
30          [1]: a ring, say RR; in the ring an ideal ker such that R[It]=RR/ker
31
32          [2]: a ring, say Kxt; the basering with additional variable t
33               containing an ideal mapI that defines the map RR-->Kxt
34EXAMPLE:  example ReesAlgebra; shows an example
35"
36{
37  // remember the data of the basering
38
39  def oldring = basering;
40  string oldchar = charstr(basering);
41  string oldvar  = varstr(basering);
42  string oldord  = ordstr(basering);
43  int n = ncols(I);
44  ideal m = maxideal(1);
45
46
47  // Create a new ring with variables for each generator of I
48
49  execute ("ring Rees = "+oldchar+",("+oldvar+",U(1.."+string(n)+")),dp");
50
51
52  // Kxt is the old ring with additional variable t
53  // Here I -> t*I, so the generators of I generate the subalgebra R[It] in Kxt
54
55  execute ("ring Kxt = "+oldchar+",("+oldvar+",t),dp");
56  ideal I = fetch(oldring,I);
57  ideal m = fetch(oldring,m);
58  int k;
59  for (k=1;k<=n;k++)
60  {
61    I[k]=t*I[k];
62  }
63
64
65  // Now we map from Rees to Kxt, identity on the original variables, and
66  // U(k) -> I[k]
67
68  ideal mapI = m,I;
69  map phi = Rees,mapI;
70  ideal zero = 0;
71  export (mapI);
72
73  // Now the Rees-Algebra is Rees/ker(phi)
74
75  setring Rees;
76  ideal ker = preimage(Kxt,phi,zero);
77  export (ker);
78
79  list result = Rees,Kxt;
80
81  return(result);
82
83}
84example
85{
86  "EXAMPLE:"; echo=2;
87  ring R = 0,(x,y),dp;
88  ideal I = x2,xy4,y5;
89  list L = ReesAlgebra(I);
90  def Rees = L[1];       // defines the ring Rees, containing the ideal ker
91  setring Rees;          // passes to the ring Rees
92  Rees;
93  ker;                   // R[It] is isomorphic to Rees/ker
94}
95
96
97////////////////////////////////////////////////////////////////////////////
98
99static
100proc ClosureRees (list L)
101"USAGE:    ClosureRees (L); L a list
102ASSUME:   L is a list containing
103          - a ring L[1], inside L[1] an ideal ker such that L[1]/ker is
104            isomorphic to the Rees Algebra R[It] of an ideal I in k[x]
105          - a ring L[2]=k[x,t], inside L[1] an ideal mapI defining the
106            map L[1] --> L[2] with image R[It]
107RETURN:   quotients of elements of k[x,t] representing generators of the
108          integral closure of R[It]. The result of ClosureRees is a list
109          images, the first size(images)-1 entries are the numerators of the
110          generators, the last one is the universal denominator
111"
112{
113  int dblvl=printlevel-voice+2;   // toggles how much data is printed
114                                  // during the procedure
115
116  def Kxt = basering;
117  def R(1) = L[1];
118  setring R(1);                   // declaration of variables used later
119  ideal ker(1)=ker;               // in STEP 2
120  poly p;
121
122  // some auxiliary variables
123
124  int i=1;         // counts the number of steps to reach the closure of R(1)
125  int check=0;     // a 'boolean' variable for several checks
126
127  /////// STEP 1: ///////////////////////////////////////////////////////////
128  // construct R(i) step by step as done in normal.lib; 2 differences:     //
129  //  - since the input is a prime ideal, no splitting will occur          //
130  //  - the intermediate rings and nonzerodivisors for J are remembered    //
131  ///////////////////////////////////////////////////////////////////////////
132
133  if (dblvl>0)
134  {
135    "// STEP 1: Compute the integral closure of R[It]";
136  }
137
138  list RS;
139
140  while (check==0)  // repeat until the closure is reached
141  {
142    // construction of HomJJ, J an ideal containing the non-normal locus
143    // of R(i)/ker(i), as done in normalizationPrimes in normal.lib for
144    // the special case that we are working with a prime ideal
145
146    if (dblvl>0)
147    {
148      "";
149      "// We are in step",i;
150      pause();
151    }
152
153    if (homog(ker(i))==1)
154    {
155      list SM=mstd(ker(i));
156    }
157    else
158    {
159      list SM=groebner(ker(i)),ker(i);
160    }
161
162    if (dblvl>0)
163    {
164      "// Standard basis of the current ideal:";
165      SM[1];
166    }
167
168    // In the first iteration, we have to compute the singular locus
169    // "from scratch". In further iterations, we can fetch it from the
170    // previous one but have to compute its radical.
171
172    if (i==1)
173    {
174      ideal J=minor(jacob(SM[2]),nvars(basering)-dim(SM[1]));
175      ideal sin=J+SM[2];
176      if (homog(sin)==1)
177      {
178        list JM=mstd(sin);
179      }
180      else
181      {
182        list JM=groebner(sin),sin;
183      }
184
186    }
187    else
188    {
190    }
191
192    if (dblvl>0)
193    {
194      "// Radical of the singular locus:";
195      J;
196    }
197
198    ideal SL=simplify(reduce(JM[2],SM[1]),2);
199    JM =J,J;
200    poly nzd=SL[1];           // universal denominator for HomJJ
201
202    if (dblvl>0)
203    {
204      "// The non-zerodivisor";
205      nzd;
206      pause();
207    }
208
209    list RR=SM[1],SM[2],JM[2],SL[1];
210    RS=HomJJ(RR);
211
212    if (RS[2]==1)             // we've reached the integral closure
213    {
214      if (dblvl>0)
215      {
216        "";
217        "// We've reached the integral closure after",i,"iterations";
218        pause();
219      }
220      check=1;
221    }
222    else                      // prepare the next iteration with new
223    {                         // ring R(i+1)
224      ideal MJ=JM[2];
225
226      def R(i+1)=RS[1];       // the data of and some variable decla-
227      setring R(i+1);         // rations in R(i+1) needed in STEP 2
228      ideal ker(i+1)=endid;
229      map phi=R(i),endphi;
230      poly p;
231
232      list JM=mstd(simplify(phi(MJ)+ker(i+1),4)); // fetch the singular locus
233      i++;
234    }
235  }
236
237  if (i==1)                     // R[It] (and thus I) was integrally
238  {                             // closed ==> we're already done
239        list result="closed";
240        return(result);
241  }
242
243
244  /////// STEP 2: ////////////////////////////////////////////////////////
245  // compute representations of the ring variables of R(i) as fractions //
246  // of elements of R(1);                                               //
247  ////////////////////////////////////////////////////////////////////////
248
249  int length=nvars(R(i));  // the number of variables of the last ring
250  int j,k,n;               // some counters
251  string mapstr;           // will be used while constructing preimages
252
253  list preimages;          // here the fractions are stored in the
254                           // form var(j)=preimages[j]/preimages[length+1]
255                           // ('=' means identification via the inclusion)
256                           // the last entry corresponds to nzd in R(i)
257
258  // For each variable (counter j) and for each intermediate ring (k):
259  // find preimages of var(j)*endphi(nzd_k-1) in R(k-1).
260  // Finally, do the same for nzd.
261
262  if (dblvl>0)
263  {
264    "// STEP 2: Compute fractions representing the ring variables of";
265    "           the last ring";
266  }
267
268  for (j=1;j<=length+1;j++)
269  {
270    setring R(i);
271
272    if (j<=length)          // do it with for ring variables...
273    {
274      p=var(j);
275    }
276    else                   // ...and finally for nzd in R(i)
277    {
278      p=1;
279    }
280
281    // get back from R(i) to R(1) step by step
282
283    for (k=i;k>1;k--)
284    {
285      // clear the fraction in the representation in R(i)
286
287      p=p*phi(nzd);
288
289      // compute the preimage of [p mod ker(k)] under phi in R(k-1):
290      // as p is an element of im(phi), there is a polynomial h such that
291      // h(vars(R(k-1)) is mapped to [p mod ker (k)], and h can be com-
292      // puted as the normal form of a w.r.t <ker(k),Z(n)-phi(k)(n)> in R(k)[Z]
293
294         // compute this normal form h...
295
296      if (j==1)    // in the first iteration: construct S(k)=R(k)[Z], fetch
297                   // endphi (the ideal defining phi) and ker(k) and construct
298                   // the ideal from above
299      {
300        execute ("ring S(k) = "+charstr(R(k))+",("+varstr(R(k))+",Z(1.."+string(ncols(endphi))+")),(dp("+string(nvars(R(k)))+"),dp("+string(ncols(endphi))+"));");
301        ideal endphi = imap(R(k),endphi);
302        ideal J = imap(R(k),ker(k));
303        for (n=1;n<=ncols(endphi);n++)
304        {
305          J=J+(Z(n)-endphi[n]);
306        }
307        J=groebner(J);
308        poly h=NF(imap(R(k),p),J);
309      }
310      else
311      {
312        setring S(k);
313        h=NF(imap(R(k),p),J);
314      }
315
316         // and compute h(vars(R(k-1)))
317
318      setring R(k-1);
319      if (j==1)  // in the first iteration: compute backmap:S(k)-->R(k-1)
320      {
321        mapstr="map backmap = S(k),";
322        for (n=1;n<=nvars(R(k));n++)
323        {
324          mapstr=mapstr+"0,";
325        }
326        execute (mapstr+"maxideal(1);");
327      }
328
329      p=NF(backmap(h),std(ker(k-1)));
330    }
331
332    // when down to R(1), store the result in the list preimages
333
334    preimages=insert(preimages,p,j-1);
335    if (dblvl>0)
336    {
337      if (j<=length)
338      {
339        "numerator of variable ",j,":",p;
340      }
341      else
342      {
343        "";
344        "and finally the 'universal' denominator:",p;
345      }
346    }
347  }
348
349  // at the end: go back to the original basering and construct gene-
350  // rators of the closure of I
351
352  setring Kxt;
353  map psi=R(1),mapI;              // from ReesAlgebra: the map Rees->Kxt
354
355  list images=psi(preimages);
356
357  if (dblvl>-1)
358  {
359    pause();
360    "";
361    "// Get back to the original basering and construct the";
362    "// generators of the closure of I";
363    "";
364    "// Back in k[x,t], the fractions, stored in the list images:";
365    for (int j=1;j<=size(images);j++)
366    {
367      if (j<size(images))
368      {
369        "numerator",j,":",images[j];
370      }
371      else
372      {
373        "";
374        "denominator:  ",images[j];
375      }
376    }
377    pause();
378  }
379  return (images);
380}
381
382
383////////////////////////////////////////////////////////////////////////////
384
385static
386proc ClosurePower(list images, list #)
387"USAGE:    ClosurePower (L [,#]); L a list, # an optional list containing an
388          integer
389ASSUME:   - L is a list containing generators of the closure of R[It] in k[x,t]
390            (the first size(L)-1 elements are the numerators, the last one
391            is the denominator)
392          - if # is given: #[1] is an integer, compute generators for the
393                           closure of I, I^2, ..., I^#[1]
394RETURN:   the integral closure of I, ... I^#[1]. If # is not given, compute
395          the closure of all powers up to the maximum degree in t occurring
396          in the closure of R[It] (so this is the last power whose closure is
397          not just the sum/product of the smaller powers). The returned
398          result is a list of elements of k[x,t] containing generators of the
399          closure of the desired powers of I. "
400{
401  int dblvl=printlevel-voice+2;   // toggles how much data is printed
402                                  // during the procedure
403
404  int j,k,d,computepow;                    // some counters
405  int pow=0;
406  int length = size(images)-1;             // the number of generators
407  poly image;
408  poly denominator = images[length+1];     // the universal denominator
409
410  if (size(#)>0)
411  {
412    pow=#[1];
413  }
414  computepow=pow;
415
416  if (dblvl>0)
417  {
418    "";
419    "// The generators of the closure of R[It]:";
420  }
421
422  intmat m[nvars(basering)-1][1];  // an intvec used for jet and maxdeg1
423  intvec tw=m,1;                   // such that t has weight 1 and all
424                                   // other variables have weight 0
425
426  // Construct the generators of the closure of R[It] as elements of k[x,t]
427  // If # is not given, determine the highest degree pow in t that occurs.
428
429  for (j=1;j<=length;j++)
430  {
431    images[j] = (images[j]/denominator); // construct the fraction
432    image = images[j];
433    if (dblvl>0)
434    {
435      "generator",j,":",image;
436    }
437
438    if (computepow==0)              // #[1] not given or ==0 => compute pow
439    {
440      if (maxdeg1(image,tw)>pow)    // from poly.lib
441      {
442        pow=maxdeg1(image,tw);
443      }
444    }
445  }
446
447  if (dblvl>0)
448  {
449    "";
450    if (computepow==0)
451    {
452      "// Compute the closure up to the given powers of I";
453    }
454    else
455    {
456     "// Compute the closure up to the maximal power of t that occured:",pow;
457    }
458  }
459
460  // Construct a list consisting of #[1] resp. pow times the zero ideal
461
462  ideal CurrentPower=0;
463  list result;
464  for (k=1;k<=pow;k++)
465  {
466    result=insert(result,CurrentPower);
467  }
468
469  // For each generator and each k, add its degree-k-coefficient to the #
470  // closure of I^k
471
472  for (j=1;j<=length;j++)
473  {
474    for (k=1;k<=pow;k++)
475    {
476      image=images[j]-jet(images[j],k-1,tw);
477      if (image<>0)
478      {
479        image=subst(image/t^k,t,0);
480        if (image<>0)
481        {
482          result[k]=result[k]+image;
483        }
484      }
485    }
486  }
487
488  if (dblvl>0)
489  {
490    "";
491    "// The 'pure' parts of degrees 1..pow:";
492    result;
493    "";
494  }
495
496  // finally, add the suitable products of generators in lower degrees
497
498  for (k=2;k<=pow;k++)
499  {
500    for (j=1;j<=(k div 2);j++)
501    {
502      result[k]=result[k]+result[j]*result[k-j];
503    }
504  }
505
506  return(result);
507}
508
509////////////////////////////////////////////////////////////////////////////
510
511proc normalI(ideal I, list #)
512"USAGE:    normalI (I [,p[,r]]); I an ideal, p and r optional integers
513RETURN:   the integral closure of I, ..., I^p, where I is an ideal in the
514          polynomial ring R=k[x(1),...x(n)]. If p is not given, or p==0,
515          compute the closure of all powers up to the maximum degree in t
516          occurring in the closure of R[It] (so this is the last power whose
517          closure is not just the sum/product of the smaller). If r
518          is given and r==1, normalI starts with a check whether I is already a
520          The result is a list containing the closure of the desired powers of
521          I as ideals of the basering.
522DISPLAY:  The procedure displays more comments for higher printlevel.
523EXAMPLE:  example normalI; shows an example
524"
525{
526  int dblvl=printlevel-voice+2;   // toggles how much data is printed
527                                  // during the procedure
528
529  def BAS=basering;               // remember the basering
530
531  // two simple cases: principal ideals and radical ideals are always
532  // integrally closed
533
534  if (size(I)<=1)        // includes the case I=(0)
535  {
536    if (dblvl>0)
537    {
538      "// Trivial case: I is a principal ideal";
539    }
540    list result=I;
541    if (size(#)>0)
542    {
543      for (int k=1;k<=#[1]-1;k++)
544      {
545        result=insert(result,I*result[k],k);
546      }
547    }
548    return(result);
549  }
550
552  if (size(#)>1)
553  {
555  }
556
558  {
559    if (dblvl>0)
560    {
561      "//Check whether I is radical";
562    }
563
565    {
566      if (dblvl>0)
567      {
568        "//Trivial case: I is a radical ideal";
569      }
570      list result=I;
571      if (size(#)>0)
572      {
573        for (int k=1;k<=#[1]-1;k++)
574        {
575          result=insert(result,I*result[k],k);
576        }
577      }
578      return(result);
579    }
580  }
581
582  // start with the computation of the Rees Algebra R[It] of I
583
584  if (dblvl>0)
585  {
587  }
588
589  list Rees = ReesAlgebra(I);
590  def R(1)=Rees[1];
591  def Kxt=Rees[2];
592  setring R(1);
593
594  if (dblvl>0)
595  {
596    R(1);
597    ker;
598    "";
599    "// Now ClosureRees computes generators for the integral closure";
600    "// of R[It] step by step";
601  }
602
603  // ClosureRees computes fractions in R[x,t] representing the generators
604  // of the closure of R[It] in k[x,t], which is the same as the closure
605  // in Q(R[It]).
606  // the first size(images)-1 entries are the numerators of the gene-
607  // rators, the last entry is the 'universal' denominator
608
609  setring Kxt;
610  list images = ClosureRees(Rees);
611
612  // ClosureRees was done after the first HomJJ-call
613  // ==> I is integrally closed, and images consists of the only entry "closed"
614
615  if ((size(images)==1) && (typeof(images[1])=="string"))
616  {
617    if (dblvl>0)
618    {
619      "//I is integrally closed!";
620    }
621
622    setring BAS;
623    list result=I;
624    if (size(#)>0)
625    {
626      for (int k=1;k<=#[1]-1;k++)
627      {
628        result=insert(result,I*result[k],k);
629      }
630    }
631    return(result);
632  }
633
634  // construct the fractions corresponding to the generators of the
635  // closure of I and its powers, depending on # (in fact, they will
636  // not be real fractions, of course). This is done in ClosurePower.
637
638  list result = ClosurePower(images,#);
639
640  // finally fetch the result to the old basering
641
642  setring BAS;
643  list result=fetch(Kxt,result);
644  return(result);
645}
646example
647{
648  "EXAMPLE:"; echo=2;
649  ring R=0,(x,y),dp;
650  ideal I = x2,xy4,y5;
651  list J = normalI(I);
652  I;
653  J;                             // J[1] is the integral closure of I
654}
655
656/*
657LIB"reesclos.lib";
658
659
660// 1.  x^i,y^i in k[x,y]
661//     geht bis i = 19 (800sec), bis i=10 wenige Sekunden,
662//     bei i = 20 ueber 1GB Hauptspeicher, in der 9. Iteration no memory
663//     (braucht 20 Iterationen)
664
665  ring r = 0,(x,y),dp;
666  int i = 6;
667  ideal I = x^i,y^i;
668  list J = normalI(I);
669  I;
670  J;
671
672
673//================================================================
674
675
676// 2. x^i,y^i,z^i in k[x,y,z]
677//    aehnlich wie 1., funktioniert aber nur bis i=5 und dauert dort
678//    >1 h
679
680
681//================================================================
682
683
684// 3. scheitert in der ersten Iteration beim Radikal
685//    Standardbasis des singulaeren Ortes: 7h (in char0),
686//    in char(p) viel schneller, obwohl kleine Koeffizienten
687//    schon bei Radikal -Test braucht er zu lang (>1h)
688
689  ring r = 0,(x,y,z),dp;
690  //ring r = 32003,(x,y,z),dp;
691
692  ideal I = x2+xy3-5z,z3+y2-xzy,x2y3z5+y3-y5;
693  list l= ReesAlgebra(I);
694  list J = normalI(I);
695  I;
696  J;
697
698*/
699
Note: See TracBrowser for help on using the repository browser.