Changeset e2fd5d in git
 Timestamp:
 Dec 24, 2010, 9:49:04 PM (12 years ago)
 Branches:
 (u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
 Children:
 b90cc8b474180d8d55d1ffc520a7e8ae0a273c30
 Parents:
 9e269a0a45fc3fa26da7b381bc6593a5b4b7df4d
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/sagbi.lib
r9e269a0 re2fd5d 11 11 SAGBI stands for 'subalgebra bases analogous to Groebner bases for ideals'. 12 12 SAGBI bases provide important tools for working with finitely presented 13 subalgebras of a polynomial ring. Note that in contrast to Groebner13 subalgebras of a polynomial ring. Note, that in contrast to Groebner 14 14 bases, SAGBI bases may be infinite. 15 15 16 16 REFERENCES: 17 Ana Bravo: Some Facts About Canonical Subalgebra Bases, 17 Ana Bravo: Some Facts About Canonical Subalgebra Bases, 18 18 MSRI Publications 51, p. 247254 19 19 … … 136 136 * 5. Compute the new algebraic relations 137 137 */ 138 int ppl = printlevelvoice+3; //variable for additional printleveldependend information 139 dbprint(ppl,"//Spoly1 initialisation and precomputation"); 138 140 def br=basering; 139 141 ideal varsBasering=maxideal(1); … … 166 168 kernNew[invars(r)+nvars(br)]=leadTermsAlgebra[i]listOfVariables[i+nvars(br)]; 167 169 } 168 // cal ulate kernel of Phi depending on method choosen 169 170 // calculate kernel of Phi depending on method choosen  171 dbprint(ppl,"//Spoly2 Groebner basis computation"); 170 172 attrib(kernOld,"isSB",1); 171 173 ideal kern=stdKernPhi(kernNew,kernOld,leadTermsAlgebra,method); 174 dbprint(ppl2,"//Spoly21 ideal kern",kern); 172 175 // calulate algebraic relations  176 dbprint(ppl,"//Spoly3 computing new algebraic relations"); 173 177 ideal algebraicRelations=nselect(kern,1..nvars(br)); 174 178 attrib(algebraicRelationsOld,"isSB",1); 175 179 ideal algebraicRelationsNew=reduce(algebraicRelations,algebraicRelationsOld); 180 /* canonicalizing: */ 181 algebraicRelationsNew=canonicalform(algebraicRelationsNew); 182 dbprint(ppl2,"//Spoly31 ideal of new algebraic relations",algebraicRelationsNew); 176 183 /* algebraicRelationsOld is a groebner basis by construction (as variable 177 184 * ordering is … … 193 200 * K[y_1,...,y_m]>K[x_1,...,x_n], x(i)>x(i), @y(i)>leadterm(algebra[i]) 194 201 */ 202 int ppl = printlevelvoice+3; //variable for additional printleveldependend information 203 dbprint(ppl,"//Spoly1 initialisation and precomputation"); 195 204 def br=basering; 196 205 int m=ncols(algebra); … … 203 212 leadCoefficients[i]=leadcoef(algebra[i]); 204 213 } 214 dbprint(ppl2,"//Spoly11 Vector of leading coefficients",leadCoefficients); 205 215 int k=1; 206 216 for (i=1;i<=n;i++) … … 215 225 //of the leadings monomials in algebra. 216 226 intmat A[n][m]=intmat(tempVec,n,m); 227 dbprint(ppl2,"//Spoly12 Matrix A",A); 217 228 //Create the preimage ring K[@y(1),...,@y(m)], where m=ncols(algebra). 218 229 string variableName=uniqueVariableName("@y"); … … 227 238 setring rNew; 228 239 //Use toric_ideal to compute the kernel 240 dbprint(ppl,"//Spoly2 call of toric_ideal"); 229 241 ideal algebraicRelations=toric_ideal(A,"ect"); 230 242 //Suitable substitution 243 dbprint(ppl,"//Spoly3 substitutions"); 231 244 ideal leadCoefficients=fetch(br,leadCoefficients); 232 245 for (i=1; i<=m; i++) … … 237 250 } 238 251 } 252 dbprint(ppl2,"//Spoly31 algebraic relations",algebraicRelations); 239 253 export algebraicRelations; 240 254 return(rNew); … … 251 265 * p=reduce(leadF,kern) in K[@y(1),...,@y(m)] <=> leadF in K[lead(algebra)] 252 266 * Otherwise some precomputation has to be done, outlined below. 267 * When using sagbiReduce,sagbi and sagbiPart the integer parRed will always be zero. Only the procedure 268 * algebraicDependence causes this procedure to be called with parRed<>0. The only difference when parRed<>0 269 * is that the reduction algorithms returns the nonzero constants it attains (instead of just returning zero as the 270 * correct remainder), as they will be expressions in parameters for an algebraic dependence. 253 271 */ 272 int ppl = printlevelvoice+3; //variable for additional printleveldependend information 273 dbprint(ppl,"//Red1 initialisation and precomputation"); 254 274 def br=basering; 255 275 int numVarsBasering=nvars(br); … … 259 279 if (numVarsBasering==nvars(r)) 260 280 { 281 dbprint(ppl1,"//Red11 Groebner basis computation"); 261 282 /* Case that ring r is the same ring as the basering. Using proc extendRing, 262 283 * stdKernPhi … … 281 302 } 282 303 kern=stdKernPhi(kern,0,leadTermsAlgebra,method); 304 dbprint(ppl2,"//Red111 Ideal kern",kern); 283 305 } 284 306 setring r; … … 290 312 intvec tempExp; 291 313 //algebraic reduction for each polynomial F[i]  314 dbprint(ppl,"//Red2 reduction, polynomial by polynomial"); 292 315 for (i=1; i<=ncols(F);i++) 293 316 { 317 dbprint(ppl1,"//Red2"+string(i)+" starting with new polynomial"); 318 dbprint(ppl2,"//Red2"+string(i)+"1 Polynomial before reduction",F[i]); 294 319 normalform=0; 295 320 while (F[i]!=0) … … 300 325 //K is always contained in the subalgebra, 301 326 //thus the remainder is zero in this case 302 if (parRed) { break; } 303 else { F[i]=0; break; } 327 if (parRed) 328 { 329 //If parRed<>0 save nonzero constants the reduction algorithms attains. 330 break; 331 } 332 else 333 { 334 F[i]=0; 335 break; 336 } 304 337 } 305 338 //note: as the ordering in br and r might not be compatible … … 338 371 F[i] = normalform; 339 372 } 373 dbprint(ppl2,"//Red2"+string(i)+"2 Polynomial after reduction",F[i]); 340 374 } 341 375 return(F); … … 389 423 * and iterations specifies the 390 424 * number of iterations. A degree boundary is not used here. 391 * Note that parRed is used for testing a special modification 392 * and can be ignored (assume parRed==0). 425 * When this method is called via the procedures sagbi and sagbiPart the integer parRed 426 * will always be zero. Only the procedure algebraicDependence calls this procedure with 427 * parRed<>0. The only difference when parRed<>0 is that the reduction algorithms returns 428 * the nonzero constants it attains (instead of just returning zero as the correct 429 * remainder), as they will be expressions in parameters for an algebraic dependence. 430 * These constants are saved in the ideal reducedParameters. 393 431 */ 394 432 { 433 int ppl = printlevelvoice+3; //variable for additional printleveldependend information 434 dbprint(ppl,"// 0 initialisation and precomputation"); 395 435 def br=basering; 396 436 int i=1; 397 437 398 438 ideal reducedParameters; 399 int numReducedParameters=1; 439 int numReducedParameters=1; //number of elements plus one in reducedParameters 400 440 int j; 401 if (parRed==0) 441 if (parRed==0) //if parRed<>0 the algebra does not contain monomials and normalisation should be avoided 402 442 { 403 443 algebra=reduceByMonomials(algebra); 404 444 algebra=simplify(simplify(algebra,3),4); 405 445 } 406 int step=1; 407 if (iterations==1) //case: infintitly many iterations 408 { 409 step=0; 410 iterations=1; 411 } 446 // canonicalizing the gen's: 447 algebra = canonicalform(algebra); 412 448 ideal P=1; 413 449 //note: P is initialized this way, so that the while loop is entered. … … 417 453 ideal spolynomialsNew; 418 454 def r=br; 419 while (size(P)>0 && i<=iterations) 420 { 421 def rNew=spolynomialsGB(algebra,r,method); 455 while (size(P)>0) 456 { 457 dbprint(ppl,"// "+string(i)+" interation of SAGBI construction algorithm"); 458 dbprint(ppl1,"// "+string(i)+"1 Computing algebraic relations"); 459 def rNew=spolynomialsGB(algebra,r,method); /* canonicalizing inside! */ 422 460 kill r; 423 461 def r=rNew; 424 462 kill rNew; 425 463 phi=r,varsBasering,algebra; 464 dbprint(ppl1,"// "+string(i)+"2 Substituting into algebraic relations"); 426 465 spolynomialsNew=simplify(phi(algebraicRelationsNew),6); 427 466 //By construction spolynomialsNew only contains the spolynomials, 428 467 //that have not already 429 468 //been calculated in the steps before. 469 dbprint(ppl1,"// "+string(i)+"3 SAGBI reduction"); 470 dbprint(ppl2,"// "+string(i)+"31 new Spolynomials before reduction",spolynomialsNew); 430 471 P=reductionGB(spolynomialsNew,algebra,r,tailreduction,method,parRed); 431 472 if (parRed) … … 452 493 P=simplify(P,6); 453 494 } 495 /* canonicalize ! */ 496 P = canonicalform(P); 497 dbprint(ppl2,"// "+string(i)+"31 new Spolynomials after reduction",P); 454 498 algebra=algebra,P; 455 499 //Note that elements and order of elements must in algebra must not be changed, 456 500 //otherwise the already calculated 457 501 //ideal in r will give wrong results. Thus it is important to use a komma here. 458 i=i+step; 459 } 460 if (step==1) 502 i=i+1; 503 if (iterations!=1 && i>iterations) //When iterations==1 the number of iterations is unlimited 504 { 505 break; 506 } 507 } 508 if (iterations!=1) 461 509 { //case that sagbiPart called this procedure 462 510 if (size(P)==0) … … 480 528 algebra=algebra,reducedParameters; 481 529 } 482 algebra=simplify(algebra,6); 530 algebra = simplify(algebra,6); 531 algebra = canonicalform(algebra); 483 532 return(algebra); 484 533 } … … 735 784 a=sagbiConstruction(algebra,1,tailreduction,method,0); 736 785 a=simplify(a,7); 737 a=interreduced(a);786 // a=interreduced(a); 738 787 return(a); 739 788 } … … 803 852 ideal a; 804 853 a=sagbiConstruction(algebra,iterations,tailreduction,method,0); 805 a=simplify(a, 7);806 a=interreduced(a);854 a=simplify(a,6); 855 // a=interreduced(a); 807 856 return(a); 808 857 } … … 821 870 822 871 823 824 // VL: finished the documentation825 // TO compare with algDependent from algebra_lib826 // TODO : remove constants from algebraic dependencies827 872 proc algebraicDependence(ideal I,int iterations) 828 873 "USAGE: algebraicDependence(I,it); I an an ideal, it is an integer 829 874 RETURN: ring 830 875 ASSUME: basering is not a qring 831 PURPOSE: In @code{it} iterations, compute algebraic dependencies between elements of I 876 PURPOSE: Returns a ring containing the ideal @code{algDep}, which contains possibly 877 @* some algebraic dependencies of the elements of I obtained through @code{it} 878 @* iterations of the SAGBI construction algorithms. See the example on how 879 @* to access these objects. 832 880 EXAMPLE: example algebraicDependence; shows an example" 833 881 { 834 assumeQring(); 882 assumeQring(); 883 int ppl = printlevelvoice+3; //variable for additional printleveldependend information 884 dbprint(ppl,"//AlgDep1 initialisation and precomputation"); 835 885 def br=basering; 836 886 int i; 837 887 I=simplify(I,2); //avoid that I contains zeros 888 889 //Create two polynomial rings, which both are extensions of the current basering. 890 //The first ring will contain the additional paramteres @c(1),...,@c(m), the second one 891 //will contain the additional variables @c(1),...,@c(m), where m=ncols(I). 838 892 string parameterName=uniqueVariableName("@c"); 839 893 list l = ringlist(basering); … … 843 897 parList[i]=string(parameterName,"(",i,")"); 844 898 } 845 l[1]=list(l[1],parList,list(list("dp",1:ncols(I)))); 899 l[1]=list(l[1],parList,list(list("dp",1:ncols(I)))); //add @c(i) to the ring as paramteres 846 900 ideal temp=0; 847 901 l[1][4]=temp; … … 862 916 for (i=1; i<=ncols(I);i++) 863 917 { 864 l2[2][i+nvars(br)]=string(extendVarName,"(",i,")"); 918 l2[2][i+nvars(br)]=string(extendVarName,"(",i,")"); //add @c(i) to the rings as variables 865 919 } 866 920 l2[3][size(l2[3])+1]=l2[3][size(l2[3])]; … … 899 953 // } 900 954 901 955 //Compute a partial SAGBI basis of the algebra generated by I[1]@c(1),...,I[m]@c(m), 956 //where the @c(n) are parameters 902 957 ideal I=fetch(br,I); 903 958 ideal algebra; … … 906 961 algebra[i]=I[i]par(i); 907 962 } 963 dbprint(ppl,"//AlgDep2 call of SAGBI construction algorithm"); 908 964 algebra=sagbiConstruction(algebra, iterations,0,0,1); 965 dbprint(ppl,"//AlgDep3 postprocessing of results"); 909 966 int j=1; 967 //If K[x_1,...,x_n] was the basering, then algebra is in K(@c(1),...,@c(m))[x_1,...x_n]. We intersect 968 //elements in algebra with K(@c(1),..,@c(n)) to get algDep. Note that @c(i) can only appear in the numerator, 969 //as the SAGBI construction algorithms just multiplies and substracts polynomials. So actually we have 970 //algDep=algebra intersect K[@c(1),...,@c(m)] 910 971 ideal algDep; 911 972 for (i=1; i<= ncols(algebra); i++) 912 973 { 913 if (leadmonom(algebra[i])==1) 974 if (leadmonom(algebra[i])==1) //leadmonom(algebra[i])==1 iff algebra[i] in K[@c(1),...,@c(m)] 914 975 { 915 976 algDep[j]=algebra[i]; … … 917 978 } 918 979 } 980 //Transfer algebraic dependencies to ring where @c(i) are not parameters, but now variables. 919 981 setring extendVarRing; 920 982 ideal algDep=imap(parameterRing,algDep); 921 983 ideal algebra=imap(parameterRing,algebra); 984 //Now get rid of constants in K that may have been added to algDep. 985 for (i=1; i<=ncols(algDep); i++) 986 { 987 if(leadmonom(algDep[i])==1) 988 { 989 algDep[i]=0; 990 } 991 } 992 algDep=simplify(algDep,2); 922 993 export algDep,algebra; 923 //print(algDep);924 994 setring br; 925 995 return(extendVarRing); … … 946 1016 static proc interreduced(ideal I) 947 1017 { 1018 /* performs subalgebra interreduction of a set of subalgebra generators */ 1019 int ppl = printlevelvoice+3; //variable for additional printleveldependend information 1020 dbprint(ppl,"//Interred1 starting interreduction"); 948 1021 ideal J,B; 949 1022 int i,j,k; … … 951 1024 for(k=1;k<=ncols(I);k++) 952 1025 { 1026 dbprint(ppl1,"//Interred1"+string(k)+" reducing next poly"); 953 1027 f=I[k]; 954 1028 I[k]=0; … … 957 1031 } 958 1032 I=simplify(I,2); 1033 dbprint(ppl,"//Interred2 interreduction completed"); 959 1034 return(I); 960 1035 } … … 1148 1223 } 1149 1224 1225 static proc canonicalform(ideal I) 1226 { 1227 /* placeholder for the canonical form of a set of gen's */ 1228 /* for the time being we agree on content(p)=1; that is coeffs with no fractions */ 1229 int i; ideal J=I; 1230 for(i=ncols(I); i>=1; i) 1231 { 1232 J[i] = canonicalform_poly(I[i]); 1233 } 1234 return(J); 1235 } 1236 1237 static proc canonicalform_poly(poly p) 1238 { 1239 /* placeholder for the canonical form of a poly */ 1240 /* for the time being we agree on content(p)=1; that is coeffs with no fractions */ 1241 number n = content(p); 1242 return( p/content(p) ); 1243 } 1244 1150 1245 /* 1151 1246 ring r= 0,(x,y),dp;
Note: See TracChangeset
for help on using the changeset viewer.