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