# source:git/Singular/LIB/reesclos.lib@9e1207

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