1 | /////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="version modstd.lib 4.2.0.0 Feb_2020 "; // $Id$ |
---|
3 | category="Commutative Algebra"; |
---|
4 | info=" |
---|
5 | LIBRARY: modstd.lib Groebner bases of ideals/modules using modular methods |
---|
6 | |
---|
7 | AUTHORS: A. Hashemi Amir.Hashemi@lip6.fr |
---|
8 | G. Pfister pfister@mathematik.uni-kl.de |
---|
9 | H. Schoenemann hannes@mathematik.uni-kl.de |
---|
10 | A. Steenpass steenpass@mathematik.uni-kl.de |
---|
11 | S. Steidel steidel@mathematik.uni-kl.de |
---|
12 | |
---|
13 | OVERVIEW: |
---|
14 | A library for computing Groebner bases of ideals/modules in the polynomial ring |
---|
15 | over the rational numbers using modular methods. |
---|
16 | |
---|
17 | REFERENCES: |
---|
18 | E. A. Arnold: Modular algorithms for computing Groebner bases. |
---|
19 | J. Symb. Comp. 35, 403-419 (2003). |
---|
20 | |
---|
21 | N. Idrees, G. Pfister, S. Steidel: Parallelization of Modular Algorithms. |
---|
22 | J. Symb. Comp. 46, 672-684 (2011). |
---|
23 | |
---|
24 | PROCEDURES: |
---|
25 | modStd(I); standard basis of I using modular methods |
---|
26 | modSyz(I); syzygy module of I using modular methods |
---|
27 | modIntersect(I,J); intersection of I and J using modular methods |
---|
28 | "; |
---|
29 | |
---|
30 | LIB "polylib.lib"; |
---|
31 | LIB "modular.lib"; |
---|
32 | |
---|
33 | proc modStd(def I, list #) |
---|
34 | "USAGE: modStd(I[, exactness]); I ideal/module, exactness int |
---|
35 | RETURN: a standard basis of I |
---|
36 | NOTE: The procedure computes a standard basis of I (over the rational |
---|
37 | numbers) by using modular methods. |
---|
38 | @* An optional parameter 'exactness' can be provided. |
---|
39 | If exactness = 1(default), the procedure computes a standard basis |
---|
40 | of I for sure; if exactness = 0, it computes a standard basis of I |
---|
41 | with high probability. |
---|
42 | SEE ALSO: modular |
---|
43 | EXAMPLE: example modStd; shows an example" |
---|
44 | { |
---|
45 | /* read optional parameter */ |
---|
46 | int exactness = 1; |
---|
47 | if (size(#) > 0) |
---|
48 | { |
---|
49 | /* For compatibility, we only test size(#) > 4. This can be changed to |
---|
50 | * size(#) > 1 in the future. */ |
---|
51 | if (size(#) > 4 || typeof(#[1]) != "int") |
---|
52 | { |
---|
53 | ERROR("wrong optional parameter"); |
---|
54 | } |
---|
55 | exactness = #[1]; |
---|
56 | } |
---|
57 | |
---|
58 | /* save options */ |
---|
59 | intvec opt = option(get); |
---|
60 | option(redSB); |
---|
61 | |
---|
62 | /* choose the right command */ |
---|
63 | string command = "groebner"; |
---|
64 | if (npars(basering) > 0) { command = "Modstd::groebner_norm"; } |
---|
65 | |
---|
66 | /* call modular() */ |
---|
67 | if (exactness) |
---|
68 | { |
---|
69 | I = modular(command, list(I), primeTest_std, |
---|
70 | deleteUnluckyPrimes_std, pTest_std, finalTest_std); |
---|
71 | } |
---|
72 | else |
---|
73 | { |
---|
74 | I = modular(command, list(I), primeTest_std, |
---|
75 | deleteUnluckyPrimes_std, pTest_std); |
---|
76 | } |
---|
77 | |
---|
78 | /* return the result */ |
---|
79 | attrib(I, "isSB", 1); |
---|
80 | option(set, opt); |
---|
81 | return(I); |
---|
82 | } |
---|
83 | example |
---|
84 | { |
---|
85 | "EXAMPLE:"; |
---|
86 | echo = 2; |
---|
87 | ring R1 = 0, (x,y,z,t), dp; |
---|
88 | ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; |
---|
89 | ideal J = modStd(I); |
---|
90 | J; |
---|
91 | I = homog(I, t); |
---|
92 | J = modStd(I); |
---|
93 | J; |
---|
94 | |
---|
95 | ring R2 = 0, (x,y,z), ds; |
---|
96 | ideal I = jacob(x5+y6+z7+xyz); |
---|
97 | ideal J = modStd(I, 0); |
---|
98 | J; |
---|
99 | |
---|
100 | ring R3 = 0, x(1..4), lp; |
---|
101 | ideal I = cyclic(4); |
---|
102 | ideal J1 = modStd(I, 1); // default |
---|
103 | ideal J2 = modStd(I, 0); |
---|
104 | size(reduce(J1, J2)); |
---|
105 | size(reduce(J2, J1)); |
---|
106 | } |
---|
107 | |
---|
108 | proc modSyz(def I) |
---|
109 | "USAGE: modSyz(I); I ideal/module |
---|
110 | RETURN: a generating set of syzygies of I |
---|
111 | NOTE: The procedure computes a the syzygy module of I (over the rational |
---|
112 | numbers) by using modular methods with high probability. |
---|
113 | The property of being a syzygy is tested. |
---|
114 | SEE ALSO: modular |
---|
115 | EXAMPLE: example modSyz; shows an example" |
---|
116 | { |
---|
117 | /* save options */ |
---|
118 | intvec opt = option(get); |
---|
119 | option(redSB); |
---|
120 | |
---|
121 | /* choose the right command */ |
---|
122 | string command = "syz"; |
---|
123 | |
---|
124 | /* call modular() */ |
---|
125 | module M = modular(command, list(I), primeTest_std, |
---|
126 | deleteUnluckyPrimes_std, pTest_syz); |
---|
127 | |
---|
128 | /* return the result */ |
---|
129 | option(set, opt); |
---|
130 | return(M); |
---|
131 | } |
---|
132 | example |
---|
133 | { |
---|
134 | "EXAMPLE:"; echo = 2; |
---|
135 | ring R1 = 0, (x,y,z,t), dp; |
---|
136 | ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; |
---|
137 | modSyz(I); |
---|
138 | simplify(syz(I),1); |
---|
139 | } |
---|
140 | |
---|
141 | proc modIntersect(def I, def J) |
---|
142 | "USAGE: modIntersect(I,J); I,J ideal/module |
---|
143 | RETURN: a generating set of the intersection of I and J |
---|
144 | NOTE: The procedure computes a the intersection of I and J |
---|
145 | (over the rational numbers) by using modular methods |
---|
146 | with high probability. |
---|
147 | No additional tests are performed. |
---|
148 | SEE ALSO: modular |
---|
149 | EXAMPLE: example modIntersect; shows an example" |
---|
150 | { |
---|
151 | /* save options */ |
---|
152 | intvec opt = option(get); |
---|
153 | option(redSB); |
---|
154 | |
---|
155 | /* choose the right command */ |
---|
156 | string command = "intersect"; |
---|
157 | |
---|
158 | /* call modular() */ |
---|
159 | def M = modular(command, list(I,J), primeTest_std, |
---|
160 | deleteUnluckyPrimes_std); |
---|
161 | |
---|
162 | /* return the result */ |
---|
163 | option(set, opt); |
---|
164 | return(M); |
---|
165 | } |
---|
166 | example |
---|
167 | { |
---|
168 | "EXAMPLE:"; echo = 2; |
---|
169 | ring R1 = 0, (x,y,z,t), dp; |
---|
170 | ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; |
---|
171 | ideal J = maxideal(2); |
---|
172 | modIntersect(I,J); |
---|
173 | simplify(intersect(I,J),1); |
---|
174 | } |
---|
175 | |
---|
176 | /* compute a normalized GB via groebner() */ |
---|
177 | static proc groebner_norm(def I) |
---|
178 | { |
---|
179 | I = simplify(groebner(I), 1); |
---|
180 | attrib(I, "isSB", 1); |
---|
181 | return(I); |
---|
182 | } |
---|
183 | |
---|
184 | /* test if the prime p is suitable for the input, i.e. it does not divide |
---|
185 | * the numerator or denominator of any of the coefficients */ |
---|
186 | static proc primeTest_std(int p, alias list args) |
---|
187 | { |
---|
188 | /* erase zero generators */ |
---|
189 | def I = simplify(args[1], 2); |
---|
190 | |
---|
191 | /* clear denominators and count the terms */ |
---|
192 | def J=I; // dummy assign, to get the type of I |
---|
193 | ideal K; |
---|
194 | int n = ncols(I); |
---|
195 | intvec sizes; |
---|
196 | number cnt; |
---|
197 | int i; |
---|
198 | for(i = n; i > 0; i--) |
---|
199 | { |
---|
200 | J[i] = cleardenom(I[i]); |
---|
201 | cnt = leadcoef(J[i])/leadcoef(I[i]); |
---|
202 | K[i] = numerator(cnt)*var(1)+denominator(cnt); |
---|
203 | } |
---|
204 | sizes = size(J[1..n]); |
---|
205 | |
---|
206 | /* change to characteristic p */ |
---|
207 | def br = basering; |
---|
208 | list lbr = ringlist(br); |
---|
209 | if (typeof(lbr[1]) == "int") { lbr[1] = p; } |
---|
210 | else { lbr[1][1] = p; } |
---|
211 | def rp = ring(lbr); |
---|
212 | setring(rp); |
---|
213 | def Jp = fetch(br, J); |
---|
214 | ideal Kp = fetch(br, K); |
---|
215 | |
---|
216 | /* test if any coefficient is missing */ |
---|
217 | if (intvec(size(Kp[1..n])) != 2:n) { setring(br); return(0); } |
---|
218 | if (intvec(size(Jp[1..n])) != sizes) { setring(br); return(0); } |
---|
219 | setring(br); |
---|
220 | return(1); |
---|
221 | } |
---|
222 | |
---|
223 | /* find entries in modresults which come from unlucky primes. |
---|
224 | * For this, sort the entries into categories depending on their leading |
---|
225 | * ideal and return the indices in all but the biggest category. */ |
---|
226 | static proc deleteUnluckyPrimes_std(alias list modresults) |
---|
227 | { |
---|
228 | int size_modresults = size(modresults); |
---|
229 | |
---|
230 | /* sort results into categories. |
---|
231 | * each category is represented by three entries: |
---|
232 | * - the corresponding leading ideal |
---|
233 | * - the number of elements |
---|
234 | * - the indices of the elements |
---|
235 | */ |
---|
236 | list cat; |
---|
237 | int size_cat; |
---|
238 | def L=modresults[1]; // dummy assign to get the type of L |
---|
239 | int i; |
---|
240 | int j; |
---|
241 | for (i = 1; i <= size_modresults; i++) |
---|
242 | { |
---|
243 | L = lead(modresults[i]); |
---|
244 | attrib(L, "isSB", 1); |
---|
245 | for (j = 1; j <= size_cat; j++) |
---|
246 | { |
---|
247 | if (size(L) == size(cat[j][1]) |
---|
248 | && size(reduce(L, cat[j][1], 5)) == 0 |
---|
249 | && size(reduce(cat[j][1], L, 5)) == 0) |
---|
250 | { |
---|
251 | cat[j][2] = cat[j][2]+1; |
---|
252 | cat[j][3][cat[j][2]] = i; |
---|
253 | break; |
---|
254 | } |
---|
255 | } |
---|
256 | if (j > size_cat) |
---|
257 | { |
---|
258 | size_cat++; |
---|
259 | cat[size_cat] = list(); |
---|
260 | cat[size_cat][1] = L; |
---|
261 | cat[size_cat][2] = 1; |
---|
262 | cat[size_cat][3] = list(i); |
---|
263 | } |
---|
264 | } |
---|
265 | |
---|
266 | /* find the biggest categories */ |
---|
267 | int cat_max = 1; |
---|
268 | int max = cat[1][2]; |
---|
269 | for (i = 2; i <= size_cat; i++) |
---|
270 | { |
---|
271 | if (cat[i][2] > max) |
---|
272 | { |
---|
273 | cat_max = i; |
---|
274 | max = cat[i][2]; |
---|
275 | } |
---|
276 | } |
---|
277 | |
---|
278 | /* return all other indices */ |
---|
279 | list unluckyIndices; |
---|
280 | for (i = 1; i <= size_cat; i++) |
---|
281 | { |
---|
282 | if (i != cat_max) { unluckyIndices = unluckyIndices + cat[i][3]; } |
---|
283 | } |
---|
284 | return(unluckyIndices); |
---|
285 | } |
---|
286 | //////////////////////////////////////////////////////////////////////////////// |
---|
287 | |
---|
288 | static proc cleardenomModule(def I) |
---|
289 | { |
---|
290 | int t=ncols(I); |
---|
291 | if(size(I)==0) |
---|
292 | { |
---|
293 | return(I); |
---|
294 | } |
---|
295 | else |
---|
296 | { |
---|
297 | for(int i=1;i<=t;i++) |
---|
298 | { |
---|
299 | I[i]=cleardenom(I[i]); |
---|
300 | } |
---|
301 | } |
---|
302 | return(I); |
---|
303 | } |
---|
304 | |
---|
305 | static proc pTest_syz(string command, alias list args, alias def result, int p) |
---|
306 | { |
---|
307 | module result_without_denom=cleardenomModule(result); |
---|
308 | return(size(module(matrix(args[1])*matrix(result_without_denom)))==0); |
---|
309 | } |
---|
310 | |
---|
311 | /* test if 'command' applied to 'args' in characteristic p is the same as |
---|
312 | 'result' mapped to characteristic p */ |
---|
313 | static proc pTest_std(string command, alias list args, alias def result, |
---|
314 | int p) |
---|
315 | { |
---|
316 | /* change to characteristic p */ |
---|
317 | def br = basering; |
---|
318 | list lbr = ringlist(br); |
---|
319 | if (typeof(lbr[1]) == "int") { lbr[1] = p; } |
---|
320 | else { lbr[1][1] = p; } |
---|
321 | def rp = ring(lbr); |
---|
322 | setring(rp); |
---|
323 | def Ip = fetch(br, args)[1]; |
---|
324 | def Gp = fetch(br, result); |
---|
325 | attrib(Gp, "isSB", 1); |
---|
326 | |
---|
327 | /* test if Ip is in Gp */ |
---|
328 | int i; |
---|
329 | for (i = ncols(Ip); i > 0; i--) |
---|
330 | { |
---|
331 | if (reduce(Ip[i], Gp, 1) != 0) |
---|
332 | { |
---|
333 | setring(br); |
---|
334 | return(0); |
---|
335 | } |
---|
336 | } |
---|
337 | |
---|
338 | /* compute command(args) */ |
---|
339 | execute("Ip = "+command+"(Ip);"); |
---|
340 | |
---|
341 | /* test if Gp is in Ip */ |
---|
342 | for (i = ncols(Gp); i > 0; i--) |
---|
343 | { |
---|
344 | if (reduce(Gp[i], Ip, 1) != 0) { setring(br); return(0); } |
---|
345 | } |
---|
346 | setring(br); |
---|
347 | return(1); |
---|
348 | } |
---|
349 | |
---|
350 | /* test if 'result' is a GB of the input ideal */ |
---|
351 | static proc finalTest_std(string command, alias list args, def result) |
---|
352 | { |
---|
353 | /* test if args[1] is in result */ |
---|
354 | attrib(result, "isSB", 1); |
---|
355 | int i; |
---|
356 | for (i = ncols(args[1]); i > 0; i--) |
---|
357 | { |
---|
358 | if (reduce(args[1][i], result, 1) != 0) { return(0); } |
---|
359 | } |
---|
360 | |
---|
361 | /* test if result is a GB */ |
---|
362 | //def G = std(result); |
---|
363 | //if (reduce_parallel(G, result)) { return(0); } |
---|
364 | //return(1); |
---|
365 | return(system("verifyGB",result)); |
---|
366 | } |
---|
367 | |
---|
368 | /* return 1, if I_reduce is _not_ in G_reduce, |
---|
369 | * 0, otherwise |
---|
370 | * (same as size(reduce(I_reduce, G_reduce))). |
---|
371 | * Uses parallelization. */ |
---|
372 | static proc reduce_parallel(def I_reduce, def G_reduce) |
---|
373 | { |
---|
374 | exportto(Modstd, I_reduce); |
---|
375 | exportto(Modstd, G_reduce); |
---|
376 | int size_I = ncols(I_reduce); |
---|
377 | int chunks = Modular::par_range(size_I); |
---|
378 | intvec range; |
---|
379 | int i; |
---|
380 | for (i = chunks; i > 0; i--) |
---|
381 | { |
---|
382 | range = Modular::par_range(size_I, i); |
---|
383 | task t(i) = "Modstd::reduce_task", list(range); |
---|
384 | } |
---|
385 | startTasks(t(1..chunks)); |
---|
386 | waitAllTasks(t(1..chunks)); |
---|
387 | int result = 0; |
---|
388 | for (i = chunks; i > 0; i--) |
---|
389 | { |
---|
390 | if (getResult(t(i))) { result = 1; break; } |
---|
391 | } |
---|
392 | kill I_reduce; |
---|
393 | kill G_reduce; |
---|
394 | return(result); |
---|
395 | } |
---|
396 | |
---|
397 | /* compute a chunk of reductions for reduce_parallel */ |
---|
398 | static proc reduce_task(intvec range) |
---|
399 | { |
---|
400 | int result = 0; |
---|
401 | int i; |
---|
402 | for (i = range[1]; i <= range[2]; i++) |
---|
403 | { |
---|
404 | if (reduce(I_reduce[i], G_reduce, 1) != 0) { result = 1; break; } |
---|
405 | } |
---|
406 | return(result); |
---|
407 | } |
---|
408 | |
---|
409 | //////////////////////////////////////////////////////////////////////////////// |
---|
410 | /* |
---|
411 | * The following procedures are kept for backward compatibility with the old |
---|
412 | * version of modstd.lib. As of now (May 2014), they are still needed in |
---|
413 | * modnormal.lib, modwalk.lib, and symodstd.lib. They can be removed here as |
---|
414 | * soon as they are not longer needed in these libraries. |
---|
415 | */ |
---|
416 | |
---|
417 | LIB "parallel.lib"; |
---|
418 | |
---|
419 | static proc mod_init() |
---|
420 | { |
---|
421 | newstruct("idealPrimeTest", "ideal Ideal"); |
---|
422 | } |
---|
423 | |
---|
424 | static proc redFork(ideal I, ideal J, int n) |
---|
425 | { |
---|
426 | attrib(J,"isSB",1); |
---|
427 | return(reduce(I,J,1)); |
---|
428 | } |
---|
429 | |
---|
430 | proc isIncluded(ideal I, ideal J, list #) |
---|
431 | "USAGE: isIncluded(I,J); I,J ideals |
---|
432 | RETURN: 1 if J includes I, |
---|
433 | @* 0 if there is an element f in I which does not reduce to 0 w.r.t. J. |
---|
434 | EXAMPLE: example isIncluded; shows an example |
---|
435 | " |
---|
436 | { |
---|
437 | def R = basering; |
---|
438 | setring R; |
---|
439 | |
---|
440 | attrib(J,"isSB",1); |
---|
441 | int i,j,k; |
---|
442 | |
---|
443 | if(size(#) > 0) |
---|
444 | { |
---|
445 | int n = #[1]; |
---|
446 | if(n >= ncols(I)) { n = ncols(I); } |
---|
447 | if(n > 1) |
---|
448 | { |
---|
449 | for(i = 1; i <= n - 1; i++) |
---|
450 | { |
---|
451 | //link l(i) = "MPtcp:fork"; |
---|
452 | link l(i) = "ssi:fork"; |
---|
453 | open(l(i)); |
---|
454 | |
---|
455 | write(l(i), quote(redFork(eval(I[ncols(I)-i]), eval(J), 1))); |
---|
456 | } |
---|
457 | |
---|
458 | int t = timer; |
---|
459 | if(reduce(I[ncols(I)], J, 1) != 0) |
---|
460 | { |
---|
461 | for(i = 1; i <= n - 1; i++) |
---|
462 | { |
---|
463 | close(l(i)); |
---|
464 | } |
---|
465 | return(0); |
---|
466 | } |
---|
467 | t = timer - t; |
---|
468 | if(t > 60) { t = 60; } |
---|
469 | int i_sleep = system("sh", "sleep "+string(t)); |
---|
470 | |
---|
471 | j = ncols(I) - n; |
---|
472 | |
---|
473 | while(j >= 0) |
---|
474 | { |
---|
475 | for(i = 1; i <= n - 1; i++) |
---|
476 | { |
---|
477 | if(status(l(i), "read", "ready")) |
---|
478 | { |
---|
479 | if(read(l(i)) != 0) |
---|
480 | { |
---|
481 | for(i = 1; i <= n - 1; i++) |
---|
482 | { |
---|
483 | close(l(i)); |
---|
484 | } |
---|
485 | return(0); |
---|
486 | } |
---|
487 | else |
---|
488 | { |
---|
489 | if(j >= 1) |
---|
490 | { |
---|
491 | write(l(i), quote(redFork(eval(I[j]), eval(J), 1))); |
---|
492 | j--; |
---|
493 | } |
---|
494 | else |
---|
495 | { |
---|
496 | k++; |
---|
497 | close(l(i)); |
---|
498 | } |
---|
499 | } |
---|
500 | } |
---|
501 | } |
---|
502 | if(k == n - 1) |
---|
503 | { |
---|
504 | j--; |
---|
505 | } |
---|
506 | i_sleep = system("sh", "sleep "+string(t)); |
---|
507 | } |
---|
508 | return(1); |
---|
509 | } |
---|
510 | } |
---|
511 | |
---|
512 | for(i = ncols(I); i >= 1; i--) |
---|
513 | { |
---|
514 | if(reduce(I[i],J,1) != 0){ return(0); } |
---|
515 | } |
---|
516 | return(1); |
---|
517 | } |
---|
518 | example |
---|
519 | { "EXAMPLE:"; echo = 2; |
---|
520 | ring r=0,(x,y,z),dp; |
---|
521 | ideal I = x+1,x+y+1; |
---|
522 | ideal J = x+1,y; |
---|
523 | isIncluded(I,J); |
---|
524 | isIncluded(J,I); |
---|
525 | isIncluded(I,J,4); |
---|
526 | |
---|
527 | ring R = 0, x(1..5), dp; |
---|
528 | ideal I1 = cyclic(4); |
---|
529 | ideal I2 = I1,x(5)^2; |
---|
530 | isIncluded(I1,I2,4); |
---|
531 | } |
---|
532 | |
---|
533 | proc deleteUnluckyPrimes(list T, list L, int ho, list #) |
---|
534 | "USAGE: deleteUnluckyPrimes(T,L,ho,#); T/L list of polys/primes, ho integer |
---|
535 | RETURN: lists T,L(,M),lT with T/L(/M) list of polys/primes(/type of #), |
---|
536 | lT ideal |
---|
537 | NOTE: - if ho = 1, the polynomials in T are homogeneous, else ho = 0, |
---|
538 | @* - lT is prevalent, i.e. the most appearing leading ideal in T |
---|
539 | EXAMPLE: example deleteUnluckyPrimes; shows an example |
---|
540 | " |
---|
541 | { |
---|
542 | ho = ((ho)||(ord_test(basering) == -1)); |
---|
543 | int j,k,c; |
---|
544 | intvec hl,hc; |
---|
545 | ideal cT,lT,cK; |
---|
546 | lT = lead(T[size(T)]); |
---|
547 | attrib(lT,"isSB",1); |
---|
548 | if(!ho) |
---|
549 | { |
---|
550 | for(j = 1; j < size(T); j++) |
---|
551 | { |
---|
552 | cT = lead(T[j]); |
---|
553 | attrib(cT,"isSB",1); |
---|
554 | if((size(reduce(cT,lT,5))!=0)||(size(reduce(lT,cT,5))!=0)) |
---|
555 | { |
---|
556 | cK = cT; |
---|
557 | c++; |
---|
558 | } |
---|
559 | } |
---|
560 | if(c > size(T) div 2){ lT = cK; } |
---|
561 | } |
---|
562 | else |
---|
563 | { |
---|
564 | hl = hilb(lT,1); |
---|
565 | for(j = 1; j < size(T); j++) |
---|
566 | { |
---|
567 | cT = lead(T[j]); |
---|
568 | attrib(cT,"isSB",1); |
---|
569 | hc = hilb(cT,1); |
---|
570 | if(hl == hc) |
---|
571 | { |
---|
572 | for(k = 1; k <= size(lT); k++) |
---|
573 | { |
---|
574 | if(lT[k] < cT[k]) { lT = cT; c++; break; } |
---|
575 | if(lT[k] > cT[k]) { c++; break; } |
---|
576 | } |
---|
577 | } |
---|
578 | else |
---|
579 | { |
---|
580 | if(hc < hl){ lT = cT; hl = hilb(lT,1); c++; } |
---|
581 | } |
---|
582 | } |
---|
583 | } |
---|
584 | |
---|
585 | int addList; |
---|
586 | if(size(#) > 0) { list M = #; addList = 1; } |
---|
587 | j = 1; |
---|
588 | attrib(lT,"isSB",1); |
---|
589 | while((j <= size(T))&&(c > 0)) |
---|
590 | { |
---|
591 | cT = lead(T[j]); |
---|
592 | attrib(cT,"isSB",1); |
---|
593 | if((size(reduce(cT,lT,5)) != 0)||(size(reduce(lT,cT,5)) != 0)) |
---|
594 | { |
---|
595 | T = delete(T,j); |
---|
596 | if(j == 1) |
---|
597 | { |
---|
598 | L = L[2..size(L)]; |
---|
599 | if(addList == 1) { M = M[2..size(M)]; } |
---|
600 | } |
---|
601 | else |
---|
602 | { |
---|
603 | if(j == size(L)) |
---|
604 | { |
---|
605 | L = L[1..size(L)-1]; |
---|
606 | if(addList == 1) { M = M[1..size(M)-1]; } |
---|
607 | } |
---|
608 | else |
---|
609 | { |
---|
610 | L = L[1..j-1],L[j+1..size(L)]; |
---|
611 | if(addList == 1) { M = M[1..j-1],M[j+1..size(M)]; } |
---|
612 | } |
---|
613 | } |
---|
614 | j--; |
---|
615 | } |
---|
616 | j++; |
---|
617 | } |
---|
618 | |
---|
619 | for(j = 1; j <= size(L); j++) |
---|
620 | { |
---|
621 | L[j] = bigint(L[j]); |
---|
622 | } |
---|
623 | |
---|
624 | if(addList == 0) { return(list(T,L,lT)); } |
---|
625 | if(addList == 1) { return(list(T,L,M,lT)); } |
---|
626 | } |
---|
627 | example |
---|
628 | { "EXAMPLE:"; echo = 2; |
---|
629 | list L = 2,3,5,7,11; |
---|
630 | ring r = 0,(y,x),Dp; |
---|
631 | ideal I1 = 2y2x,y6; |
---|
632 | ideal I2 = yx2,y3x,x5,y6; |
---|
633 | ideal I3 = y2x,x3y,x5,y6; |
---|
634 | ideal I4 = y2x,11x3y,x5; |
---|
635 | ideal I5 = y2x,yx3,x5,7y6; |
---|
636 | list T = I1,I2,I3,I4,I5; |
---|
637 | deleteUnluckyPrimes(T,L,1); |
---|
638 | list P = poly(x),poly(x2),poly(x3),poly(x4),poly(x5); |
---|
639 | deleteUnluckyPrimes(T,L,1,P); |
---|
640 | } |
---|
641 | |
---|
642 | proc primeTest(def II, bigint p) |
---|
643 | { |
---|
644 | if(typeof(II) == "string") |
---|
645 | { |
---|
646 | ideal I = `II`.Ideal; |
---|
647 | } |
---|
648 | else |
---|
649 | { |
---|
650 | ideal I = II; |
---|
651 | } |
---|
652 | |
---|
653 | I = simplify(I, 2); // erase zero generators |
---|
654 | |
---|
655 | int i,j; |
---|
656 | poly f; |
---|
657 | number cnt; |
---|
658 | for(i = 1; i <= size(I); i++) |
---|
659 | { |
---|
660 | f = cleardenom(I[i]); |
---|
661 | if(f == 0) { return(0); } |
---|
662 | cnt = leadcoef(I[i])/leadcoef(f); |
---|
663 | if((bigint(numerator(cnt)) mod p) == 0) { return(0); } |
---|
664 | if((bigint(denominator(cnt)) mod p) == 0) { return(0); } |
---|
665 | for(j = size(f); j > 0; j--) |
---|
666 | { |
---|
667 | if((bigint(leadcoef(f[j])) mod p) == 0) { return(0); } |
---|
668 | } |
---|
669 | } |
---|
670 | return(1); |
---|
671 | } |
---|
672 | |
---|
673 | proc primeList(ideal I, int n, list #) |
---|
674 | "USAGE: primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal, |
---|
675 | n integer |
---|
676 | RETURN: the intvec of n greatest primes <= 2147483647 (resp. n greatest primes |
---|
677 | < L[size(L)] union with L) such that none of these primes divides any |
---|
678 | coefficient occurring in I |
---|
679 | NOTE: The number of cores to use can be defined by ncores, default is 1. |
---|
680 | EXAMPLE: example primeList; shows an example |
---|
681 | " |
---|
682 | { |
---|
683 | intvec L; |
---|
684 | int i,p; |
---|
685 | int ncores = 1; |
---|
686 | |
---|
687 | //----------------- Initialize optional parameter ncores --------------------- |
---|
688 | if(size(#) > 0) |
---|
689 | { |
---|
690 | if(size(#) == 1) |
---|
691 | { |
---|
692 | if(typeof(#[1]) == "int") |
---|
693 | { |
---|
694 | ncores = #[1]; |
---|
695 | # = list(); |
---|
696 | } |
---|
697 | } |
---|
698 | else |
---|
699 | { |
---|
700 | ncores = #[2]; |
---|
701 | } |
---|
702 | } |
---|
703 | |
---|
704 | if(size(#) == 0) |
---|
705 | { |
---|
706 | p = 2147483647; |
---|
707 | while(!primeTest(I,p)) |
---|
708 | { |
---|
709 | p = prime(p-1); |
---|
710 | if(p == 2) { ERROR("no more primes"); } |
---|
711 | } |
---|
712 | L[1] = p; |
---|
713 | } |
---|
714 | else |
---|
715 | { |
---|
716 | L = #[1]; |
---|
717 | p = prime(L[size(L)]-1); |
---|
718 | while(!primeTest(I,p)) |
---|
719 | { |
---|
720 | p = prime(p-1); |
---|
721 | if(p == 2) { ERROR("no more primes"); } |
---|
722 | } |
---|
723 | L[size(L)+1] = p; |
---|
724 | } |
---|
725 | if(p == 2) { ERROR("no more primes"); } |
---|
726 | if(ncores == 1) |
---|
727 | { |
---|
728 | for(i = 2; i <= n; i++) |
---|
729 | { |
---|
730 | p = prime(p-1); |
---|
731 | while(!primeTest(I,p)) |
---|
732 | { |
---|
733 | p = prime(p-1); |
---|
734 | if(p == 2) { ERROR("no more primes"); } |
---|
735 | } |
---|
736 | L[size(L)+1] = p; |
---|
737 | } |
---|
738 | } |
---|
739 | else |
---|
740 | { |
---|
741 | int neededSize = size(L)+n-1;; |
---|
742 | list parallelResults; |
---|
743 | list arguments; |
---|
744 | int neededPrimes = neededSize-size(L); |
---|
745 | idealPrimeTest Id; |
---|
746 | Id.Ideal = I; |
---|
747 | export(Id); |
---|
748 | while(neededPrimes > 0) |
---|
749 | { |
---|
750 | arguments = list(); |
---|
751 | for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0)) |
---|
752 | *ncores; i > 0; i--) |
---|
753 | { |
---|
754 | p = prime(p-1); |
---|
755 | if(p == 2) { ERROR("no more primes"); } |
---|
756 | arguments[i] = list("Id", p); |
---|
757 | } |
---|
758 | parallelResults = parallelWaitAll("primeTest", arguments, 0, ncores); |
---|
759 | for(i = size(arguments); i > 0; i--) |
---|
760 | { |
---|
761 | if(parallelResults[i]) |
---|
762 | { |
---|
763 | L[size(L)+1] = arguments[i][2]; |
---|
764 | } |
---|
765 | } |
---|
766 | neededPrimes = neededSize-size(L); |
---|
767 | } |
---|
768 | kill Id; |
---|
769 | if(size(L) > neededSize) |
---|
770 | { |
---|
771 | L = L[1..neededSize]; |
---|
772 | } |
---|
773 | } |
---|
774 | return(L); |
---|
775 | } |
---|
776 | example |
---|
777 | { "EXAMPLE:"; echo = 2; |
---|
778 | ring r = 0,(x,y,z),dp; |
---|
779 | ideal I = 2147483647x+y, z-181; |
---|
780 | intvec L = primeList(I,10); |
---|
781 | size(L); |
---|
782 | L[1]; |
---|
783 | L[size(L)]; |
---|
784 | L = primeList(I,5,L); |
---|
785 | size(L); |
---|
786 | L[size(L)]; |
---|
787 | } |
---|
788 | |
---|
789 | proc modStdL(def I, list #) |
---|
790 | "USAGE: modStdL(I[, exactness]); I ideal/module, exactness int |
---|
791 | RETURN: a standard basis of I |
---|
792 | NOTE: The procedure computes a standard basis of I (over the rational |
---|
793 | numbers) by using modular methods via an external Singular. |
---|
794 | @* An optional parameter 'exactness' can be provided. |
---|
795 | If exactness = 1(default), the procedure computes a standard basis |
---|
796 | of I for sure; if exactness = 0, it computes a standard basis of I |
---|
797 | with high probability. |
---|
798 | SEE ALSO: modular, modStd |
---|
799 | EXAMPLE: example modStdL; shows an example" |
---|
800 | { |
---|
801 | link l="ssi:tcp localhost:"+system("Singular"); |
---|
802 | write(l,quote(option(noloadLib))); // suppress "loaded..." |
---|
803 | read(l); //dummy: return value of option |
---|
804 | write(l,quote(load("modstd.lib","with"))); // load library |
---|
805 | read(l); //dummy: return value of load |
---|
806 | if (size(#)==0) |
---|
807 | { |
---|
808 | write(l,quote(modStd(eval(I)))); |
---|
809 | } |
---|
810 | else |
---|
811 | { |
---|
812 | write(l,quote(modStd(eval(I),eval(#[1])))); |
---|
813 | } |
---|
814 | return(read(l)); |
---|
815 | } |
---|
816 | example |
---|
817 | { |
---|
818 | "EXAMPLE:"; |
---|
819 | echo = 2; |
---|
820 | ring R1 = 0, (x,y,z,t), dp; |
---|
821 | ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; |
---|
822 | ideal J = modStdL(I); |
---|
823 | J; |
---|
824 | I = homog(I, t); |
---|
825 | J = modStdL(I); |
---|
826 | J; |
---|
827 | |
---|
828 | ring R2 = 0, (x,y,z), ds; |
---|
829 | ideal I = jacob(x5+y6+z7+xyz); |
---|
830 | ideal J = modStdL(I, 0); |
---|
831 | J; |
---|
832 | |
---|
833 | ring R3 = 0, x(1..4), lp; |
---|
834 | ideal I = cyclic(4); |
---|
835 | ideal J1 = modStdL(I, 1); // default |
---|
836 | ideal J2 = modStdL(I, 0); |
---|
837 | size(reduce(J1, J2)); |
---|
838 | size(reduce(J2, J1)); |
---|
839 | } |
---|
840 | |
---|
841 | ////////////////////////////// further examples //////////////////////////////// |
---|
842 | |
---|
843 | /* |
---|
844 | ring r = 0, (x,y,z), lp; |
---|
845 | poly s1 = 5x3y2z+3y3x2z+7xy2z2; |
---|
846 | poly s2 = 3xy2z2+x5+11y2z2; |
---|
847 | poly s3 = 4xyz+7x3+12y3+1; |
---|
848 | poly s4 = 3x3-4y3+yz2; |
---|
849 | ideal i = s1, s2, s3, s4; |
---|
850 | |
---|
851 | ring r = 0, (x,y,z), lp; |
---|
852 | poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7; |
---|
853 | poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y; |
---|
854 | poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z; |
---|
855 | poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y; |
---|
856 | ideal i = s1, s2, s3, s4; |
---|
857 | |
---|
858 | ring r = 0, (x,y,z), lp; |
---|
859 | poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz; |
---|
860 | poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4; |
---|
861 | poly s3 = 8x3 + 12y3 + xz2 + 3; |
---|
862 | poly s4 = 7x2y4 + 18xy3z2 + y3z3; |
---|
863 | ideal i = s1, s2, s3, s4; |
---|
864 | |
---|
865 | int n = 6; |
---|
866 | ring r = 0,(x(1..n)),lp; |
---|
867 | ideal i = cyclic(n); |
---|
868 | ring s = 0, (x(1..n),t), lp; |
---|
869 | ideal i = imap(r,i); |
---|
870 | i = homog(i,t); |
---|
871 | |
---|
872 | ring r = 0, (x(1..4),s), (dp(4),dp); |
---|
873 | poly s1 = 1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4); |
---|
874 | poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4); |
---|
875 | poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4); |
---|
876 | poly s4 = x(3) + s^4*x(1)*x(2) + s^19*x(1)*x(3)*x(4) +s^24*x(2)*x(3)*x(4); |
---|
877 | poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4); |
---|
878 | ideal i = s1, s2, s3, s4, s5; |
---|
879 | |
---|
880 | ring r = 0, (x,y,z), ds; |
---|
881 | int a = 16; |
---|
882 | int b = 15; |
---|
883 | int c = 4; |
---|
884 | int t = 1; |
---|
885 | poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3 |
---|
886 | +x^(c-2)*y^c*(y2+t*x)^2; |
---|
887 | ideal i = jacob(f); |
---|
888 | |
---|
889 | ring r = 0, (x,y,z), ds; |
---|
890 | int a = 25; |
---|
891 | int b = 25; |
---|
892 | int c = 5; |
---|
893 | int t = 1; |
---|
894 | poly f = x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3 |
---|
895 | +x^(c-2)*y^c*(y2+t*x)^2; |
---|
896 | ideal i = jacob(f),f; |
---|
897 | |
---|
898 | ring r = 0, (x,y,z), ds; |
---|
899 | int a = 10; |
---|
900 | poly f = xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a; |
---|
901 | ideal i = jacob(f); |
---|
902 | |
---|
903 | ring r = 0, (x,y,z), ds; |
---|
904 | int a = 6; |
---|
905 | int b = 8; |
---|
906 | int c = 10; |
---|
907 | int alpha = 5; |
---|
908 | int beta = 5; |
---|
909 | int t = 1; |
---|
910 | poly f = x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3) |
---|
911 | +x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2; |
---|
912 | ideal i = jacob(f); |
---|
913 | */ |
---|
914 | |
---|