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