Changeset e801fe in git
- Timestamp:
- Feb 16, 1998, 1:07:04 PM (26 years ago)
- Branches:
- (u'spielwiese', '2fa36c576e6a4ddbb1093b43c7f8e9835e17e52a')
- Children:
- 1e7d8d9fbdaf1094cc869ca3a6be98b9e1fbfc9a
- Parents:
- 437e2c21158af6efe6402cc1da62ee40760c03e0
- Location:
- Singular/LIB
- Files:
-
- 3 edited
-
invar.lib (modified) (13 diffs)
-
normal.lib (modified) (55 diffs)
-
primdec.lib (modified) (108 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/invar.lib
r437e2c2 re801fe 1 // $Id: invar.lib,v 1. 3 1997-09-18 09:58:24 Singular Exp $1 // $Id: invar.lib,v 1.4 1998-02-16 12:07:01 pfister Exp $ 2 2 /////////////////////////////////////////////////////// 3 3 // invar.lib … … 7 7 ////////////////////////////////////////////////////// 8 8 9 LIBRARY: invar.lib PROCEDURE FOR COMPUTING INVARIANTS UNDER(C,+)-ACTIONS9 LIBRARY: invar.lib PROCEDURES FOR COMPUTING INVARIANTS OF (C,+)-ACTIONS 10 10 11 11 invariantRing(matrix m,poly p,poly q,int choose) … … 142 142 143 143 144 proc reduction(poly p, ideal mo, list #)145 USAGE: reduction(p, mo,q); p poly, moideal, q (optional) monomial144 proc reduction(poly p, ideal dom, list #) 145 USAGE: reduction(p,dom,q); p poly, dom ideal, q (optional) monomial 146 146 RETURN: poly= (p-H(f1,...,fr))/q^a, if Lt(p)=H(Lt(f1),...,Lt(fr)) for 147 147 some polynomial H in r variables over the base field, 148 148 a maximal such that q^a devides p-H(f1,...,fr), 149 mo=(f1,...,fr)149 dom =(f1,...,fr) 150 150 NOTE: 151 151 EXAMPLE: example reduction; shows an example 152 152 { 153 153 int i; 154 int z=size( mo);154 int z=size(dom); 155 155 def bsr=basering; 156 156 … … 165 165 execute "ring s="+charstr(basering)+",("+varstr(basering)+",@(0..z)),dp;"; 166 166 167 //costructes the ideal mo=(p-@(0),mo[1]-@(1),...,mo[z]-@(z))168 ideal mo=imap(bsr,mo);167 //costructes the ideal dom=(p-@(0),dom[1]-@(1),...,dom[z]-@(z)) 168 ideal dom=imap(bsr,dom); 169 169 for (i=1;i<=z;i++) 170 170 { 171 mo[i]=lead(mo[i])-var(nvars(bsr)+i+1);172 } 173 mo=lead(imap(bsr,p))-@(0),mo;171 dom[i]=lead(dom[i])-var(nvars(bsr)+i+1); 172 } 173 dom=lead(imap(bsr,p))-@(0),dom; 174 174 175 175 //eliminates the variables of the basering bsr 176 //i.e. computes mointersected with K[@(0),...,@(z)]177 ideal kern=eliminate( mo,imap(bsr,v));176 //i.e. computes dom intersected with K[@(0),...,@(z)] 177 ideal kern=eliminate(dom,imap(bsr,v)); 178 178 179 179 // test wether @(0)-h(@(1),...,@(z)) is in ker for some poly h … … 186 186 { 187 187 h=(1/h)*kern[i]; 188 // defines the map psi : s ---> bsr defined by @(i) ---> p, mo[i]188 // defines the map psi : s ---> bsr defined by @(i) ---> p,dom[i] 189 189 setring bsr; 190 map psi=s,maxideal(1),p, mo;190 map psi=s,maxideal(1),p,dom; 191 191 poly re=psi(h); 192 192 … … 210 210 ring q=0,(x,y,z,u,v,w),dp; 211 211 poly p=x2yz-x2v; 212 ideal mo=x-w,u2w+1,yz-v;213 reduction(p, mo);214 reduction(p, mo,w);215 } 216 217 /////////////////////////////////////////////////////////////////////////////// 218 219 proc completeReduction(poly p, ideal mo, list #)220 USAGE: completeReduction(p, mo,q); p poly, moideal,212 ideal dom =x-w,u2w+1,yz-v; 213 reduction(p,dom); 214 reduction(p,dom,w); 215 } 216 217 /////////////////////////////////////////////////////////////////////////////// 218 219 proc completeReduction(poly p, ideal dom, list #) 220 USAGE: completeReduction(p,dom,q); p poly, dom ideal, 221 221 q (optional) monomial 222 RETURN: poly= the polynomial p reduced with movia the procedure222 RETURN: poly= the polynomial p reduced with dom via the procedure 223 223 reduction as long as possible 224 224 NOTE: … … 226 226 { 227 227 poly p1=p; 228 poly p2=reduction(p, mo,#);228 poly p2=reduction(p,dom,#); 229 229 while (p1!=p2) 230 230 { 231 231 p1=p2; 232 p2=reduction(p1, mo,#);232 p2=reduction(p1,dom,#); 233 233 } 234 234 return(p2); … … 238 238 ring q=0,(x,y,z,u,v,w),dp; 239 239 poly p=x2yz-x2v; 240 ideal mo=x-w,u2w+1,yz-v;241 completeReduction(p, mo);242 completeReduction(p, mo,w);243 } 244 /////////////////////////////////////////////////////////////////////////////// 245 246 proc inSubring(poly p, ideal mo)247 USAGE: inSubring(p, mo); p poly, moideal248 RETURN: list= 1,string(@(0)-h(@(1),...,@(size( mo)))) :if p = h(mo[1],...,mo[size(mo)])249 0,string(h(@(0),...,@(size( mo)))) :if there is only a nonlinear relation250 h(p, mo[1],...,mo[size(mo)])=0.240 ideal dom =x-w,u2w+1,yz-v; 241 completeReduction(p,dom); 242 completeReduction(p,dom,w); 243 } 244 /////////////////////////////////////////////////////////////////////////////// 245 246 proc inSubring(poly p, ideal dom) 247 USAGE: inSubring(p,dom); p poly, dom ideal 248 RETURN: list= 1,string(@(0)-h(@(1),...,@(size(dom)))) :if p = h(dom[1],...,dom[size(dom)]) 249 0,string(h(@(0),...,@(size(dom)))) :if there is only a nonlinear relation 250 h(p,dom[1],...,dom[size(dom)])=0. 251 251 NOTE: 252 252 EXAMPLE: example inSubring; shows an example 253 253 { 254 int z=size( mo);254 int z=size(dom); 255 255 int i; 256 256 def gnir=basering; … … 262 262 } 263 263 string eli=string(mile); 264 // the intersection of ideal nett=(p-@(0), mo[1]-@(1),...)264 // the intersection of ideal nett=(p-@(0),dom[1]-@(1),...) 265 265 // with the ring k[@(0),...,@(n)] is computed, the result is ker 266 266 execute "ring r1="+charstr(basering)+",("+varstr(basering)+",@(0..z)),dp;"; 267 ideal nett=imap(gnir, mo);267 ideal nett=imap(gnir,dom); 268 268 poly p; 269 269 for (i=1;i<=z;i++) … … 299 299 ring q=0,(x,y,z,u,v,w),dp; 300 300 poly p=xyzu2w-1yzu2w2+u4w2-1xu2vw+u2vw2+xyz-1yzw+2u2w-1xv+vw+2; 301 ideal mo=x-w,u2w+1,yz-v;302 inSubring(p, mo);301 ideal dom =x-w,u2w+1,yz-v; 302 inSubring(p,dom); 303 303 } 304 304 … … 478 478 " "; 479 479 karl; 480 pause;480 // pause; 481 481 " "; 482 482 } … … 517 517 " "; 518 518 karl; 519 pause;519 // pause; 520 520 " "; 521 521 } … … 643 643 644 644 } 645 /////////////////////////////////////////////////////////////////////////////// -
Singular/LIB/normal.lib
r437e2c2 re801fe 1 // $Id: normal.lib,v 1.1 1997-11-05 18:16:55 Singular Exp $2 1 /////////////////////////////////////////////////////////////////////////////// 3 2 // normal.lib … … 10 9 11 10 normal(ideal I) 12 // computes a set of rings such that their product is the 13 // normalization of the reduced basering/I 14 15 LIB "sing.lib"; 11 // computes a set of rings such that their product is the 12 // normalization of the reduced basering/I 13 14 LIB "sing.lib"; 16 15 LIB "primdec.lib"; 17 16 LIB "elim.lib"; … … 27 26 p = nonzero divisor of R 28 27 RETURN: 1 if R = R:J, 0 if not 29 EXAMPLE: example isR_HomJR; shows an example 28 EXAMPLE: example isR_HomJR; shows an example 30 29 { 31 30 int n, ii; … … 43 42 n=1; 44 43 for (ii=1; ii<=size(f); ii++ ) 45 { 46 if ( reduce(f[ii],lp) != 0) 44 { 45 if ( reduce(f[ii],lp) != 0) 47 46 { n = 0; break; } 48 47 } … … 73 72 NOTE: This is useful when enlarging P but keeping the weights of the old 74 73 variables 75 EXAMPLE: example extraweight; shows an example 74 EXAMPLE: example extraweight; shows an example 76 75 { 77 76 int ii,q,fi,fo,fia; … … 88 87 os = osP[fi+1,find(osP,")",fi)-fi-1]; 89 88 if( find(os,",") ) 90 { 89 { 91 90 execute "nw = "+os+";"; 92 91 if( size(nw) > ii ) … … 101 100 { 102 101 execute "q = "+os+";"; 103 if( q > ii ) 104 { 102 if( q > ii ) 103 { 105 104 nw = 0; nw[q-ii] = 0; 106 105 nw = nw + 1; //creates an intvec 1,...,1 of length q-ii … … 114 113 } 115 114 //-------------- adjust weight vector to length = nvars(P) ------------------- 116 if( fo > 1 ) 115 if( fo > 1 ) 117 116 { // case when weights were found 118 rw = rw[2..size(rw)]; 119 if( size(rw) > nvars(P) ) 120 { 121 rw = rw[1..nvars(P)]; 122 } 123 if( size(rw) < nvars(P) ) 124 { 125 nw=0; nw[nvars(P)-size(rw)]=0; nw=nw+1; rw=rw,nw; 126 } 127 } 128 else 117 rw = rw[2..size(rw)]; 118 if( size(rw) > nvars(P) ) 119 { 120 rw = rw[1..nvars(P)]; 121 } 122 if( size(rw) < nvars(P) ) 123 { 124 nw=0; nw[nvars(P)-size(rw)]=0; nw=nw+1; rw=rw,nw; 125 } 126 } 127 else 129 128 { // case when no weights were found 130 rw[nvars(P)]= 0; rw=rw+1; 129 rw[nvars(P)]= 0; rw=rw+1; 131 130 } 132 131 return(rw); … … 158 157 R is the quotient ring of P modulo the standard basis SBid 159 158 RETURN: a list of two objects 160 _[1]: a polynomial ring, containing two ideals, 'endid' and 'endphi' 159 _[1]: a polynomial ring, containing two ideals, 'endid' and 'endphi' 161 160 s.t. _[1]/endid = Hom_R(J,J) and 162 161 endphi describes the canonical map R -> Hom_R(J,J) 163 162 _[2]: an integer which is 1 if phi is an isomorphism, 0 if not 164 EXAMPLE: example HomJJ; shows an example 163 EXAMPLE: example HomJJ; shows an example 165 164 { 166 165 //---------- initialisation --------------------------------------------------- … … 180 179 181 180 //---- set attributes for special cases where algorithm can be simplified ----- 182 if( homo==1 ) 183 { 181 if( homo==1 ) 182 { 184 183 rw = extraweight(P); 185 184 } … … 213 212 214 213 //---------- computation of p*Hom(J,J) as R-ideal ----------------------------- 215 f = quotient(p*J,J); 216 if(y==1) 217 { 214 f = quotient(p*J,J); 215 if(y==1) 216 { 218 217 "the module Hom(rad(J),rad(J)) presented by the values on"; 219 218 "the non-zerodivisor"; … … 238 237 rf = interred(reduce(f,f2)); // represents p*Hom(J,J)/p*R = Hom(J,J)/R 239 238 if ( size(rf) == 0 ) 240 { 239 { 241 240 if ( homog(f) && find(ordstr(basering),"s")==0 ) 242 241 { 243 ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp); 242 ring newR1 = char(P),(X(1..nvars(P))),(a(rw),dp); 244 243 } 245 244 else 246 245 { 247 ring newR1 = char(P),(X(1..nvars(P))),dp; 246 ring newR1 = char(P),(X(1..nvars(P))),dp; 248 247 } 249 248 ideal endphi = maxideal(1); 250 ideal endid = fetch(P,id); 249 ideal endid = fetch(P,id); 251 250 attrib(endid,"isCohenMacaulay",isCo); 252 251 attrib(endid,"isPrim",isPr); … … 301 300 q = size(f); 302 301 syzf = syz(f); 303 302 304 303 if ( homo==1 ) 305 304 { … … 307 306 for ( ii=2; ii<=q; ii++ ) 308 307 { 309 rw = rw, deg(f[ii])-deg(f[1]); 310 rw1 = rw1, deg(f[ii])-deg(f[1]); 311 } 312 ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp); 308 rw = rw, deg(f[ii])-deg(f[1]); 309 rw1 = rw1, deg(f[ii])-deg(f[1]); 310 } 311 ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),(a(rw1),dp); 313 312 } 314 313 else 315 314 { 316 ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp; 315 ring newR1 = char(R),(X(1..nvars(R)),T(1..q)),dp; 317 316 } 318 317 … … 320 319 ideal SBid = psi1(SBid); 321 320 attrib(SBid,"isSB",1); 322 321 323 322 qring newR = std(SBid); 324 323 map psi = R,ideal(X(1..nvars(R))); … … 331 330 332 331 //---------- computation of Hom(J,J) as ring ---------------------------------- 333 // determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1), 332 // determine kernel of: R[T1,...,Tq] -> J:J >-> R[1/p]=R[t]/(t*p-1), 334 333 // Ti -> fi/p -> t*fi (p=f1=f[1]), to get ring structure. This is of course 335 334 // the same as the kernel of R[T1,...,Tq] -> pJ:J >-> R, Ti -> fi. … … 339 338 pf = f[1]*f; 340 339 T = matrix(ideal(T(1..q)),1,q); 341 Lin = ideal(T*syzf); 340 Lin = ideal(T*syzf); 342 341 if(y==1) 343 342 { … … 374 373 else 375 374 { 376 ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp; 375 ring newRing = char(R),(X(1..nvars(R)),T(2..q)),dp; 377 376 } 378 377 … … 381 380 382 381 map phi = basering,maxideal(1); 383 list Le = elimpart(endid); 382 list Le = elimpart(endid); 384 383 //this proc and the next loop try to 385 384 q = size(Le[2]); //substitute as many variables as possible 386 rw1 = 0; 385 rw1 = 0; 387 386 rw1[nvars(basering)] = 0; 388 387 rw1 = rw1+1; … … 395 394 kill ps; 396 395 397 for( ii=1; ii<=size(rw1); ii++ ) 398 { 399 if( Le[4][ii]==0 ) 400 { 396 for( ii=1; ii<=size(rw1); ii++ ) 397 { 398 if( Le[4][ii]==0 ) 399 { 401 400 rw1[ii]=0; //look for substituted vars 402 401 } … … 404 403 Le=elimpart(endid); 405 404 q = q + size(Le[2]); 406 } 405 } 407 406 endphi = phi(endphi); 408 407 … … 412 411 413 412 if (homo==1 && nvars(newRing)-q >1 && size(endid) >0 ) 414 { 413 { 415 414 jj=1; 416 415 for( ii=2; ii<size(rw1); ii++) 417 { 416 { 418 417 jj++; 419 if( rw1[ii]==0 ) 420 { 418 if( rw1[ii]==0 ) 419 { 421 420 rw=rw[1..jj-1],rw[jj+1..size(rw)]; 422 jj=jj-1; 421 jj=jj-1; 423 422 } 424 423 } … … 430 429 else 431 430 { 432 ring lastRing = char(R),(T(1..nvars(newRing)-q)),dp; 433 } 434 435 ideal lastmap; 431 ring lastRing = char(R),(T(1..nvars(newRing)-q)),dp; 432 } 433 434 ideal lastmap; 436 435 q = 1; 437 436 for(ii=1; ii<=size(rw1); ii++ ) … … 439 438 if ( rw1[ii]==1 ) { lastmap[ii] = T(q); q=q+1; } 440 439 if ( rw1[ii]==0 ) { lastmap[ii] = 0; } 441 } 440 } 442 441 map phi = newRing,lastmap; 443 442 ideal endid = phi(endid); … … 490 489 list L = HomJJ(Li); 491 490 def end = L[1]; // defines ring L[1], containing ideals endid and endphi 492 setring end; // makes end the basering 491 setring end; // makes end the basering 493 492 end; 494 493 endid; // end/endid is isomorphic to End(r/id) as ring … … 496 495 psi; 497 496 498 ring r = 32003,(x,y,z),dp; 497 ring r = 32003,(x,y,z),dp; 499 498 ideal id = x2-xy-xz+yz; 500 499 ideal J =y-z,x-z; … … 503 502 list L = HomJJ(Li,0); 504 503 def end = L[1]; // defines ring L[1], containing ideals endid and endphi 505 setring end; // makes end the basering 504 setring end; // makes end the basering 506 505 end; 507 506 endid; // end/endid is isomorphic to End(r/id) as ring … … 528 527 if( typeof(attrib(id,"isEquidimensional"))=="int" ) 529 528 { 530 if(attrib(id,"isEquidimensional")==1) 529 if(attrib(id,"isEquidimensional")==1) 531 530 { 532 attrib(prim[1],"isEquidimensional",1); 531 attrib(prim[1],"isEquidimensional",1); 533 532 } 534 533 } … … 539 538 if( typeof(attrib(id,"isCompleteIntersection"))=="int" ) 540 539 { 541 if(attrib(id,"isCompleteIntersection")==1) 540 if(attrib(id,"isCompleteIntersection")==1) 542 541 { 543 attrib(prim[1],"isCompleteIntersection",1); 542 attrib(prim[1],"isCompleteIntersection",1); 544 543 } 545 544 } … … 548 547 attrib(prim[1],"isCompleteIntersection",0); 549 548 } 550 549 551 550 if( typeof(attrib(id,"isPrim"))=="int" ) 552 551 { … … 559 558 if( typeof(attrib(id,"isIsolatedSingularity"))=="int" ) 560 559 { 561 if(attrib(id,"isIsolatedSingularity")==1) 560 if(attrib(id,"isIsolatedSingularity")==1) 562 561 {attrib(prim[1],"isIsolatedSingularity",1); } 563 562 } 564 563 else 565 564 { 566 attrib(prim[1],"isIsolatedSingularity",0); 565 attrib(prim[1],"isIsolatedSingularity",0); 567 566 } 568 567 if( typeof(attrib(id,"isCohenMacaulay"))=="int" ) 569 568 { 570 if(attrib(id,"isCohenMacaulay")==1) 569 if(attrib(id,"isCohenMacaulay")==1) 571 570 { attrib(prim[1],"isCohenMacaulay",1); } 572 571 } 573 572 else 574 573 { 575 attrib(prim[1],"isCohenMacaulay",0); 574 attrib(prim[1],"isCohenMacaulay",0); 576 575 } 577 576 if( typeof(attrib(id,"isRegInCodim2"))=="int" ) … … 582 581 else 583 582 { 584 attrib(prim[1],"isRegInCodim2",0); 583 attrib(prim[1],"isRegInCodim2",0); 585 584 } 586 585 return(normalizationPrimes(prim[1],maxideal(1))); … … 616 615 if( typeof(attrib(id,"isEquidimensional"))=="int" ) 617 616 { 618 if(attrib(id,"isEquidimensional")==1) 617 if(attrib(id,"isEquidimensional")==1) 619 618 { 620 attrib(prim[i],"isEquidimensional",1); 619 attrib(prim[i],"isEquidimensional",1); 621 620 } 622 621 } … … 624 623 { 625 624 attrib(prim[i],"isEquidimensional",0); 626 } 625 } 627 626 if( typeof(attrib(id,"isIsolatedSingularity"))=="int" ) 628 627 { 629 if(attrib(id,"isIsolatedSingularity")==1) 628 if(attrib(id,"isIsolatedSingularity")==1) 630 629 {attrib(prim[i],"isIsolatedSingularity",1); } 631 630 } 632 631 else 633 632 { 634 attrib(prim[i],"isIsolatedSingularity",0); 633 attrib(prim[i],"isIsolatedSingularity",0); 635 634 } 636 635 637 636 keepresult=normalizationPrimes(prim[i],maxideal(1)); 638 637 for(j=1;j<=size(keepresult);j++) … … 688 687 ideal PP=fetch(BAS,ihp); 689 688 export PP; 690 export KK; 689 export KK; 691 690 result=newR7; 692 691 setring BAS; … … 859 858 // if((dim(SM[1]))==depth) 860 859 // { 861 // attrib(SM[2],"isCohenMacaulay",1); 860 // attrib(SM[2],"isCohenMacaulay",1); 862 861 // "it is CohenMacaulay"; 863 // } 862 // } 864 863 // } 865 864 866 865 //compute the singular locus+lower dimensional components 867 866 if(((attrib(SM[2],"isIsolatedSingularity")==0)||(homog(SM[2])==0)) … … 946 945 " "; 947 946 maxideal(1); 948 " "; 947 " "; 949 948 " "; 950 949 } … … 963 962 // export SB,MB; 964 963 // result=BAS; 965 // return(result); 964 // return(result); 966 965 // } 967 966 // timer-ti; … … 972 971 if(RR[2]==0) 973 972 { 974 def newR=RR[1]; 973 def newR=RR[1]; 975 974 setring newR; 976 975 map psi=BAS,endphi; … … 980 979 // timer-ti; 981 980 setring BAS; 982 return(tluser); 981 return(tluser); 983 982 } 984 983 MB=SM[2]; … … 992 991 setring BAS; 993 992 return(result); 994 993 995 994 } 996 995 else … … 1026 1025 } 1027 1026 } 1028 1027 1029 1028 //test for non-normality 1030 1029 //Hom(I,I)<>R … … 1052 1051 // export SB,MB; 1053 1052 // result=BAS; 1054 // return(result); 1053 // return(result); 1055 1054 // } 1056 1055 // timer-ti; … … 1058 1057 // list RR=SM[1],SM[2],JM[2],SL[1]; 1059 1058 // ti=timer; 1060 list RS; 1059 list RS; 1061 1060 // list RS=HomJJ(RR); 1062 1061 // timer-ti; … … 1068 1067 // list tluser=normalizationPrimes(SM); 1069 1068 // setring BAS; 1070 // return(tluser); 1069 // return(tluser); 1071 1070 // } 1072 1071 … … 1077 1076 } 1078 1077 // ti=timer; 1079 1078 1080 1079 1081 1080 if((attrib(JM[2],"isRad")==0)&&(attrib(SM[2],"isEquidimensional")==0)) … … 1083 1082 //J=radical(JM[2]); 1084 1083 J=radical(SM[2]+ideal(SL[1])); 1085 1084 1086 1085 // evtl. test auf J=SM[2]+ideal(SL[1]) dann schon normal 1087 1086 } … … 1103 1102 // timer-ti; 1104 1103 1105 JM=J,J; 1104 JM=J,J; 1106 1105 1107 1106 //evtl. fuer SL[1] anderen Nichtnullteiler aus J waehlen … … 1117 1116 // keepresult1=insert(keepresult1,keepresult2[lauf]); 1118 1117 // } 1119 // return(keepresult1); 1118 // return(keepresult1); 1120 1119 // } 1121 1120 RR=SM[1],SM[2],JM[2],SL[1]; … … 1140 1139 export KK; 1141 1140 setring BAS; 1142 // return(RS[1]); 1141 // return(RS[1]); 1143 1142 return(lastR); 1144 1143 } … … 1153 1152 // normalizationPrimes(endid); 1154 1153 setring BAS; 1155 return(tluser); 1154 return(tluser); 1156 1155 } 1157 1156 else … … 1162 1161 +ordstr(basering)+");"; 1163 1162 if(y==1) 1164 { 1163 { 1165 1164 "zero-divisor found"; 1166 1165 } … … 1180 1179 } 1181 1180 attrib(vid,"isCompleteIntersection",0); 1182 1181 1183 1182 keepresult1=normalizationPrimes(vid,ihp); 1184 1183 … … 1189 1188 execute "ring newR2="+charstr(basering)+",("+varstr(basering)+"),(" 1190 1189 +ordstr(basering)+");"; 1191 1190 1192 1191 ideal vid=fetch(BAS,new2); 1193 1192 ideal ihp=fetch(BAS,ihp); … … 1213 1212 keepresult1=insert(keepresult1,keepresult2[lauf]); 1214 1213 } 1215 return(keepresult1); 1216 } 1214 return(keepresult1); 1215 } 1217 1216 } 1218 1217 example 1219 1218 { "EXAMPLE:";echo = 2; 1219 LIB"normal.lib"; 1220 1220 //Huneke 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1221 ring qr=31991,(a,b,c,d,e),dp; 1222 ideal i= 1223 5abcde-a5-b5-c5-d5-e5, 1224 ab3c+bc3d+a3be+cd3e+ade3, 1225 a2bc2+b2cd2+a2d2e+ab2e2+c2de2, 1226 abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5, 1227 ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5, 1228 a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6, 1229 a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4, 1230 b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5; 1231 1232 list pr=normal(i); 1233 def r1=pr[1]; 1234 setring r1; 1235 KK; 1236 1236 } 1237 1238 1239 1240 ///////////////////////////////////////////////////////////////////////////// 1241 1242 LIB"normal.lib"; 1243 int aa=timer;list pr=normal(i);timer-aa; 1244 1245 1246 1247 //Huneke 1248 ring qr=31991,(a,b,c,d,e),dp; 1249 ideal i= 1250 5abcde-a5-b5-c5-d5-e5, 1251 ab3c+bc3d+a3be+cd3e+ade3, 1252 a2bc2+b2cd2+a2d2e+ab2e2+c2de2, 1253 abc5-b4c2d-2a2b2cde+ac3d2e-a4de2+bcd2e3+abe5, 1254 ab2c4-b5cd-a2b3de+2abc2d2e+ad4e2-a2bce3-cde5, 1255 a3b2cd-bc2d4+ab2c3e-b5de-d6e+3abcd2e2-a2be4-de6, 1256 a4b2c-abc2d3-ab5e-b3c2de-ad5e+2a2bcde2+cd2e4, 1257 b6c+bc6+a2b4e-3ab2c2de+c4d2e-a3cde2-abd3e2+bce5; 1258 1259 1260 //Vasconcelos 1261 ring r=32003,(x,y,z,w,t),dp; 1262 ideal i= 1263 x2+zw, 1264 y3+xwt, 1265 xw3+z3t+ywt2, 1266 y2w4-xy2z2t-w3t3; 1267 1268 //Theo1 1269 ring r=32003,(x,y,z),wp(2,3,6); 1270 ideal i=zy2-zx3-x6; 1271 1272 //Theo2 1273 ring r=32003,(x,y,z),wp(3,4,12); 1274 ideal i=z*(y3-x4)+x8; 1275 1276 //Theo2a 1277 ring r=32003,(T(1..4)),wp(3,4,12,17); 1278 ideal i= 1279 T(1)^8-T(1)^4*T(3)+T(2)^3*T(3), 1280 T(1)^4*T(2)^2-T(2)^2*T(3)+T(1)*T(4), 1281 T(1)^7+T(1)^3*T(2)^3-T(1)^3*T(3)+T(2)*T(4), 1282 T(1)^6*T(2)*T(3)+T(1)^2*T(2)^4*T(3)+T(1)^3*T(2)^2*T(4)-T(1)^2*T(2)*T(3)^2+T(4)^2; 1283 1284 //Theo3 1285 ring r=32003,(x,y,z),wp(3,5,15); 1286 ideal i=z*(y3-x5)+x10; 1287 1288 1289 //Theo4 1290 ring r=32003,(x,y,z),dp; 1291 ideal i=(x-y)*(x-z)*(y-z); 1292 1293 //Theo5 1294 ring r=32003,(x,y,z),wp(2,1,2); 1295 ideal i=z3-xy4; 1296 1297 //Theo6 1298 ring r=32003,(x,y,z),dp; 1299 ideal i=x2y2+x2z2+y2z2; 1300 1301 ring r=32003,(a,b,c,d,e,f),dp; 1302 ideal i= 1303 bf, 1304 af, 1305 bd, 1306 ad; 1307 1308 //Beispiel, wo vorher Primaerzerlegung schneller 1309 //ist CM 1310 //Sturmfels 1311 ring r=32003,(b,s,t,u,v,w,x,y,z),dp; 1312 ideal i= 1313 bv+su, 1314 bw+tu, 1315 sw+tv, 1316 by+sx, 1317 bz+tx, 1318 sz+ty, 1319 uy+vx, 1320 uz+wx, 1321 vz+wy, 1322 bvz; 1323 1324 //J S/Y 1325 ring r=32003,(x,y,z,t),dp; 1326 ideal i= 1327 x2z+xzt, 1328 xyz, 1329 xy2-xyt, 1330 x2y+xyt; 1331 1332 //St_S/Y 1333 ring r=32003,(b,s,t,u,v,w,x,y,z),dp; 1334 ideal i= 1335 wy-vz, 1336 vx-uy, 1337 tv-sw, 1338 su-bv, 1339 tuy-bvz; 1340 1341 //dauert laenger 1342 //Horrocks: 1343 ring r=32003,(a,b,c,d,e,f),dp; 1344 ideal i= 1345 adef-16000be2f+16001cef2, 1346 ad2f+8002bdef+8001cdf2, 1347 abdf-16000b2ef+16001bcf2, 1348 a2df+8002abef+8001acf2, 1349 ad2e-8000bde2-7999cdef, 1350 acde-16000bce2+16001c2ef, 1351 a2de-8000abe2-7999acef, 1352 acd2+8002bcde+8001c2df, 1353 abd2-8000b2de-7999bcdf, 1354 a2d2+9603abde-10800b2e2-9601acdf+800bcef+11601c2f2, 1355 abde-8000b2e2-acdf-16001bcef-8001c2f2, 1356 abcd-16000b2ce+16001bc2f, 1357 a2cd+8002abce+8001ac2f, 1358 a2bd-8000ab2e-7999abcf, 1359 ab3f-3bdf3, 1360 a2b2f-2adf3-16000bef3+16001cf4, 1361 a3bf+4aef3, 1362 ac3e-10668cde3, 1363 a2c2e+10667ade3+16001be4+5334ce3f, 1364 a3ce+10669ae3f, 1365 bc3d+8001cd3e, 1366 ac3d+8000bc3e+16001cd2e2+8001c4f, 1367 b2c2d+16001ad4+4000bd3e+12001cd3f, 1368 b2c2e-10668bc3f-10667cd2ef, 1369 abc2e-cde2f, 1370 b3cd-8000bd3f, 1371 b3ce-10668b2c2f-10667bd2ef, 1372 abc2f-cdef2, 1373 a2bce-16000be3f+16001ce2f2, 1374 ab3d-8000b4e-8001b3cf+16000bd2f2, 1375 ab2cf-bdef2, 1376 a2bcf-16000be2f2+16001cef3, 1377 a4d-8000a3be+8001a3cf-2ae2f2; 1378 1379 1380 ring r=32003,(b,s,t,u,v,w,x,y,z),dp; 1381 1382 ideal k= 1383 wy-vz, 1384 vx-uy, 1385 tv-sw, 1386 su-bv, 1387 tuy-bvz; 1388 ideal j=x2y2+x2z2+y2z2; 1389 ideal i=mstd(intersect(j,k))[2]; 1390 1391 //22 1392 ring r=32003,(b,s,t,u,v,w,x,y,z),dp; 1393 ideal i= 1394 wx2y3-vx2y2z+wx2yz2+wy3z2-vx2z3-vy2z3, 1395 vx3y2-ux2y3+vx3z2-ux2yz2+vxy2z2-uy3z2, 1396 tvx2y2-swx2y2+tvx2z2-swx2z2+tvy2z2-swy2z2, 1397 sux2y2-bvx2y2+sux2z2-bvx2z2+suy2z2-bvy2z2, 1398 tux2y3-bvx2y2z+tux2yz2+tuy3z2-bvx2z3-bvy2z3; 1399 1400 1401 //riemenschneider 1402 //33 1403 //normal+primary 3 1404 //primary 9 1405 //radical 1 1406 //minAssPrimes 2 1407 ring r=32000,(p,q,s,t,u,v,w,x,y,z),wp(1,1,1,1,1,1,2,1,1,1); 1408 ideal i= 1409 xz, 1410 vx, 1411 ux, 1412 su, 1413 qu, 1414 txy, 1415 stx, 1416 qtx, 1417 uv2z-uwz, 1418 uv3-uvw, 1419 puv2-puw; 1420 1421 -
Singular/LIB/primdec.lib
r437e2c2 re801fe 1 // $Id: primdec.lib,v 1. 8 1998-01-27 18:51:23 Singular Exp $1 // $Id: primdec.lib,v 1.9 1998-02-16 12:07:04 pfister Exp $ 2 2 /////////////////////////////////////////////////////// 3 3 // primdec.lib … … 5 5 // the ideas of Gianni,Trager,Zacharias 6 6 // written by Gerhard Pfister 7 // 8 // algorithms for primary decomposition based on 9 // the ideas of Shimoyama/Yokoyama 10 // written by Wolfram Decker and Hans Schoenemann 7 11 ////////////////////////////////////////////////////// 8 12 … … 12 16 //computes the minimal associated primes of the ideal I 13 17 14 decomp(ideal I)18 primdecGTZ (ideal I) 15 19 // Computes a complete primary decomposition via Gianni,Trager,Zacharias 16 20 … … 24 28 //computes the radicals of the equidimensional parts of I 25 29 30 min_ass_prim_charsets (ideal I, int choose) 31 // minimal associated primes via the characteristic set 32 // package written by Michael Messollen. 33 // The integer choose must be either 0 or 1. 34 // If choose=0, the given ordering of the variables is used. 35 // If choose=1, the system tries to find an "optimal ordering", 36 // which in some cases may considerably speed up the algorithm 37 38 // You may also may want to try one of the algorithms for 39 // minimal associated primes in the the library 40 // primdec.lib, written by Gerhard Pfister. 41 // These algorithms are variants of the algorithm 42 // of Gianni-Trager-Zacharias 43 44 primdecSY (ideal I, int choose) 45 46 // Computes a complete primary decomposition via 47 // a variant of the pseudoprimary approach of 48 // Shimoyama-Yokoyama. 49 // The integer choose must be either 0, 1, 2 or 3. 50 // If choose=0, min_ass_prim_charsets with the given 51 // ordering of the variables is used. 52 // If choose=1, min_ass_prim_charsets with the "optimized" 53 // ordering of the variables is used. 54 // If choose=2, minAssPrimes from primdec.lib is used 55 // If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used 56 57 LIB "general.lib"; 58 LIB "elim.lib"; 59 LIB "poly.lib"; 26 60 LIB "random.lib"; 27 LIB "elim.lib";28 61 /////////////////////////////////////////////////////////////////////////////// 29 62 … … 62 95 if(ordstr(basering)[1,2]!="dp") 63 96 { 64 execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"), dp;";97 execute "ring @Phelp=("+charstr(@P)+"),("+varstr(@P)+"),(C,dp);"; 65 98 ideal inew=std(imap(@P,id)); 66 99 ideal @h=imap(@P,h); … … 175 208 ideal fac=tsil[2]; 176 209 poly f=tsil[3]; 177 210 178 211 ideal star=quotient(laedi,f); 179 180 212 action=1; 181 213 … … 190 222 { 191 223 g=1; 224 verg=laedi; 225 192 226 for(j=1;j<=size(fac);j++) 193 227 { … … 198 232 } 199 233 verg=quotient(laedi,g); 234 200 235 if(specialIdealsEqual(verg,star)==1) 201 236 { … … 211 246 } 212 247 } 213 l=star,fac,f; 248 l=star,fac,f; 214 249 return(l); 215 250 } 216 217 251 218 252 //////////////////////////////////////////////////////////////////////////////// 219 253 proc testFactor(list act,poly p) 220 254 { 255 poly keep=p; 256 221 257 int i; 222 258 poly q=act[1][1]^act[2][1]; … … 230 266 { 231 267 "ERROR IN FACTOR"; 268 basering; 269 232 270 act; 271 keep; 272 pause; 273 233 274 p; 234 275 q; … … 274 315 return(@l); 275 316 } 276 @l=factorize(p,2); 277 if(npars(basering)>0) 278 { 317 // @l=factorize(p,2); 318 @l=factorize(p); 319 // if(npars(basering)>0) 320 // { 279 321 for(@j=1;@j<=size(@l[1]);@j++) 280 322 { … … 320 362 } 321 363 } 322 }364 // } 323 365 return(@l); 324 366 } … … 381 423 } 382 424 attrib(k2,"isSB",1); 383 if(size(reduce(k1,k2 ))==0)425 if(size(reduce(k1,k2,1))==0) 384 426 { 385 427 return(1); … … 514 556 m=m-1; 515 557 } 516 //check whether i[m] =(c*var(n)+h)^e modulo prm for some 558 //check whether i[m] =(c*var(n)+h)^e modulo prm for some 517 559 //h in K[var(n+1),...,var(nvars(basering))], c in K 518 560 //if not (0) is returned, else var(n)+h is added to prm … … 524 566 i[m]=poly(e)^e*leadcoef(i[m])^(e-1)*i[m]; 525 567 526 if (reduce(i[m]-t^e,prm ) !=0)568 if (reduce(i[m]-t^e,prm,1) !=0) 527 569 { 528 570 return(ideal(0)); … … 606 648 { 607 649 i++; 608 if((size(ser)>0)&&(size(reduce(ser,l[2*i-1] ))==0))650 if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0)) 609 651 { 610 652 l[2*i-1]=ideal(1); … … 736 778 for(i=1;i<=size(l)/2;i++) 737 779 { 738 if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1)) 780 attrib(l[2*i-1],"isSB",1); 781 782 if((size(ser)>0)&&(size(reduce(ser,l[2*i-1],1))==0)&&(deg(l[2*i-1][1])>0)) 783 { 784 "Achtung in split"; 785 786 l[2*i-1]=ideal(1); 787 l[2*i]=ideal(1); 788 } 789 if((size(l[2*i])==0)&&(specialIdealsEqual(keepprime[i],l[2*i-1])!=1)) 739 790 { 740 791 keepprime[i]=std(keepprime[i]); … … 782 833 def @P = basering; 783 834 int nva = nvars(basering); 784 int @k,@s,@n,@k1 ;785 list primary,lres, act,@lh,@wh;786 map phi,psi ;787 ideal jmap, helpprim,@qh,@qht;835 int @k,@s,@n,@k1,zz; 836 list primary,lres,lres1,act,@lh,@wh; 837 map phi,psi,phi1,psi1; 838 ideal jmap,jmap1,jmap2,helpprim,@qh,@qht,ser1; 788 839 intvec @vh,@hilb; 789 840 string @ri; … … 796 847 return(primary); 797 848 } 798 j=interred(j); 849 850 j=interred(j); 799 851 attrib(j,"isSB",1); 800 852 if(vdim(j)==deg(j[1])) 801 { 802 if((size(ser)>0)&&(size(reduce(ser,j))==0)) 803 { 804 primary[1]=ideal(1); 805 primary[2]=ideal(1); 806 return(primary); 807 } 853 { 808 854 act=factor(j[1]); 809 855 for(@k=1;@k<=size(act[1]);@k++) … … 822 868 @qh[1]=act[1][@k]; 823 869 primary[2*@k]=interred(@qh); 824 } 825 return(primary); 870 attrib( primary[2*@k-1],"isSB",1); 871 872 if((size(ser)>0)&&(size(reduce(ser,primary[2*@k-1],1))==0)) 873 { 874 primary[2*@k-1]=ideal(1); 875 primary[2*@k]=ideal(1); 876 } 877 } 878 return(primary); 826 879 } 827 880 … … 829 882 { 830 883 primary[1]=j; 831 if((size(ser)>0)&&(size(reduce(ser,j ))==0))884 if((size(ser)>0)&&(size(reduce(ser,j,1))==0)) 832 885 { 833 886 primary[1]=ideal(1); … … 894 947 } 895 948 else 896 { 949 { 897 950 primary[1]=j; 898 951 if((size(#)>0)&&(act[2][1]>1)) … … 913 966 if(size(#)==0) 914 967 { 968 915 969 primary=splitPrimary(primary,ser,@wr,act); 970 916 971 } 917 972 … … 922 977 923 978 @k=0; 924 int zz;925 979 while(@k<(size(primary)/2)) 926 980 { … … 937 991 } 938 992 } 993 939 994 @k=0; 995 ideal keep; 940 996 while(@k<(size(primary)/2)) 941 997 { … … 943 999 if (size(primary[2*@k])==0) 944 1000 { 945 // "the univariat polynomials";946 // int qwe=timer;947 // system("finduni",primary[2*@k-1]);948 1001 949 1002 jmap=randomLast(100); 950 // timer-qwe; 951 // primary[2*@k-1]; 952 // pause; 1003 jmap1=maxideal(1); 1004 jmap2=maxideal(1); 1005 @qht=primary[2*@k-1]; 1006 953 1007 for(@n=2;@n<=size(primary[2*@k-1]);@n++) 954 1008 { 955 1009 if(deg(lead(primary[2*@k-1][@n]))==1) 956 1010 { 1011 for(zz=1;zz<=nva;zz++) 1012 { 1013 if(lead(primary[2*@k-1][@n])/var(zz)!=0) 1014 { 1015 jmap1[zz]=-1/leadcoef(primary[2*@k-1][@n])*primary[2*@k-1][@n] 1016 +2/leadcoef(primary[2*@k-1][@n])*lead(primary[2*@k-1][@n]); 1017 jmap2[zz]=primary[2*@k-1][@n]; 1018 @qht[@n]=var(zz); 1019 1020 } 1021 } 957 1022 jmap[nva]=subst(jmap[nva],lead(primary[2*@k-1][@n]),0); 958 1023 } 959 1024 } 1025 if(size(subst(jmap[nva],var(1),0)-var(nva))!=0) 1026 { 1027 jmap[nva]=subst(jmap[nva],var(1),0); 1028 } 1029 phi1=@P,jmap1; 960 1030 phi=@P,jmap; 961 jmap[nva]=-(jmap[nva]-2*var(nva)); 1031 1032 for(@n=1;@n<=nva;@n++) 1033 { 1034 jmap[@n]=-(jmap[@n]-2*var(@n)); 1035 } 962 1036 psi=@P,jmap; 963 @qht=primary[2*@k-1]; 1037 psi1=@P,jmap2; 1038 964 1039 @qh=phi(@qht); 1040 965 1041 if(npars(@P)>0) 966 1042 { 967 1043 @ri= "ring @Phelp =" 968 +string(char(@P))+",("+varstr(@P)+","+parstr(@P)+",@t),dp;"; 1044 +string(char(@P))+", 1045 ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; 969 1046 } 970 1047 else 971 1048 { 972 1049 @ri= "ring @Phelp =" 973 +string(char(@P))+",("+varstr(@P)+",@t), dp;";1050 +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; 974 1051 } 975 1052 execute(@ri); … … 979 1056 @hilb=hilb(@qh1,1); 980 1057 @ri= "ring @Phelp1 =" 981 +string(char(@P))+",("+varstr(@Phelp)+"), lp;";1058 +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; 982 1059 execute(@ri); 983 1060 ideal @qh=homog(imap(@P,@qh),@t); … … 989 1066 kill @Phelp1; 990 1067 @qh=clearSB(@qh); 991 attrib(@qh,"isSB",1); 992 993 @lh=zero_decomp (@qh,psi(ser),@wr); 994 1068 attrib(@qh,"isSB",1); 1069 ser1=phi1(ser); 1070 @lh=zero_decomp (@qh,phi(ser1),@wr); 1071 // @lh=zero_decomp (@qh,psi(ser),@wr); 1072 1073 995 1074 kill lres; 996 1075 list lres; … … 999 1078 helpprim=@lh[2]; 1000 1079 lres[1]=primary[2*@k-1]; 1001 lres[2]=psi(helpprim); 1002 if(size(reduce(lres[2],lres[1]))==0) 1080 ser1=psi(helpprim); 1081 lres[2]=psi1(ser1); 1082 if(size(reduce(lres[2],lres[1],1))==0) 1003 1083 { 1004 1084 primary[2*@k]=primary[2*@k-1]; … … 1014 1094 { 1015 1095 @f=act[1][@n]^act[2][@n]; 1016 lres[2*@n-1]=interred(primary[2*@k-1]+psi(@f)); 1096 ser1=psi(@f); 1097 lres[2*@n-1]=interred(primary[2*@k-1]+psi1(ser1)); 1017 1098 helpprim=@lh[2*@n]; 1018 lres[2*@n]=psi(helpprim); 1099 ser1=psi(helpprim); 1100 lres[2*@n]=psi1(ser1); 1019 1101 } 1020 1102 } 1021 1103 else 1022 1104 { 1023 lres=psi(@lh); 1105 lres1=psi(@lh); 1106 lres=psi1(lres1); 1024 1107 } 1025 1108 } … … 1027 1110 { 1028 1111 @ri= "ring @Phelp =" 1029 +string(char(@P))+",("+varstr(@P)+","+parstr(@P)+",@t),dp;"; 1112 +string(char(@P))+", 1113 ("+varstr(@P)+","+parstr(@P)+",@t),(C,dp);"; 1030 1114 } 1031 1115 else 1032 1116 { 1033 1117 @ri= "ring @Phelp =" 1034 +string(char(@P))+",("+varstr(@P)+",@t), dp;";1118 +string(char(@P))+",("+varstr(@P)+",@t),(C,dp);"; 1035 1119 } 1036 1120 execute(@ri); … … 1069 1153 } 1070 1154 @ri= "ring @Phelp1 =" 1071 +string(char(@P))+",("+varstr(@Phelp)+"), lp;";1155 +string(char(@P))+",("+varstr(@Phelp)+"),(C,lp);"; 1072 1156 execute(@ri); 1073 1157 list @lr=imap(@Phelp,@lr); … … 1156 1240 return(1) 1157 1241 } 1158 p= gcd(p,i[k]);1242 p=GCD(p,i[k]); 1159 1243 if(deg(p)==0) 1160 1244 { … … 1204 1288 def bsr= basering; 1205 1289 string @id=string(h); 1206 execute "ring @r=0,("+@pa+","+varstr(bsr)+"), dp;";1290 execute "ring @r=0,("+@pa+","+varstr(bsr)+"),(C,dp);"; 1207 1291 execute "ideal @i="+@id+";"; 1208 1292 poly @p=lcmP(@i); … … 1250 1334 if(deg(i[k])!=0) 1251 1335 { 1252 q= gcd(p,i[k]);1336 q=GCD(p,i[k]); 1253 1337 if(deg(q)==0) 1254 1338 { … … 1514 1598 @jh=@jh+var(@n); 1515 1599 } 1516 quotring=quotring+"), lp;";1600 quotring=quotring+"),(C,lp);"; 1517 1601 1518 1602 return(quotring); … … 1577 1661 #[1]=1; 1578 1662 def @P=basering; 1663 list qr=simplifyIdeal(i); 1664 map phi=@P,qr[2]; 1665 i=qr[1]; 1666 1579 1667 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(" 1580 1668 +ordstr(basering)+");"; … … 1593 1681 setring @P; 1594 1682 list @res=imap(gnir,tluser); 1595 return( @res);1683 return(phi(@res)); 1596 1684 } 1597 1685 list @res,empty; … … 1606 1694 setring @P; 1607 1695 list @res=maxideal(1); 1608 return( @res);1696 return(phi(@res)); 1609 1697 } 1610 1698 if(dim(@pr[1])>1) 1611 1699 { 1612 1700 setring @P; 1613 kill gnir; 1614 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;"; 1701 // kill gnir; 1702 execute "ring gnir1 = ("+charstr(basering)+"), 1703 ("+varstr(basering)+"),(C,lp);"; 1615 1704 ideal i=fetch(@P,i); 1616 1705 list @pr=facstd(i); 1617 ideal ser; 1618 } 1619 } 1620 option(noredSB); 1706 // ideal ser; 1707 setring gnir; 1708 @pr=fetch(gnir1,@pr); 1709 kill gnir1; 1710 } 1711 } 1712 option(noredSB); 1621 1713 int j,k,odim,ndim,count; 1622 1714 attrib(@pr[1],"isSB",1); … … 1662 1754 for(j=1;j<=size(@pr);j++) 1663 1755 { 1664 @res[j]=decomp(@pr[j],2); 1665 // @res[j]=decomp(@pr[j],2,@pr[j],ser); 1666 // for(k=1;k<=size(@res[j]);k++) 1667 // { 1668 // ser=intersect1(ser,@res[j][k]); 1669 // } 1756 //@pr[j]; 1757 //pause; 1758 @res[j]=decomp(@pr[j],2); 1759 // @res[j]=decomp(@pr[j],2,@pr[j],ser); 1760 // for(k=1;k<=size(@res[j]);k++) 1761 // { 1762 // ser=intersect(ser,@res[j][k]); 1763 // } 1670 1764 } 1671 1765 } … … 1674 1768 setring @P; 1675 1769 list @res=imap(gnir,@res); 1676 return( @res);1770 return(phi(@res)); 1677 1771 } 1678 1772 example … … 1696 1790 def P=basering; 1697 1791 1698 execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"), lp;";1792 execute "ring ir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);"; 1699 1793 list l=fetch(P,li); 1700 1794 list @erg; … … 1735 1829 if(size(reduce(i1,i2,1))==0) 1736 1830 { 1737 if(size(reduce(@erg[i],@erg[k] ))==0)1831 if(size(reduce(@erg[i],@erg[k],1))==0) 1738 1832 { 1739 1833 @erg[k]=ideal(1); … … 1743 1837 if(size(reduce(i2,i1,1))==0) 1744 1838 { 1745 if(size(reduce(@erg[k],@erg[i] ))==0)1839 if(size(reduce(@erg[k],@erg[i],1))==0) 1746 1840 { 1747 1841 break; … … 1780 1874 return(radical(i,1)); 1781 1875 } 1876 1782 1877 /////////////////////////////////////////////////////////////////////////////// 1783 proc decomp 1878 proc decomp(ideal i,list #) 1784 1879 USAGE: decomp(i); i ideal (for primary decomposition) (resp. 1785 1880 decomp(i,1); (for the minimal associated primes) ) … … 1793 1888 list primary,indep,ltras; 1794 1889 intvec @vh,isat; 1795 int @wr,@k,@n,@m,@n1,@n2,@n3,homo ;1890 int @wr,@k,@n,@m,@n1,@n2,@n3,homo,seri,keepdi; 1796 1891 ideal peek=i; 1797 1892 ideal ser,tras; 1798 1893 1799 int @aa=timer;1800 1801 homo=homog(i);1802 1894 if(size(#)>0) 1803 1895 { … … 1807 1899 if(size(#)>1) 1808 1900 { 1901 seri=1; 1809 1902 peek=#[2]; 1810 1903 ser=#[3]; … … 1813 1906 else 1814 1907 { 1815 peek=#[1]; 1816 ser=#[2]; 1817 } 1818 } 1819 1908 seri=1; 1909 peek=#[1]; 1910 ser=#[2]; 1911 } 1912 } 1913 1914 homo=homog(i); 1915 1820 1916 if(homo==1) 1821 1917 { 1822 ltras=mstd(i); 1823 attrib(ltras[1],"isSB",1); 1824 tras=ltras[1]; 1825 if(dim(tras)==0) 1826 { 1918 if(attrib(i,"isSB")!=1) 1919 { 1920 ltras=mstd(i); 1921 attrib(ltras[1],"isSB",1); 1922 } 1923 else 1924 { 1925 ltras=i,i; 1926 } 1927 tras=ltras[1]; 1928 if(dim(tras)==0) 1929 { 1827 1930 primary[1]=ltras[2]; 1828 1931 primary[2]=maxideal(1); … … 1836 1939 return(primary); 1837 1940 } 1838 intvec @hilb=hilb(tras,1); 1941 intvec @hilb=hilb(tras,1); 1942 intvec keephilb=@hilb; 1839 1943 } 1840 1944 … … 1853 1957 //---------------------------------------------------------------- 1854 1958 1855 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;"; 1856 1959 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),(C,lp);"; 1857 1960 option(redSB); 1858 ideal ser=fetch(@P,ser); 1859 ideal peek=std(fetch(@P,peek)); 1860 homo=homog(peek); 1961 1962 ideal ser=fetch(@P,ser); 1861 1963 1862 1964 if(homo==1) 1863 1965 { 1864 if( ordstr(@P)[1,2]!="lp")1865 { 1866 1966 if((ordstr(@P)[1]!="(C,lp)")&&(ordstr(@P)[3]!="(C,lp)")) 1967 { 1968 ideal @j=std(fetch(@P,i),@hilb); 1867 1969 } 1868 1970 else … … 1876 1978 ideal @j=std(fetch(@P,i)); 1877 1979 } 1878 1980 option(noredSB); 1981 if(seri==1) 1982 { 1983 ideal peek=fetch(@P,peek); 1984 attrib(peek,"isSB",1); 1985 } 1986 else 1987 { 1988 ideal peek=@j; 1989 } 1990 if(size(ser)==0) 1991 { 1992 ideal fried; 1993 @n=size(@j); 1994 for(@k=1;@k<=@n;@k++) 1995 { 1996 if(deg(lead(@j[@k]))==1) 1997 { 1998 fried[size(fried)+1]=@j[@k]; 1999 @j[@k]=0; 2000 } 2001 } 2002 if(size(fried)>0) 2003 { 2004 @j=simplify(@j,2); 2005 attrib(@j,"isSB",1); 2006 list pr=decomp(@j); 2007 for(@k=1;@k<=size(pr);@k++) 2008 { 2009 @j=pr[@k]+fried; 2010 pr[@k]=@j; 2011 } 2012 setring @P; 2013 return(fetch(gnir,pr)); 2014 } 2015 } 2016 1879 2017 //---------------------------------------------------------------- 1880 2018 //j is the ring … … 1883 2021 if (dim(@j)==-1) 1884 2022 { 1885 setring @P; 1886 option(noredSB); 1887 return(ideal(0)); 2023 setring @P; 2024 return(ideal(0)); 1888 2025 } 1889 2026 … … 1911 2048 } 1912 2049 setring @P; 1913 option(noredSB);1914 2050 primary=fetch(gnir,gprimary); 1915 2051 1916 2052 return(primary); 1917 2053 } 1918 2054 1919 2055 //------------------------------------------------------------------ 1920 2056 //the zero-dimensional case … … 1923 2059 if (dim(@j)==0) 1924 2060 { 1925 list gprimary= zero_decomp(@j,ser,@wr);1926 1927 setring @P;1928 option(noredSB);1929 1930 1931 1932 1933 1934 1935 2061 option(redSB); 2062 list gprimary= zero_decomp(@j,ser,@wr); 2063 option(noredSB); 2064 setring @P; 2065 primary=fetch(gnir,gprimary); 2066 if(size(ser)>0) 2067 { 2068 primary=cleanPrimary(primary); 2069 } 2070 return(primary); 2071 } 1936 2072 1937 2073 poly @gs,@gh,@p; … … 1941 2077 int jdim=dim(@j); 1942 2078 list fett; 1943 int lauf,di ;2079 int lauf,di,newtest; 1944 2080 //------------------------------------------------------------------ 1945 2081 //search for a maximal independent set indep,i.e. … … 1970 2106 indep=maxIndependSet(@j); 1971 2107 } 1972 2108 1973 2109 ideal jkeep=@j; 1974 2110 … … 1979 2115 else 1980 2116 { 1981 execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),dp;"; 1982 } 1983 ideal jwork=std(imap(gnir,@j)); 2117 execute "ring @Phelp=("+charstr(gnir)+"),("+varstr(gnir)+"),(C,dp);"; 2118 } 2119 2120 if(homo==1) 2121 { 2122 if((ordstr(@P)[3]=="d")||(ordstr(@P)[1]=="d")||(ordstr(@P)[1]=="w") 2123 ||(ordstr(@P)[3]=="w")) 2124 { 2125 ideal jwork=imap(@P,tras); 2126 attrib(jwork,"isSB",1); 2127 } 2128 else 2129 { 2130 ideal jwork=std(imap(gnir,@j),@hilb); 2131 } 2132 2133 } 2134 else 2135 { 2136 ideal jwork=std(imap(gnir,@j)); 2137 } 2138 list hquprimary; 1984 2139 poly @p,@q; 1985 ideal @h,fac ;2140 ideal @h,fac,ser; 1986 2141 di=dim(jwork); 2142 keepdi=di; 2143 1987 2144 setring gnir; 1988 2145 for(@m=1;@m<=size(indep);@m++) … … 1990 2147 isat=0; 1991 2148 @n2=0; 2149 option(redSB); 1992 2150 if((indep[@m][1]==varstr(basering))&&(@m==1)) 1993 2151 //this is the good case, nothing to do, just to have the same notations … … 1998 2156 ideal @j=fetch(gnir,@j); 1999 2157 attrib(@j,"isSB",1); 2158 ideal ser=fetch(gnir,ser); 2159 2000 2160 } 2001 2161 else … … 2013 2173 ideal @j=std(phi(@j)); 2014 2174 } 2015 } 2175 ideal ser=phi(ser); 2176 2177 } 2178 option(noredSB); 2016 2179 if((deg(@j[1])==0)||(dim(@j)<jdim)) 2017 2180 { … … 2050 2213 // @j considered in the quotientring 2051 2214 ideal @j=imap(gnir1,@j); 2052 ideal ser=imap(gnir ,ser);2215 ideal ser=imap(gnir1,ser); 2053 2216 2054 2217 kill gnir1; … … 2069 2232 } 2070 2233 //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..] 2071 2234 option(redSB); 2072 2235 list uprimary= zero_decomp(@j,ser,@wr); 2073 2236 option(noredSB); 2074 2237 } 2075 2238 else 2076 2239 { 2077 2240 list uprimary; 2241 uprimary[1]=ideal(1); 2078 2242 uprimary[2]=ideal(1); 2079 uprimary[1]=ideal(1);2080 2243 } 2081 2244 … … 2130 2293 //here the intersection with the polynomialring 2131 2294 //mentioned above is really computed 2132 2133 for(@n=@n3/2+1;@n<=@n2/2;@n++) 2134 { 2295 for(@n=@n3/2+1;@n<=@n2/2;@n++) 2296 { 2135 2297 if(specialIdealsEqual(quprimary[2*@n-1],quprimary[2*@n])) 2136 2298 { … … 2147 2309 } 2148 2310 } 2311 2149 2312 if(size(@h)>0) 2150 2313 { … … 2155 2318 @h=imap(gnir,@h); 2156 2319 if(@wr!=1) 2157 // if(@wr==0) 2158 { 2159 @q=minSat(jwork,@h)[2]; 2320 { 2321 @q=minSat(jwork,@h)[2]; 2160 2322 } 2161 2323 else … … 2177 2339 } 2178 2340 jwork=std(jwork,@q); 2179 if(dim(jwork)<di) 2341 keepdi=dim(jwork); 2342 if(keepdi<di) 2180 2343 { 2181 2344 setring gnir; … … 2195 2358 { 2196 2359 @j=ideal(1); 2360 quprimary[1]=ideal(1); 2197 2361 quprimary[2]=ideal(1); 2198 quprimary[1]=ideal(1); 2199 } 2200 2362 } 2363 if((size(quprimary)==0)) 2364 { 2365 keepdi=di-1; 2366 } 2201 2367 //--------------------------------------------------------------- 2202 2368 //notice that j=sat(j,gh) intersected with (j,gh^n) … … 2205 2371 if((deg(@j[1])!=0)&&(@wr!=1)) 2206 2372 { 2207 int uq=size(quprimary); 2208 if(uq>0) 2209 { 2373 if(size(quprimary)>0) 2374 { 2375 setring @Phelp; 2376 ser=imap(gnir,ser); 2377 hquprimary=imap(gnir,quprimary); 2210 2378 if(@wr==0) 2211 2379 { 2212 ideal htest=quprimary[1]; 2213 2214 for (@n1=2;@n1<=size(quprimary)/2;@n1++) 2380 ideal htest=hquprimary[1]; 2381 for (@n1=2;@n1<=size(hquprimary)/2;@n1++) 2215 2382 { 2216 htest=intersect(htest, quprimary[2*@n1-1]);2383 htest=intersect(htest,hquprimary[2*@n1-1]); 2217 2384 } 2218 2385 } 2219 2386 else 2220 2387 { 2221 ideal htest= quprimary[2];2222 2223 for (@n1=2;@n1<=size( quprimary)/2;@n1++)2388 ideal htest=hquprimary[2]; 2389 2390 for (@n1=2;@n1<=size(hquprimary)/2;@n1++) 2224 2391 { 2225 htest=intersect(htest, quprimary[2*@n1]);2392 htest=intersect(htest,hquprimary[2*@n1]); 2226 2393 } 2227 2394 } 2395 2228 2396 if(size(ser)>0) 2229 { 2230 htest=intersect(htest,ser); 2231 } 2232 ser=std(htest); 2233 } 2234 //we are not ready yet 2235 if (specialIdealsEqual(ser,peek)!=1) 2236 { 2397 { 2398 ser=intersect(htest,ser); 2399 } 2400 else 2401 { 2402 ser=htest; 2403 } 2404 setring gnir; 2405 ser=imap(@Phelp,ser); 2406 } 2407 if(size(reduce(ser,peek,1))!=0) 2408 { 2237 2409 for(@m=1;@m<=size(restindep);@m++) 2238 2410 { 2411 // if(restindep[@m][3]>=keepdi) 2412 // { 2239 2413 isat=0; 2240 2414 @n2=0; 2415 option(redSB); 2416 2241 2417 if(restindep[@m][1]==varstr(basering)) 2242 2418 //this is the good case, nothing to do, just to have the same notations … … 2256 2432 if(homo==1) 2257 2433 { 2258 ideal @j=std(phi(jkeep), @hilb);2434 ideal @j=std(phi(jkeep),keephilb); 2259 2435 } 2260 2436 else … … 2262 2438 ideal @j=std(phi(jkeep)); 2263 2439 } 2440 ideal ser=phi(ser); 2264 2441 } 2265 2442 option(noredSB); 2443 2266 2444 for (lauf=1;lauf<=size(@j);lauf++) 2267 2445 { … … 2295 2473 // @j considered in the quotientring 2296 2474 ideal @j=imap(gnir1,@j); 2297 ideal ser=imap(gnir ,ser);2475 ideal ser=imap(gnir1,ser); 2298 2476 2299 2477 kill gnir1; … … 2312 2490 } 2313 2491 //the primary decomposition of j*K(var(nnp+1),..,var(nva))[..the rest..] 2314 2315 list uprimary= zero_decomp(@j,ser,@wr); 2316 2492 2493 option(redSB); 2494 list uprimary= zero_decomp(@j,ser,@wr); 2495 option(noredSB); 2496 2497 2317 2498 //we need the intersection of the ideals in the list quprimary with the 2318 2499 //polynomialring, i.e. let q=(f1,...,fr) in the quotientring such an ideal … … 2349 2530 @n2=size(quprimary); 2350 2531 @n3=@n2; 2351 2532 2352 2533 for(@n1=1;@n1<=size(collectprimary)/2;@n1++) 2353 2534 { … … 2362 2543 } 2363 2544 } 2364 2545 2546 2365 2547 //here the intersection with the polynomialring 2366 2548 //mentioned above is really computed … … 2380 2562 } 2381 2563 quprimary[2*@n]=sat2(quprimary[2*@n],lnew[2*@n])[1]; 2382 } 2383 if(@wr==0) 2384 { 2385 ser=std(intersect(ser,quprimary[2*@n-1])); 2386 } 2387 else 2388 { 2389 ser=std(intersect(ser,quprimary[2*@n])); 2390 } 2564 } 2391 2565 } 2392 } 2393 //we are not ready yet 2394 if (specialIdealsEqual(ser,peek)!=1) 2566 if(@n2>=@n3+2) 2567 { 2568 setring @Phelp; 2569 ser=imap(gnir,ser); 2570 hquprimary=imap(gnir,quprimary); 2571 for(@n=@n3/2+1;@n<=@n2/2;@n++) 2572 { 2573 if(@wr==0) 2574 { 2575 ser=intersect(ser,hquprimary[2*@n-1]); 2576 } 2577 else 2578 { 2579 ser=intersect(ser,hquprimary[2*@n]); 2580 } 2581 } 2582 setring gnir; 2583 ser=imap(@Phelp,ser); 2584 } 2585 2586 // } 2587 } 2588 if(size(reduce(ser,peek,1))!=0) 2395 2589 { 2396 2590 if(@wr>0) … … 2404 2598 // here we collect now both results primary(sat(j,gh)) 2405 2599 // and primary(j,gh^n) 2406 2407 2600 @n=size(quprimary); 2408 2601 for (@k=1;@k<=size(htprimary);@k++) … … 2412 2605 } 2413 2606 } 2607 2414 2608 } 2415 2609 //------------------------------------------------------------ … … 2417 2611 //the final result: primary 2418 2612 //------------------------------------------------------------ 2613 2419 2614 setring @P; 2420 2615 primary=imap(gnir,quprimary); 2421 2422 option(noredSB);2423 2616 return(primary); 2424 2617 } … … 2438 2631 2439 2632 /////////////////////////////////////////////////////////////////////////////// 2633 proc primdecGTZ(ideal i) 2634 { 2635 return(decomp(i)); 2636 } 2637 /////////////////////////////////////////////////////////////////////////////// 2638 proc primdecSY(ideal i,int j) 2639 { 2640 return(prim_dec(i,j)); 2641 } 2642 /////////////////////////////////////////////////////////////////////////////// 2643 2440 2644 proc radical(ideal i,list #) 2441 2645 { … … 2444 2648 if(size(#)>0) 2445 2649 { 2446 il=#[1]; 2650 il=#[1]; 2447 2651 } 2448 2652 ideal re=1; 2449 2653 option(redSB); 2654 list qr=simplifyIdeal(i); 2655 map phi=@P,qr[2]; 2656 i=qr[1]; 2657 2450 2658 list pr=facstd(i); 2451 2659 … … 2456 2664 { 2457 2665 ideal @res=maxideal(1); 2458 return( @res);2666 return(phi(@res)); 2459 2667 } 2460 2668 if(dim(pr[1])>1) 2461 2669 { 2462 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+"),lp;"; 2670 execute "ring gnir = ("+charstr(basering)+"), 2671 ("+varstr(basering)+"),(C,lp);"; 2463 2672 ideal i=fetch(@P,i); 2464 2673 list @pr=facstd(i); … … 2469 2678 option(noredSB); 2470 2679 int s=size(pr); 2680 2471 2681 if(s==1) 2472 2682 { 2473 return(radicalEHV(i,ideal(1),il)); 2683 i=radicalEHV(i,ideal(1),il); 2684 return(phi(i)); 2474 2685 } 2475 2686 intvec pos; 2476 2687 pos[s]=0; 2477 2478 2688 if(il==1) 2479 2689 { … … 2482 2692 int odim=dim(pr[1]); 2483 2693 int count=1; 2484 2694 2485 2695 for(j=2;j<=s;j++) 2486 2696 { … … 2502 2712 } 2503 2713 } 2504 2505 2714 for(j=1;j<=s;j++) 2506 2715 { 2507 2716 if(pos[j]==0) 2508 2717 { 2509 re=intersect 1(re,radicalEHV(pr[s+1-j],re,il));2510 } 2511 } 2512 return( re);2718 re=intersect(re,radicalEHV(pr[s+1-j],re,il)); 2719 } 2720 } 2721 return(phi(re)); 2513 2722 } 2514 2723 … … 2516 2725 { 2517 2726 def R=basering; 2518 execute "ring gnir = ("+charstr(basering)+"),("+varstr(basering)+",@t),dp;"; 2727 execute "ring gnir = ("+charstr(basering)+"), 2728 ("+varstr(basering)+",@t),(C,dp);"; 2519 2729 ideal i=var(nvars(basering))*imap(R,i)+(var(nvars(basering))-1)*imap(R,j); 2520 2730 ideal j=eliminate(i,var(nvars(basering))); … … 2523 2733 return(phi(j)); 2524 2734 } 2525 2735 2526 2736 2527 2737 /////////////////////////////////////////////////////////////////////////////// … … 2544 2754 return(ideal(0)); 2545 2755 } 2546 2756 2547 2757 def @P = basering; 2548 2758 list indep,allindep,restindep,fett,@mu; … … 2552 2762 ideal rad=ideal(1); 2553 2763 ideal te=ser; 2764 2554 2765 poly @p,@q; 2555 2766 string @va,quotring; … … 2563 2774 @j1=m[2]; 2564 2775 int jdim=dim(@j); 2565 if(size(reduce(ser,@j ))==0)2776 if(size(reduce(ser,@j,1))==0) 2566 2777 { 2567 2778 return(ser); … … 2595 2806 { 2596 2807 fac=factorize(@j[1],1); 2597 poly@p=1;2808 @p=1; 2598 2809 for(@k=1;@k<=size(fac);@k++) 2599 2810 { … … 2603 2814 2604 2815 return(ideal(@p)); 2605 } 2816 } 2606 2817 //------------------------------------------------------------------ 2607 2818 //the case of a complete intersection … … 2619 2830 if (jdim==0) 2620 2831 { 2621 @j1=system("finduni",@j); 2832 @j1=system("finduni",@j); 2622 2833 for(@k=1;@k<=size(@j1);@k++) 2623 2834 { … … 2634 2845 return(@j); 2635 2846 } 2636 2847 2637 2848 //------------------------------------------------------------------ 2638 2849 //search for a maximal independent set indep,i.e. … … 2643 2854 2644 2855 indep=maxIndependSet(@j); 2645 2856 2646 2857 di=dim(@j); 2647 2858 … … 2772 2983 { 2773 2984 fac=ideal(0); 2774 for(lauf=1;lauf<=ncols(@h);lauf++) 2985 for(lauf=1;lauf<=ncols(@h);lauf++) 2775 2986 { 2776 2987 if(deg(@h[lauf])>0) … … 2786 2997 } 2787 2998 2788 2999 2789 3000 @mu=mstd(quotient(@j+ideal(@q),rad)); 2790 3001 @j=@mu[1]; 2791 3002 attrib(@j,"isSB",1); 2792 3003 2793 3004 } 2794 3005 if((deg(rad[1])>0)&&(deg(collectrad[1])>0)) 2795 3006 { 2796 int xyz=timer;2797 "bei collecterad";2798 3007 rad=intersect(rad,collectrad); 2799 timer-xyz;2800 3008 } 2801 3009 else … … 2808 3016 2809 3017 te=simplify(reduce(te*rad,@j),2); 2810 3018 2811 3019 if((dim(@j)<di)||(size(te)==0)) 2812 3020 { … … 2822 3030 { 2823 3031 return(rad); 2824 } 2825 ideal tes=radicalKL(@mu,rad,@wr); 2826 int sml=timer; 2827 "bei rad"; 2828 rad=intersect(rad,tes); 2829 timer-sml; 2830 // rad=intersect(rad,radicalKL(@mu,ideal(1),@wr)); 3032 } 3033 // rad=intersect(rad,radicalKL(@mu,rad,@wr)); 3034 rad=intersect(rad,radicalKL(@mu,ideal(1),@wr)); 2831 3035 2832 3036 … … 2849 3053 il=#[1]; 2850 3054 } 3055 2851 3056 option(redSB); 2852 3057 list m=mstd(i); 2853 3058 I=m[2]; 2854 3059 option(noredSB); 2855 if(size(reduce(re,m[1] ))==0)3060 if(size(reduce(re,m[1],1))==0) 2856 3061 { 2857 3062 return(re); 2858 3063 } 2859 3064 int cod=nvars(basering)-dim(m[1]); 2860 if( nvars(basering)<9)3065 if((nvars(basering)<=5)&&(size(m[2])<=5)) 2861 3066 { 2862 3067 if(cod==size(m[2])) 2863 3068 { 2864 2865 2866 } 2867 3069 J=minor(jacob(I),cod); 3070 return(quotient(I,J)); 3071 } 3072 2868 3073 for(l=1;l<=cod;l++) 2869 3074 { … … 2877 3082 radI1=quotient(radI0,L); 2878 3083 2879 if(size(reduce(radI1,m[1] ))==0)3084 if(size(reduce(radI1,m[1],1))==0) 2880 3085 { 2881 3086 return(I); … … 2883 3088 if(il==1) 2884 3089 { 3090 2885 3091 return(radI1); 2886 3092 } … … 2892 3098 return(radI1); 2893 3099 } 2894 return(intersect(radI1,radicalEHV(I2,re,il))); 3100 return(intersect(radI1,radicalEHV(I2,re,il))); 2895 3101 } 2896 3102 } … … 2901 3107 { 2902 3108 M=prune(M); //to obtain a small embedding 2903 return(quotient(M,freemodule(nrows(M)))); 3109 ideal ann=quotient1(M,freemodule(nrows(M))); 3110 return(ann); 2904 3111 } 2905 3112 … … 2919 3126 } 2920 3127 return(ann); 2921 } 2922 3128 } 3129 2923 3130 //computes all equidimensional parts of the ideal i 2924 3131 proc prepareAss(ideal i) … … 2937 3144 { 2938 3145 list re=mres(i,0); //fehler in sres 2939 } 3146 } 2940 3147 for(e=cod;e<=nvars(basering);e++) 2941 3148 { 2942 3149 ann=AnnExt_R(e,re); 3150 2943 3151 if(nvars(basering)-dim(std(ann))==e) 2944 3152 { … … 2947 3155 } 2948 3156 return(er); 2949 } 3157 } 2950 3158 2951 3159 //computes the annihilator of Ext^n(R/i,R) with given resolution re … … 2958 3166 module k=res(f,2)[2]; //the kernel 2959 3167 matrix g=transpose(re[n]); //the image of Hom(_,R) 3168 2960 3169 ideal ann=quotient(g,k); //the anihilator 2961 3170 } … … 2964 3173 ideal ann=Ann(transpose(re[n])); 2965 3174 } 2966 return(ann); 2967 } 2968 3175 return(ann); 3176 } 3177 3178 proc quotient1(module a,module b) 3179 { 3180 int i; 3181 ideal re=quotient(a,module(b[1])); 3182 for(i=2;i<=size(b);i++) 3183 { 3184 re=intersect1(re,quotient(a,module(b[i]))); 3185 } 3186 return(re); 3187 } 3188 3189 3190 3191 proc analyze(list pr) 3192 { 3193 int ii,jj; 3194 for(ii=1;ii<=size(pr)/2;ii++) 3195 { 3196 dim(std(pr[2*ii])); 3197 idealsEqual(pr[2*ii-1],pr[2*ii]); 3198 "==========================="; 3199 } 3200 3201 for(ii=size(pr)/2;ii>1;ii--) 3202 { 3203 for(jj=1;jj<ii;jj++) 3204 { 3205 if(size(reduce(pr[2*jj],std(pr[2*ii],1)))==0) 3206 { 3207 "eingebette Komponente"; 3208 jj; 3209 ii; 3210 } 3211 } 3212 } 3213 } 3214 3215 3216 proc simplifyIdeal(ideal i) 3217 { 3218 def r=basering; 3219 3220 int j,k; 3221 map phi; 3222 poly p; 3223 3224 ideal iwork=i; 3225 ideal imap1=maxideal(1); 3226 ideal imap2=maxideal(1); 3227 3228 3229 for(j=1;j<=nvars(basering);j++) 3230 { 3231 for(k=1;k<=size(i);k++) 3232 { 3233 if(deg(iwork[k]/var(j))==0) 3234 { 3235 p=-1/leadcoef(iwork[k]/var(j))*iwork[k]; 3236 imap1[j]=p+2*var(j); 3237 phi=r,imap1; 3238 iwork=phi(iwork); 3239 iwork=subst(iwork,var(j),0); 3240 iwork[k]=var(j); 3241 imap1=maxideal(1); 3242 imap2[j]=-p; 3243 break; 3244 } 3245 } 3246 } 3247 return(iwork,imap2); 3248 } 3249 3250 3251 /////////////////////////////////////////////////////// 3252 // ini_mod 3253 // input: a polynomial p 3254 // output: the initial term of p as needed 3255 // in the context of characteristic sets 3256 ////////////////////////////////////////////////////// 3257 3258 proc ini_mod(poly p) 3259 { 3260 if (p==0) 3261 { 3262 return(0); 3263 } 3264 int n; matrix m; 3265 for( n=nvars(basering); n>0; n=n-1) 3266 { 3267 m=coef(p,var(n)); 3268 if(m[1,1]!=1) 3269 { 3270 p=m[2,1]; 3271 break; 3272 } 3273 } 3274 if(deg(p)==0) 3275 { 3276 p=0; 3277 } 3278 return(p); 3279 } 3280 /////////////////////////////////////////////////////// 3281 // min_ass_prim_charsets 3282 // input: generators of an ideal PS and an integer cho 3283 // If cho=0, the given ordering of the variables is used. 3284 // Otherwise, the system tries to find an "optimal ordering", 3285 // which in some cases may considerably speed up the algorithm 3286 // output: the minimal associated primes of PS 3287 // algorithm: via characteriostic sets 3288 ////////////////////////////////////////////////////// 3289 3290 3291 proc min_ass_prim_charsets (ideal PS, int cho) 3292 { 3293 if((cho<0) and (cho>1)) 3294 { 3295 "ERROR: <int> must be 0 or 1" 3296 return(); 3297 } 3298 if(system("version")>933) 3299 { 3300 option(notWarnSB); 3301 } 3302 if(cho==0) 3303 { 3304 return(min_ass_prim_charsets0(PS)); 3305 } 3306 else 3307 { 3308 return(min_ass_prim_charsets1(PS)); 3309 } 3310 } 3311 /////////////////////////////////////////////////////// 3312 // min_ass_prim_charsets0 3313 // input: generators of an ideal PS 3314 // output: the minimal associated primes of PS 3315 // algorithm: via characteristic sets 3316 // the given ordering of the variables is used 3317 ////////////////////////////////////////////////////// 3318 3319 3320 proc min_ass_prim_charsets0 (ideal PS) 3321 { 3322 3323 matrix m=char_series(PS); // We compute an irreducible 3324 // characteristic series 3325 int i,j,k; 3326 list PSI; 3327 list PHI; // the ideals given by the characteristic series 3328 for(i=nrows(m);i>=1; i--) 3329 { 3330 PHI[i]=ideal(m[i,1..ncols(m)]); 3331 } 3332 // We compute the radical of each ideal in PHI 3333 ideal I,JS,II; 3334 int sizeJS, sizeII; 3335 for(i=size(PHI);i>=1; i--) 3336 { 3337 I=0; 3338 for(j=size(PHI[i]);j>0;j=j-1) 3339 { 3340 I=I+ini_mod(PHI[i][j]); 3341 } 3342 JS=std(PHI[i]); 3343 sizeJS=size(JS); 3344 for(j=size(I);j>0;j=j-1) 3345 { 3346 II=0; 3347 sizeII=0; 3348 k=0; 3349 while(k<=sizeII) // successive saturation 3350 { 3351 option(returnSB); 3352 II=quotient(JS,I[j]); 3353 option(noreturnSB); 3354 //std 3355 // II=std(II); 3356 sizeII=size(II); 3357 if(sizeII==sizeJS) 3358 { 3359 for(k=1;k<=sizeII;k++) 3360 { 3361 if(leadexp(II[k])!=leadexp(JS[k])) break; 3362 } 3363 } 3364 JS=II; 3365 sizeJS=sizeII; 3366 } 3367 } 3368 PSI=insert(PSI,JS); 3369 } 3370 int sizePSI=size(PSI); 3371 // We eliminate redundant ideals 3372 for(i=1;i<sizePSI;i++) 3373 { 3374 for(j=i+1;j<=sizePSI;j++) 3375 { 3376 if(size(PSI[i])!=0) 3377 { 3378 if(size(PSI[j])!=0) 3379 { 3380 if(size(NF(PSI[i],PSI[j],1))==0) 3381 { 3382 PSI[j]=ideal(0); 3383 } 3384 else 3385 { 3386 if(size(NF(PSI[j],PSI[i],1))==0) 3387 { 3388 PSI[i]=ideal(0); 3389 } 3390 } 3391 } 3392 } 3393 } 3394 } 3395 for(i=sizePSI;i>=1;i--) 3396 { 3397 if(size(PSI[i])==0) 3398 { 3399 PSI=delete(PSI,i); 3400 } 3401 } 3402 return (PSI); 3403 } 3404 3405 /////////////////////////////////////////////////////// 3406 // min_ass_prim_charsets1 3407 // input: generators of an ideal PS 3408 // output: the minimal associated primes of PS 3409 // algorithm: via characteristic sets 3410 // input: generators of an ideal PS and an integer i 3411 // The system tries to find an "optimal ordering" of 3412 // the variables 3413 ////////////////////////////////////////////////////// 3414 3415 3416 proc min_ass_prim_charsets1 (ideal PS) 3417 { 3418 def oldring=basering; 3419 string n=system("neworder",PS); 3420 execute "ring r="+charstr(oldring)+",("+n+"),dp;"; 3421 ideal PS=imap(oldring,PS); 3422 matrix m=char_series(PS); // We compute an irreducible 3423 // characteristic series 3424 int i,j,k; 3425 ideal I; 3426 list PSI; 3427 list PHI; // the ideals given by the characteristic series 3428 list ITPHI; // their initial terms 3429 for(i=nrows(m);i>=1; i--) 3430 { 3431 PHI[i]=ideal(m[i,1..ncols(m)]); 3432 I=0; 3433 for(j=size(PHI[i]);j>0;j=j-1) 3434 { 3435 I=I,ini_mod(PHI[i][j]); 3436 } 3437 I=I[2..ncols(I)]; 3438 ITPHI[i]=I; 3439 } 3440 setring oldring; 3441 matrix m=imap(r,m); 3442 list PHI=imap(r,PHI); 3443 list ITPHI=imap(r,ITPHI); 3444 // We compute the radical of each ideal in PHI 3445 ideal I,JS,II; 3446 int sizeJS, sizeII; 3447 for(i=size(PHI);i>=1; i--) 3448 { 3449 I=0; 3450 for(j=size(PHI[i]);j>0;j=j-1) 3451 { 3452 I=I+ITPHI[i][j]; 3453 } 3454 JS=std(PHI[i]); 3455 sizeJS=size(JS); 3456 for(j=size(I);j>0;j=j-1) 3457 { 3458 II=0; 3459 sizeII=0; 3460 k=0; 3461 while(k<=sizeII) // successive iteration 3462 { 3463 option(returnSB); 3464 II=quotient(JS,I[j]); 3465 option(noreturnSB); 3466 //std 3467 // II=std(II); 3468 sizeII=size(II); 3469 if(sizeII==sizeJS) 3470 { 3471 for(k=1;k<=sizeII;k++) 3472 { 3473 if(leadexp(II[k])!=leadexp(JS[k])) break; 3474 } 3475 } 3476 JS=II; 3477 sizeJS=sizeII; 3478 } 3479 } 3480 PSI=insert(PSI,JS); 3481 } 3482 int sizePSI=size(PSI); 3483 // We eliminate redundant ideals 3484 for(i=1;i<sizePSI;i++) 3485 { 3486 for(j=i+1;j<=sizePSI;j++) 3487 { 3488 if(size(PSI[i])!=0) 3489 { 3490 if(size(PSI[j])!=0) 3491 { 3492 if(size(NF(PSI[i],PSI[j],1))==0) 3493 { 3494 PSI[j]=ideal(0); 3495 } 3496 else 3497 { 3498 if(size(NF(PSI[j],PSI[i],1))==0) 3499 { 3500 PSI[i]=ideal(0); 3501 } 3502 } 3503 } 3504 } 3505 } 3506 } 3507 for(i=sizePSI;i>=1;i--) 3508 { 3509 if(size(PSI[i])==0) 3510 { 3511 PSI=delete(PSI,i); 3512 } 3513 } 3514 return (PSI); 3515 } 3516 3517 3518 ///////////////////////////////////////////////////// 3519 // proc prim_dec 3520 // input: generators of an ideal I and an integer choose 3521 // If choose=0, min_ass_prim_charsets with the given 3522 // ordering of the variables is used. 3523 // If choose=1, min_ass_prim_charsets with the "optimized" 3524 // ordering of the variables is used. 3525 // If choose=2, minAssPrimes from primdec.lib is used 3526 // If choose=3, minAssPrimes+factorizing Buchberger from primdec.lib is used 3527 // output: a primary decomposition of I, i.e., a list 3528 // of pairs consisting of a standard basis of a primary component 3529 // of I and a standard basis of the corresponding associated prime. 3530 // To compute the minimal associated primes of a given ideal 3531 // min_ass_prim_l is called, i.e., the minimal associated primes 3532 // are computed via characteristic sets. 3533 // In the homogeneous case, the performance of the procedure 3534 // will be improved if I is already given by a minimal set of 3535 // generators. Apply minbase if necessary. 3536 ////////////////////////////////////////////////////////// 3537 3538 3539 proc prim_dec(ideal I, int choose) 3540 { 3541 if((choose<0) or (choose>3)) 3542 { 3543 "ERROR: <int> must be 0 or 1 or 2 or 3"; 3544 return(); 3545 } 3546 if(system("version")>933) 3547 { 3548 option(notWarnSB); 3549 } 3550 ideal H=1; // The intersection of the primary components 3551 list U; // the leaves of the decomposition tree, i.e., 3552 // pairs consisting of a primary component of I 3553 // and the corresponding associated prime 3554 list W; // the non-leaf vertices in the decomposition tree. 3555 // every entry has 6 components: 3556 // 1- the vertex itself , i.e., a standard bais of the 3557 // given ideal I (type 1), or a standard basis of a 3558 // pseudo-primary component arising from 3559 // pseudo-primary decomposition (type 2), or a 3560 // standard basis of a remaining component arising from 3561 // pseudo-primary decomposition or extraction (type 3) 3562 // 2- the type of the vertex as indicated above 3563 // 3- the weighted_tree_depth of the vertex 3564 // 4- the tester of the vertex 3565 // 5- a standard basis of the associated prime 3566 // of a vertex of type 2, or 0 otherwise 3567 // 6- a list of pairs consisting of a standard 3568 // basis of a minimal associated prime ideal 3569 // of the father of the vertex and the 3570 // irreducible factors of the "minimal 3571 // divisor" of the seperator or extractor 3572 // corresponding to the prime ideal 3573 // as computed by the procedure minsat, 3574 // if the vertex is of type 3, or 3575 // the empty list otherwise 3576 ideal SI=std(I); 3577 int ncolsSI=ncols(SI); 3578 int ncolsH=1; 3579 W[1]=list(I,1,0,poly(1),ideal(0),list()); // The root of the tree 3580 int weighted_tree_depth; 3581 int i,j; 3582 int check; 3583 list V; // current vertex 3584 list VV; // new vertex 3585 list QQ; 3586 list WI; 3587 ideal Qi,SQ,SRest,fac; 3588 poly tester; 3589 3590 while(1) 3591 { 3592 i=1; 3593 while(1) 3594 { 3595 while(i<=size(W)) // find vertex V of smallest weighted tree-depth 3596 { 3597 if (W[i][3]<=weighted_tree_depth) break; 3598 i++; 3599 } 3600 if (i<=size(W)) break; 3601 i=1; 3602 weighted_tree_depth++; 3603 } 3604 V=W[i]; 3605 W=delete(W,i); // delete V from W 3606 3607 // now proceed by type of vertex V 3608 3609 if (V[2]==2) // extraction needed 3610 { 3611 SQ,SRest,fac=extraction(V[1],V[5]); 3612 // standard basis of primary component, 3613 // standard basis of remaining component, 3614 // irreducible factors of 3615 // the "minimal divisor" of the extractor 3616 // as computed by the procedure minsat, 3617 check=0; 3618 for(j=1;j<=ncolsH;j++) 3619 { 3620 if (NF(H[j],SQ,1)!=0) // Q is not redundant 3621 { 3622 check=1; 3623 break; 3624 } 3625 } 3626 if(check==1) // Q is not redundant 3627 { 3628 QQ=list(); 3629 QQ[1]=list(SQ,V[5]); // primary component, associated prime, 3630 // i.e., standard bases thereof 3631 U=U+QQ; 3632 H=intersect(H,SQ); 3633 H=std(H); 3634 ncolsH=ncols(H); 3635 check=0; 3636 if(ncolsH==ncolsSI) 3637 { 3638 for(j=1;j<=ncolsSI;j++) 3639 { 3640 if(leadexp(H[j])!=leadexp(SI[j])) 3641 { 3642 check=1; 3643 break; 3644 } 3645 } 3646 } 3647 else 3648 { 3649 check=1; 3650 } 3651 if(check==0) // H==I => U is a primary decomposition 3652 { 3653 return(U); 3654 } 3655 } 3656 if (SRest[1]!=1) // the remaining component is not 3657 // the whole ring 3658 { 3659 if (rad_con(V[4],SRest)==0) // the new vertex is not the 3660 // root of a redundant subtree 3661 { 3662 VV[1]=SRest; // remaining component 3663 VV[2]=3; // pseudoprimdec_special 3664 VV[3]=V[3]+1; // weighted depth 3665 VV[4]=V[4]; // the tester did not change 3666 VV[5]=ideal(0); 3667 VV[6]=list(list(V[5],fac)); 3668 W=insert(W,VV,size(W)); 3669 } 3670 } 3671 } 3672 else 3673 { 3674 if (V[2]==3) // pseudo_prim_dec_special is needed 3675 { 3676 QQ,SRest=pseudo_prim_dec_special_charsets(V[1],V[6],choose); 3677 // QQ = quadruples: 3678 // standard basis of pseudo-primary component, 3679 // standard basis of corresponding prime, 3680 // seperator, irreducible factors of 3681 // the "minimal divisor" of the seperator 3682 // as computed by the procedure minsat, 3683 // SRest=standard basis of remaining component 3684 } 3685 else // V is the root, pseudo_prim_dec is needed 3686 { 3687 QQ,SRest=pseudo_prim_dec_charsets(I,SI,choose); 3688 // QQ = quadruples: 3689 // standard basis of pseudo-primary component, 3690 // standard basis of corresponding prime, 3691 // seperator, irreducible factors of 3692 // the "minimal divisor" of the seperator 3693 // as computed by the procedure minsat, 3694 // SRest=standard basis of remaining component 3695 3696 } 3697 //check 3698 for(i=size(QQ);i>=1;i--) 3699 //for(i=1;i<=size(QQ);i++) 3700 { 3701 tester=QQ[i][3]*V[4]; 3702 Qi=QQ[i][2]; 3703 if(NF(tester,Qi,1)!=0) // the new vertex is not the 3704 // root of a redundant subtree 3705 { 3706 VV[1]=QQ[i][1]; 3707 VV[2]=2; 3708 VV[3]=V[3]+1; 3709 VV[4]=tester; // the new tester as computed above 3710 VV[5]=Qi; // QQ[i][2]; 3711 VV[6]=list(); 3712 W=insert(W,VV,size(W)); 3713 } 3714 } 3715 if (SRest[1]!=1) // the remaining component is not 3716 // the whole ring 3717 { 3718 if (rad_con(V[4],SRest)==0) // the vertex is not the root 3719 // of a redundant subtree 3720 { 3721 VV[1]=SRest; 3722 VV[2]=3; 3723 VV[3]=V[3]+2; 3724 VV[4]=V[4]; // the tester did not change 3725 VV[5]=ideal(0); 3726 WI=list(); 3727 for(i=1;i<=size(QQ);i++) 3728 { 3729 WI=insert(WI,list(QQ[i][2],QQ[i][4])); 3730 } 3731 VV[6]=WI; 3732 W=insert(W,VV,size(W)); 3733 } 3734 } 3735 } 3736 } 3737 } 3738 3739 ////////////////////////////////////////////////////////////////////////// 3740 // proc pseudo_prim_dec_charsets 3741 // input: Generators of an arbitrary ideal I, a standard basis SI of I, 3742 // and an integer choo 3743 // If choo=0, min_ass_prim_charsets with the given 3744 // ordering of the variables is used. 3745 // If choo=1, min_ass_prim_charsets with the "optimized" 3746 // ordering of the variables is used. 3747 // If choo=2, minAssPrimes from primdec.lib is used 3748 // If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used 3749 // output: a pseudo primary decomposition of I, i.e., a list 3750 // of pseudo primary components together with a standard basis of the 3751 // remaining component. Each pseudo primary component is 3752 // represented by a quadrupel: A standard basis of the component, 3753 // a standard basis of the corresponding associated prime, the 3754 // seperator of the component, and the irreducible factors of the 3755 // "minimal divisor" of the seperator as computed by the procedure minsat, 3756 // calls proc pseudo_prim_dec_i 3757 ////////////////////////////////////////////////////////////////////////// 3758 3759 3760 proc pseudo_prim_dec_charsets (ideal I, ideal SI, int choo) 3761 { 3762 list L; // The list of minimal associated primes, 3763 // each one given by a standard basis 3764 if((choo==0) or (choo==1)) 3765 { 3766 L=min_ass_prim_charsets(I,choo); 3767 } 3768 else 3769 { 3770 if(choo==2) 3771 { 3772 L=minAssPrimes(I); 3773 } 3774 else 3775 { 3776 L=minAssPrimes(I,1); 3777 } 3778 for(int i=size(L);i>=1;i=i-1) 3779 { 3780 L[i]=std(L[i]); 3781 } 3782 } 3783 return (pseudo_prim_dec_i(SI,L)); 3784 } 3785 3786 //////////////////////////////////////////////////////////////// 3787 // proc pseudo_prim_dec_special_charsets 3788 // input: a standard basis of an ideal I whose radical is the 3789 // intersection of the radicals of ideals generated by one prime ideal 3790 // P_i together with one polynomial f_i, the list V6 must be the list of 3791 // pairs (standard basis of P_i, irreducible factors of f_i), 3792 // and an integer choo 3793 // If choo=0, min_ass_prim_charsets with the given 3794 // ordering of the variables is used. 3795 // If choo=1, min_ass_prim_charsets with the "optimized" 3796 // ordering of the variables is used. 3797 // If choo=2, minAssPrimes from primdec.lib is used 3798 // If choo=3, minAssPrimes+factorizing Buchberger from primdec.lib is used 3799 // output: a pseudo primary decomposition of I, i.e., a list 3800 // of pseudo primary components together with a standard basis of the 3801 // remaining component. Each pseudo primary component is 3802 // represented by a quadrupel: A standard basis of the component, 3803 // a standard basis of the corresponding associated prime, the 3804 // seperator of the component, and the irreducible factors of the 3805 // "minimal divisor" of the seperator as computed by the procedure minsat, 3806 // calls proc pseudo_prim_dec_i 3807 //////////////////////////////////////////////////////////////// 3808 3809 3810 proc pseudo_prim_dec_special_charsets (ideal SI,list V6, int choo) 3811 { 3812 int i,j,l; 3813 list m; 3814 list L; 3815 int sizeL; 3816 ideal P,SP; ideal fac; 3817 int dimSP; 3818 for(l=size(V6);l>=1;l--) // creates a list of associated primes 3819 // of I, possibly redundant 3820 { 3821 P=V6[l][1]; 3822 fac=V6[l][2]; 3823 for(i=ncols(fac);i>=1;i--) 3824 { 3825 SP=P+fac[i]; 3826 SP=std(SP); 3827 if(SP[1]!=1) 3828 { 3829 if((choo==0) or (choo==1)) 3830 { 3831 m=min_ass_prim_charsets(SP,choo); // a list of SB 3832 } 3833 else 3834 { 3835 if(choo==2) 3836 { 3837 m=minAssPrimes(SP); 3838 } 3839 else 3840 { 3841 m=minAssPrimes(SP,1); 3842 } 3843 for(j=size(m);j>=1;j=j-1) 3844 { 3845 m[j]=std(m[j]); 3846 } 3847 } 3848 dimSP=dim(SP); 3849 for(j=size(m);j>=1; j--) 3850 { 3851 if(dim(m[j])==dimSP) 3852 { 3853 L=insert(L,m[j],size(L)); 3854 } 3855 } 3856 } 3857 } 3858 } 3859 sizeL=size(L); 3860 for(i=1;i<sizeL;i++) // get rid of redundant primes 3861 { 3862 for(j=i+1;j<=sizeL;j++) 3863 { 3864 if(size(L[i])!=0) 3865 { 3866 if(size(L[j])!=0) 3867 { 3868 if(size(NF(L[i],L[j],1))==0) 3869 { 3870 L[j]=ideal(0); 3871 } 3872 else 3873 { 3874 if(size(NF(L[j],L[i],1))==0) 3875 { 3876 L[i]=ideal(0); 3877 } 3878 } 3879 } 3880 } 3881 } 3882 } 3883 for(i=sizeL;i>=1;i--) 3884 { 3885 if(size(L[i])==0) 3886 { 3887 L=delete(L,i); 3888 } 3889 } 3890 return (pseudo_prim_dec_i(SI,L)); 3891 } 3892 3893 3894 //////////////////////////////////////////////////////////////// 3895 // proc pseudo_prim_dec_i 3896 // input: A standard basis of an arbitrary ideal I, and standard bases 3897 // of the minimal associated primes of I 3898 // output: a pseudo primary decomposition of I, i.e., a list 3899 // of pseudo primary components together with a standard basis of the 3900 // remaining component. Each pseudo primary component is 3901 // represented by a quadrupel: A standard basis of the component Q_i, 3902 // a standard basis of the corresponding associated prime P_i, the 3903 // seperator of the component, and the irreducible factors of the 3904 // "minimal divisor" of the seperator as computed by the procedure minsat, 3905 //////////////////////////////////////////////////////////////// 3906 3907 3908 proc pseudo_prim_dec_i (ideal SI, list L) 3909 { 3910 list Q; 3911 if (size(L)==1) // one minimal associated prime only 3912 // the ideal is already pseudo primary 3913 { 3914 Q=SI,L[1],1; 3915 list QQ; 3916 QQ[1]=Q; 3917 return (QQ,ideal(1)); 3918 } 3919 3920 poly f0,f,g; 3921 ideal fac; 3922 int i,j,k,l; 3923 ideal SQi; 3924 ideal I'=SI; 3925 list QP; 3926 int sizeL=size(L); 3927 for(i=1;i<=sizeL;i++) 3928 { 3929 fac=0; 3930 for(j=1;j<=sizeL;j++) // compute the seperator sep_i 3931 // of the i-th component 3932 { 3933 if (i!=j) // search g not in L[i], but L[j] 3934 { 3935 for(k=1;k<=ncols(L[j]);k++) 3936 { 3937 if(NF(L[j][k],L[i],1)!=0) 3938 { 3939 break; 3940 } 3941 } 3942 fac=fac+L[j][k]; 3943 } 3944 } 3945 // delete superfluous polynomials 3946 fac=simplify(fac,8); 3947 // saturation 3948 SQi,f0,f,fac=minsat_ppd(SI,fac); 3949 I'=I',f; 3950 QP=SQi,L[i],f0,fac; 3951 // the quadrupel: 3952 // a standard basis of Q_i, 3953 // a standard basis of P_i, 3954 // sep_i, 3955 // irreducible factors of 3956 // the "minimal divisor" of the seperator 3957 // as computed by the procedure minsat, 3958 Q[i]=QP; 3959 } 3960 I'=std(I'); 3961 return (Q, I'); 3962 // I' = remaining component 3963 } 3964 3965 3966 //////////////////////////////////////////////////////////////// 3967 // proc extraction 3968 // input: A standard basis of a pseudo primary ideal I, and a standard 3969 // basis of the unique minimal associated prime P of I 3970 // output: an extraction of I, i.e., a standard basis of the primary 3971 // component Q of I with associated prime P, a standard basis of the 3972 // remaining component, and the irreducible factors of the 3973 // "minimal divisor" of the extractor as computed by the procedure minsat 3974 //////////////////////////////////////////////////////////////// 3975 3976 3977 proc extraction (ideal SI, ideal SP) 3978 { 3979 list indsets=system("indsetall",SP,0); 3980 poly f; 3981 if(size(indsets)!=0) //check, whether dim P != 0 3982 { 3983 intvec v; // a maximal independent set of variables 3984 // modulo P 3985 string U; // the independent variables 3986 string A; // the dependent variables 3987 int j,k; 3988 int a; // the size of A 3989 int degf; 3990 ideal g; 3991 list polys; 3992 int sizepolys; 3993 list newpoly; 3994 def R=basering; 3995 //intvec hv=hilb(SI,1); 3996 for (k=1;k<=size(indsets);k++) 3997 { 3998 v=indsets[k]; 3999 for (j=1;j<=nvars(R);j++) 4000 { 4001 if (v[j]==1) 4002 { 4003 U=U+varstr(j)+","; 4004 } 4005 else 4006 { 4007 A=A+varstr(j)+","; 4008 a++; 4009 } 4010 } 4011 4012 U[size(U)]=")"; // we compute the extractor of I (w.r.t. U) 4013 execute "ring RAU="+charstr(basering)+",("+A+U+",(dp("+string(a)+"),dp);"; 4014 ideal I=imap(R,SI); 4015 //I=std(I,hv); // the standard basis in (R[U])[A] 4016 I=std(I); // the standard basis in (R[U])[A] 4017 A[size(A)]=")"; 4018 execute "ring Rloc=("+charstr(basering)+","+U+",("+A+",dp;"; 4019 ideal I=imap(RAU,I); 4020 //"std in lokalisierung:"+newline,I; 4021 ideal h; 4022 for(j=ncols(I);j>=1;j--) 4023 { 4024 h[j]=leadcoef(I[j]); // consider I in (R(U))[A] 4025 } 4026 setring R; 4027 g=imap(Rloc,h); 4028 kill RAU,Rloc; 4029 U=""; 4030 A=""; 4031 a=0; 4032 f=lcm(g); 4033 newpoly[1]=f; 4034 polys=polys+newpoly; 4035 newpoly=list(); 4036 } 4037 f=polys[1]; 4038 degf=deg(f); 4039 sizepolys=size(polys); 4040 for (k=2;k<=sizepolys;k++) 4041 { 4042 if (deg(polys[k])<degf) 4043 { 4044 f=polys[k]; 4045 degf=deg(f); 4046 } 4047 } 4048 } 4049 else 4050 { 4051 f=1; 4052 } 4053 poly f0,h0; ideal SQ; ideal fac; 4054 if(f!=1) 4055 { 4056 SQ,f0,h0,fac=minsat(SI,f); 4057 return(SQ,std(SI+h0),fac); 4058 // the tripel 4059 // a standard basis of Q, 4060 // a standard basis of remaining component, 4061 // irreducible factors of 4062 // the "minimal divisor" of the extractor 4063 // as computed by the procedure minsat 4064 } 4065 else 4066 { 4067 return(SI,ideal(1),ideal(1)); 4068 } 4069 } 4070 4071 ///////////////////////////////////////////////////// 4072 // proc minsat 4073 // input: a standard basis of an ideal I and a polynomial p 4074 // output: a standard basis IS of the saturation of I w.r. to p, 4075 // the maximal squarefree factor f0 of p, 4076 // the "minimal divisor" f of f0 such that the saturation of 4077 // I w.r. to f equals the saturation of I w.r. to f0 (which is IS), 4078 // the irreducible factors of f 4079 ////////////////////////////////////////////////////////// 4080 4081 4082 proc minsat(ideal SI, poly p) 4083 { 4084 ideal fac=factorize(p,1); //the irreducible factors of p 4085 fac=sort(fac)[1]; 4086 int i,k; 4087 poly f0=1; 4088 for(i=ncols(fac);i>=1;i--) 4089 { 4090 f0=f0*fac[i]; 4091 } 4092 poly f=1; 4093 ideal iold; 4094 list quotM; 4095 quotM[1]=SI; 4096 quotM[2]=fac; 4097 quotM[3]=f0; 4098 // we deal seperately with the first quotient; 4099 // factors, which do not contribute to this one, 4100 // are omitted 4101 iold=quotM[1]; 4102 quotM=minquot(quotM); 4103 fac=quotM[2]; 4104 if(quotM[3]==1) 4105 { 4106 return(quotM[1],f0,f,fac); 4107 } 4108 while(special_ideals_equal(iold,quotM[1])==0) 4109 { 4110 f=f*quotM[3]; 4111 iold=quotM[1]; 4112 quotM=minquot(quotM); 4113 } 4114 return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f) 4115 } 4116 4117 ///////////////////////////////////////////////////// 4118 // proc minsat_ppd 4119 // input: a standard basis of an ideal I and a polynomial p 4120 // output: a standard basis IS of the saturation of I w.r. to p, 4121 // the maximal squarefree factor f0 of p, 4122 // the "minimal divisor" f of f0 such that the saturation of 4123 // I w.r. to f equals the saturation of I w.r. to f0 (which is IS), 4124 // the irreducible factors of f 4125 ////////////////////////////////////////////////////////// 4126 4127 4128 proc minsat_ppd(ideal SI, ideal fac) 4129 { 4130 fac=sort(fac)[1]; 4131 int i,k; 4132 poly f0=1; 4133 for(i=ncols(fac);i>=1;i--) 4134 { 4135 f0=f0*fac[i]; 4136 } 4137 poly f=1; 4138 ideal iold; 4139 list quotM; 4140 quotM[1]=SI; 4141 quotM[2]=fac; 4142 quotM[3]=f0; 4143 // we deal seperately with the first quotient; 4144 // factors, which do not contribute to this one, 4145 // are omitted 4146 iold=quotM[1]; 4147 quotM=minquot(quotM); 4148 fac=quotM[2]; 4149 if(quotM[3]==1) 4150 { 4151 return(quotM[1],f0,f,fac); 4152 } 4153 while(special_ideals_equal(iold,quotM[1])==0) 4154 { 4155 f=f*quotM[3]; 4156 iold=quotM[1]; 4157 quotM=minquot(quotM); 4158 k++; 4159 } 4160 return(quotM[1],f0,f,fac); // the quadrupel ((I:p),f0,f, irr. factors of f) 4161 } 4162 ///////////////////////////////////////////////////////////////// 4163 // proc minquot 4164 // input: a list with 3 components: a standard basis 4165 // of an ideal I, a set of irreducible polynomials, and 4166 // there product f0 4167 // output: a standard basis of the ideal (I:f0), the irreducible 4168 // factors of the "minimal divisor" f of f0 with (I:f0) = (I:f), 4169 // the "minimal divisor" f 4170 ///////////////////////////////////////////////////////////////// 4171 4172 proc minquot(list tsil) 4173 { 4174 int i,j,k,action; 4175 ideal verg; 4176 list l; 4177 poly g; 4178 ideal laedi=tsil[1]; 4179 ideal fac=tsil[2]; 4180 poly f=tsil[3]; 4181 4182 //std 4183 // ideal star=quotient(laedi,f); 4184 // star=std(star); 4185 option(returnSB); 4186 ideal star=quotient(laedi,f); 4187 option(noreturnSB); 4188 if(special_ideals_equal(laedi,star)==1) 4189 { 4190 return(laedi,ideal(1),1); 4191 } 4192 action=1; 4193 while(action==1) 4194 { 4195 if(size(fac)==1) 4196 { 4197 action=0; 4198 break; 4199 } 4200 for(i=1;i<=size(fac);i++) 4201 { 4202 g=1; 4203 for(j=1;j<=size(fac);j++) 4204 { 4205 if(i!=j) 4206 { 4207 g=g*fac[j]; 4208 } 4209 } 4210 //std 4211 // verg=quotient(laedi,g); 4212 // verg=std(verg); 4213 option(returnSB); 4214 verg=quotient(laedi,g); 4215 option(noreturnSB); 4216 if(special_ideals_equal(verg,star)==1) 4217 { 4218 f=g; 4219 fac[i]=0; 4220 fac=simplify(fac,2); 4221 break; 4222 } 4223 if(i==size(fac)) 4224 { 4225 action=0; 4226 } 4227 } 4228 } 4229 l=star,fac,f; 4230 return(l); 4231 } 4232 ///////////////////////////////////////////////// 4233 // proc special_ideals_equal 4234 // input: standard bases of ideal k1 and k2 such that 4235 // k1 is contained in k2, or k2 is contained ink1 4236 // output: 1, if k1 equals k2, 0 otherwise 4237 ////////////////////////////////////////////////// 4238 4239 proc special_ideals_equal( ideal k1, ideal k2) 4240 { 4241 int j; 4242 if(size(k1)==size(k2)) 4243 { 4244 for(j=1;j<=size(k1);j++) 4245 { 4246 if(leadexp(k1[j])!=leadexp(k2[j])) 4247 { 4248 return(0); 4249 } 4250 } 4251 return(1); 4252 } 4253 return(0); 4254 } 4255 4256 4257
Note: See TracChangeset
for help on using the changeset viewer.
