[380a17b] | 1 | /////////////////////////////////////////////////////////////////////////////// |
---|
[6391eb] | 2 | version="version modwalk.lib 4.1.2.0 Feb_2019 "; // $Id$ |
---|
[e92c6a] | 3 | category = "Commutative Algebra"; |
---|
| 4 | info=" |
---|
| 5 | LIBRARY: modwalk.lib Groebner basis convertion |
---|
| 6 | |
---|
| 7 | AUTHORS: S. Oberfranz oberfran@mathematik.uni-kl.de |
---|
| 8 | |
---|
| 9 | OVERVIEW: |
---|
| 10 | |
---|
| 11 | A library for converting Groebner bases of an ideal in the polynomial |
---|
| 12 | ring over the rational numbers using modular methods. The procedures are |
---|
| 13 | inspired by the following paper: |
---|
| 14 | Elizabeth A. Arnold: Modular algorithms for computing Groebner bases. |
---|
| 15 | Journal of Symbolic Computation 35, 403-419 (2003). |
---|
| 16 | |
---|
| 17 | PROCEDURES: |
---|
[daf501] | 18 | |
---|
| 19 | modWalk(I,#); standard basis conversion of I by Groebner Walk using modular methods |
---|
[b7405b] | 20 | modrWalk(I,radius,#); standard basis conversion of I by Random Walk using modular methods |
---|
[daf501] | 21 | modfWalk(I,#); standard basis conversion of I by Fractal Walk using modular methods |
---|
| 22 | modfrWalk(I,radius,#); standard basis conversion of I by Random Fractal Walk using modular methods |
---|
[4698cf8] | 23 | |
---|
| 24 | KEYWORDS: walk, groebner;Groebnerwalk |
---|
| 25 | SEE ALSO: grwalk_lib;swalk_lib;rwalk_lib |
---|
[e92c6a] | 26 | "; |
---|
| 27 | |
---|
| 28 | LIB "rwalk.lib"; |
---|
| 29 | LIB "grwalk.lib"; |
---|
[daf501] | 30 | LIB "modular.lib"; |
---|
[468f1c] | 31 | LIB "ring.lib"; |
---|
[daf501] | 32 | |
---|
| 33 | proc modWalk(ideal I, list #) |
---|
[b7405b] | 34 | "USAGE: modWalk(I, [, v, w]); I ideal, v intvec or string, w intvec |
---|
| 35 | If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp). |
---|
| 36 | If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with |
---|
| 37 | respect to dp or Dp, respectively. |
---|
| 38 | If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise, |
---|
| 39 | the output will be a standard basis with respect to lp. |
---|
| 40 | If no optional argument is given, I is assumed to be a standard basis with respect to dp |
---|
| 41 | and a standard basis with respect to lp will be computed. |
---|
[daf501] | 42 | RETURN: a standard basis of I |
---|
| 43 | NOTE: The procedure computes a standard basis of I (over the rational |
---|
| 44 | numbers) by using modular methods. |
---|
| 45 | SEE ALSO: modular |
---|
| 46 | EXAMPLE: example modWalk; shows an example" |
---|
[e92c6a] | 47 | { |
---|
[daf501] | 48 | /* save options */ |
---|
| 49 | intvec opt = option(get); |
---|
| 50 | option(redSB); |
---|
[75f460] | 51 | |
---|
[daf501] | 52 | /* call modular() */ |
---|
| 53 | if (size(#) > 0) { |
---|
[c448cd] | 54 | I = modular("gwalk", list(I,#), primeTest_std, |
---|
| 55 | deleteUnluckyPrimes_std, pTest_std, finalTest_std); |
---|
[daf501] | 56 | } |
---|
| 57 | else { |
---|
[c448cd] | 58 | I = modular("gwalk", list(I), primeTest_std, |
---|
| 59 | deleteUnluckyPrimes_std, pTest_std, finalTest_std); |
---|
[ffd4a9] | 60 | } |
---|
[9cf0188] | 61 | |
---|
[daf501] | 62 | /* return the result */ |
---|
| 63 | attrib(I, "isSB", 1); |
---|
| 64 | option(set, opt); |
---|
| 65 | return(I); |
---|
[e92c6a] | 66 | } |
---|
| 67 | example |
---|
| 68 | { |
---|
[daf501] | 69 | "EXAMPLE:"; |
---|
| 70 | echo = 2; |
---|
| 71 | ring R1 = 0, (x,y,z,t), dp; |
---|
| 72 | ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; |
---|
| 73 | I = std(I); |
---|
| 74 | ring R2 = 0, (x,y,z,t), lp; |
---|
| 75 | ideal I = fetch(R1, I); |
---|
| 76 | ideal J = modWalk(I); |
---|
| 77 | J; |
---|
[3308a0] | 78 | ring S1 = 0, (a,b,c,d), Dp; |
---|
| 79 | ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c; |
---|
| 80 | I = std(I); |
---|
| 81 | ring S2 = 0, (c,d,b,a), lp; |
---|
| 82 | ideal I = fetch(S1,I); |
---|
| 83 | // I is assumed to be a Dp-Groebner basis. |
---|
| 84 | // We compute a lp-Groebner basis. |
---|
| 85 | ideal J = modWalk(I,"Dp"); |
---|
| 86 | J; |
---|
| 87 | intvec w = 3,2,1,2; |
---|
| 88 | ring S3 = 0, (c,d,b,a), (a(w),lp); |
---|
| 89 | ideal I = fetch(S1,I); |
---|
| 90 | // I is assumed to be a Dp-Groebner basis. |
---|
| 91 | // We compute a (a(w),lp)-Groebner basis. |
---|
| 92 | ideal J = modWalk(I,"Dp",w); |
---|
| 93 | J; |
---|
[e92c6a] | 94 | } |
---|
| 95 | |
---|
[b7405b] | 96 | proc modrWalk(ideal I, int radius, list #) |
---|
| 97 | "USAGE: modrWalk(I, radius[, v, w]); |
---|
| 98 | I ideal, radius int, pertdeg int, v intvec or string, w intvec |
---|
| 99 | If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp). |
---|
| 100 | If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with |
---|
| 101 | respect to dp or Dp, respectively. |
---|
| 102 | If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise, |
---|
| 103 | the output will be a standard basis with respect to lp. |
---|
| 104 | If no optional argument is given, I is assumed to be a standard basis with respect to dp |
---|
| 105 | and a standard basis with respect to lp will be computed. |
---|
[daf501] | 106 | RETURN: a standard basis of I |
---|
| 107 | NOTE: The procedure computes a standard basis of I (over the rational |
---|
| 108 | numbers) by using modular methods. |
---|
| 109 | SEE ALSO: modular |
---|
| 110 | EXAMPLE: example modrWalk; shows an example" |
---|
[e92c6a] | 111 | { |
---|
[daf501] | 112 | /* save options */ |
---|
| 113 | intvec opt = option(get); |
---|
| 114 | option(redSB); |
---|
[e92c6a] | 115 | |
---|
[daf501] | 116 | /* call modular() */ |
---|
| 117 | if (size(#) > 0) { |
---|
[b7405b] | 118 | I = modular("rwalk", list(I,radius,1,#), primeTest_std, |
---|
[c448cd] | 119 | deleteUnluckyPrimes_std, pTest_std, finalTest_std); |
---|
[e92c6a] | 120 | } |
---|
[daf501] | 121 | else { |
---|
[b7405b] | 122 | I = modular("rwalk", list(I,radius,1), primeTest_std, deleteUnluckyPrimes_std, |
---|
[3308a0] | 123 | pTest_std,finalTest_std); |
---|
[e92c6a] | 124 | } |
---|
| 125 | |
---|
[daf501] | 126 | /* return the result */ |
---|
| 127 | attrib(I, "isSB", 1); |
---|
| 128 | option(set, opt); |
---|
| 129 | return(I); |
---|
[e92c6a] | 130 | } |
---|
| 131 | example |
---|
| 132 | { |
---|
[daf501] | 133 | "EXAMPLE:"; |
---|
| 134 | echo = 2; |
---|
| 135 | ring R1 = 0, (x,y,z,t), dp; |
---|
| 136 | ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; |
---|
| 137 | I = std(I); |
---|
| 138 | ring R2 = 0, (x,y,z,t), lp; |
---|
| 139 | ideal I = fetch(R1, I); |
---|
| 140 | int radius = 2; |
---|
[b7405b] | 141 | ideal J = modrWalk(I,radius); |
---|
[daf501] | 142 | J; |
---|
[3308a0] | 143 | ring S1 = 0, (a,b,c,d), Dp; |
---|
| 144 | ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c; |
---|
| 145 | I = std(I); |
---|
| 146 | ring S2 = 0, (c,d,b,a), lp; |
---|
| 147 | ideal I = fetch(S1,I); |
---|
| 148 | // I is assumed to be a Dp-Groebner basis. |
---|
| 149 | // We compute a lp-Groebner basis. |
---|
| 150 | ideal J = modrWalk(I,radius,"Dp"); |
---|
| 151 | J; |
---|
| 152 | intvec w = 3,2,1,2; |
---|
| 153 | ring S3 = 0, (c,d,b,a), (a(w),lp); |
---|
| 154 | ideal I = fetch(S1,I); |
---|
| 155 | // I is assumed to be a Dp-Groebner basis. |
---|
| 156 | // We compute a (a(w),lp)-Groebner basis. |
---|
| 157 | ideal J = modrWalk(I,radius,"Dp",w); |
---|
| 158 | J; |
---|
[e92c6a] | 159 | } |
---|
| 160 | |
---|
[daf501] | 161 | proc modfWalk(ideal I, list #) |
---|
[b7405b] | 162 | "USAGE: modfWalk(I, [, v, w]); I ideal, v intvec or string, w intvec |
---|
| 163 | If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp). |
---|
| 164 | If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with |
---|
| 165 | respect to dp or Dp, respectively. |
---|
| 166 | If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise, |
---|
| 167 | the output will be a standard basis with respect to lp. |
---|
| 168 | If no optional argument is given, I is assumed to be a standard basis with respect to dp |
---|
| 169 | and a standard basis with respect to lp will be computed. |
---|
[daf501] | 170 | RETURN: a standard basis of I |
---|
| 171 | NOTE: The procedure computes a standard basis of I (over the rational |
---|
| 172 | numbers) by using modular methods. |
---|
| 173 | SEE ALSO: modular |
---|
| 174 | EXAMPLE: example modfWalk; shows an example" |
---|
[e92c6a] | 175 | { |
---|
[daf501] | 176 | /* save options */ |
---|
| 177 | intvec opt = option(get); |
---|
| 178 | option(redSB); |
---|
| 179 | |
---|
| 180 | /* call modular() */ |
---|
| 181 | if (size(#) > 0) { |
---|
[c448cd] | 182 | I = modular("fwalk", list(I,#), primeTest_std, |
---|
| 183 | deleteUnluckyPrimes_std, pTest_std, finalTest_std); |
---|
[e92c6a] | 184 | } |
---|
[daf501] | 185 | else { |
---|
[c448cd] | 186 | I = modular("fwalk", list(I), primeTest_std, |
---|
| 187 | deleteUnluckyPrimes_std, pTest_std, finalTest_std); |
---|
[e92c6a] | 188 | } |
---|
| 189 | |
---|
[daf501] | 190 | /* return the result */ |
---|
| 191 | attrib(I, "isSB", 1); |
---|
| 192 | option(set, opt); |
---|
| 193 | return(I); |
---|
| 194 | } |
---|
| 195 | example |
---|
[e92c6a] | 196 | { |
---|
[daf501] | 197 | "EXAMPLE:"; |
---|
| 198 | echo = 2; |
---|
| 199 | ring R1 = 0, (x,y,z,t), dp; |
---|
| 200 | ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; |
---|
| 201 | I = std(I); |
---|
| 202 | ring R2 = 0, (x,y,z,t), lp; |
---|
| 203 | ideal I = fetch(R1, I); |
---|
| 204 | ideal J = modfWalk(I); |
---|
| 205 | J; |
---|
[3308a0] | 206 | ring S1 = 0, (a,b,c,d), Dp; |
---|
| 207 | ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c; |
---|
| 208 | I = std(I); |
---|
| 209 | ring S2 = 0, (c,d,b,a), lp; |
---|
| 210 | ideal I = fetch(S1,I); |
---|
| 211 | // I is assumed to be a Dp-Groebner basis. |
---|
| 212 | // We compute a lp-Groebner basis. |
---|
| 213 | ideal J = modfWalk(I,"Dp"); |
---|
| 214 | J; |
---|
| 215 | intvec w = 3,2,1,2; |
---|
| 216 | ring S3 = 0, (c,d,b,a), (a(w),lp); |
---|
| 217 | ideal I = fetch(S1,I); |
---|
| 218 | // I is assumed to be a Dp-Groebner basis. |
---|
| 219 | // We compute a (a(w),lp)-Groebner basis. |
---|
| 220 | ideal J = modfWalk(I,"Dp",w); |
---|
| 221 | J; |
---|
[daf501] | 222 | } |
---|
[e92c6a] | 223 | |
---|
[daf501] | 224 | proc modfrWalk(ideal I, int radius, list #) |
---|
[b7405b] | 225 | "USAGE: modfrWalk(I, radius [, v, w]); I ideal, radius int, v intvec or string, w intvec |
---|
| 226 | If v intvec, then I is assumed to be a standard basis with respect to (a(v),lp). |
---|
| 227 | If v string, then either v="dp" or v="Dp" and I is assumed to be a standard basis with |
---|
| 228 | respect to dp or Dp, respectively. |
---|
| 229 | If w is given, then a standard basis with respect to (a(w),lp) will be computed. Otherwise, |
---|
| 230 | the output will be a standard basis with respect to lp. |
---|
| 231 | If no optional argument is given, I is assumed to be a standard basis with respect to dp |
---|
| 232 | and a standard basis with respect to lp will be computed. |
---|
[daf501] | 233 | RETURN: a standard basis of I |
---|
| 234 | NOTE: The procedure computes a standard basis of I (over the rational |
---|
| 235 | numbers) by using modular methods. |
---|
| 236 | SEE ALSO: modular |
---|
| 237 | EXAMPLE: example modfrWalk; shows an example" |
---|
| 238 | { |
---|
| 239 | /* save options */ |
---|
| 240 | intvec opt = option(get); |
---|
| 241 | option(redSB); |
---|
| 242 | |
---|
| 243 | /* call modular() */ |
---|
| 244 | if (size(#) > 0) { |
---|
[c448cd] | 245 | I = modular("frandwalk", list(I,radius,#), primeTest_std, |
---|
| 246 | deleteUnluckyPrimes_std, pTest_std, finalTest_std); |
---|
[e92c6a] | 247 | } |
---|
[daf501] | 248 | else { |
---|
[c448cd] | 249 | I = modular("frandwalk", list(I,radius), primeTest_std, |
---|
| 250 | deleteUnluckyPrimes_std, pTest_std, finalTest_std); |
---|
[e92c6a] | 251 | } |
---|
| 252 | |
---|
[daf501] | 253 | /* return the result */ |
---|
| 254 | attrib(I, "isSB", 1); |
---|
| 255 | option(set, opt); |
---|
| 256 | return(I); |
---|
[b62e23] | 257 | } |
---|
| 258 | example |
---|
[daf501] | 259 | { |
---|
| 260 | "EXAMPLE:"; |
---|
| 261 | echo = 2; |
---|
| 262 | ring R1 = 0, (x,y,z,t), dp; |
---|
| 263 | ideal I = 3x3+x2+1, 11y5+y3+2, 5z4+z2+4; |
---|
| 264 | I = std(I); |
---|
| 265 | ring R2 = 0, (x,y,z,t), lp; |
---|
| 266 | ideal I = fetch(R1, I); |
---|
| 267 | int radius = 2; |
---|
| 268 | ideal J = modfrWalk(I,radius); |
---|
| 269 | J; |
---|
[3308a0] | 270 | ring S1 = 0, (a,b,c,d), Dp; |
---|
| 271 | ideal I = 5b2, ac2+9d3+3a2+5b, 2a2c+7abd+bcd+4a2, 2ad2+6b2d+7c3+8ad+4c; |
---|
| 272 | I = std(I); |
---|
| 273 | ring S2 = 0, (c,d,b,a), lp; |
---|
| 274 | ideal I = fetch(S1,I); |
---|
| 275 | // I is assumed to be a Dp-Groebner basis. |
---|
| 276 | // We compute a lp-Groebner basis. |
---|
| 277 | ideal J = modfrWalk(I,radius,"Dp"); |
---|
| 278 | J; |
---|
| 279 | intvec w = 3,2,1,2; |
---|
| 280 | ring S3 = 0, (c,d,b,a), (a(w),lp); |
---|
| 281 | ideal I = fetch(S1,I); |
---|
| 282 | // I is assumed to be a Dp-Groebner basis. |
---|
| 283 | // We compute a (a(w),lp)-Groebner basis. |
---|
| 284 | ideal J = modfrWalk(I,radius,"Dp",w); |
---|
| 285 | J; |
---|
[c448cd] | 286 | } |
---|
| 287 | |
---|
| 288 | /* test if the prime p is suitable for the input, i.e. it does not divide |
---|
| 289 | * the numerator or denominator of any of the coefficients */ |
---|
| 290 | static proc primeTest_std(int p, alias list args) |
---|
| 291 | { |
---|
| 292 | /* erase zero generators */ |
---|
| 293 | ideal I = simplify(args[1], 2); |
---|
| 294 | |
---|
| 295 | /* clear denominators and count the terms */ |
---|
| 296 | ideal J; |
---|
| 297 | ideal K; |
---|
| 298 | int n = ncols(I); |
---|
| 299 | intvec sizes; |
---|
| 300 | number cnt; |
---|
| 301 | int i; |
---|
| 302 | for(i = n; i > 0; i--) { |
---|
| 303 | J[i] = cleardenom(I[i]); |
---|
| 304 | cnt = leadcoef(J[i])/leadcoef(I[i]); |
---|
| 305 | K[i] = numerator(cnt)*var(1)+denominator(cnt); |
---|
| 306 | } |
---|
| 307 | sizes = size(J[1..n]); |
---|
| 308 | |
---|
| 309 | /* change to characteristic p */ |
---|
| 310 | def br = basering; |
---|
| 311 | list lbr = ringlist(br); |
---|
| 312 | if (typeof(lbr[1]) == "int") { |
---|
| 313 | lbr[1] = p; |
---|
| 314 | } |
---|
| 315 | else { |
---|
| 316 | lbr[1][1] = p; |
---|
| 317 | } |
---|
| 318 | def rp = ring(lbr); |
---|
| 319 | setring(rp); |
---|
| 320 | ideal Jp = fetch(br, J); |
---|
| 321 | ideal Kp = fetch(br, K); |
---|
| 322 | |
---|
| 323 | /* test if any coefficient is missing */ |
---|
| 324 | if (intvec(size(Kp[1..n])) != 2:n) { |
---|
| 325 | setring(br); |
---|
| 326 | return(0); |
---|
| 327 | } |
---|
| 328 | if (intvec(size(Jp[1..n])) != sizes) { |
---|
| 329 | setring(br); |
---|
| 330 | return(0); |
---|
| 331 | } |
---|
| 332 | setring(br); |
---|
| 333 | return(1); |
---|
| 334 | } |
---|
| 335 | |
---|
| 336 | /* find entries in modresults which come from unlucky primes. |
---|
| 337 | * For this, sort the entries into categories depending on their leading |
---|
| 338 | * ideal and return the indices in all but the biggest category. */ |
---|
| 339 | static proc deleteUnluckyPrimes_std(alias list modresults) |
---|
| 340 | { |
---|
| 341 | int size_modresults = size(modresults); |
---|
| 342 | |
---|
| 343 | /* sort results into categories. |
---|
| 344 | * each category is represented by three entries: |
---|
| 345 | * - the corresponding leading ideal |
---|
| 346 | * - the number of elements |
---|
| 347 | * - the indices of the elements |
---|
| 348 | */ |
---|
| 349 | list cat; |
---|
| 350 | int size_cat; |
---|
| 351 | ideal L; |
---|
| 352 | int i; |
---|
| 353 | int j; |
---|
| 354 | for (i = 1; i <= size_modresults; i++) { |
---|
| 355 | L = lead(modresults[i]); |
---|
| 356 | attrib(L, "isSB", 1); |
---|
| 357 | for (j = 1; j <= size_cat; j++) { |
---|
| 358 | if (size(L) == size(cat[j][1]) |
---|
[6518eba] | 359 | && size(reduce(L, cat[j][1],5)) == 0 |
---|
| 360 | && size(reduce(cat[j][1], L,5)) == 0) |
---|
| 361 | { |
---|
[c448cd] | 362 | cat[j][2] = cat[j][2]+1; |
---|
| 363 | cat[j][3][cat[j][2]] = i; |
---|
| 364 | break; |
---|
| 365 | } |
---|
| 366 | } |
---|
| 367 | if (j > size_cat) { |
---|
| 368 | size_cat++; |
---|
| 369 | cat[size_cat] = list(); |
---|
| 370 | cat[size_cat][1] = L; |
---|
| 371 | cat[size_cat][2] = 1; |
---|
| 372 | cat[size_cat][3] = list(i); |
---|
| 373 | } |
---|
| 374 | } |
---|
| 375 | |
---|
| 376 | /* find the biggest categories */ |
---|
| 377 | int cat_max = 1; |
---|
| 378 | int max = cat[1][2]; |
---|
| 379 | for (i = 2; i <= size_cat; i++) { |
---|
| 380 | if (cat[i][2] > max) { |
---|
| 381 | cat_max = i; |
---|
| 382 | max = cat[i][2]; |
---|
| 383 | } |
---|
| 384 | } |
---|
| 385 | |
---|
| 386 | /* return all other indices */ |
---|
| 387 | list unluckyIndices; |
---|
| 388 | for (i = 1; i <= size_cat; i++) { |
---|
| 389 | if (i != cat_max) { |
---|
| 390 | unluckyIndices = unluckyIndices + cat[i][3]; |
---|
| 391 | } |
---|
| 392 | } |
---|
| 393 | return(unluckyIndices); |
---|
| 394 | } |
---|
| 395 | |
---|
| 396 | /* test if 'command' applied to 'args' in characteristic p is the same as |
---|
| 397 | 'result' mapped to characteristic p */ |
---|
| 398 | static proc pTest_std(string command, alias list args, alias ideal result, |
---|
| 399 | int p) |
---|
| 400 | { |
---|
| 401 | /* change to characteristic p */ |
---|
| 402 | def br = basering; |
---|
| 403 | list lbr = ringlist(br); |
---|
| 404 | if (typeof(lbr[1]) == "int") { |
---|
| 405 | lbr[1] = p; |
---|
| 406 | } |
---|
| 407 | else { |
---|
| 408 | lbr[1][1] = p; |
---|
| 409 | } |
---|
| 410 | def rp = ring(lbr); |
---|
| 411 | setring(rp); |
---|
| 412 | ideal Ip = fetch(br, args)[1]; |
---|
| 413 | list Arg = fetch(br, args); |
---|
| 414 | string exstr; |
---|
| 415 | ideal 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 | if (reduce(Ip[i], Gp, 1) != 0) { |
---|
| 422 | setring(br); |
---|
| 423 | return(0); |
---|
| 424 | } |
---|
| 425 | } |
---|
| 426 | |
---|
| 427 | /* compute command(args) */ |
---|
| 428 | exstr = "Ip = "+command+" (Ip"; |
---|
[bc55ff] | 429 | |
---|
[c448cd] | 430 | for(i=2; i<=size(Arg); i++) { |
---|
| 431 | exstr = exstr+",Arg["+string(eval(i))+"]"; |
---|
| 432 | } |
---|
| 433 | exstr = exstr+");"; |
---|
[bc55ff] | 434 | |
---|
[c448cd] | 435 | execute(exstr); |
---|
| 436 | |
---|
| 437 | /* test if Gp is in Ip */ |
---|
| 438 | for (i = ncols(Gp); i > 0; i--) { |
---|
| 439 | if (reduce(Gp[i], Ip, 1) != 0) { |
---|
| 440 | setring(br); |
---|
| 441 | return(0); |
---|
| 442 | } |
---|
| 443 | } |
---|
| 444 | |
---|
| 445 | setring(br); |
---|
| 446 | return(1); |
---|
| 447 | } |
---|
| 448 | |
---|
| 449 | /* test if 'result' is a GB of the input ideal */ |
---|
| 450 | static proc finalTest_std(string command, alias list args, ideal result) |
---|
| 451 | { |
---|
| 452 | /* test if args[1] is in result */ |
---|
| 453 | attrib(result, "isSB", 1); |
---|
| 454 | int i; |
---|
| 455 | for (i = ncols(args[1]); i > 0; i--) { |
---|
| 456 | if (reduce(args[1][i], result, 1) != 0) { |
---|
| 457 | return(0); |
---|
| 458 | } |
---|
| 459 | } |
---|
| 460 | |
---|
[b7405b] | 461 | /* test if result is in args[1]. */ |
---|
| 462 | /* args[1] is given by a Groebner basis. Thus we may */ |
---|
| 463 | /* reduce the result with respect to args[1]. */ |
---|
| 464 | int n=nvars(basering); |
---|
| 465 | string ord_str = "dp"; |
---|
| 466 | |
---|
| 467 | for(i=2; i<=size(args); i++) |
---|
| 468 | { |
---|
| 469 | if(typeof(args[i]) == "list") { |
---|
| 470 | if(typeof(args[i][1]) == "intvec") { |
---|
| 471 | ord_str = "(a("+string(args[i][1])+"),lp("+string(n) + "),C)"; |
---|
| 472 | break; |
---|
| 473 | } |
---|
| 474 | if(typeof(args[i][1]) == "string") { |
---|
| 475 | if(args[i][1] == "Dp") { |
---|
| 476 | ord_str = "Dp"; |
---|
| 477 | } |
---|
| 478 | break; |
---|
| 479 | } |
---|
| 480 | } |
---|
| 481 | } |
---|
| 482 | ideal xI = args[1]; |
---|
| 483 | ring xR = basering; |
---|
[67b671] | 484 | ring yR = create_ring(ringlist(xR)[1], "("+varstr(xR)+")", ord_str, "no_minpoly"); |
---|
[b7405b] | 485 | ideal yI = fetch(xR,xI); |
---|
| 486 | ideal yresult = fetch(xR,result); |
---|
| 487 | attrib(yI, "isSB", 1); |
---|
| 488 | for(i=size(yresult); i>0; i--) |
---|
| 489 | { |
---|
| 490 | if(reduce(yresult[i],yI) != 0) |
---|
| 491 | { |
---|
| 492 | return(0); |
---|
| 493 | } |
---|
| 494 | } |
---|
| 495 | setring xR; |
---|
| 496 | kill yR; |
---|
| 497 | |
---|
| 498 | /* test if result is a Groebner basis */ |
---|
[60c3fb] | 499 | link l1="ssi:fork"; |
---|
| 500 | open(l1); |
---|
| 501 | link l2="ssi:fork"; |
---|
| 502 | open(l2); |
---|
| 503 | list l=list(l1,l2); |
---|
| 504 | write(l1,quote(TestSBred(result))); |
---|
| 505 | write(l2,quote(TestSBstd(result))); |
---|
| 506 | i=waitfirst(l); |
---|
| 507 | if(i==1) { |
---|
| 508 | i=read(l1); |
---|
| 509 | } |
---|
| 510 | else { |
---|
| 511 | i=read(l2); |
---|
| 512 | } |
---|
| 513 | close(l1); |
---|
| 514 | close(l2); |
---|
| 515 | return(i); |
---|
[c448cd] | 516 | } |
---|
| 517 | |
---|
| 518 | /* return 1, if I_reduce is _not_ in G_reduce, |
---|
| 519 | * 0, otherwise |
---|
| 520 | * (same as size(reduce(I_reduce, G_reduce))). |
---|
| 521 | * Uses parallelization. */ |
---|
| 522 | static proc reduce_parallel(def I_reduce, def G_reduce) |
---|
| 523 | { |
---|
[b7405b] | 524 | exportto(Modwalk, I_reduce); |
---|
| 525 | exportto(Modwalk, G_reduce); |
---|
[c448cd] | 526 | int size_I = ncols(I_reduce); |
---|
| 527 | int chunks = Modular::par_range(size_I); |
---|
| 528 | intvec range; |
---|
| 529 | int i; |
---|
| 530 | for (i = chunks; i > 0; i--) { |
---|
| 531 | range = Modular::par_range(size_I, i); |
---|
[b7405b] | 532 | task t(i) = "Modwalk::reduce_task", list(range); |
---|
[c448cd] | 533 | } |
---|
| 534 | startTasks(t(1..chunks)); |
---|
| 535 | waitAllTasks(t(1..chunks)); |
---|
| 536 | int result = 0; |
---|
| 537 | for (i = chunks; i > 0; i--) { |
---|
| 538 | if (getResult(t(i))) { |
---|
| 539 | result = 1; |
---|
| 540 | break; |
---|
| 541 | } |
---|
| 542 | } |
---|
| 543 | kill I_reduce; |
---|
| 544 | kill G_reduce; |
---|
| 545 | return(result); |
---|
| 546 | } |
---|
| 547 | |
---|
| 548 | /* compute a chunk of reductions for reduce_parallel */ |
---|
| 549 | static proc reduce_task(intvec range) |
---|
| 550 | { |
---|
| 551 | int result = 0; |
---|
| 552 | int i; |
---|
| 553 | for (i = range[1]; i <= range[2]; i++) { |
---|
| 554 | if (reduce(I_reduce[i], G_reduce, 1) != 0) { |
---|
| 555 | result = 1; |
---|
| 556 | break; |
---|
| 557 | } |
---|
| 558 | } |
---|
| 559 | return(result); |
---|
[60c3fb] | 560 | } |
---|
| 561 | |
---|
| 562 | /* test if result is a GB with std*/ |
---|
| 563 | static proc TestSBstd(ideal result) |
---|
| 564 | { |
---|
| 565 | ideal G = std(result); |
---|
| 566 | if(reduce_parallel(G,result)) { |
---|
| 567 | return(0); |
---|
| 568 | } |
---|
| 569 | return(1); |
---|
| 570 | } |
---|
| 571 | |
---|
| 572 | /* test if result is a GB by reducing s-polynomials*/ |
---|
| 573 | static proc TestSBred(ideal result) |
---|
| 574 | { |
---|
| 575 | int i,j; |
---|
| 576 | for(i=1; i<=size(result); i++) { |
---|
| 577 | for(j=i; j<=size(result); j++) { |
---|
| 578 | if(reduce(sPolynomial(result[i],result[j]),result)!=0) { |
---|
| 579 | return(0); |
---|
| 580 | } |
---|
| 581 | } |
---|
| 582 | } |
---|
| 583 | return(1); |
---|
| 584 | } |
---|
| 585 | |
---|
| 586 | /* compute s-polynomial of f and g */ |
---|
| 587 | static proc sPolynomial(poly f,poly g) |
---|
| 588 | { |
---|
| 589 | int i; |
---|
| 590 | poly lcmp = 1; |
---|
| 591 | |
---|
| 592 | intvec lexpf = leadexp(f); |
---|
| 593 | intvec lexpg = leadexp(g); |
---|
| 594 | |
---|
| 595 | for(i=1; i<=nvars(basering); i++) { |
---|
| 596 | if(lexpf[i]>=lexpg[i]) { |
---|
| 597 | lcmp=lcmp*var(i)**lexpf[i]; |
---|
| 598 | } |
---|
| 599 | else { |
---|
| 600 | lcmp=lcmp*var(i)**lexpg[i]; |
---|
| 601 | } |
---|
| 602 | } |
---|
| 603 | |
---|
| 604 | poly fmult=lcmp/leadmonom(f); |
---|
| 605 | poly gmult=lcmp/leadmonom(g); |
---|
| 606 | poly result=leadcoef(g)*fmult*f-leadcoef(f)*gmult*g; |
---|
| 607 | |
---|
| 608 | return(result); |
---|
| 609 | } |
---|