source: git/Singular/LIB/locnormal.lib @ b0732eb

jengelh-datetimespielwiese
Last change on this file since b0732eb was b0732eb, checked in by Hans Schoenemann <hannes@…>, 11 years ago
add: new libs from master
  • Property mode set to 100644
File size: 13.8 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2version="$Id: normalTools.lib,v 1.0 2010/05/19 Exp$";
3category="Commutative Algebra";
4info="
5LIBRARY: locnormal.lib   Normalization of affine domains using local methods
6AUTHORS:  J. Boehm        boehm@mathematik.uni-kl.de
7          W. Decker       decker@mathematik.uni-kl.de
8          S. Laplagne     slaplagn@dm.uba.ar
9          G. Pfister      pfister@mathematik.uni-kl.de
10          S. Steidel      steidel@mathematik.uni-kl.de
11          A. Steenpass    steenpass@mathematik.uni-kl.de
12
13OVERVIEW:
14
15Suppose A is an affine domain over a perfect field.@*
16This library implements a local-to-global strategy for finding the normalization
17of A. Following [1], the idea is to stratify the singular locus of A, apply the
18normalization algorithm given in [2] locally at each stratum, and put the local
19results together. This approach is inherently parallel.@*
20Furthermore we allow for the optional modular computation of the local results
21as provided by modnormal.lib. See again [1] for details.
22
23REFERENCES:
24
25[1] Janko Boehm, Wolfram Decker, Santiago Laplagne, Gerhard Pfister, Stefan Steidel,
26Andreas Steenpass: Parallel algorithms for normalization, http://arxiv.org/abs/1110.4299, 2011.
27
28[2] Gert-Martin Greuel, Santiago Laplagne, Frank Seelisch: Normalization of Rings,
29Journal of Symbolic Computation 9 (2010), p. 887-901
30
31KEYWORDS:
32normalization; local methods; modular methods
33
34SEE ALSO: normal_lib, modnormal_lib
35
36PROCEDURES:
37locNormal(I, [...]);  normalization of R/I using local methods
38
39";
40
41LIB "normal.lib";
42LIB "sing.lib";
43LIB "modstd.lib";
44
45static proc changeDenom(ideal U1, poly c1, poly c2, ideal I){
46// Given a ring in the form 1/c1 * U, it computes a new ideal U2 such that the
47// ring is 1/c2 * U2.
48// The base ring is R, but the computations are to be done in R / I.
49  int a;      // counter
50  def R = basering;
51  qring Q = groebner(I);
52  ideal U1 = fetch(R, U1);
53  poly c1 = fetch(R, c1);
54  poly c2 = fetch(R, c2);
55  ideal U2 = changeDenomQ(U1, c1, c2);
56  setring R;
57  ideal U2 = fetch(Q, U2);
58  return(U2);
59}
60
61///////////////////////////////////////////////////////////////////////////////
62
63static proc changeDenomQ(ideal U1, poly c1, poly c2){
64// Given a ring in the form 1/c1 * U, it computes a new U2 st the ring
65// is 1/c2 * U2.
66// The base ring is already a quotient ring R / I.
67  int a;      // counter
68  ideal U2;
69  poly p;
70  for(a = 1; a <= ncols(U1); a++){
71    p = lift(c1, c2*U1[a])[1,1];
72    U2[a] = p;
73  }
74  return(U2);
75}
76
77/////////////////////////////////////////////////////////////////////////////////
78
79proc locNormal(ideal I, list #)
80"USAGE:  locNormal(I [,options]); I = prime ideal, options = list of options. @*
81         Optional parameters in list options (can be entered in any order):@*
82         modular: use a modular approach for the local computations. The number of primes is
83                  increased one at a time, starting with 2 primes, until the result stabelizes.@*
84         noVerificication: if the modular approach is used, the result will not be verified.
85ASSUME:  I is a prime ideal (the algorithm will also work for radical ideals as long as the
86         normal command does not detect that the ideal under consideration is not prime).
87RETURN:  a list of an ideal U and a universal denominator d such that U/d is the normalization.
88REMARKS: We use the local-to-global algorithm given in [1] to compute the normalization of
89         A = R/I, where R is the basering.@*
90         The idea is to stratify the singular locus of A, apply the normalization algorithm
91         given in [2] locally at each stratum, and put the local results together.@*
92         If the option modular is given, the result is returned as a probabilistic result
93         or verified, depending on whether the option noVerificication is used or not.@*
94         The normalization of A is represented as an R-module by returning a list of U and d,
95         where U is an ideal of A and d is an element of A such that U/d is the normalization of A.
96         In fact, U and d are returned as an ideal and a polynomial of the base ring R.
97
98         References:
99
100         [1] Janko Boehm, Wolfram Decker, Santiago Laplagne, Gerhard Pfister, Stefan Steidel,
101             Andreas Steenpass: Parallel algorithms for normalization, http://arxiv.org/abs/1110.4299, 2011.@*
102         [2] Gert-Martin Greuel, Santiago Laplagne, Frank Seelisch: Normalization of Rings,
103             Journal of Symbolic Computation 9 (2010), p. 887-901
104KEYWORDS: normalization; local methods; modular methods.
105SEE ALSO: normal_lib, modnormal_lib.
106EXAMPLE: example locNormal; shows an example
107"
108{
109// Computes the normalization by localizing at the different components of the singularity.
110  int i;
111  int totalLocalTime;
112  int dbg = printlevel - voice + 2;
113  def R = basering;
114
115  int totalTime = timer;
116  intvec LTimer;
117  int t;
118  int printTimings=0;
119
120  int locmod;
121  for ( i=1; i <= size(#); i++ )
122  {
123    if ( typeof(#[i]) == "string" )
124    {
125      if (#[i]=="modular") { locmod = 1;}
126      if (#[i]=="printTimings") { printTimings = 1;}
127    }
128  }
129
130  // Computes the Singular Locus
131  list IM = mstd(I);
132  I = IM[1];
133  int d = dim(I);
134  ideal IMin = IM[2];
135  qring Q = I;   // We work in the quotient by the groebner base of the ideal I
136  option("redSB");
137  option("returnSB");
138  ideal I = fetch(R, I);
139  attrib(I, "isSB", 1);
140  ideal IMin = fetch(R, IMin);
141  dbprint(dbg, "Computing the jacobian ideal...");
142  ideal J = minor(jacob(IMin), nvars(basering) - d, I);
143  t=timer;
144  J = modStd(J);
145  if (printTimings==1) {"Time for modStd Jacobian "+string(timer-t);}
146
147  setring R;
148  ideal J = fetch(Q, J);
149  // We compute a universal denominator
150  poly condu = getSmallest(J);
151  J = J, I;
152  if(dbg >= 2){
153    "Conductor: ", condu;
154    "The original singular locus is";
155    groebner(J);
156    if(dbg >= 2){pause();}
157    "";
158  }
159  t=timer;
160  list pd = locIdeals(J);
161  dbprint(dbg,pd);
162  if (printTimings==1) {
163    "Number of maximal components to localize at: ", size(pd);
164    "";
165  }
166
167  ideal U;
168  ideal resT;
169  ideal resu;
170  poly denomOld;
171  poly denom;
172  totalLocalTime = timer;
173  int maxLocalTime;
174  list Lnor;
175  list parallelArguments;
176  for(i = 1; i <= size(pd); i++){
177    parallelArguments[i] = list(pd[i], I, condu, i, locmod, printTimings, #);
178  }
179  list parallelResults = parallelWaitAll("locNormal_parallelTask",
180    parallelArguments);
181  for(i = 1; i <= size(pd); i++){
182    // We sum the result to the previous results.
183    resu = resu, parallelResults[i][1];
184    Lnor[i] = parallelResults[i][1];
185    if(parallelResults[i][2] > maxLocalTime) {
186      maxLocalTime = parallelResults[i][2];
187    }
188    LTimer[i] = parallelResults[i][2];
189  }
190  if (printTimings==1) {
191    "List of local times: "; LTimer;
192    "Maximal local time: "+string(maxLocalTime);
193  }
194  totalLocalTime = timer - totalLocalTime;
195  if (printTimings==1) {
196    "Total time local computations: "+string(totalLocalTime);
197  }
198
199  t=timer;
200  resu = modStd(resu);
201  if (printTimings==1) {
202     "Time for combining the local results, modStd "+string(timer-t);
203  }
204
205  totalTime = timer - totalTime;
206  if (printTimings==1) {
207    "Total time locNormal: "+string(totalTime);
208    "Simulated parallel time: "+string(totalTime + maxLocalTime - totalLocalTime);
209  }
210  return(list(resu, condu));
211}
212
213example
214{ "EXAMPLE:";
215ring R = 0,(x,y,z),dp;
216int k = 4;
217poly f = (x^(k+1)+y^(k+1)+z^(k+1))^2-4*(x^(k+1)*y^(k+1)+y^(k+1)*z^(k+1)+z^(k+1)*x^(k+1));
218f = subst(f,z,3x-2y+1);
219ring S = 0,(x,y),dp;
220poly f = imap(R,f);
221ideal i = f;
222list L = locNormal(i);
223}
224
225proc locNormal_parallelTask(ideal pdi, ideal I, poly condu, int i, int locmod,
226  int printTimings, list #)
227{
228  pdi=pdi;
229  int t = timer;
230  list opt = list(list("inputJ", pdi)) + #;
231  if (printTimings==1) {
232    "Local component ",i," of degree "+string(deg(pdi))+" and dimension "
233      +string(dim(pdi));
234  }
235  list n;
236  ideal norT;
237  poly denomT;
238  if (locmod==1) {
239     n = modNormal(I,1, opt);
240     norT = n[1];
241     denomT = n[2];
242  } else {
243     n = normal(I,opt);
244     if(size(n[2]) > 1){
245       ERROR("Input was not prime...");
246     }
247     norT = n[2][1];
248     denomT = norT[size(norT)];
249  }
250  t = timer-t;
251  // We compute the normalization of I localized at a component of the Singular
252  // Locus
253  if (printTimings==1) {
254    "Output of normalization of component ", i, ": "; norT;
255    "";
256  }
257  ideal nor = changeDenom(norT, denomT, condu, I);
258  return(list(nor, t));
259}
260
261static proc locComps(list l)
262{
263  int d = size(l);
264  int i;
265  int j;
266  intvec m = 1:d;   // 1 = maximal ideal
267  ideal IT;
268
269  // Check for maximal ideals
270  for(i = 1; i<d; i++)
271  {
272    for(j = i+1; j <= d; j++)
273    {
274      if(subset(l[i], l[j]))
275      {
276        m[i] = 0;
277        break;
278      }
279    }
280  }
281  list outL;
282  for(i = 1; i<= d; i++)
283  {
284    if(m[i] == 1){
285      // Maximal ideal
286      IT = l[i];
287      for(j = 1; j <= d; j++){
288        if(j != i)
289        {
290          if(subset(l[j], l[i]))
291          {
292            IT = intersect(IT, l[j]);
293          }
294        }
295      }
296      outL = insert(outL, IT);
297    }
298  }
299  return(outL);
300}
301
302// I C J ??
303static proc subset(ideal I, ideal J)
304{
305  J = groebner(J);
306  return(size(reduce(I, J)) == 0);
307}
308
309
310// Computes the different localizations of the radical of I at all the points of the space.
311static proc locIdeals(ideal I){
312  int i, j;
313  I = groebner(I);
314
315  // Minimal associated primes of I
316  list l = minAssGTZ(I);
317  //"Total number of components of the Singular Locus: ", size(l);
318
319  int s = size(l);
320  int d = dim(I);
321  if (d==0) {return(l)};
322  intvec m = (1:d);
323  // 1 = maximal ideal
324  ideal IT;
325
326  // inters will contain all the different intersections of components of I.
327  // It is a list of list. The j-th list contains the intersections of j components of I.
328  list inters;
329  list compIndex;   // Indicate the index of the last component in the corresponding intersection
330                    // This is used to intersect it with the remaining components.
331  for(i = 1; i<= d; i++)
332  {
333    l[i] = groebner(l[i]);
334  }
335
336  // We add all the components to the list of intersections
337  inters[1] = l;
338  compIndex[1] = 1..s;
339
340  ideal J;
341  int e;
342  int a;
343
344  // Intersections of two or more components
345  for(e = 1; e <= d; e++)
346  {
347    inters[e+1]  = list();
348    a = 1;
349    for(i = 1; j <= size(inters[e]); i++)
350    {
351      for(j = compIndex[e][i] + 1; j <= s; j++)
352      {
353        J = l[j] + inters[e][i];
354        J = groebner(J);
355        if(J[1] != 1)
356        {
357          inters[e+1] = inters[e+1]+list(J);
358          if(a == 1)
359          {
360            compIndex[e+1] = intvec(j);
361          } else
362          {
363            compIndex[e+1][a] = j;
364          }
365          a++;
366        }
367      }
368    }
369  }
370
371  list ids;
372  for(e = 1; e <= d+1; e++)
373  {
374    ids = ids + inters[e];
375  }
376  return(locComps(ids));
377}
378
379///////////////////////////////////////////////////////////////////////////
380//
381//                            EXAMPLES
382//
383///////////////////////////////////////////////////////////////////////////
384/*
385// plane curves
386
387ring r24 = 0,(x,y,z),dp;
388int k = 2;
389poly f = (x^(k+1)+y^(k+1)+z^(k+1))^2-4*(x^(k+1)*y^(k+1)+y^(k+1)*z^(k+1)+z^(k+1)*x^(k+1));
390f = subst(f,z,2x-y+1);
391ring s24 = 0,(x,y),dp;
392poly f = imap(r24,f);
393ideal i = f;
394
395locNormal(i);
396//modNormal(i,1);
397
398
399ring r24 = 0,(x,y,z),dp;
400int k = 3;
401poly f = (x^(k+1)+y^(k+1)+z^(k+1))^2-4*(x^(k+1)*y^(k+1)+y^(k+1)*z^(k+1)+z^(k+1)*x^(k+1));
402f = subst(f,z,2x-y+1);
403ring s24 = 0,(x,y),dp;
404poly f = imap(r24,f);
405ideal i = f;
406
407locNormal(i);
408//modNormal(i,1,"noVerification");
409
410
411ring r24 = 0,(x,y,z),dp;
412int k = 4;
413poly f = (x^(k+1)+y^(k+1)+z^(k+1))^2-4*(x^(k+1)*y^(k+1)+y^(k+1)*z^(k+1)+z^(k+1)*x^(k+1));
414f = subst(f,z,2x-y+1);
415ring s24 = 0,(x,y),dp;
416poly f = imap(r24,f);
417ideal i = f;
418
419locNormal(i);
420//modNormal(i,1,"noVerification");
421
422
423ring r24 = 0,(x,y,z),dp;
424int k = 5;
425poly f = (x^(k+1)+y^(k+1)+z^(k+1))^2-4*(x^(k+1)*y^(k+1)+y^(k+1)*z^(k+1)+z^(k+1)*x^(k+1));
426f = subst(f,z,2x-y+1);
427ring s24 = 0,(x,y),dp;
428poly f = imap(r24,f);
429ideal i = f;
430
431locNormal(i);
432
433
434
435ring s24 = 0,(x,y),dp;
436int a=7;
437ideal i = ((x-1)^a-y^3)*((x+1)^a-y^3)*((x)^a-y^3)*((x-2)^a-y^3)*((x+2)^a-y^3)+y^15;
438
439locNormal(i);
440//modNormal(i,1);
441
442
443ring s24 = 0,(x,y),dp;
444int a=7;
445ideal i = ((x-1)^a-y^3)*((x+1)^a-y^3)*((x)^a-y^3)*((x-2)^a-y^3)*((x+2)^a-y^3)+y^15;
446
447locNormal(i);
448//modNormal(i,1);
449
450ring s24 = 0,(x,y),dp;
451int a=7;
452ideal i = ((x-1)^a-y^3)*((x+1)^a-y^3)*((x)^a-y^3)*((x-2)^a-y^3)*((x+2)^a-y^3)+y^15;
453
454locNormal(i);
455//modNormal(i,1,"noVerification");
456
457
458
459
460ring r=0,(x,y),dp;
461ideal i=9127158539954x10+3212722859346x8y2+228715574724x6y4-34263110700x4y6
462-5431439286x2y8-201803238y10-134266087241x8-15052058268x6y2+12024807786x4y4
463+506101284x2y6-202172841y8+761328152x6-128361096x4y2+47970216x2y4-6697080y6
464-2042158x4+660492x2y2-84366y4+2494x2-474y2-1;
465
466locNormal(i);
467//modNormal(i,1);
468
469
470// surfaces in A3
471
472
473ring r7 = 0,(x,y,t),dp;
474int a=11;
475ideal i = x*y*(x-y)*(x+y)*(y-1)*t+(x^a-y^2)*(x^10-(y-1)^2);
476locNormal(i);
477//modNormal(i,1,"noVerification");
478
479ring r7 = 0,(x,y,t),dp;
480int a=12;
481ideal i = x*y*(x-y)*(x+y)*(y-1)*t+(x^a-y^2)*(x^10-(y-1)^2);
482locNormal(i);
483//modNormal(i,1,"noVerification");
484
485
486ring r7 = 0,(x,y,t),dp;
487int a=13;
488ideal i = x*y*(x-y)*(x+y)*(y-1)*t+(x^a-y^2)*(x^10-(y-1)^2);
489
490locNormal(i);
491modNormal(i,1,"noVerification");
492
493
494ring r22 = 0,(x,y,z),dp;
495ideal i = z2-(y2-1234x3)^2*(15791x2-y3)*(1231y2-x2*(x+158))*(1357y5-3x11);
496
497locNormal(i);
498//modNormal(i,1,"noVerification");
499
500
501ring r22 = 0,(x,y,z),dp;
502ideal i = z2-(y2-1234x3)^3*(15791x2-y3)*(1231y2-x2*(x+158))*(1357y5-3x11);
503
504locNormal(i);
505//modNormal(i,1,"noVerification");
506
507
508ring r23 = 0,(x,y,z),dp;
509ideal i = z5-((13x-17y)*(5x2-7y3)*(3x3-2y2)*(19y2-23x2*(x+29)))^2;
510
511locNormal(i);
512//modNormal(i,1,"noVerification");
513
514
515// curve in A3
516
517ring r23 = 0,(x,y,z),dp;
518ideal i = z3-(19y2-23x2*(x+29))^2,x3-(11y2-13z2*(z+1));
519
520locNormal(i);
521//modNormal(i,1,"noVerification");
522
523
524ring r23 = 0,(x,y,z),dp;
525ideal i = z3-(19y2-23x2*(x+29))^2,x3-(11y2-13z2*(z+1))^2;
526
527locNormal(i);
528//modNormal(i,1,"noVerification");
529
530// surface in A4
531
532ring r23 = 0,(x,y,z,w),dp;
533ideal i = z2-(y3-123456w2)*(15791x2-y3)^2, w*z-(1231y2-x*(111x+158));
534
535
536locNormal(i);
537//modNormal(i,1,"noVerification");
538*/
539
Note: See TracBrowser for help on using the repository browser.