# source:git/Singular/LIB/pointid.lib@7de8e4

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