[e27e47] | 1 | ///////////////////////////////////////////////////////////////////// |
---|
| 2 | |
---|
[341696] | 3 | version="$Id$"; |
---|
[e27e47] | 4 | category="Commutative Algebra"; |
---|
| 5 | |
---|
| 6 | info=" |
---|
| 7 | LIBRARY: EHV.lib PROCEDURES FOR PRIMARY DECOMPOSITION OF IDEALS |
---|
| 8 | AUTHORS: Kai Dehmann, dehmann@mathematik.uni-kl.de; |
---|
| 9 | |
---|
| 10 | OVERVIEW: |
---|
| 11 | Algorithms for primary decomposition and radical-computation |
---|
| 12 | based on the ideas of Eisenbud, Huneke, and Vasconcelos. |
---|
| 13 | |
---|
| 14 | PROCEDURES: |
---|
| 15 | equiMaxEHV(I); equidimensional part of I |
---|
| 16 | removeComponent(I,e); intersection of the primary components |
---|
| 17 | of I of dimension >= e |
---|
| 18 | AssOfDim(I,e); an ideal such that the associated primes |
---|
| 19 | are exactly the associated primes of I |
---|
| 20 | having dimension e |
---|
| 21 | equiRadEHV(I [,Strategy]); equidimensional radical of I |
---|
| 22 | radEHV(I [,Strategy]); radical of I |
---|
| 23 | IntAssOfDim1(I,e); intersection of the associated primes of I |
---|
| 24 | having dimension e |
---|
| 25 | IntAssOfDim2(I,e); another way of computing the intersection |
---|
| 26 | of the associated primes of I |
---|
| 27 | having dimension e |
---|
| 28 | decompEHV(I); decomposition of a zero-dimensional |
---|
| 29 | radical ideal I |
---|
| 30 | AssEHV(I [,Strategy]); associated primes of I |
---|
| 31 | minAssEHV(I [,Strategy]); minimal associated primes of I |
---|
| 32 | localize(I,P,l); the contraction of the ideal generated by I |
---|
| 33 | in the localization w.r.t P |
---|
| 34 | componentEHV(I,P,L [,Strategy]); a P-primary component for I |
---|
| 35 | primdecEHV(I [,Strategy]); a minimal primary decomposition of I |
---|
| 36 | compareLists(L, K); procedure for comparing the output of |
---|
| 37 | primary decomposition algorithms (checks |
---|
| 38 | if the computed associated primes coincide) |
---|
| 39 | "; |
---|
| 40 | |
---|
| 41 | LIB "ring.lib"; |
---|
| 42 | LIB "general.lib"; |
---|
| 43 | LIB "elim.lib"; |
---|
| 44 | LIB "poly.lib"; |
---|
| 45 | LIB "random.lib"; |
---|
| 46 | LIB "inout.lib"; |
---|
| 47 | LIB "matrix.lib"; |
---|
| 48 | LIB "algebra.lib"; |
---|
[a7fcfab] | 49 | LIB "normal.lib"; |
---|
[e27e47] | 50 | |
---|
| 51 | |
---|
| 52 | ///////////////////////////////////////////////////////////////////// |
---|
| 53 | // // |
---|
| 54 | // G E N E R A L A L G O R I T H M S // |
---|
| 55 | // // |
---|
| 56 | ///////////////////////////////////////////////////////////////////// |
---|
| 57 | |
---|
| 58 | ///////////////////////////////////////////////////////////////////// |
---|
| 59 | static proc AnnExtEHV(int n,list re) |
---|
| 60 | "USAGE: AnnExtEHV(n,re); n integer, re resolution |
---|
| 61 | RETURN: ideal, the annihilator of Ext^n(R/I,R) with given |
---|
| 62 | resolution re of I" |
---|
| 63 | { |
---|
| 64 | if(printlevel > 2){"Entering AnnExtEHV.";} |
---|
| 65 | |
---|
| 66 | if(n < 0) |
---|
| 67 | { |
---|
| 68 | ideal ann = ideal(1); |
---|
| 69 | if(printlevel > 2){"Leaving AnnExtEHV.";} |
---|
| 70 | return(ann); |
---|
| 71 | } |
---|
| 72 | int l = size(re); |
---|
| 73 | |
---|
| 74 | if(n < l) |
---|
| 75 | { |
---|
| 76 | matrix f = transpose(re[n+1]); |
---|
| 77 | if(n == 0) |
---|
| 78 | { |
---|
| 79 | matrix g = 0*gen(ncols(f)); |
---|
| 80 | } |
---|
| 81 | else |
---|
| 82 | { |
---|
| 83 | matrix g = transpose(re[n]); |
---|
| 84 | } |
---|
| 85 | module k = syz(f); |
---|
| 86 | ideal ann = quotient1(g,k); |
---|
| 87 | if(printlevel > 2){"Leaving AnnExtEHV.";} |
---|
| 88 | return(ann); |
---|
| 89 | } |
---|
| 90 | |
---|
| 91 | if(n == l) |
---|
| 92 | { |
---|
| 93 | ideal ann = Ann(transpose(re[n])); |
---|
| 94 | if(printlevel > 2){"Leaving AnnExtEHV.";} |
---|
| 95 | return(ann); |
---|
| 96 | } |
---|
| 97 | |
---|
| 98 | ideal ann = ideal(1); |
---|
| 99 | if(printlevel > 2){"Leaving AnnExtEHV.";} |
---|
| 100 | return(ann); |
---|
| 101 | } |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | ///////////////////////////////////////////////////////////////////// |
---|
| 105 | //static |
---|
| 106 | proc isSubset(ideal I,ideal J) |
---|
| 107 | "USAGE: isSubset(I,J); I, J ideals |
---|
| 108 | RETURN: integer, 1 if I is a subset of J and 0 otherwise |
---|
| 109 | NOTE: if J is not a standard basis the result may be wrong |
---|
| 110 | " |
---|
| 111 | { |
---|
| 112 | int s = size(I); |
---|
| 113 | for(int i=1; i<=s; i++) |
---|
| 114 | { |
---|
| 115 | if(reduce(I[i],J)!=0) |
---|
| 116 | { |
---|
| 117 | return(0); |
---|
| 118 | } |
---|
| 119 | } |
---|
| 120 | return(1); |
---|
| 121 | } |
---|
| 122 | |
---|
| 123 | |
---|
| 124 | ///////////////////////////////////////////////////////////////////// |
---|
| 125 | // // |
---|
| 126 | // T H E E Q U I D I M E N S I O N A L P A R T // |
---|
| 127 | // // |
---|
| 128 | ///////////////////////////////////////////////////////////////////// |
---|
| 129 | |
---|
| 130 | ///////////////////////////////////////////////////////////////////// |
---|
| 131 | proc equiMaxEHV(ideal I) |
---|
| 132 | "USAGE: equiMaxEHV(I); I ideal |
---|
| 133 | RETURN: ideal, the equidimensional part of I. |
---|
| 134 | NOTE: Uses algorithm of Eisenbud, Huneke, and Vasconcelos. |
---|
| 135 | EXAMPLE: example equiMaxEHV; shows an example |
---|
| 136 | " |
---|
| 137 | { |
---|
| 138 | if(printlevel > 2){"Entering equiMaxEHV.";} |
---|
| 139 | if(ord_test(basering)!=1) |
---|
| 140 | { |
---|
| 141 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 142 | } |
---|
| 143 | ideal J = groebner(I); |
---|
| 144 | int cod = nvars(basering)-dim(J); |
---|
| 145 | |
---|
| 146 | //If I is the entire ring... |
---|
| 147 | if(cod > nvars(basering)) |
---|
| 148 | { |
---|
| 149 | //...then return the ideal generated by 1. |
---|
| 150 | return(ideal(1)); |
---|
| 151 | } |
---|
| 152 | |
---|
| 153 | //Compute a resolution of I. |
---|
| 154 | if(printlevel > 2){"Computing resolution.";} |
---|
| 155 | if(homog(I)==1) |
---|
| 156 | { |
---|
| 157 | list re = sres(J,cod+1); |
---|
| 158 | re = minres(re); |
---|
| 159 | } |
---|
| 160 | else |
---|
| 161 | { |
---|
| 162 | list re = mres(I,cod+1); |
---|
| 163 | } |
---|
| 164 | if(printlevel > 2){"Finished computing resolution.";} |
---|
| 165 | |
---|
| 166 | //Compute the annihilator of the cod-th EXT-module. |
---|
| 167 | ideal ann = AnnExtEHV(cod,re); |
---|
| 168 | attrib(ann,"isEquidimensional",1); |
---|
| 169 | if(printlevel > 2){"Leaving equiMaxEHV.";} |
---|
| 170 | return(ann); |
---|
| 171 | } |
---|
| 172 | |
---|
| 173 | example |
---|
| 174 | { |
---|
| 175 | "EXAMPLE:"; |
---|
| 176 | echo = 2; |
---|
| 177 | ring r = 0,(x,y,z),dp; |
---|
| 178 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
| 179 | equiMaxEHV(I); |
---|
| 180 | } |
---|
| 181 | |
---|
| 182 | |
---|
| 183 | ///////////////////////////////////////////////////////////////////// |
---|
| 184 | proc removeComponent(ideal I, int e) |
---|
| 185 | "USAGE: removeComponent(I,e); I ideal, e integer |
---|
| 186 | RETURN: ideal, the intersection of the primary components |
---|
| 187 | of I of dimension >= e |
---|
| 188 | EXAMPLE: example removeComponent; shows an example" |
---|
| 189 | { |
---|
| 190 | if(ord_test(basering)!=1) |
---|
| 191 | { |
---|
| 192 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 193 | } |
---|
| 194 | |
---|
| 195 | ideal J = groebner(I); |
---|
| 196 | |
---|
| 197 | //Compute a resolution of I |
---|
| 198 | if(homog(I)==1) |
---|
| 199 | { |
---|
| 200 | list re = sres(J,0); |
---|
| 201 | re = minres(re); |
---|
| 202 | } |
---|
| 203 | else |
---|
| 204 | { |
---|
| 205 | list re = mres(I,0); |
---|
| 206 | } |
---|
| 207 | |
---|
| 208 | int f = nvars(basering); |
---|
| 209 | int cod; |
---|
| 210 | ideal ann; |
---|
| 211 | int g = nvars(basering) - e; |
---|
| 212 | while(f > g) |
---|
| 213 | { |
---|
| 214 | ann = AnnExtEHV(f,re); |
---|
| 215 | cod = nvars(basering) - dim(groebner(ann)); |
---|
| 216 | if( cod == f ) |
---|
| 217 | { |
---|
| 218 | I = quotient(I,ann); |
---|
| 219 | } |
---|
| 220 | f = f-1; |
---|
| 221 | } |
---|
| 222 | return(I); |
---|
| 223 | } |
---|
| 224 | |
---|
| 225 | example |
---|
| 226 | { |
---|
| 227 | "EXAMPLE:"; |
---|
| 228 | echo = 2; |
---|
| 229 | ring r = 0,(x,y,z),dp; |
---|
| 230 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
| 231 | removeComponent(I,1); |
---|
| 232 | } |
---|
| 233 | |
---|
| 234 | |
---|
| 235 | ///////////////////////////////////////////////////////////////////// |
---|
| 236 | proc AssOfDim(ideal I, int e) |
---|
| 237 | "USAGE: AssOfDim(I,e); I ideal, e integer |
---|
| 238 | RETURN: ideal, such that the associated primes are exactly |
---|
| 239 | the associated primes of I having dimension e |
---|
| 240 | EXAMPLE: example AssOfDim; shows an example" |
---|
| 241 | { |
---|
| 242 | if(ord_test(basering)!=1) |
---|
| 243 | { |
---|
| 244 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 245 | } |
---|
| 246 | int g = nvars(basering) - e; |
---|
| 247 | |
---|
| 248 | //Compute a resolution of I. |
---|
| 249 | ideal J = std(I); |
---|
| 250 | if(homog(I)==1) |
---|
| 251 | { |
---|
| 252 | list re = sres(J,g+1); |
---|
| 253 | re = minres(re); |
---|
| 254 | } |
---|
| 255 | else |
---|
| 256 | { |
---|
| 257 | list re = mres(I,g+1); |
---|
| 258 | } |
---|
| 259 | |
---|
| 260 | ideal ann = AnnExtEHV(g,re); |
---|
| 261 | int cod = nvars(basering) - dim(std(ann)); |
---|
| 262 | |
---|
| 263 | //If the codimension of I_g:=Ann(Ext^g(R/I,R)) equals g... |
---|
| 264 | if(cod == g) |
---|
| 265 | { |
---|
| 266 | //...then return the equidimensional part of I_g... |
---|
| 267 | ann = equiMaxEHV(ann); |
---|
| 268 | attrib(ann,"isEquidimensional",1); |
---|
| 269 | return(ann); |
---|
| 270 | } |
---|
| 271 | //...otherwise... |
---|
| 272 | else |
---|
| 273 | { |
---|
| 274 | //...I has no associated primes of dimension e. |
---|
| 275 | return(ideal(1)); |
---|
| 276 | } |
---|
| 277 | } |
---|
| 278 | |
---|
| 279 | example |
---|
| 280 | { |
---|
| 281 | "EXAMPLE:"; |
---|
| 282 | echo = 2; |
---|
| 283 | ring r = 0,(x,y,z),dp; |
---|
| 284 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
| 285 | AssOfDim(I,1); |
---|
| 286 | } |
---|
| 287 | |
---|
| 288 | ///////////////////////////////////////////////////////////////////// |
---|
| 289 | // // |
---|
| 290 | // T H E R A D I C A L // |
---|
| 291 | // // |
---|
| 292 | ///////////////////////////////////////////////////////////////////// |
---|
| 293 | |
---|
| 294 | ///////////////////////////////////////////////////////////////////// |
---|
| 295 | static proc aJacob(ideal I, int a) |
---|
| 296 | "USAGE: aJacob(I,a); I ideal, a integer |
---|
| 297 | RETURN: ideal, the ath-Jacobian ideal of I" |
---|
| 298 | { |
---|
| 299 | matrix M = jacob(I); |
---|
| 300 | int n = nvars(basering); |
---|
| 301 | if(n-a <= 0) |
---|
| 302 | { |
---|
| 303 | return(ideal(1)); |
---|
| 304 | } |
---|
| 305 | if(n-a > nrows(M) or n-a > ncols(M)) |
---|
| 306 | { |
---|
| 307 | return(ideal(0)); |
---|
| 308 | } |
---|
| 309 | ideal J = minor(M,n-a); |
---|
| 310 | return(J); |
---|
| 311 | } |
---|
| 312 | |
---|
| 313 | |
---|
| 314 | ///////////////////////////////////////////////////////////////////// |
---|
| 315 | proc equiRadEHV(ideal I, list #) |
---|
| 316 | "USAGE: equiRadEHV(I [,Strategy]); I ideal, Strategy list |
---|
| 317 | RETURN: ideal, the equidimensional radical of I, |
---|
| 318 | i.e. the intersection of the minimal associated primes of I |
---|
| 319 | having the same dimension as I |
---|
| 320 | NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos, |
---|
| 321 | Works only in characteristic 0 or p large. |
---|
| 322 | The (optional) second argument determines the strategy used: |
---|
| 323 | Strategy[1] > strategy for the equidimensional part |
---|
| 324 | = 0 : uses equiMaxEHV |
---|
| 325 | = 1 : uses equidimMax |
---|
| 326 | Strategy[2] > strategy for the radical |
---|
| 327 | = 0 : combination of strategy 1 and 2 |
---|
| 328 | = 1 : computation of the radical just with the |
---|
| 329 | help of regular sequences |
---|
| 330 | = 2 : does not try to find a regular sequence |
---|
| 331 | Strategy[3] > strategy for the computation of ideal quotients |
---|
| 332 | = n : uses quot(.,.,n) for the ideal quotient computations |
---|
| 333 | If no second argument is given then Strategy=(0,0,0) is used. |
---|
| 334 | EXAMPLE: example equiRadEHV; shows an example" |
---|
| 335 | { |
---|
| 336 | if(printlevel > 2){"Entering equiRadEHV.";} |
---|
| 337 | if(ord_test(basering)!=1) |
---|
| 338 | { |
---|
| 339 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 340 | } |
---|
| 341 | if((char(basering)<100)&&(char(basering)!=0)) |
---|
| 342 | { |
---|
| 343 | "WARNING: The characteristic is too small, the result may be wrong"; |
---|
| 344 | } |
---|
| 345 | |
---|
| 346 | //Define the Strategy to be used. |
---|
| 347 | if(size(#) > 0) |
---|
| 348 | { |
---|
| 349 | if(#[1]!=1) |
---|
| 350 | { |
---|
| 351 | int equStr = 0; |
---|
| 352 | } |
---|
| 353 | else |
---|
| 354 | { |
---|
| 355 | int equStr = 1; |
---|
| 356 | } |
---|
| 357 | if(size(#) > 1) |
---|
| 358 | { |
---|
| 359 | if(#[2]!=1 and #[2]!=2) |
---|
| 360 | { |
---|
| 361 | int strategy = 0; |
---|
| 362 | } |
---|
| 363 | else |
---|
| 364 | { |
---|
| 365 | int strategy = #[2]; |
---|
| 366 | } |
---|
| 367 | if(size(#) > 2) |
---|
| 368 | { |
---|
| 369 | int quoStr = #[3]; |
---|
| 370 | } |
---|
| 371 | else |
---|
| 372 | { |
---|
| 373 | int quoStr = 0; |
---|
| 374 | } |
---|
| 375 | } |
---|
| 376 | else |
---|
| 377 | { |
---|
| 378 | int strategy = 0; |
---|
| 379 | int quoStr = 0; |
---|
| 380 | } |
---|
| 381 | } |
---|
| 382 | else |
---|
| 383 | { |
---|
| 384 | int equStr = 0; |
---|
| 385 | int strategy = 0; |
---|
| 386 | int quoStr = 0; |
---|
| 387 | } |
---|
| 388 | |
---|
| 389 | ideal J,I0,radI0,L,radI1,I2,radI2; |
---|
| 390 | int l,n; |
---|
| 391 | intvec op = option(get); |
---|
| 392 | matrix M; |
---|
| 393 | |
---|
| 394 | option(redSB); |
---|
| 395 | list m = mstd(I); |
---|
| 396 | option(set,op); |
---|
| 397 | |
---|
| 398 | int d = dim(m[1]); |
---|
| 399 | if(d==-1) |
---|
| 400 | { |
---|
| 401 | return(ideal(1)); |
---|
| 402 | } |
---|
| 403 | |
---|
| 404 | if(strategy != 2) |
---|
| 405 | { |
---|
| 406 | ///////////////////////////////////////////// |
---|
| 407 | // Computing the equidimensional radical // |
---|
| 408 | // via regular sequenves // |
---|
| 409 | ///////////////////////////////////////////// |
---|
| 410 | |
---|
| 411 | if(printlevel > 2){"Trying to find a regular sequence.";} |
---|
| 412 | int cod = nvars(basering)-d; |
---|
| 413 | |
---|
| 414 | //Complete intersection case: |
---|
| 415 | if(cod==size(m[2])) |
---|
| 416 | { |
---|
| 417 | J = aJacob(m[2],d); |
---|
| 418 | if(printlevel > 2){"Leaving equiRadEHV.";} |
---|
| 419 | return(quot1(m[2],J,quoStr)); |
---|
| 420 | } |
---|
| 421 | |
---|
| 422 | //First codim elements of I are a complete intersection: |
---|
| 423 | for(l=1; l<=cod; l++) |
---|
| 424 | { |
---|
| 425 | I0[l] = m[2][l]; |
---|
| 426 | } |
---|
| 427 | n = dim(groebner(I0))+cod-nvars(basering); |
---|
| 428 | |
---|
| 429 | //Last codim elements of I are a complete intersection: |
---|
| 430 | if(n!=0) |
---|
| 431 | { |
---|
| 432 | for(l=1; l<=cod; l++) |
---|
| 433 | { |
---|
| 434 | I0[l] = m[2][size(m[2])-l+1]; |
---|
| 435 | } |
---|
| 436 | n = dim(groebner(I0))+cod-nvars(basering); |
---|
| 437 | } |
---|
| 438 | |
---|
| 439 | //Taking a generic linear combination of the input: |
---|
| 440 | if(n!=0) |
---|
| 441 | { |
---|
| 442 | M = transpose(sparsetriag(size(m[2]),cod,95,1)); |
---|
| 443 | I0 = ideal(M*transpose(m[2])); |
---|
| 444 | n = dim(groebner(I0))+cod-nvars(basering); |
---|
| 445 | } |
---|
| 446 | |
---|
| 447 | //Taking a more generic linear combination of the input: |
---|
| 448 | if(n!=0) |
---|
| 449 | { |
---|
| 450 | while(strategy == 1 and n!=0) |
---|
| 451 | { |
---|
| 452 | M = transpose(sparsetriag(size(m[2]),cod,0,100)); |
---|
| 453 | I0 = ideal(M*transpose(m[2])); |
---|
| 454 | n = dim(groebner(I0))+cod-nvars(basering); |
---|
| 455 | } |
---|
| 456 | } |
---|
| 457 | |
---|
| 458 | if(n==0) |
---|
| 459 | { |
---|
| 460 | J = aJacob(I0,d); |
---|
| 461 | if(printlevel > 2){"1st quotient.";} |
---|
| 462 | radI0 = quot1(I0,J,quoStr); |
---|
| 463 | if(printlevel > 2){"2nd quotient.";} |
---|
| 464 | L = quot1(radI0,m[2],quoStr); |
---|
| 465 | if(printlevel > 2){"3rd quotient.";} |
---|
| 466 | radI1 = quot1(radI0,L,quoStr); |
---|
| 467 | attrib(radI1,"isEquidimensional",1); |
---|
| 468 | attrib(radI1,"isRadical",1); |
---|
| 469 | if(printlevel > 2){"Leaving equiRadEHV.";} |
---|
| 470 | return(radI1); |
---|
| 471 | } |
---|
| 472 | } |
---|
| 473 | |
---|
| 474 | //////////////////////////////////////////////////// |
---|
| 475 | // Computing the equidimensional radical directly // |
---|
| 476 | //////////////////////////////////////////////////// |
---|
| 477 | |
---|
| 478 | if(printlevel > 2){"Computing the equidimensional radical directly";} |
---|
| 479 | |
---|
| 480 | //Compute the equidimensional part depending on the chosen strategy |
---|
| 481 | if(equStr == 0) |
---|
| 482 | { |
---|
| 483 | I = equiMaxEHV(I); |
---|
| 484 | } |
---|
| 485 | if(equStr == 1) |
---|
| 486 | { |
---|
| 487 | I = equidimMax(I); |
---|
| 488 | } |
---|
| 489 | int a = nvars(basering)-1; |
---|
| 490 | |
---|
| 491 | while(a > d) |
---|
| 492 | { |
---|
| 493 | if(printlevel > 2){"While-Loop: "+string(a);} |
---|
| 494 | J = aJacob(I,a); |
---|
| 495 | while(dim(groebner(J+I))==d) |
---|
| 496 | { |
---|
| 497 | if(printlevel > 2){"Quotient-Computation.";} |
---|
| 498 | I = quot1(I,J,quoStr); |
---|
| 499 | if(printlevel > 2){"Computing the a-th Jacobian";} |
---|
| 500 | J = aJacob(I,a); |
---|
| 501 | } |
---|
| 502 | a = a-1; |
---|
| 503 | } |
---|
| 504 | if(printlevel > 2){"We left While-Loop.";} |
---|
| 505 | if(printlevel > 2){"Computing the a-th Jacobian";} |
---|
| 506 | J = aJacob(I,d); |
---|
| 507 | if(printlevel > 2){"Quotient-Computation.";} |
---|
| 508 | I = quot1(I,J,quoStr); |
---|
| 509 | attrib(I,"isEquidimensional",1); |
---|
| 510 | attrib(I,"isRadical",1); |
---|
| 511 | if(printlevel > 2){"Leaving equiRadEHV.";} |
---|
| 512 | return(I); |
---|
| 513 | } |
---|
| 514 | |
---|
| 515 | example |
---|
| 516 | { |
---|
| 517 | "EXAMPLE:"; |
---|
| 518 | echo = 2; |
---|
| 519 | ring r = 0,(x,y,z),dp; |
---|
| 520 | poly p = z2+1; |
---|
| 521 | poly q = z3+2; |
---|
| 522 | ideal i = p*q^2,y-z2; |
---|
| 523 | ideal pr= equiRadEHV(i); |
---|
| 524 | pr; |
---|
| 525 | } |
---|
| 526 | |
---|
| 527 | ///////////////////////////////////////////////////////////////////// |
---|
| 528 | proc radEHV(ideal I, list #) |
---|
| 529 | "USAGE: radEHV(I [,Strategy]); ideal I, Strategy list |
---|
| 530 | RETURN: ideal, the radical of I |
---|
| 531 | NOTE: uses the algorithm of Eisenbud/Huneke/Vasconcelos |
---|
| 532 | Works only in characteristic 0 or p large. |
---|
| 533 | The (optional) second argument determines the strategy used: |
---|
| 534 | Strategy[1] > strategy for the equidimensional part |
---|
| 535 | = 0 : uses equiMaxEHV |
---|
| 536 | = 1 : uses equidimMax |
---|
| 537 | Strategy[2] > strategy for the radical |
---|
| 538 | = 0 : combination of strategy 1 and 2 |
---|
| 539 | = 1 : computation of the radical just with the |
---|
| 540 | help of regular sequences |
---|
| 541 | = 2 : does not try to find a regular sequence |
---|
| 542 | Strategy[3] > strategy for the computation of ideal quotients |
---|
| 543 | = n : uses quot(.,.,n) for the ideal quotient computations |
---|
| 544 | If no second argument is given then Strategy=(0,0,0) is used. |
---|
| 545 | EXAMPLE: example radEHV; shows an example" |
---|
| 546 | { |
---|
| 547 | if(printlevel > 2){"Entering radEHV.";} |
---|
| 548 | if(ord_test(basering)!=1) |
---|
| 549 | { |
---|
| 550 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 551 | } |
---|
| 552 | |
---|
| 553 | //Compute the equidimensional radical J of I. |
---|
| 554 | ideal J = equiRadEHV(I,#); |
---|
| 555 | |
---|
| 556 | //If I is the entire ring... |
---|
| 557 | if(deg(J[1]) <= 0) |
---|
| 558 | { |
---|
| 559 | //...then return the ideal generated by 1... |
---|
| 560 | return(ideal(1)); |
---|
| 561 | } |
---|
| 562 | |
---|
| 563 | //...else remove the maximal dimensional components and |
---|
| 564 | //compute the radical K of the lower dimensional components ... |
---|
| 565 | ideal K = radEHV(sat(I,J)[1],#); |
---|
| 566 | |
---|
| 567 | //..and intersect it with J. |
---|
| 568 | K = intersect(J,K); |
---|
| 569 | attrib(K,"isRadical",1); |
---|
| 570 | return(K); |
---|
| 571 | } |
---|
| 572 | |
---|
| 573 | example |
---|
| 574 | { |
---|
| 575 | "EXAMPLE:"; |
---|
| 576 | echo = 2; |
---|
| 577 | ring r = 0,(x,y,z),dp; |
---|
| 578 | poly p = z2+1; |
---|
| 579 | poly q = z3+2; |
---|
| 580 | ideal i = p*q^2,y-z2; |
---|
| 581 | ideal pr= radical(i); |
---|
| 582 | pr; |
---|
| 583 | } |
---|
| 584 | |
---|
| 585 | ///////////////////////////////////////////////////////////////////// |
---|
| 586 | proc IntAssOfDim1(ideal I, int e) |
---|
| 587 | "USAGE: IntAssOfDim1(I,e); I idea, e integer |
---|
| 588 | RETURN: ideal, the intersection of the associated primes of I having dimension e |
---|
| 589 | EXAMPLE: example IntAssOfDim1; shows an example" |
---|
| 590 | { |
---|
| 591 | if(ord_test(basering)!=1) |
---|
| 592 | { |
---|
| 593 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 594 | } |
---|
| 595 | int g = nvars(basering) - e; |
---|
| 596 | |
---|
| 597 | //Compute a resolution of I. |
---|
| 598 | ideal J = groebner(I); |
---|
| 599 | if(homog(I)==1) |
---|
| 600 | { |
---|
| 601 | list re = sres(J,g+1); |
---|
| 602 | re = minres(re); |
---|
| 603 | } |
---|
| 604 | else |
---|
| 605 | { |
---|
| 606 | list re = mres(I,g+1); |
---|
| 607 | } |
---|
| 608 | |
---|
| 609 | ideal ann = AnnExtEHV(g,re); |
---|
| 610 | int cod = nvars(basering) - dim(groebner(ann)); |
---|
| 611 | //If the codimension of I_g:=Ann(Ext^g(R/I,I)) equals g... |
---|
| 612 | if(cod == g) |
---|
| 613 | { |
---|
| 614 | //...then return the equidimensional radical of I_g... |
---|
| 615 | ann = equiRadEHV(ann); |
---|
| 616 | attrib(ann,"isEquidimensional",1); |
---|
| 617 | attrib(ann,"isRadical",1); |
---|
| 618 | return(ann); |
---|
| 619 | } |
---|
| 620 | else |
---|
| 621 | { |
---|
| 622 | //...else I has no associated primes of dimension e. |
---|
| 623 | return(ideal(1)); |
---|
| 624 | } |
---|
| 625 | } |
---|
| 626 | |
---|
| 627 | example |
---|
| 628 | { |
---|
| 629 | "EXAMPLE:"; |
---|
| 630 | echo = 2; |
---|
| 631 | ring r = 0,(x,y,z),dp; |
---|
| 632 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
| 633 | IntAssOfDim1(I,1); |
---|
| 634 | } |
---|
| 635 | |
---|
| 636 | ///////////////////////////////////////////////////////////////////// |
---|
| 637 | proc IntAssOfDim2(ideal I, int e) |
---|
| 638 | "USAGE: IntAssOfDim2(I,e); I ideal, e integer |
---|
| 639 | RETURN: ideal, the intersection of the associated primes of I having dimension e |
---|
| 640 | EXAMPLE: example IntAssOfDim2; shows an example" |
---|
| 641 | { |
---|
| 642 | if(ord_test(basering)!=1) |
---|
| 643 | { |
---|
| 644 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 645 | } |
---|
| 646 | ideal I1 = removeComponent(I,e); |
---|
| 647 | ideal I2 = removeComponent(I,e+1); |
---|
| 648 | ideal b = quotient(I1,I2); |
---|
| 649 | b = equiRadEHV(b); |
---|
| 650 | return(b); |
---|
| 651 | } |
---|
| 652 | |
---|
| 653 | example |
---|
| 654 | { |
---|
| 655 | "EXAMPLE:"; |
---|
| 656 | echo = 2; |
---|
| 657 | ring r = 0,(x,y,z),dp; |
---|
| 658 | ideal I = intersect(ideal(z),ideal(x,y),ideal(x2,z2),ideal(x5,y5,z5)); |
---|
| 659 | IntAssOfDim2(I,1); |
---|
| 660 | } |
---|
| 661 | |
---|
| 662 | |
---|
| 663 | ///////////////////////////////////////////////////////////////////// |
---|
| 664 | // // |
---|
| 665 | // Z E R O - D I M E N S I O N A L D E C O M P O S I T I O N // |
---|
| 666 | // // |
---|
| 667 | ///////////////////////////////////////////////////////////////////// |
---|
| 668 | |
---|
| 669 | ///////////////////////////////////////////////////////////////////// |
---|
| 670 | proc decompEHV(ideal I) |
---|
| 671 | "USAGE: decompEHV(I); I zero-dimensional radical ideal |
---|
| 672 | RETURN: list, the associated primes of I |
---|
| 673 | EXAMPLE: example decompEHV; shows an example" |
---|
| 674 | { |
---|
| 675 | if(printlevel > 2){"Entering decompEHV.";} |
---|
| 676 | if(ord_test(basering)!=1) |
---|
| 677 | { |
---|
| 678 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 679 | } |
---|
| 680 | |
---|
| 681 | int e,m,vd,nfact; |
---|
| 682 | list l,k,L; |
---|
| 683 | poly f,h; |
---|
| 684 | ideal J; |
---|
| 685 | def base = basering; |
---|
| 686 | |
---|
| 687 | if( attrib(I,"isSB")!=1 ) |
---|
| 688 | { |
---|
| 689 | I = groebner(I); |
---|
| 690 | } |
---|
| 691 | |
---|
| 692 | while(1) |
---|
| 693 | { |
---|
| 694 | //Choose a random polynomial f from R. |
---|
| 695 | e = random(0,100); |
---|
| 696 | f = sparsepoly(e); |
---|
| 697 | |
---|
| 698 | //Check if f lies not in I. |
---|
| 699 | if(reduce(f,I)!=0) |
---|
| 700 | { |
---|
| 701 | J = quotient1(I,f); |
---|
| 702 | |
---|
| 703 | //If f is a zerodivisor modulo I... |
---|
| 704 | if( isSubset(J,I) == 0 ) |
---|
| 705 | { |
---|
| 706 | //...then use recursion... |
---|
| 707 | if(printlevel > 2){"We found a zerodivisor -- recursion";} |
---|
| 708 | l = decompEHV(J) + decompEHV(I+f); |
---|
| 709 | return(l); |
---|
| 710 | } |
---|
| 711 | if(printlevel > 2){"We found a non-zero-divisor.";} |
---|
| 712 | |
---|
| 713 | //...else compute the vectorspace dimension vd of I and... |
---|
| 714 | vd = vdim(I); |
---|
| 715 | |
---|
| 716 | //...compute m minimal such that 1,f,f^2,...,f^m are linearly dependent. |
---|
| 717 | qring Q = I; |
---|
| 718 | poly g = fetch(base,f); |
---|
| 719 | k = algDependent(g); |
---|
| 720 | def R = k[2]; |
---|
| 721 | setring R; |
---|
| 722 | if(size(ker)!=1) |
---|
| 723 | { |
---|
| 724 | //Calculate a generator for ker. |
---|
| 725 | ker = mstd(ker)[2]; |
---|
| 726 | } |
---|
| 727 | poly g = ker[1]; |
---|
| 728 | m = deg(g); |
---|
| 729 | |
---|
| 730 | //If m and vd coincide... |
---|
| 731 | if(m==vd) |
---|
| 732 | { |
---|
| 733 | if(printlevel > 2){"We have a good candidate.";} |
---|
| 734 | //...then factorize g. |
---|
| 735 | if(printlevel > 2){"Factorizing.";} |
---|
| 736 | L = factorize(g,2); //returns non-constant factors and multiplicities |
---|
| 737 | nfact = size(L[1]); |
---|
| 738 | |
---|
| 739 | //If g is irreducible... |
---|
| 740 | if(nfact==1 and L[2][1]==1) |
---|
| 741 | { |
---|
| 742 | if(printlevel > 2){"The element is irreducible.";} |
---|
| 743 | setring base; |
---|
| 744 | //..then I is a maximal ideal... |
---|
| 745 | l[1] = I; |
---|
| 746 | kill R,Q; |
---|
| 747 | return(l); |
---|
| 748 | } |
---|
| 749 | //...else... |
---|
| 750 | else |
---|
| 751 | { |
---|
| 752 | |
---|
| 753 | if(printlevel > 2){"The element is not irreducible -- recursion.";} |
---|
| 754 | //...take a non-trivial factor g1 of g ... |
---|
| 755 | poly g1 = L[1][1]; |
---|
| 756 | //..and insert f in g1... |
---|
| 757 | execute("ring newR = (" + charstr(R) + "),(y(1)),(" + ordstr(R) + ");"); |
---|
| 758 | poly g2 = imap(R,g1); |
---|
| 759 | setring base; |
---|
| 760 | h = fetch(newR,g2); |
---|
| 761 | h = subst(h,var(1),f); |
---|
| 762 | //...and use recursion... |
---|
| 763 | kill R,Q,newR; |
---|
| 764 | l = l + decompEHV(quotient1(I,h)); |
---|
| 765 | l = l + decompEHV(I+h); |
---|
| 766 | return(l); |
---|
| 767 | } |
---|
| 768 | } |
---|
| 769 | setring base; |
---|
| 770 | kill R,Q; |
---|
| 771 | } |
---|
| 772 | } |
---|
| 773 | } |
---|
| 774 | |
---|
| 775 | example |
---|
| 776 | { |
---|
| 777 | "EXAMPLE:"; |
---|
| 778 | echo = 2; |
---|
| 779 | ring r = 32003,(x,y),dp; |
---|
| 780 | ideal i = x2+y2-10,x2+xy+2y2-16; |
---|
| 781 | decompEHV(i); |
---|
| 782 | } |
---|
| 783 | |
---|
| 784 | |
---|
| 785 | ///////////////////////////////////////////////////////////////////// |
---|
| 786 | // // |
---|
| 787 | // A S S O C I A T E D P R I M E S // |
---|
| 788 | // // |
---|
| 789 | ///////////////////////////////////////////////////////////////////// |
---|
| 790 | |
---|
| 791 | |
---|
| 792 | ///////////////////////////////////////////////////////////////////// |
---|
| 793 | static proc idempotent(ideal I) |
---|
| 794 | "USAGE: idempotent(I); ideal I (weighted) homogeneous radical ideal, |
---|
| 795 | I intersected K[x(1),...,x(k)] zero-dimensional |
---|
| 796 | where deg(x(i))=0 for all i <= k and deg(x(i))>0 for all i>k. |
---|
| 797 | RETURN: a list of ideals I(1),...,I(t) such that |
---|
| 798 | K[x]/I = K[x]/I(1) x ... x K[x]/I(t)" |
---|
| 799 | { |
---|
| 800 | if(printlevel > 2){"Entering idempotent.";} |
---|
| 801 | if(ord_test(basering)!=1) |
---|
| 802 | { |
---|
| 803 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 804 | } |
---|
| 805 | int n = nvars(basering); |
---|
| 806 | poly f,g; |
---|
| 807 | string j; |
---|
| 808 | int i,k,splits; |
---|
| 809 | list l; |
---|
| 810 | |
---|
| 811 | //Collect all variables of degree 0. |
---|
| 812 | for(i=1; i<=n; i++) |
---|
| 813 | { |
---|
| 814 | if(deg(var(i)) > 0) |
---|
| 815 | { |
---|
| 816 | f = f*var(i); |
---|
| 817 | } |
---|
| 818 | else |
---|
| 819 | { |
---|
| 820 | splits = 1; |
---|
| 821 | j = j + string(var(i)) + "," ; |
---|
| 822 | } |
---|
| 823 | } |
---|
| 824 | //If there are no variables of degree 0 |
---|
| 825 | //then there are no idempotents and we are done... |
---|
| 826 | if(splits == 0) |
---|
| 827 | { |
---|
| 828 | l[1]=I; |
---|
| 829 | return(l); |
---|
| 830 | } |
---|
| 831 | |
---|
| 832 | //...else compute J = I intersected K[x(1),...,x(k)]... |
---|
| 833 | ideal J = eliminate(I,f); |
---|
| 834 | def base = basering; |
---|
| 835 | j = j[1,size(j)-1]; |
---|
| 836 | j = "ring @r = (" + charstr(basering) + "),(" + j + "),(" + ordstr(basering) + ");"; |
---|
| 837 | execute(j); |
---|
| 838 | ideal J = imap(base,J); |
---|
| 839 | |
---|
| 840 | //...and compute the associated primes of the zero-dimensional ideal J. |
---|
| 841 | list L = decompEHV(J); |
---|
| 842 | int s = size(L); |
---|
| 843 | ideal K,Z; |
---|
| 844 | poly g; |
---|
| 845 | //For each associated prime ideal P_i of J... |
---|
| 846 | for(i=1; i<=s; i++) |
---|
| 847 | { |
---|
| 848 | K = ideal(1); |
---|
| 849 | //...comnpute the intersection K of the other associated prime ideals... |
---|
| 850 | for(k=1; k<=s; k++) |
---|
| 851 | { |
---|
| 852 | if(i!=k) |
---|
| 853 | { |
---|
| 854 | K = intersect(K,L[k]); |
---|
| 855 | } |
---|
| 856 | } |
---|
| 857 | |
---|
| 858 | //...and find an element that lies in K but not in P_i... |
---|
| 859 | g = randomid(K,1)[1]; |
---|
| 860 | Z = L[i]; |
---|
| 861 | Z = std(Z); |
---|
| 862 | while(reduce(g,Z)==0) |
---|
| 863 | { |
---|
| 864 | g = randomid(K,1)[1]; |
---|
| 865 | } |
---|
| 866 | setring base; |
---|
| 867 | g = imap(@r,g); |
---|
| 868 | //...and compute the corresponding ideal I(i) |
---|
| 869 | l[i] = quotient(I,g); |
---|
| 870 | setring @r; |
---|
| 871 | |
---|
| 872 | } |
---|
| 873 | setring base; |
---|
| 874 | return(l); |
---|
| 875 | } |
---|
| 876 | |
---|
| 877 | |
---|
| 878 | ///////////////////////////////////////////////////////////////////// |
---|
| 879 | static proc equiAssEHV(ideal I) |
---|
| 880 | "USAGE: equiAssEHV(I); I equidimensional, radical, and homogeneous ideal |
---|
| 881 | RETURN: a list, the associated prime ideals of I" |
---|
| 882 | { |
---|
| 883 | if(printlevel > 2){"Entering equiAssEHV.";} |
---|
| 884 | if(ord_test(basering)!=1) |
---|
| 885 | { |
---|
| 886 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 887 | } |
---|
| 888 | list L; |
---|
| 889 | def base = basering; |
---|
| 890 | int n = nvars(basering); |
---|
| 891 | int i,j; |
---|
| 892 | |
---|
| 893 | //Compute the normalization of I. |
---|
| 894 | if(printlevel > 2){"Entering Normalization.";} |
---|
[6af4b1] | 895 | list norOut = normal(I, "noDeco"); |
---|
| 896 | list K = norOut[1]; |
---|
[e27e47] | 897 | if(printlevel > 2){"Leaving Normalisation.";} |
---|
| 898 | |
---|
| 899 | //The normalization algorithm returns k factors. |
---|
| 900 | int k = size(K); |
---|
| 901 | if(printlevel > 1){"Normalization algorithm splits ideal in " + string(k) + " factors.";} |
---|
| 902 | |
---|
| 903 | //Examine each factor of the normalization. |
---|
| 904 | def P; |
---|
| 905 | for(i=1; i<=k; i++) |
---|
| 906 | { |
---|
| 907 | P = K[i]; |
---|
| 908 | setring P; |
---|
| 909 | |
---|
| 910 | //Use procedure idempotent to split the i-th factor of |
---|
| 911 | //the normalization in a product of integral domains. |
---|
| 912 | if(printlevel > 1){"Examining " + string(i) +". factor.";} |
---|
| 913 | list l = idempotent(norid); |
---|
| 914 | int leng = size(l); |
---|
| 915 | if(printlevel > 1){"Idempotent algorithm splits factor " + string(i) + " in " + string(leng) + " factors.";} |
---|
| 916 | |
---|
| 917 | //Intersect the minimal primes corresponding to the |
---|
| 918 | //integral domains obtained from idempotent with the groundring, |
---|
| 919 | //i.e. compute the preimages w.r.t. the corresponding normalization map |
---|
| 920 | ideal J; |
---|
| 921 | for(j=1; j<=leng; j++) |
---|
| 922 | { |
---|
| 923 | J = l[j]; |
---|
| 924 | setring base; |
---|
| 925 | L[size(L)+1] = preimage(P,normap,J); |
---|
| 926 | setring P; |
---|
| 927 | } |
---|
| 928 | kill l, leng, J; |
---|
| 929 | setring base; |
---|
| 930 | } |
---|
| 931 | |
---|
| 932 | return(L); |
---|
| 933 | } |
---|
| 934 | |
---|
| 935 | |
---|
| 936 | ///////////////////////////////////////////////////////////////////// |
---|
| 937 | proc AssEHV(ideal I, list #) |
---|
| 938 | "USAGE: AssEHV(I [,Strategy]); I Ideal, Strategy list |
---|
| 939 | RETURN: a list, the associated prime ideals of I |
---|
| 940 | NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos. |
---|
| 941 | The (optional) second argument determines the strategy used: |
---|
| 942 | Strategy[1] > strategy for the equidimensional part |
---|
| 943 | = 0 : uses equiMaxEHV |
---|
| 944 | = 1 : uses equidimMax |
---|
| 945 | Strategy[2] > strategy for the equidimensional radical |
---|
| 946 | = 0 : uses equiRadEHV |
---|
| 947 | = 1 : uses equiRadical |
---|
| 948 | Strategy[3] > strategy for equiRadEHV |
---|
| 949 | = 0 : combination of strategy 1 and 2 |
---|
| 950 | = 1 : computation of the radical just with the |
---|
| 951 | help of regular sequences |
---|
| 952 | = 2 : does not try to find a regular sequence |
---|
| 953 | Strategy[4] > strategy for the computation of ideal quotients |
---|
| 954 | = n : uses quot(.,.,n) for the ideal quotient computations |
---|
| 955 | If no second argument is given, Strategy=(0,0,0,0) is used. |
---|
| 956 | EXAMPLE: example AssEHV; shows an example" |
---|
| 957 | { |
---|
| 958 | if(ord_test(basering)!=1) |
---|
| 959 | { |
---|
| 960 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 961 | } |
---|
| 962 | if(printlevel > 2){"Entering AssEHV";} |
---|
| 963 | |
---|
| 964 | //Specify the strategy to be used. |
---|
| 965 | if(size(#)==0) |
---|
| 966 | { |
---|
| 967 | # = 0,0,0,0; |
---|
| 968 | } |
---|
| 969 | if(size(#)==1) |
---|
| 970 | { |
---|
| 971 | # = #[1],0,0,0; |
---|
| 972 | } |
---|
| 973 | if(size(#)==2) |
---|
| 974 | { |
---|
| 975 | # = #[1],#[2],0,0; |
---|
| 976 | } |
---|
| 977 | if(size(#)==3) |
---|
| 978 | { |
---|
| 979 | # = #[1],#[2],#[3],0; |
---|
| 980 | } |
---|
| 981 | |
---|
| 982 | list L; |
---|
| 983 | ideal K; |
---|
| 984 | def base = basering; |
---|
| 985 | int m,j; |
---|
| 986 | ideal J = groebner(I); |
---|
| 987 | int n = nvars(basering); |
---|
| 988 | int d = dim(J); |
---|
| 989 | |
---|
| 990 | //Compute a resolution of I. |
---|
| 991 | if(printlevel > 2){"Computing resolution.";} |
---|
| 992 | if(homog(I)==1) |
---|
| 993 | { |
---|
| 994 | list re = sres(J,0); |
---|
| 995 | re = minres(re); |
---|
| 996 | } |
---|
| 997 | else |
---|
| 998 | { |
---|
| 999 | list re = mres(I,0); |
---|
| 1000 | } |
---|
| 1001 | ideal ann; |
---|
| 1002 | int cod; |
---|
| 1003 | |
---|
| 1004 | //For 0<=i<= dim(I) compute the intersection of the i-dimensional associated primes of I |
---|
| 1005 | for(int i=0; i<=d; i++) |
---|
| 1006 | { |
---|
| 1007 | if(printlevel > 1){"Are there components of dimension " + string(i) + "?";} |
---|
| 1008 | ann = AnnExtEHV(n-i,re); |
---|
| 1009 | cod = n - dim(groebner(ann)); |
---|
| 1010 | |
---|
| 1011 | //If there are associated primes of dimension i... |
---|
| 1012 | if(cod == n-i) |
---|
| 1013 | { |
---|
| 1014 | if(printlevel > 1){"Yes. There are components of dimension " + string(i) + ".";} |
---|
| 1015 | //...then compute the intersection K of all associated primes of I of dimension i |
---|
| 1016 | if(#[2]==0) |
---|
| 1017 | { |
---|
| 1018 | if(size(#) > 3) |
---|
| 1019 | { |
---|
| 1020 | K = equiRadEHV(ann,#[1],#[3],#[4]); |
---|
| 1021 | } |
---|
| 1022 | if(size(#) > 2) |
---|
| 1023 | { |
---|
| 1024 | K = equiRadEHV(ann,#[1],#[3]); |
---|
| 1025 | } |
---|
| 1026 | else |
---|
| 1027 | { |
---|
| 1028 | K = equiRadEHV(ann,#[1]); |
---|
| 1029 | } |
---|
| 1030 | } |
---|
| 1031 | if(#[2]==1) |
---|
| 1032 | { |
---|
| 1033 | K = equiRadical(ann); |
---|
| 1034 | } |
---|
| 1035 | attrib(K,"isEquidimensional",1); |
---|
| 1036 | attrib(K,"isRadical",1); |
---|
| 1037 | |
---|
| 1038 | //If K is already homogeneous then use equiAssEHV to recover the associated primes of K... |
---|
| 1039 | if(homog(K)==1) |
---|
| 1040 | { |
---|
| 1041 | if(printlevel > 2){"Input already homogeneous.";} |
---|
| 1042 | L = L + equiAssEHV(K); |
---|
| 1043 | } |
---|
| 1044 | |
---|
| 1045 | //...else... |
---|
| 1046 | else |
---|
| 1047 | { |
---|
| 1048 | //...homogenize K w.r.t. t,... |
---|
| 1049 | if(printlevel > 2){"Input not homogeneous; must homogenize.";} |
---|
| 1050 | changeord("homoR","dp"); |
---|
| 1051 | ideal homoJ = fetch(base,K); |
---|
| 1052 | homoJ = groebner(homoJ); |
---|
| 1053 | execute("ring newR = (" + charstr(base) + "),(x(1..n),t),dp;"); |
---|
| 1054 | ideal homoK = fetch(homoR,homoJ); |
---|
| 1055 | homoK = homog(homoK,t); |
---|
| 1056 | attrib(homoK,"isEquidimensional",1); |
---|
| 1057 | attrib(homoK,"isRadical",1); |
---|
| 1058 | |
---|
| 1059 | //...compute the associated primes of the homogenization using equiAssEHV,... |
---|
| 1060 | list l = equiAssEHV(homoK); |
---|
| 1061 | |
---|
| 1062 | //...and set t=1 in the generators of the associated primes just computed. |
---|
| 1063 | ideal Z; |
---|
| 1064 | for(j=1; j<=size(l); j++) |
---|
| 1065 | { |
---|
| 1066 | Z = subst(l[j],t,1); |
---|
| 1067 | setring base; |
---|
| 1068 | L[size(L)+1] = fetch(newR,Z); |
---|
| 1069 | setring newR; |
---|
| 1070 | } |
---|
| 1071 | setring base; |
---|
| 1072 | kill homoR; |
---|
| 1073 | kill newR; |
---|
| 1074 | } |
---|
| 1075 | } |
---|
| 1076 | else |
---|
| 1077 | { |
---|
| 1078 | if(printlevel > 1){"No. There are no components of dimension " + string(i) + ".";} |
---|
| 1079 | } |
---|
| 1080 | } |
---|
| 1081 | return(L); |
---|
| 1082 | } |
---|
| 1083 | |
---|
| 1084 | example |
---|
| 1085 | { |
---|
| 1086 | "EXAMPLE:"; |
---|
| 1087 | echo = 2; |
---|
| 1088 | ring r = 0,(x,y,z),dp; |
---|
| 1089 | poly p = z2+1; |
---|
| 1090 | poly q = z3+2; |
---|
| 1091 | ideal i = p*q^2,y-z2; |
---|
| 1092 | list pr = AssEHV(i); |
---|
| 1093 | pr; |
---|
| 1094 | } |
---|
| 1095 | |
---|
| 1096 | ///////////////////////////////////////////////////////////////////// |
---|
| 1097 | proc minAssEHV(ideal I, list #) |
---|
| 1098 | "USAGE: minAssEHV(I [,Strategy]); I ideal, Strategy list |
---|
| 1099 | RETURN: a list, the minimal associated prime ideals of I |
---|
| 1100 | NOTE: Uses the algorithm of Eisenbud/Huneke/Vasconcelos. |
---|
| 1101 | The (optional) second argument determines the strategy used: |
---|
| 1102 | Strategy[1] > strategy for the equidimensional part |
---|
| 1103 | = 0 : uses equiMaxEHV |
---|
| 1104 | = 1 : uses equidimMax |
---|
| 1105 | Strategy[2] > strategy for the equidimensional radical |
---|
| 1106 | = 0 : uses equiRadEHV, resp. radicalEHV |
---|
| 1107 | = 1 : uses equiRadical, resp. radical |
---|
| 1108 | Strategy[3] > strategy for equiRadEHV |
---|
| 1109 | = 0 : combination of strategy 1 and 2 |
---|
| 1110 | = 1 : computation of the radical just with the |
---|
| 1111 | help of regular sequences |
---|
| 1112 | = 2 : does not try to find a regular sequence |
---|
| 1113 | Strategy[4] > strategy for the computation of ideal quotients |
---|
| 1114 | = n : uses quot(.,.,n) for the ideal quotient computations |
---|
| 1115 | If no second argument is given, Strategy=(0,0,0,0) is used. |
---|
| 1116 | EXAMPLE: example minAssEHV; shows an example" |
---|
| 1117 | { |
---|
| 1118 | if(ord_test(basering)!=1) |
---|
| 1119 | {ERROR("// Not implemented for this ordering, please change to global ordering.");} |
---|
| 1120 | |
---|
| 1121 | //Specify the strategy to be used. |
---|
| 1122 | if(size(#)==0) |
---|
| 1123 | { |
---|
| 1124 | # = 0,0,0,0; |
---|
| 1125 | } |
---|
| 1126 | if(size(#)==1) |
---|
| 1127 | { |
---|
| 1128 | # = #[1],0,0,0; |
---|
| 1129 | } |
---|
| 1130 | if(size(#)==2) |
---|
| 1131 | { |
---|
| 1132 | # = #[1],#[2],0,0; |
---|
| 1133 | } |
---|
| 1134 | if(size(#)==3) |
---|
| 1135 | { |
---|
| 1136 | # = #[1],#[2],#[3],0; |
---|
| 1137 | } |
---|
| 1138 | |
---|
| 1139 | //Compute the radical of I. |
---|
| 1140 | if(#[2]==0) |
---|
| 1141 | { |
---|
| 1142 | I = radEHV(I,#); |
---|
| 1143 | } |
---|
| 1144 | if(#[2]==1) |
---|
| 1145 | { |
---|
| 1146 | I = radical(I); |
---|
| 1147 | } |
---|
| 1148 | |
---|
| 1149 | //Compute the associated primes of the radical. |
---|
| 1150 | return(AssEHV(I,#)); |
---|
| 1151 | } |
---|
| 1152 | |
---|
| 1153 | example |
---|
| 1154 | { |
---|
| 1155 | "EXAMPLE:"; |
---|
| 1156 | echo = 2; |
---|
| 1157 | ring r = 0,(x,y,z),dp; |
---|
| 1158 | poly p = z2+1; |
---|
| 1159 | poly q = z3+2; |
---|
| 1160 | ideal i = p*q^2,y-z2; |
---|
| 1161 | list pr = minAssEHV(i); |
---|
| 1162 | pr; |
---|
| 1163 | } |
---|
| 1164 | |
---|
| 1165 | |
---|
| 1166 | ///////////////////////////////////////////////////////////////////// |
---|
| 1167 | // // |
---|
| 1168 | // P R I M A R Y D E C O M P O S I T I O N // |
---|
| 1169 | // // |
---|
| 1170 | ///////////////////////////////////////////////////////////////////// |
---|
| 1171 | |
---|
| 1172 | |
---|
| 1173 | ///////////////////////////////////////////////////////////////////// |
---|
| 1174 | proc localize(ideal I, ideal P, list l) |
---|
| 1175 | "USAGE: localize(I,P,l); I ideal, P an associated prime ideal of I, |
---|
| 1176 | l list of all associated primes of I |
---|
| 1177 | RETURN: ideal, the contraction of the ideal generated by I |
---|
| 1178 | in the localization w.r.t P |
---|
| 1179 | EXAMPLE: example localize; shows an example" |
---|
| 1180 | { |
---|
| 1181 | if(ord_test(basering)!=1) |
---|
| 1182 | { |
---|
| 1183 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 1184 | } |
---|
| 1185 | |
---|
| 1186 | ideal Intersection = ideal(1); |
---|
| 1187 | int s = size(l); |
---|
| 1188 | if(attrib(P,"isSB")!=1) |
---|
| 1189 | { |
---|
| 1190 | P = groebner(P); |
---|
| 1191 | } |
---|
| 1192 | |
---|
| 1193 | //Compute the intersection of all associated primes of I that are not contained in P. |
---|
| 1194 | for(int i=1; i<=s; i++) |
---|
| 1195 | { |
---|
| 1196 | if(isSubset(l[i],P)!=1) |
---|
| 1197 | { |
---|
| 1198 | Intersection = intersect(Intersection,l[i]); |
---|
| 1199 | } |
---|
| 1200 | } |
---|
| 1201 | Intersection = groebner(Intersection); |
---|
| 1202 | |
---|
| 1203 | //If the intersection is the entire ring... |
---|
| 1204 | if(reduce(1,Intersection)==0) |
---|
| 1205 | { |
---|
| 1206 | //...then return I... |
---|
| 1207 | return(I); |
---|
| 1208 | } |
---|
| 1209 | //...else try to find an element f that lies in the intersection but outside P... |
---|
| 1210 | poly f = 0; |
---|
| 1211 | while(reduce(f,P) == 0) |
---|
| 1212 | { |
---|
| 1213 | f = randomid(Intersection,1)[1]; |
---|
| 1214 | } |
---|
| 1215 | |
---|
| 1216 | //...and saturate I w.r.t. f. |
---|
| 1217 | I = sat(I,f)[1]; |
---|
| 1218 | return(I); |
---|
| 1219 | } |
---|
| 1220 | |
---|
| 1221 | ///////////////////////////////////////////////////////////////////// |
---|
| 1222 | proc componentEHV(ideal I, ideal P, list L, list #) |
---|
| 1223 | "USAGE: componentEHV(I,P,L [,Strategy]); |
---|
| 1224 | I ideal, P associated prime of I, |
---|
| 1225 | L list of all associated primes of I, Strategy list |
---|
| 1226 | RETURN: ideal, a P-primary component Q for I |
---|
| 1227 | NOTE: The (optional) second argument determines the strategy used: |
---|
| 1228 | Strategy[1] > strategy for equidimensional part |
---|
| 1229 | = 0 : uses equiMaxEHV |
---|
| 1230 | = 1 : uses equidimMax |
---|
| 1231 | If no second argument is given then Strategy=0 is used. |
---|
| 1232 | EXAMPLE: example componentEHV; shows an example" |
---|
| 1233 | { |
---|
| 1234 | if(printlevel > 2){"Entering componentEHV.";} |
---|
| 1235 | if(ord_test(basering)!=1) |
---|
| 1236 | { |
---|
| 1237 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 1238 | } |
---|
| 1239 | |
---|
| 1240 | //If no strategy is specified use standard strategy. |
---|
| 1241 | if(size(#)==0) |
---|
| 1242 | { |
---|
| 1243 | # = 0; |
---|
| 1244 | } |
---|
| 1245 | |
---|
| 1246 | ideal T = P; |
---|
| 1247 | ideal Q; |
---|
| 1248 | |
---|
| 1249 | //Compute the localization of I at P... |
---|
| 1250 | ideal IP = groebner(localize(I,P,L)); |
---|
| 1251 | |
---|
| 1252 | //...and compute the saturation of the localization w.r.t. P. |
---|
| 1253 | ideal IP2 = sat(IP,P)[1]; |
---|
| 1254 | |
---|
| 1255 | //As long as we have not found a primary component... |
---|
| 1256 | int isPrimaryComponent = 0; |
---|
| 1257 | while(isPrimaryComponent!=1) |
---|
| 1258 | { |
---|
| 1259 | //...compute the equidimensional part Q of I+P^n... |
---|
| 1260 | if(#[1]==0) |
---|
| 1261 | { |
---|
| 1262 | Q = equiMaxEHV(I+T); |
---|
| 1263 | } |
---|
| 1264 | if(#[1]==1) |
---|
| 1265 | { |
---|
| 1266 | Q = equidimMax(I+T); |
---|
| 1267 | } |
---|
| 1268 | //...and check if it is a primary component for P. |
---|
| 1269 | if(isSubset(intersect(IP2,Q),IP)==1) |
---|
| 1270 | { |
---|
| 1271 | isPrimaryComponent = 1; |
---|
| 1272 | } |
---|
| 1273 | else |
---|
| 1274 | { |
---|
| 1275 | T = T*P; |
---|
| 1276 | } |
---|
| 1277 | } |
---|
| 1278 | if(printlevel > 2){"Leaving componentEHV.";} |
---|
| 1279 | return(Q); |
---|
| 1280 | } |
---|
| 1281 | |
---|
| 1282 | example |
---|
| 1283 | { |
---|
| 1284 | "EXAMPLE:"; |
---|
| 1285 | echo = 2; |
---|
| 1286 | ring r = 0,(x,y,z),dp; |
---|
| 1287 | poly p = z2+1; |
---|
| 1288 | poly q = z3+2; |
---|
| 1289 | ideal i = p*q^2,y-z2; |
---|
| 1290 | list pr = AssEHV(i); |
---|
| 1291 | componentEHV(i,pr[1],pr); |
---|
| 1292 | } |
---|
| 1293 | |
---|
| 1294 | |
---|
| 1295 | ///////////////////////////////////////////////////////////////////// |
---|
| 1296 | proc primdecEHV(ideal I, list #) |
---|
| 1297 | "USAGE: primdecEHV(I [,Strategy]); I ideal, Strategy list |
---|
| 1298 | RETURN: a list pr of primary ideals and their associated primes: |
---|
| 1299 | pr[i][1] the i-th primary component, |
---|
| 1300 | pr[i][2] the i-th prime component. |
---|
| 1301 | NOTE: Algorithm of Eisenbud/Huneke/Vasconcelos. |
---|
| 1302 | The (optional) second argument determines the strategy used: |
---|
| 1303 | Strategy[1] > strategy for equidimensional part |
---|
| 1304 | = 0 : uses equiMaxEHV |
---|
| 1305 | = 1 : uses equidimMax |
---|
| 1306 | Strategy[2] > strategy for equidimensional radical |
---|
| 1307 | = 0 : uses equiRadEHV, resp. radicalEHV |
---|
| 1308 | = 1 : uses equiRadical, resp. radical |
---|
| 1309 | Strategy[3] > strategy for equiRadEHV |
---|
| 1310 | = 0 : combination of strategy 1 and 2 |
---|
| 1311 | = 1 : computation of the radical just with the |
---|
| 1312 | help of regular sequences |
---|
| 1313 | = 2 : does not try to find a regular sequence |
---|
| 1314 | Strategy[4] > strategy for the computation of ideal quotients |
---|
| 1315 | = n : uses quot(.,.,n) for the ideal quotient computations |
---|
| 1316 | If no second argument is given then Strategy=(0,0,0,0) is used. |
---|
| 1317 | EXAMPLE: example primdecEHV; shows an example" |
---|
| 1318 | { |
---|
| 1319 | if(printlevel > 2){"Entering primdecEHV.";} |
---|
| 1320 | if(ord_test(basering)!=1) |
---|
| 1321 | { |
---|
| 1322 | ERROR("// Not implemented for this ordering, please change to global ordering."); |
---|
| 1323 | } |
---|
| 1324 | list L,K; |
---|
| 1325 | |
---|
| 1326 | //Specify the strategy to be used. |
---|
| 1327 | if(size(#)==0) |
---|
| 1328 | { |
---|
| 1329 | # = 0,0,0,0; |
---|
| 1330 | } |
---|
| 1331 | if(size(#)==1) |
---|
| 1332 | { |
---|
| 1333 | # = #[1],0,0,0; |
---|
| 1334 | } |
---|
| 1335 | if(size(#)==2) |
---|
| 1336 | { |
---|
| 1337 | # = #[1],#[2],0,0; |
---|
| 1338 | } |
---|
| 1339 | if(size(#)==3) |
---|
| 1340 | { |
---|
| 1341 | # = #[1],#[2],#[3],0; |
---|
| 1342 | } |
---|
| 1343 | |
---|
| 1344 | //Compute the associated primes of I... |
---|
| 1345 | L = AssEHV(I,#); |
---|
| 1346 | if(printlevel > 0){"We have " + string(size(L)) + " prime components.";} |
---|
| 1347 | |
---|
| 1348 | //...and compute for each associated prime of I a corresponding primary component. |
---|
| 1349 | int l = size(L); |
---|
| 1350 | for(int i=1; i<=l; i++) |
---|
| 1351 | { |
---|
| 1352 | K[i] = list(); |
---|
| 1353 | K[i][2] = L[i]; |
---|
| 1354 | K[i][1] = componentEHV(I,L[i],L,#[1]); |
---|
| 1355 | } |
---|
| 1356 | if(printlevel > 2){"Leaving primdecEHV.";} |
---|
| 1357 | return(K); |
---|
| 1358 | |
---|
| 1359 | } |
---|
| 1360 | |
---|
| 1361 | example |
---|
| 1362 | { |
---|
| 1363 | "EXAMPLE:"; |
---|
| 1364 | echo = 2; |
---|
| 1365 | ring r = 0,(x,y,z),dp; |
---|
| 1366 | poly p = z2+1; |
---|
| 1367 | poly q = z3+2; |
---|
| 1368 | ideal i = p*q^2,y-z2; |
---|
| 1369 | list pr = primdecEHV(i); |
---|
| 1370 | pr; |
---|
| 1371 | } |
---|
| 1372 | |
---|
| 1373 | |
---|
| 1374 | ///////////////////////////////////////////////////////////////////// |
---|
| 1375 | // // |
---|
| 1376 | // A L G O R I T H M S F O R T E S T I N G // |
---|
| 1377 | // // |
---|
| 1378 | ///////////////////////////////////////////////////////////////////// |
---|
| 1379 | |
---|
| 1380 | |
---|
| 1381 | ///////////////////////////////////////////////////////////////////// |
---|
| 1382 | proc compareLists(list L, list K) |
---|
| 1383 | "USAGE: checkLists(L,K); L,K list of ideals |
---|
| 1384 | RETURN: integer, 1 if the lists are the same up to ordering and 0 otherwise |
---|
| 1385 | EXAMPLE: example checkLists; shows an example" |
---|
| 1386 | { |
---|
| 1387 | int s1 = size(L); |
---|
| 1388 | int s2 = size(K); |
---|
| 1389 | if(s1!=s2) |
---|
| 1390 | { |
---|
| 1391 | return(0); |
---|
| 1392 | } |
---|
| 1393 | list L1, K1; |
---|
| 1394 | int i,j,t; |
---|
| 1395 | list N; |
---|
| 1396 | for(i=1; i<=s1; i++) |
---|
| 1397 | { |
---|
| 1398 | L1[i]=std(L[i][2]); |
---|
| 1399 | K1[i]=std(K[i][2]); |
---|
| 1400 | } |
---|
| 1401 | for(i=1; i<=s1; i++) |
---|
| 1402 | { |
---|
| 1403 | for(j=1; j<=s1; j++) |
---|
| 1404 | { |
---|
| 1405 | if(isSubset(L1[i],K1[j])) |
---|
| 1406 | { |
---|
| 1407 | if(isSubset(K1[j],L1[i])) |
---|
| 1408 | { |
---|
| 1409 | for(t=1; t<=size(N); t++) |
---|
| 1410 | { |
---|
| 1411 | if(N[t]==j) |
---|
| 1412 | { |
---|
| 1413 | return(0); |
---|
| 1414 | } |
---|
| 1415 | } |
---|
| 1416 | N[size(N)+1]=j; |
---|
| 1417 | } |
---|
| 1418 | |
---|
| 1419 | } |
---|
| 1420 | } |
---|
| 1421 | } |
---|
| 1422 | return(1); |
---|
| 1423 | } |
---|
| 1424 | |
---|
| 1425 | example |
---|
| 1426 | { |
---|
| 1427 | "EXAMPLE:"; |
---|
| 1428 | echo = 2; |
---|
| 1429 | ring r = 0,(x,y),dp; |
---|
| 1430 | ideal i = x2,xy; |
---|
| 1431 | list L1 = primdecGTZ(i); |
---|
| 1432 | list L2 = primdecEHV(i); |
---|
| 1433 | compareLists(L1,L2); |
---|
| 1434 | } |
---|
| 1435 | |
---|