[380a17b] | 1 | /////////////////////////////////////////////////////////////////////////////// |
---|
[8e9e13] | 2 | version="version modstd.lib 4.2.0.0 Feb_2020 "; // $Id$ |
---|
[a2e51b] | 3 | category="Commutative Algebra"; |
---|
[d68d30] | 4 | info=" |
---|
[8ac0543] | 5 | LIBRARY: modstd.lib Groebner bases of ideals/modules using modular methods |
---|
[d68d30] | 6 | |
---|
| 7 | AUTHORS: A. Hashemi Amir.Hashemi@lip6.fr |
---|
[a2e51b] | 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 |
---|
[d68d30] | 12 | |
---|
| 13 | OVERVIEW: |
---|
[8ac0543] | 14 | A library for computing Groebner bases of ideals/modules in the polynomial ring |
---|
| 15 | over the rational numbers using modular methods. |
---|
[d68d30] | 16 | |
---|
[a2e51b] | 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). |
---|
[d68d30] | 23 | |
---|
| 24 | PROCEDURES: |
---|
[a2e51b] | 25 | modStd(I); standard basis of I using modular methods |
---|
[d596350] | 26 | modSyz(I); syzygy module of I using modular methods |
---|
[2206a04] | 27 | modIntersect(I,J); intersection of I and J using modular methods |
---|
[d68d30] | 28 | "; |
---|
| 29 | |
---|
[4689d0] | 30 | LIB "polylib.lib"; |
---|
[a2e51b] | 31 | LIB "modular.lib"; |
---|
[a7050d] | 32 | |
---|
[61ed0b] | 33 | proc modStd(def I, list #) |
---|
| 34 | "USAGE: modStd(I[, exactness]); I ideal/module, exactness int |
---|
[a2e51b] | 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. |
---|
[b87290] | 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 |
---|
[a2e51b] | 41 | with high probability. |
---|
| 42 | SEE ALSO: modular |
---|
| 43 | EXAMPLE: example modStd; shows an example" |
---|
[d68d30] | 44 | { |
---|
[a2e51b] | 45 | /* read optional parameter */ |
---|
| 46 | int exactness = 1; |
---|
[61ed0b] | 47 | if (size(#) > 0) |
---|
| 48 | { |
---|
[a2e51b] | 49 | /* For compatibility, we only test size(#) > 4. This can be changed to |
---|
| 50 | * size(#) > 1 in the future. */ |
---|
[61ed0b] | 51 | if (size(#) > 4 || typeof(#[1]) != "int") |
---|
| 52 | { |
---|
[a2e51b] | 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"; |
---|
[61ed0b] | 64 | if (npars(basering) > 0) { command = "Modstd::groebner_norm"; } |
---|
[a2e51b] | 65 | |
---|
| 66 | /* call modular() */ |
---|
[61ed0b] | 67 | if (exactness) |
---|
| 68 | { |
---|
[510427] | 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 | { |
---|
[a2e51b] | 76 | I = modular(command, list(I), primeTest_std, |
---|
| 77 | deleteUnluckyPrimes_std, pTest_std, finalTest_std); |
---|
[510427] | 78 | } |
---|
[a2e51b] | 79 | } |
---|
[61ed0b] | 80 | else |
---|
| 81 | { |
---|
[a2e51b] | 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); |
---|
[d68d30] | 90 | } |
---|
| 91 | example |
---|
| 92 | { |
---|
[a2e51b] | 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)); |
---|
[d68d30] | 114 | } |
---|
| 115 | |
---|
[d596350] | 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. |
---|
[5f5a82] | 121 | The property of being a syzygy is tested. |
---|
[d596350] | 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 | |
---|
[2206a04] | 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 |
---|
[5f5a82] | 154 | with high probability. |
---|
| 155 | No additional tests are performed. |
---|
[2206a04] | 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 | |
---|
[a2e51b] | 184 | /* compute a normalized GB via groebner() */ |
---|
[8ac0543] | 185 | static proc groebner_norm(def I) |
---|
[d68d30] | 186 | { |
---|
[a2e51b] | 187 | I = simplify(groebner(I), 1); |
---|
| 188 | attrib(I, "isSB", 1); |
---|
| 189 | return(I); |
---|
[d68d30] | 190 | } |
---|
| 191 | |
---|
[a2e51b] | 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) |
---|
[d68d30] | 195 | { |
---|
[a2e51b] | 196 | /* erase zero generators */ |
---|
[61ed0b] | 197 | def I = simplify(args[1], 2); |
---|
[a2e51b] | 198 | |
---|
| 199 | /* clear denominators and count the terms */ |
---|
[61ed0b] | 200 | def J=I; // dummy assign, to get the type of I |
---|
[a2e51b] | 201 | ideal K; |
---|
| 202 | int n = ncols(I); |
---|
| 203 | intvec sizes; |
---|
| 204 | number cnt; |
---|
| 205 | int i; |
---|
[61ed0b] | 206 | for(i = n; i > 0; i--) |
---|
| 207 | { |
---|
[a2e51b] | 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); |
---|
[61ed0b] | 217 | if (typeof(lbr[1]) == "int") { lbr[1] = p; } |
---|
| 218 | else { lbr[1][1] = p; } |
---|
[a2e51b] | 219 | def rp = ring(lbr); |
---|
| 220 | setring(rp); |
---|
[61ed0b] | 221 | def Jp = fetch(br, J); |
---|
[a2e51b] | 222 | ideal Kp = fetch(br, K); |
---|
| 223 | |
---|
| 224 | /* test if any coefficient is missing */ |
---|
[61ed0b] | 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); } |
---|
[a2e51b] | 227 | setring(br); |
---|
| 228 | return(1); |
---|
[d68d30] | 229 | } |
---|
| 230 | |
---|
[a2e51b] | 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) |
---|
[d68d30] | 235 | { |
---|
[a2e51b] | 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; |
---|
[d596350] | 246 | def L=modresults[1]; // dummy assign to get the type of L |
---|
[a2e51b] | 247 | int i; |
---|
| 248 | int j; |
---|
[61ed0b] | 249 | for (i = 1; i <= size_modresults; i++) |
---|
| 250 | { |
---|
[a2e51b] | 251 | L = lead(modresults[i]); |
---|
| 252 | attrib(L, "isSB", 1); |
---|
[61ed0b] | 253 | for (j = 1; j <= size_cat; j++) |
---|
| 254 | { |
---|
[a2e51b] | 255 | if (size(L) == size(cat[j][1]) |
---|
[6518eba] | 256 | && size(reduce(L, cat[j][1], 5)) == 0 |
---|
| 257 | && size(reduce(cat[j][1], L, 5)) == 0) |
---|
[61ed0b] | 258 | { |
---|
[a2e51b] | 259 | cat[j][2] = cat[j][2]+1; |
---|
| 260 | cat[j][3][cat[j][2]] = i; |
---|
| 261 | break; |
---|
[7f30e2] | 262 | } |
---|
[a2e51b] | 263 | } |
---|
[61ed0b] | 264 | if (j > size_cat) |
---|
| 265 | { |
---|
[a2e51b] | 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]; |
---|
[61ed0b] | 277 | for (i = 2; i <= size_cat; i++) |
---|
| 278 | { |
---|
| 279 | if (cat[i][2] > max) |
---|
| 280 | { |
---|
[a2e51b] | 281 | cat_max = i; |
---|
| 282 | max = cat[i][2]; |
---|
| 283 | } |
---|
| 284 | } |
---|
| 285 | |
---|
| 286 | /* return all other indices */ |
---|
| 287 | list unluckyIndices; |
---|
[61ed0b] | 288 | for (i = 1; i <= size_cat; i++) |
---|
| 289 | { |
---|
| 290 | if (i != cat_max) { unluckyIndices = unluckyIndices + cat[i][3]; } |
---|
[a2e51b] | 291 | } |
---|
| 292 | return(unluckyIndices); |
---|
[d68d30] | 293 | } |
---|
[5f5a82] | 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 | } |
---|
[d68d30] | 312 | |
---|
[d596350] | 313 | static proc pTest_syz(string command, alias list args, alias def result, int p) |
---|
| 314 | { |
---|
[5f5a82] | 315 | module result_without_denom=cleardenomModule(result); |
---|
| 316 | return(size(module(matrix(args[1])*matrix(result_without_denom)))==0); |
---|
[d596350] | 317 | } |
---|
| 318 | |
---|
[a2e51b] | 319 | /* test if 'command' applied to 'args' in characteristic p is the same as |
---|
| 320 | 'result' mapped to characteristic p */ |
---|
[61ed0b] | 321 | static proc pTest_std(string command, alias list args, alias def result, |
---|
[9761b8] | 322 | int p) |
---|
[d68d30] | 323 | { |
---|
[a2e51b] | 324 | /* change to characteristic p */ |
---|
| 325 | def br = basering; |
---|
| 326 | list lbr = ringlist(br); |
---|
[61ed0b] | 327 | if (typeof(lbr[1]) == "int") { lbr[1] = p; } |
---|
| 328 | else { lbr[1][1] = p; } |
---|
[a2e51b] | 329 | def rp = ring(lbr); |
---|
| 330 | setring(rp); |
---|
[61ed0b] | 331 | def Ip = fetch(br, args)[1]; |
---|
| 332 | def Gp = fetch(br, result); |
---|
[a2e51b] | 333 | attrib(Gp, "isSB", 1); |
---|
| 334 | |
---|
| 335 | /* test if Ip is in Gp */ |
---|
| 336 | int i; |
---|
[61ed0b] | 337 | for (i = ncols(Ip); i > 0; i--) |
---|
| 338 | { |
---|
| 339 | if (reduce(Ip[i], Gp, 1) != 0) |
---|
| 340 | { |
---|
[a2e51b] | 341 | setring(br); |
---|
| 342 | return(0); |
---|
| 343 | } |
---|
| 344 | } |
---|
[d68d30] | 345 | |
---|
[a2e51b] | 346 | /* compute command(args) */ |
---|
| 347 | execute("Ip = "+command+"(Ip);"); |
---|
[d68d30] | 348 | |
---|
[a2e51b] | 349 | /* test if Gp is in Ip */ |
---|
[61ed0b] | 350 | for (i = ncols(Gp); i > 0; i--) |
---|
| 351 | { |
---|
| 352 | if (reduce(Gp[i], Ip, 1) != 0) { setring(br); return(0); } |
---|
[a2e51b] | 353 | } |
---|
| 354 | setring(br); |
---|
| 355 | return(1); |
---|
[d68d30] | 356 | } |
---|
| 357 | |
---|
[510427] | 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) |
---|
[d68d30] | 360 | { |
---|
[a2e51b] | 361 | /* test if args[1] is in result */ |
---|
| 362 | attrib(result, "isSB", 1); |
---|
| 363 | int i; |
---|
[61ed0b] | 364 | for (i = ncols(args[1]); i > 0; i--) |
---|
| 365 | { |
---|
| 366 | if (reduce(args[1][i], result, 1) != 0) { return(0); } |
---|
[a2e51b] | 367 | } |
---|
| 368 | |
---|
| 369 | /* test if result is a GB */ |
---|
[ed9ea8] | 370 | //def G = std(result); |
---|
| 371 | //if (reduce_parallel(G, result)) { return(0); } |
---|
| 372 | //return(1); |
---|
| 373 | return(system("verifyGB",result)); |
---|
[d68d30] | 374 | } |
---|
| 375 | |
---|
[510427] | 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 | |
---|
[a2e51b] | 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. */ |
---|
[4454fb7] | 397 | static proc reduce_parallel(def I_reduce, def G_reduce) |
---|
[d68d30] | 398 | { |
---|
[a2e51b] | 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; |
---|
[61ed0b] | 405 | for (i = chunks; i > 0; i--) |
---|
| 406 | { |
---|
[a2e51b] | 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; |
---|
[61ed0b] | 413 | for (i = chunks; i > 0; i--) |
---|
| 414 | { |
---|
| 415 | if (getResult(t(i))) { result = 1; break; } |
---|
[a2e51b] | 416 | } |
---|
| 417 | kill I_reduce; |
---|
| 418 | kill G_reduce; |
---|
| 419 | return(result); |
---|
[d68d30] | 420 | } |
---|
| 421 | |
---|
[a2e51b] | 422 | /* compute a chunk of reductions for reduce_parallel */ |
---|
| 423 | static proc reduce_task(intvec range) |
---|
[d68d30] | 424 | { |
---|
[a2e51b] | 425 | int result = 0; |
---|
| 426 | int i; |
---|
[61ed0b] | 427 | for (i = range[1]; i <= range[2]; i++) |
---|
| 428 | { |
---|
| 429 | if (reduce(I_reduce[i], G_reduce, 1) != 0) { result = 1; break; } |
---|
[a2e51b] | 430 | } |
---|
| 431 | return(result); |
---|
[d68d30] | 432 | } |
---|
| 433 | |
---|
[90c8bb] | 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 |
---|
[d1cc8c3] | 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. |
---|
[90c8bb] | 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); |
---|
[6518eba] | 579 | if((size(reduce(cT,lT,5))!=0)||(size(reduce(lT,cT,5))!=0)) |
---|
[90c8bb] | 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); |
---|
[6518eba] | 618 | if((size(reduce(cT,lT,5)) != 0)||(size(reduce(lT,cT,5)) != 0)) |
---|
[90c8bb] | 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); |
---|
[754c382] | 688 | if((bigint(numerator(cnt)) mod p) == 0) { return(0); } |
---|
| 689 | if((bigint(denominator(cnt)) mod p) == 0) { return(0); } |
---|
[90c8bb] | 690 | for(j = size(f); j > 0; j--) |
---|
| 691 | { |
---|
[754c382] | 692 | if((bigint(leadcoef(f[j])) mod p) == 0) { return(0); } |
---|
[90c8bb] | 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 |
---|
[3a0213] | 703 | coefficient occurring in I |
---|
[90c8bb] | 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 | } |
---|
[d68d30] | 813 | |
---|
[8e9e13] | 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 | |
---|
[d68d30] | 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; |
---|
[4cc3cac] | 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; |
---|
[d68d30] | 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; |
---|
[4cc3cac] | 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; |
---|
[d68d30] | 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; |
---|
[4cc3cac] | 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; |
---|
[d68d30] | 937 | ideal i = jacob(f); |
---|
| 938 | */ |
---|
| 939 | |
---|