# source:git/Singular/LIB/pointid.lib@bb9471

spielwiese
Last change on this file since bb9471 was bb9471, checked in by Christian Eder, 14 years ago
Stefan Steidels library for computing groebner basis with the algorithm of Mora ---------------------------------------------------------------------- git-svn-id: file:///usr/local/Singular/svn/trunk@11623 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 16.5 KB
Line
1////////////////////////////////////////////////////////////////////////////////
2info="
3LIBRARY:  pointid.lib     Procedures for computing a factorized lex GB of
4                          the vanishing ideal of a set of points via the
5                          Axis-of-Evil Theorem (M.G. Marinari, T. Mora)
6
7AUTHOR:   Stefan Steidel, steidel@mathematik.uni-kl.de
8
9OVERVIEW:
10    The algorithm of Cerlienco-Mureddu [Marinari M.G., Mora T., A remark on a
11    remark by Macaulay or Enhancing Lazard Structural Theorem. Bull. of the
12    Iranian Math. Soc., 29 (2003), 103-145] associates to each ordered set of
13    points A:={a1,...,as} in K^n, ai:=(ai1,...,ain)
14@*      - a set of monomials N and
15@*      - a bijection phi: A --> N.
16    Here I(A):={f in K[x(1),...,x(n)]|f(ai)=0, for all 1<=i<=s} denotes the
17    vanishing ideal of A and N = Mon(x(1),...,x(n)) \ {LM(f)|f in I(A)} is the
18    set of monomials not in the leading ideal of I(A) (w.r.t. the lexicograph.
19    ordering with x(n)>...>x(1)). N is also called the set of non-monomials of
20    I(A). NOTE: #A = #N and N is a monomial basis of K[x(1..n)]/I(A).
21    In particular, this allows to deduce the set of corner-monomials, i.e. the
22    minimal basis M:={m1,...,mr}, m1<...<mr, of its associated monomial ideal
23    M(I(A)), such that M(I(A))= {k*mi|k in Mon(x(1),...,x(n)), mi in M} and
24    (by interpolation) the unique reduced lexicographical Groebner basis
25    G := {f1,...,fr} such that LM(fi)=mi for each i, i.e. I(A)=<G>.
26    Moreover, a variation of this algorithm allows to deduce a canonical linear
27    factorization of each element of such a Groebner basis in the sense ot the
28    Axis-of-Evil Theorem by M.G. Marinari and T. Mora. More precisely:
29    A combinatorial algorithm and interpolation allow to deduce polynomials
30@*
31@*                y_mdi = x(m) - g_mdi(x(1),...,x(m-1)),
32@*
33    i=1,...,r, m=1,...,n, d in a finite index-set F, satisfying
34@*
35@*                 fi = (product of y_mdi) modulo (f1,...,f(i-1))
36@*
37    where the product runs over all m=1,...,n and d in F.
38
39PROCEDURES:
40 nonMonomials(id);   non-monomials of the vanishing ideal id of a set of points
41 cornerMonomials(N); corner-monomials of the set of non-monomials N
42 facGBIdeal(id);     GB G of id and linear factors of each element of G
43";
44
45LIB "poly.lib";
46
47////////////////////////////////////////////////////////////////////////////////
48
49static proc subst1(id, int m)
50{
51//id = poly/ideal/list, substitute the first m variables occuring in id by 1
52
53  int i,j;
54  def I = id;
55  if(typeof(I) == "list")
56  {
57    for(j = 1; j <= size(I); j++)
58    {
59      for(i = 1; i <= m; i++)
60      {
61        I[j] = subst(I[j],var(i),1);
62      }
63    }
64    return(I);
65  }
66  else
67  {
68    for(i = 1; i <= m; i++)
69    {
70      I = subst(I,var(i),1);
71    }
72    return(I);
73  }
74}
75
76////////////////////////////////////////////////////////////////////////////////
77
78proc nonMonomials(id)
79"USAGE:  nonMonomials(id); id = <list of vectors> or <list of lists> or <module>
80         or <matrix>.
81         Let A= {a1,...,as} be a set of points in K^n, ai:=(ai1,...,ain), then
82         A can be given as
83@*         - a list of vectors (the ai are vectors) or
84@*         - a list of lists (the ai are lists of numbers) or
85@*         - a module s.t. the ai are generators or
86@*         - a matrix s.t. the ai are columns
87ASSUME:  basering must have ordering rp, i.e. of the form 0,x(1..n),rp;
88         (the first entry of a point belongs to the lex-smallest variable, etc.)
89RETURN:  ideal, the non-monomials of the vanishing ideal I(A) of A
90PURPOSE: compute the set of non-monomials Mon(x(1),...,x(n)) \ {LM(f)|f in I(A)}
91         of the vanishing ideal I(A) of the given set of points A in K^n, where
92         K[x(1),...,x(n)] is equipped with the lexicographical ordering induced
93         by x(1)<...<x(n) by using the algorithm of Cerlienco-Mureddu
94EXAMPLE: example nonMonomials; shows an example
95"
96{
97  list A;
98  int i,j;
99  if(typeof(id) == "list")
100  {
101    for(i = 1; i <= size(id); i++)
102    {
103      if(typeof(id[i]) == "list")
104      {
105        vector a;
106        for(j = 1; j <= size(id[i]); j++)
107        {
108          a = a+id[i][j]*gen(j);
109        }
110        A[size(A)+1] = a;
111        kill a;
112      }
113      if(typeof(id[i]) == "vector")
114      {
115        A[size(A)+1] = id[i];
116      }
117    }
118  }
119  else
120  {
121    if(typeof(id) == "module")
122    {
123      for(i = 1; i <= size(id); i++)
124      {
125        A[size(A)+1] = id[i];
126      }
127    }
128    else
129    {
130      if(typeof(id) == "matrix")
131      {
132        for(i = 1; i <= ncols(id); i++)
133        {
134          A[size(A)+1] = id[i];
135        }
136      }
137      else
138      {
139        ERROR("Wrong type of input!!");
140      }
141    }
142  }
143
144  int n = nvars(basering);
145  int s;
146  int m,d;
147  ideal N = 1;
148  ideal N1,N2;
149  list Y,D,W,Z;
150  poly my,t;
151  for(s = 2; s <= size(A); s++)
152  {
153
154//-- compute m = max{ j | ex. i<s: pr_j(a_i) = pr_j(a_s)} ---------------------
155//-- compute d = |{ a_i | i<s: pr_m(a_i) = pr_m(a_s), phi(a_i) in T[1,m+1]}| --
156    m = 0;
157    Y = A[1..s-1];
158    N2 = N[1..s-1];
159    while(size(Y) > 0)      //assume all points different (m <= size(basering))
160    {
161      D = Y;
162      N1 = N2;
163      Y = list();
164      N2 = ideal();
165      m++;
166      for(i = 1; i <= size(D); i++)
167      {
168        if(A[s][m] == D[i][m])
169        {
170          Y[size(Y)+1] = D[i];
171          N2[size(N2)+1] = N1[i];
172        }
173      }
174    }
175    m = m - 1;
176    d = size(D);
177    if(m < n-1)
178    {
179      for(i = 1; i <= size(N1); i++)
180      {
181        if(N1[i] >= var(m+2))
182        {
183          d = d - 1;
184        }
185      }
186    }
187
188//------- compute t = my * x(m+1)^d where my is obtained recursively --------
189    if(m == 0)
190    {
191      my = 1;
192    }
193    else
194    {
195      W = list();
196      Z = list();
197      for(i = 2; i <= s-1; i++)
198      {
199        if((leadexp(N[i])[m+1] == d) && (coeffs(N[i],var(m+1))[nrows(coeffs(N[i],var(m+1))),1] < var(m+1)))
200        {
201          W[size(W)+1] = A[i];
202          Z[size(Z)+1] = A[i][1..m];
203        }
204      }
205      W[size(W)+1] = A[s];
206      Z[size(Z)+1] = A[s][1..m];
207
208      my = nonMonomials(Z)[size(Z)];
209    }
210
211    t = my * var(m+1)^d;
212
213//---------------------------- t is added to N --------------------------------
214    N[size(N)+1] = t;
215  }
216  return(N);
217}
218example
219{ "EXAMPLE:"; echo = 2;
220   ring R1 = 0,x(1..3),rp;
221   vector a1 = [4,0,0];
222   vector a2 = [2,1,4];
223   vector a3 = [2,4,0];
224   vector a4 = [3,0,1];
225   vector a5 = [2,1,3];
226   vector a6 = [1,3,4];
227   vector a7 = [2,4,3];
228   vector a8 = [2,4,2];
229   vector a9 = [1,0,2];
230   list A = a1,a2,a3,a4,a5,a6,a7,a8,a9;
231   nonMonomials(A);
232
233   matrix MAT[9][3] = 4,0,0,2,1,4,2,4,0,3,0,1,2,1,3,1,3,4,2,4,3,2,4,2,1,0,2;
234   MAT = transpose(MAT);
235   print(MAT);
236   nonMonomials(MAT);
237
238   module MOD = gen(3),gen(2)-2*gen(3),2*gen(1)+2*gen(3),2*gen(2)-2*gen(3),gen(1)+3*gen(3),gen(1)+gen(2)+3*gen(3),gen(1)+gen(2)+gen(3);
239   print(MOD);
240   nonMonomials(MOD);
241
242   ring R2 = 0,x(1..2),rp;
243   list l1 = 0,0;
244   list l2 = 0,1;
245   list l3 = 2,0;
246   list l4 = 0,2;
247   list l5 = 1,0;
248   list l6 = 1,1;
249   list L = l1,l2,l3,l4,l5,l6;
250   nonMonomials(L);
251}
252
253////////////////////////////////////////////////////////////////////////////////
254
255proc cornerMonomials(ideal N)
256"USAGE:  cornerMonomials(N); N ideal
257ASSUME:  N is given by monomials satisfying the condition that if a monomial is
258         in N then any of its factors is in N (N is then called an order ideal)
259RETURN:  ideal, the corner-monomials of the order ideal N @*
260         The corner-monomials are the leading monomials of an ideal I s.t. N is
261         a basis of basering/I.
262NOTE:    In our applications I is the vanishing ideal of a finte set of points.
263EXAMPLE: example cornerMonomials; shows an example
264"
265{
266  int n = nvars(basering);
267  int i,j,h;
268  list C;
269  poly m,p;
270  int Z;
271  int vars;
272  intvec v;
273  ideal M;
274
275//-------------------- Test: 1 in N ?, if no, return <1> ----------------------
276  for(i = 1; i <= size(N); i++)
277  {
278    if(1 == N[i])
279    {
280      h = 1;
281      break;
282    }
283  }
284  if(h == 0)
285  {
286    return(ideal(1));
287  }
288
289//----------------------------- compute the set M -----------------------------
290  for(i = 1; i <= n; i++)
291  {
292    C[size(C)+1] = list(var(i),1);
293  }
294  while(size(C) > 0)
295  {
296    m = C[1][1];
297    Z = C[1][2];
298    C = delete(C,1);
299    vars = 0;
301    for(i = 1; i <= n; i++)
302    {
303      vars = vars + (v[i] != 0);
304    }
305
306    if(vars == Z)
307    {
308      h = 0;
309      for(i = 1; i <= size(N); i++)
310      {
311        if(m == N[i])
312        {
313          h = 1;
314          break;
315        }
316      }
317
318      if(h == 0)
319      {
320        M[size(M)+1] = m;
321      }
322      else
323      {
324        for(i = 1; i <= n; i++)
325        {
326          p = m * var(i);
327          for(j = 1; j <= size(C); j++)
328          {
329            if(p <= C[j][1] || j == size(C))
330            {
331              if(p == C[j][1])
332              {
333                C[j][2] = C[j][2] + 1;
334              }
335              else
336              {
337                if(p < C[j][1])
338                {
339                  C = insert(C,list(p,1),j-1);
340                }
341                else
342                {
343                  C[size(C)+1] = list(p,1);
344                }
345              }
346              break;
347            }
348          }
349        }
350      }
351    }
352  }
353  return(M);
354}
355example
356{ "EXAMPLE:"; echo = 2;
357   ring R = 0,x(1..3),rp;
358   poly n1 = 1;
359   poly n2 = x(1);
360   poly n3 = x(2);
361   poly n4 = x(1)^2;
362   poly n5 = x(3);
363   poly n6 = x(1)^3;
364   poly n7 = x(2)*x(3);
365   poly n8 = x(3)^2;
366   poly n9 = x(1)*x(2);
367   ideal N = n1,n2,n3,n4,n5,n6,n7,n8,n9;
368   cornerMonomials(N);
369}
370
371////////////////////////////////////////////////////////////////////////////////
372
373proc facGBIdeal(id)
374"USAGE:  facGBIdeal(id); id = <list of vectors> or <list of lists> or <module>
375         or <matrix>.
376         Let A= {a1,...,as} be a set of points in K^n, ai:=(ai1,...,ain), then
377         A can be given as
378@*         - a list of vectors (the ai are vectors) or
379@*         - a list of lists (the ai are lists of numbers) or
380@*         - a module s.t. the ai are generators or
381@*         - a matrix s.t. the ai are columns
382ASSUME:  basering must have ordering rp, i.e. of the form 0,x(1..n),rp;
383         (the first entry of a point belongs to the lex-smallest variable, etc.)
384RETURN:  a list where the first entry contains the Groebner basis G of I(A)
385         and the second entry contains the linear factors of each element of G
386NOTE:    combinatorial algorithm due to the Axis-of-Evil Theorem of M.G.
387         Marinari, T. Mora
388EXAMPLE: example facGBIdeal; shows an example
389"
390{
391  list A;
392  int i,j;
393  if(typeof(id) == "list")
394  {
395    for(i = 1; i <= size(id); i++)
396    {
397      if(typeof(id[i]) == "list")
398      {
399        vector a;
400        for(j = 1; j <= size(id[i]); j++)
401        {
402          a = a+id[i][j]*gen(j);
403        }
404        A[size(A)+1] = a;
405        kill a;
406      }
407      if(typeof(id[i]) == "vector")
408      {
409        A[size(A)+1] = id[i];
410      }
411    }
412  }
413  else
414  {
415    if(typeof(id) == "module")
416    {
417      for(i = 1; i <= size(id); i++)
418      {
419        A[size(A)+1] = id[i];
420      }
421    }
422    else
423    {
424      if(typeof(id) == "matrix")
425      {
426        for(i = 1; i <= ncols(id); i++)
427        {
428          A[size(A)+1] = id[i];
429        }
430      }
431      else
432      {
433        ERROR("Wrong type of input!!");
434      }
435    }
436  }
437
438  int n = nvars(basering);
439  def S = basering;
440  def R;
441
442  ideal N = nonMonomials(A);
443  ideal eN;
444  ideal M = cornerMonomials(N);
445  poly my, emy;
446
447  int d,k,l,m;
448
449  int d1;
450  poly y(1);
451
452  list N2,D,H;
453  poly z,h;
454
455  int dm;
456  list Am,Z;
457  ideal E;
458  list V,eV;
459  poly p;
460  poly y(2..n),y;
461  poly xi;
462
463  poly f;
464  ideal G1;          // stores the elements of G, i.e. G1 = G the GB of I(A)
465  ideal Y;           // stores the linear factors of GB-elements in each slope
466  list G2;           // contains the linear factors of each element of G
467
468  for(j = 1; j <= size(M); j++)
469  {
470    my = M[j];
471    emy = subst(my,var(1),1);         // auxiliary polynomial
472    eN = subst(N,var(1),1);           // auxiliary monomial ideal
473    Y = ideal();
474
476    y(1) = 1;
477    i = 0;
478    k = 1;
479    while(i < d1)
480    {
481//---------- searching for phi^{-1}(x_1^i*x_2^d_2*...*x_n^d_n) ----------------
482      while(my*var(1)^i/var(1)^d1 != N[k])
483      {
484        k++;
485      }
486      y(1) = y(1)*(var(1)-A[k][1]);
487      Y[size(Y)+1] = cleardenom(var(1)-A[k][1]);
488      i++;
489    }
490    f = y(1);                             // gamma_1my
491
492//--------------- Recursion over number of variables --------------------------
493    z = 1;                                // zeta_mmy
494    for(m = 2; m <= n; m++)
495    {
496      z = z * y(m-1);
497
498      D = list();
499      H = list();
500      for(i = 1; i <= size(A); i++)
501      {
502        h = z;
503        for(k = 1; k <= n; k++)
504        {
505          h = subst(h,var(k),A[i][k]);
506        }
507        if(h != 0)
508        {
509          D[size(D)+1] = A[i];
510          H[size(H)+1] = i;
511        }
512      }
513
514      if(size(D) == 0)
515      {
516        break;
517      }
518
520      while(dm == 0)
521      {
522        m = m + 1;
524      }
525
526      N2 = list();                        // N2 = N_m
527      emy = subst(emy,var(m),1);
528      eN = subst(eN,var(m),1);
529      for(i = 1; i <= size(N); i++)
530      {
531        if((emy == eN[i]) && (my > N[i]))
532        {
533          N2[size(N2)+1] = N[i];
534        }
535      }
536
537      y(m) = 1;
538      xi = z;
539      for(d = 1; d <= dm; d++)
540      {
541        Am = list();
542        Z = list();                       // Z = pr_m(Am)
543
544//------- V contains all ny*x_m^{d_m-d}*x_m+1^d_m+1*...+x_n^d_n in N_m --------
545        eV = subst1(N2,m-1);
546        V = list();
547        for(i = 1; i <= size(eV); i++)
548        {
549          if(eV[i] == subst1(my,m-1)/var(m)^d)
550          {
551            V[size(V)+1] = eV[i];
552          }
553        }
554
555//------- A_m = phi^{-1}(V) intersect D_md-1 ----------------------------------
556        for(i = 1; i <= size(D); i++)
557        {
558          p = N[H[i]];
559          p = subst1(p,m-1);
560          for(l = 1; l <= size(V); l++)
561          {
562            if(p == V[l])
563            {
564              Am[size(Am)+1] = D[i];
565              Z[size(Z)+1] = D[i][1..m];
566              break;
567            }
568          }
569        }
570
571        E = nonMonomials(Z);
572
573        R = extendring(size(E), "c(", "lp");
574        setring R;
575        ideal E = imap(S,E);
576        list Am = imap(S,Am);
577        poly g = var(size(E)+m);
578        for(i = 1; i <= size(E); i++)
579        {
580          g = g + c(i)*E[i];
581        }
582
583        ideal I = ideal();
584        poly h;
585        for (i = 1; i <= size(Am); i++)
586        {
587          h = g;
588          for(k = 1; k <= n; k++)
589          {
590            h = subst(h,var(size(E)+k),Am[i][k]);
591          }
592          I[size(I)+1] = h;
593        }
594
595        ideal sI = std(I);
596        g = reduce(g,sI);
597
598        setring S;
599        y = imap(R,g);
600        Y[size(Y)+1] = cleardenom(y);
601        xi = xi * y;
602
603        D = list();
604        H = list();
605        for(i = 1; i <= size(A); i++)
606        {
607          h = xi;
608          for(k = 1; k <= n; k++)
609          {
610            h = subst(h,var(k),A[i][k]);
611          }
612          if(h != 0)
613          {
614            D[size(D)+1] = A[i];
615            H[size(H)+1] = i;
616          }
617        }
618
619        y(m) = y(m) * y;
620
621        if(size(D) == 0)
622        {
623          break;
624        }
625      }
626
627      f = f * y(m);
628    }
629    G1[size(G1)+1] = f;
630    G2[size(G2)+1] = Y;
631  }
632  return(list(G1,G2));
633}
634example
635{ "EXAMPLE:"; echo = 2;
636   ring R = 0,x(1..3),rp;
637   vector a1 = [4,0,0];
638   vector a2 = [2,1,4];
639   vector a3 = [2,4,0];
640   vector a4 = [3,0,1];
641   vector a5 = [2,1,3];
642   vector a6 = [1,3,4];
643   vector a7 = [2,4,3];
644   vector a8 = [2,4,2];
645   vector a9 = [1,0,2];
646   list A = a1,a2,a3,a4,a5,a6,a7,a8,a9;
647   facGBIdeal(A);
648
649   matrix MAT[9][3] = 4,0,0,2,1,4,2,4,0,3,0,1,2,1,3,1,3,4,2,4,3,2,4,2,1,0,2;
650   MAT = transpose(MAT);
651   print(MAT);
652   facGBIdeal(MAT);
653
654   module MOD = gen(3),gen(2)-2*gen(3),2*gen(1)+2*gen(3),2*gen(2)-2*gen(3),gen(1)+3*gen(3),gen(1)+gen(2)+3*gen(3),gen(1)+gen(2)+gen(3);
655   print(MOD);
656   facGBIdeal(MOD);
657
658   list l1 = 0,0,1;
659   list l2 = 0,1,-2;
660   list l3 = 2,0,2;
661   list l4 = 0,2,-2;
662   list l5 = 1,0,3;
663   list l6 = 1,1,3;
664   list L = l1,l2,l3,l4,l5,l6;
665   facGBIdeal(L);
666}
Note: See TracBrowser for help on using the repository browser.