Changeset 12603f in git
- Timestamp:
- Mar 8, 2010, 10:11:46 AM (13 years ago)
- Branches:
- (u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
- Children:
- 56b0c82650ffab78c72a5629796b5fe751482448
- Parents:
- 8e865cd8f45ac33fd409c1996172450700565c7b
- Location:
- Singular/LIB
- Files:
-
- 1 edited
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/swalk.lib
r8e865c r12603f 1 ///////////////////////////////////////////////////////////////2 3 1 version="$Id$"; 4 2 category="Commutative Algebra"; 5 6 3 info=" 7 LIBRARY: swalk.lib Sagbi Walk Conversion Algorithm 8 AUTHOR: Junaid Alam Khan 4 LIBRARY: swalk.lib Sagbi Walk Conversion Algorithm 5 AUTHOR: Junaid Alam Khan junaidalamkhan@gmail.com 6 7 OVERVIEW: 8 A library for computing the Sagbi basis of subalgebra through Sagbi walk 9 algorithm. 10 11 THEORY: The concept of SAGBI ( Subalgebra Analog to Groebner Basis for Ideals) 12 is defined in [L. Robbiano, M. Sweedler: Subalgebra Bases, volume 42, volume 13 1430 of Lectures Note in Mathematics series, Springer-Verlag (1988),61-87]. 14 The Sagbi Walk algorithm is the subalgebra analogue to the Groebner Walk 15 algorithm which has been proposed in [S. Collart, M. Kalkbrener and D.Mall: 16 Converting bases with the Grobner Walk. J. Symbolic Computation 24 (1997), 17 465-469]. 9 18 10 19 PROCEDURE: 11 swalk(ideal[,intvec]); Sagbi basis of subalgebra via Sagbiwalk algorithm20 swalk(ideal[,intvec]); Sagbi basis of subalgebra via Sagbi walk algorithm 12 21 "; 13 22 14 LIB "algebra.lib";15 23 LIB "sagbi.lib"; 16 24 ////////////////////////////////////////////////////////////////////////////// 17 25 proc swalk(ideal Gox, list #) 18 "USAGE: swalk(i); ideal; 19 swalk( i,v,w); i ideal, v,w integer vectors; 20 RETURN: If list #= v,w (resp empty) then it compute the sagbi basis of 21 the subalgebra defined by the generators of ideal, 22 calculated via the Sagbi walk algorithm from the ordering 23 (a(v),lp) (resp dp) tothe ordering (a(w),lp) ( resp lp). 26 "USAGE: swalk(i[,v,w]); i ideal, v,w int vectors 27 RETURN: The sagbi basis of the subalgebra defined by the generators of i, 28 calculated via the Sagbi walk algorithm from the ordering dp to lp 29 if v,w are not given (resp. from the ordering (a(v),lp) to the 30 ordering (a(w),lp) if v and w are given). 24 31 EXAMPLE: example swalk; shows an example 25 32 " … … 33 40 option(redSB); 34 41 def xR = basering; 35 execute("ring ostR = ("+charstr(xR)+"),("+varstr(xR)+"),"+ord_str+";"); 42 list rl=ringlist(xR); 43 rl[3][1][1]="dp"; 44 def ostR=ring(rl); 45 setring ostR; 36 46 def new_ring = basering; 37 47 ideal Gnew = fetch(xR, Gox); 38 48 Gnew=sagbi(Gnew,1); 39 49 Gnew=interreduceSd(Gnew); 40 vector curr_weight= Ctv(icurr_weight);41 vector target_weight= Ctv(itarget_weight);50 vector curr_weight=changeTypeInt(icurr_weight); 51 vector target_weight=changeTypeInt(itarget_weight); 42 52 ideal Gold; 43 53 list l; … … 53 63 if(curr_weight==target_weight){n=1;} 54 64 else { 55 l= Df(Gold);65 l=collectDiffExpo(Gold); 56 66 ulast=last(curr_weight, target_weight, l); 57 vector new_weight= newvec(curr_weight,target_weight,ulast);67 vector new_weight=(1-ulast)*curr_weight+ulast*target_weight; 58 68 vector w=cleardenom(new_weight); 59 v= CT(w);69 v=changeType(w); 60 70 list p= ringlist(old_ring); 61 71 p[3][1][2]= v; … … 77 87 example 78 88 { 79 "EXAMPLE:"; 80 //** compute a Sagbi basis of I w.r.t. lp. 89 "EXAMPLE:";echo = 2; 81 90 ring r = 0,(x,y), lp; 82 91 ideal I =x2,y2,xy+y,2xy2+y3; 83 intvec v=1,1; 84 intvec w=1,0; 85 swalk(I,v,w); 86 } 87 88 proc inprod(vector v,vector w) 92 swalk(I); 93 } 94 95 ////////////////////////////////////////////////////////////////////////////// 96 static proc inprod(vector v,vector w) 89 97 "USAGE: inprod(v,w); v,w vectors 90 98 RETURN: inner product of vector v and w … … 108 116 } 109 117 110 111 proc df(poly f)112 "USAGE: d f(f); f polynomial113 RETURN: a list of integers vectors which are the differ nce of exponent vector114 of leading monomial of f with the exponent vector of115 of othermonmials in f.116 EXAMPLE: example d f; shows an example118 ////////////////////////////////////////////////////////////////////////////// 119 static proc diffExpo(poly f) 120 "USAGE: diffExpo(f); f polynomial 121 RETURN: a list of integers vectors which are the difference of exponent 122 vector of leading monomial of f with the exponent vector of of other 123 monmials in f. 124 EXAMPLE: example diffExpo; shows an example 117 125 " 118 126 { … … 120 128 int i; 121 129 intvec v; 122 for(i= 2;i<=size(f);i++)130 for(i=size(f);i>=2;i--) 123 131 { 124 132 v=leadexp(f[1])-leadexp(f[i]); 125 l[ size(l)+1]=v;133 l[i-1]=v; 126 134 } 127 135 return(l); … … 131 139 ring r=0,(x,y,z),lp; 132 140 poly f = xy+z2 ; 133 d f(f);134 } 135 136 137 proc Df( ideal i)138 "USAGE: Df(i); i ideal139 RETURN: a list which contain df(f), for all generators f of ideal i140 EXAMPLE: example Df; shows an example141 diffExpo(f); 142 } 143 144 ////////////////////////////////////////////////////////////////////////////// 145 static proc collectDiffExpo( ideal i) 146 "USAGE: collectDiffExpo(i); i ideal 147 RETURN: a list which contains diffExpo(f), for all generators f of ideal i 148 EXAMPLE: example collectDiffExpo; shows an example 141 149 " 142 150 { 143 151 list l; 144 152 int j; 145 for(j= 1;j<=size(i);j++)153 for(j=ncols(i); j>=1;j--) 146 154 { 147 l[ size(l)+1]=df(i[j]);155 l[j]=diffExpo(i[j]); 148 156 } 149 return(l);150 157 return(l); 158 } 151 159 example 152 160 { "EXAMPLE:"; echo = 2; 153 161 ring r=0,(x,y,z),lp; 154 162 ideal I = xy+z2,y3+x2y2; 155 Df(I); 156 } 157 158 159 proc newvec(vector v, vector w, number u) 160 "USAGE: newvec(v,w,u); v,w vectors, u number 161 RETURN: the vector (1-u)v+(u)w. 162 EXAMPLE: example newvec; shows an example 163 " 164 { 165 vector wnew ; 166 wnew=(1-u)*v+(u)*w ; 167 return(wnew); 168 } 169 example 170 { "EXAMPLE:"; echo = 2; 171 ring r=0,(x,y,z),lp; 172 vector v=[0,0,1]; 173 vector w=[1,0,0]; 174 number u=2/3; 175 newvec(v,w,u); 176 } 177 178 proc CT(vector v) 179 "USAGE: CT(v); v vector 180 RETURN: change the type of vector v into integer vector. 181 EXAMPLE: example CT; shows an example 163 collectDiffExpo(I); 164 } 165 166 ////////////////////////////////////////////////////////////////////////////// 167 static proc changeType(vector v) 168 "USAGE: changeType(v); v vector 169 RETURN: change the type of vector 170 v into integer vector. 171 172 EXAMPLE: example changeType; shows an example 182 173 " 183 174 { … … 194 185 ring r=0,(x,y,z),lp; 195 186 vector v = [2,1,3]; 196 CT(v); 197 } 198 199 200 proc Ctv( intvec v) 201 "USAGE: Ctv(v); v integer vector 187 changeType(v); 188 } 189 ////////////////////////////////////////////////////////////////////////////// 190 static proc changeTypeInt( intvec v) 191 "USAGE: changeTypeInt(v); v integer vector 202 192 RETURN: change the type of integer vector v into vector. 203 EXAMPLE: example Ctv; shows an example204 " 205 { 206 vector w;207 int j ;208 for(j=1;j<=size(v);j++)209 {210 w=w+v[j]*gen(j);211 }212 return(w);193 EXAMPLE: example changeTypeInt; shows an example 194 " 195 { 196 vector w; 197 int j ; 198 for(j=1;j<=size(v);j++) 199 { 200 w=w+v[j]*gen(j); 201 } 202 return(w); 213 203 } 214 204 example … … 216 206 ring r=0,(x,y,z),lp; 217 207 intvec v = 4,2,3; 218 Ctv(v); 219 } 220 221 proc last( vector c, vector t,list l) 208 changeTypeInt(v); 209 } 210 211 ////////////////////////////////////////////////////////////////////////////// 212 static proc last( vector c, vector t,list l) 222 213 "USAGE: last(c,t,l); c,t vectors, l list 223 RETURN: a parametric value which corresponds to vector lies along 224 the path between c and t using list l of integer225 vectors. This vector is thelast vector on old Sagbi cone214 RETURN: a parametric value which corresponds to vector lies along the path 215 between c and t using list l of integer vectors. This vector is the 216 last vector on old Sagbi cone 226 217 EXAMPLE: example last; shows an example 227 218 " … … 232 223 vector v; 233 224 for(i=1;i<=size(l);i++) 234 225 { 235 226 for(j=1;j<=size(l[i]);j++) 236 227 { 237 228 v=0; 238 229 for(k=1;k<=size(l[i][j]);k++) 239 230 { 240 231 v=v+l[i][j][k]*gen(k); 241 232 } 242 233 poly n= inprod(c,v); 243 234 poly q= inprod(t,v); … … 246 237 number z=a-b; 247 238 if(b<0) 248 239 { 249 240 u=a/z; 250 241 if(u<ul) {ul=u;} 251 242 } 252 243 kill a,b,z,n,q ; 253 254 255 return(ul);244 } 245 } 246 return(ul); 256 247 } 257 248 example … … 261 252 vector w=[1,0,0]; 262 253 ideal i=z2+xy,x2y2+y3; 263 list l= Df(i);254 list l=collectDiffExpo(i); 264 255 last(v,w,l) 265 256 } 266 257 267 proc initialForm(poly P) 258 ////////////////////////////////////////////////////////////////////////////// 259 static proc initialForm(poly P) 268 260 "USAGE: initialForm(P); P polynomial 269 261 RETURN: sum of monomials of P with maximum w-degree … … 275 267 int i=1; 276 268 while(deg(P[i])==deg(P[1])) 277 269 { 278 270 q=q+P[i]; 279 271 i++; 280 272 if(i>size(P)) {break;} 281 282 273 } 274 return(q); 283 275 } 284 276 example … … 289 281 } 290 282 291 proc Initial(ideal I) 283 ////////////////////////////////////////////////////////////////////////////// 284 static proc Initial(ideal I) 292 285 "USAGE: Initial(I); I ideal 293 RETURN: an ideal which is generate by the InitialForm of the generators of I. 286 RETURN: an ideal which is generate by the InitialForm 287 of the generators of I. 294 288 EXAMPLE: example Initial; shows an example 295 289 " … … 297 291 ideal J; 298 292 int i; 299 for(i=1;i<= size(I);i++)300 293 for(i=1;i<=ncols(I);i++) 294 { 301 295 J[i]=initialForm(I[i]); 302 296 } 303 297 return(J); 304 298 } … … 310 304 } 311 305 312 proc Lift(ideal In,ideal InG,ideal Gold) 313 "USAGE: Lift(In, InG, Gold); In,InG, Gold ideals 314 Gold given by Sagbi basis {g_1,...,g_t} 315 In given by tne initial forms In(g_1),...,In(g_t) 316 InG={h_1,...,h_s} a Sagbi basis of In 317 RETURN: P_j, a polynomial in K[y_1,..,y_t] such that 318 h_j=P_j(In(g_1),...,In_(g_t)) and return f_j=P_j(g_1,...,g_t) 306 ////////////////////////////////////////////////////////////////////////////// 307 static proc Lift(ideal In,ideal InG,ideal Gold) 308 "USAGE: Lift(In, InG, Gold); In, InG, Gold ideals; 309 Gold given by Sagbi basis {g_1,...,g_t}, 310 In given by tne initial forms In(g_1),...,In(g_t), 311 InG = {h_1,...,h_s} a Sagbi basis of In 312 RETURN: P_j, a polynomial in K[y_1,..,y_t] such that h_j = 313 P_j(In(g_1),...,In_(g_t)) 314 and return f_j = P_j(g_1,...,g_t) 319 315 EXAMPLE: example Lift; shows an example 320 316 " … … 322 318 int i; 323 319 ideal J; 324 for(i=1;i<= size(InG);i++)320 for(i=1;i<=ncols(InG);i++) 325 321 { 326 322 poly f=InG[i]; … … 344 340 } 345 341 346 347 proc Convert( ideal Gold )342 ////////////////////////////////////////////////////////////////////////////// 343 static proc Convert( ideal Gold ) 348 344 "USAGE: Convert(I); I ideal 349 345 RETURN: Convert old Sagbi basis into new Sagbi basis … … 363 359 } 364 360 365 proc interreduceSd(ideal I) 361 ////////////////////////////////////////////////////////////////////////////// 362 static proc interreduceSd(ideal I) 366 363 "USAGE: interreduceSd(I); I ideal 367 364 RETURN: interreduceSd the set of generators if I with … … 374 371 int i,j,k; 375 372 poly f; 376 for(k=1;k<= size(I);k++)377 373 for(k=1;k<=ncols(I);k++) 374 {l[k]=I[k];} 378 375 for(j=1;j<=size(l);j++) 379 376 { 380 377 f=l[j]; 381 378 M=delete(l,j); 382 379 for(i=1;i<=size(M);i++) 383 380 { B[i]=M[i];} 384 381 f=sagbiNF(f,B,1); 385 382 J=J+f; 386 383 } 387 384 return(J); 388 385 } … … 394 391 } 395 392 393 ////////////////////////////////////////////////////////////////////////////// 396 394 static proc OrderStringalp(string Wpal,list #) 397 395 { … … 478 476 } 479 477 478 ////////////////////////////////////////////////////////////////////////////// 480 479 static proc OrderStringalp_NP(string Wpal,list #) 481 480 { … … 489 488 if(size(#) == 1) 490 489 { 491 if(typeof(#[1]) == "intvec") { 490 if(typeof(#[1]) == "intvec") 491 { 492 492 curr_weight = #[1]; 493 493 494 if(Wpal == "al"){ 494 if(Wpal == "al") 495 { 495 496 order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)"; 496 497 } 497 else { 498 else 499 { 498 500 order_str = "(Wp("+string(#[1])+"),C)"; 499 501 } 500 502 } 501 503 else { 502 if(typeof(#[1]) == "int"){ 504 if(typeof(#[1]) == "int") 505 { 503 506 nP = #[1]; 504 507 } 505 else { 508 else 509 { 506 510 print("// ** the input must be \"(ideal, int)\" or "); 507 511 print("// ** \"(ideal, intvec)\""); … … 514 518 if(size(#) == 2) 515 519 { 516 if(typeof(#[1]) == "intvec" and typeof(#[2]) == "int"){ 520 if(typeof(#[1]) == "intvec" and typeof(#[2]) == "int") 521 { 517 522 curr_weight = #[1]; 518 523 519 if(Wpal == "al"){ 524 if(Wpal == "al") 525 { 520 526 order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)"; 521 527 } 522 else { 528 else 529 { 523 530 order_str = "(Wp("+string(#[1])+"),C)"; 524 531 } 525 532 } 526 else{ 527 if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec"){ 533 else 534 { 535 if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec") 536 { 528 537 curr_weight = #[1]; 529 538 target_weight = #[2]; 530 539 531 if(Wpal == "al"){ 540 if(Wpal == "al") 541 { 532 542 order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)"; 533 543 } 534 else { 544 else 545 { 535 546 order_str = "(Wp("+string(#[1])+"),C)"; 536 547 } 537 548 } 538 else{ 549 else 550 { 539 551 print("// ** the input must be \"(ideal,intvec,int)\" or "); 540 552 print("// ** \"(ideal,intvec,intvec)\""); … … 544 556 } 545 557 else { 546 if(size(#) == 3) { 558 if(size(#) == 3) 559 { 547 560 if(typeof(#[1]) == "intvec" and typeof(#[2]) == "intvec" and 548 561 typeof(#[3]) == "int") … … 551 564 target_weight = #[2]; 552 565 nP = #[3]; 553 if(Wpal == "al"){ 566 if(Wpal == "al") 567 { 554 568 order_str = "(a("+string(#[1])+"),lp("+string(n) + "),C)"; 555 569 } 556 else { 570 else 571 { 557 572 order_str = "(Wp("+string(#[1])+"),C)"; 558 573 } 559 574 } 560 else{ 575 else 576 { 561 577 print("// ** the input must be \"(ideal,intvec,intvec,int)\""); 562 578 print("// ** and a lex. GB will be computed from \"dp\" to \"lp\""); … … 564 580 } 565 581 } 566 else{ 582 else 583 { 567 584 print("// ** The given input is wrong"); 568 585 print("// ** and a lex. GB will be computed from \"dp\" to \"lp\""); … … 579 596 } 580 597 598 ///////////////////////////////////////////////////////////////////////////////* 599 Further examples 600 ring r=0,(x,y,z),lp; 601 602 ideal I=x2y4, y4z2, xy4z+y2z, xy6z2+y10z5; 603 604 ideal Q=x2y4, y4z2, xy4z+y2z, xy6z2+y14z7; 605 606 ideal J=x2y4, y4z2, xy4z+y2z, xy6z2+y18z9; 607 608 ideal K=x2,y2,xy+y,2xy2+y5,z3+x; 609 610 ideal L=x2+y,y2+z,x+z2; 611 */ 612 613
Note: See TracChangeset
for help on using the changeset viewer.