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