1 | //////////////////////////////////////////////////////////////////////////////// |
---|
2 | info=" |
---|
3 | LIBRARY: 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 | |
---|
7 | AUTHOR: Stefan Steidel, steidel@mathematik.uni-kl.de |
---|
8 | |
---|
9 | OVERVIEW: |
---|
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 | |
---|
39 | PROCEDURES: |
---|
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 | |
---|
45 | LIB "poly.lib"; |
---|
46 | |
---|
47 | //////////////////////////////////////////////////////////////////////////////// |
---|
48 | |
---|
49 | static 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 | |
---|
78 | proc 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 |
---|
87 | ASSUME: 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.) |
---|
89 | RETURN: ideal, the non-monomials of the vanishing ideal I(A) of A |
---|
90 | PURPOSE: 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 |
---|
94 | EXAMPLE: 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 | } |
---|
218 | example |
---|
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 | |
---|
255 | proc cornerMonomials(ideal N) |
---|
256 | "USAGE: cornerMonomials(N); N ideal |
---|
257 | ASSUME: 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) |
---|
259 | RETURN: 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. |
---|
262 | NOTE: In our applications I is the vanishing ideal of a finte set of points. |
---|
263 | EXAMPLE: 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; |
---|
300 | v = leadexp(m); |
---|
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 | } |
---|
355 | example |
---|
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 | |
---|
373 | proc 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 |
---|
382 | ASSUME: 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.) |
---|
384 | RETURN: 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 |
---|
386 | NOTE: combinatorial algorithm due to the Axis-of-Evil Theorem of M.G. |
---|
387 | Marinari, T. Mora |
---|
388 | EXAMPLE: 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 | |
---|
475 | d1 = leadexp(my)[1]; |
---|
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 | |
---|
519 | dm = leadexp(my)[m]; |
---|
520 | while(dm == 0) |
---|
521 | { |
---|
522 | m = m + 1; |
---|
523 | dm = leadexp(my)[m]; |
---|
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 | } |
---|
634 | example |
---|
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 | } |
---|