Changeset 565e86 in git
 Timestamp:
 Dec 23, 2008, 10:39:31 PM (14 years ago)
 Branches:
 (u'spielwiese', '91fdef05f09f54b8d58d92a472e9c4a43aa4656f')
 Children:
 199ae72040d35b09813089e33464df761b5aad78
 Parents:
 6045c31371ebacbf1acfc2f082225f05016f06b2
 Location:
 Singular/LIB
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/bfct.lib
r6045c3 r565e86 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: bfct.lib,v 1. 7 20081209 16:50:21 levandov Exp $";2 version="$Id: bfct.lib,v 1.8 20081223 21:39:31 levandov Exp $"; 3 3 category="Noncommutative"; 4 4 info=" 5 LIBRARY: bfct.lib M. Noro's Algorithm forBernsteinSato polynomial5 LIBRARY: bfct.lib Algorithms for bfunctions and BernsteinSato polynomial 6 6 AUTHORS: Daniel Andres, daniel.andres@math.rwthaachen.de 7 7 @* Viktor Levandovskyy, levandov@math.rwthaachen.de 8 8 9 9 THEORY: Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R, 10 @* one is interested in the global BernsteinSato polynomial b(s) in K[s], 11 @* defined to be the monic polynomial, satisfying a functional identity 12 @* L * f^{s+1} = b(s) f^s, for some operator L in D[s]. 13 @* Here, D stands for an nth Weyl algebra K<x_1,...,x_n,d_1,...,d_n  d_j x_j = x_j d_j +1> 14 @* One is interested in the following data: 15 @* global BernsteinSato polynomial in K[s] and 16 @* the list of all roots of b(s), which are known to be rational, with their multiplicities. 10 @* one is interested in the global bFunction (also known as BernsteinSato 11 @* polynomial) b(s) in K[s], defined to be the monic polynomial, satisfying a functional 12 @* identity L * F^{s+1} = b(s) F^s, for some operator L in D[s]. Here, D stands for an 13 @* nth Weyl algebra K<x_1,...,x_n,d_1,...,d_n  d_j x_j = x_j d_j +1>. 14 @* One is interested in the following data: 15 @*  BernsteinSato polynomial b(s) in K[s], 16 @*  the list of its roots, which are known to be rational 17 @*  the multiplicities of the roots. 18 @* References: Saito, Strurmfels, Takayama: Groebner Deformations of Hypergeometric 19 @* Differential Equations (2000), Noro: An Efficient Modular Algorithm for Computing 20 @* the Global bfunction, (2002). 21 17 22 18 23 MAIN PROCEDURES: 19 24 20 bfct(f[,s,t,v]); compute the global BernsteinSato polynomial of a given poly21 bfct syz(f[,r,s,t,u,v]); compute the global BernsteinSato polynomial of a given poly22 bfct ann(f[,s]); compute the global BernsteinSato polynomial of a given poly23 bfct onestep(f[,s,t]); compute the global BernsteinSato polynomial of a given poly24 bfct ideal(I,w[,s,t]); compute the global bfunction of a given ideal w.r.t. a given weight25 p intersect(f,I); compute the intersection of the ideal generated by a given poly with a given ideal26 p intersectsyz(f,I[,p,s,t]); compute the intersection of the ideal generated by a given poly with a given ideal27 lin reduce(f,I[,s]); reducea poly by linear reductions only28 ncsolve(I[,s]); find and compute a linear dependency of the elements of an ideal 25 bfct(f[,s,t,v]); computes the BernsteinSato polynomial of poly f 26 bfctSyz(f[,r,s,t,u,v]); computes the BernsteinSato polynomial of poly f 27 bfctAnn(f[,s]); computes the BernsteinSato polynomial of poly f 28 bfctOneGB(f[,s,t]); computes the BernsteinSato polynomial of poly f 29 bfctIdeal(I,w[,s,t]); computes the global bfunction of ideal I w.r.t. the weight w 30 pIntersect(f,I); intersection of the ideal I with the subalgebra K[f] for a poly f 31 pIntersectSyz(f,I[,p,s,t]); intersection of the ideal I with the subalgebra K[f] for a poly f 32 linReduce(f,I[,s]); reduces a poly by linear reductions only 33 linSyzSolve(I[,s]); computes a linear dependency of the elements of ideal I 29 34 30 35 AUXILIARY PROCEDURES: 31 36 32 ispositive(v); check whether all entries of an intvec are positive 33 isin(l,i); check whether an element is a member of a list 34 scalarprod(v,w); compute the standard scalar product of two intvecs 35 vec2poly(v[,i]); convert a coefficient vector to a poly 37 allPositive(v); checks whether all entries of an intvec are positive 38 scalarProd(v,w); computes the standard scalar product of two intvecs 39 vec2poly(v[,i]); constructs an univariate poly with given coefficients 36 40 37 41 SEE ALSO: dmod_lib, dmodapp_lib, gmssing_lib … … 40 44 41 45 LIB "qhmoduli.lib"; // for Max 42 LIB "dmod.lib"; // for SannfsBFCT etc 43 LIB "dmodapp.lib"; // for initialideal etc 44 45 46 LIB "dmod.lib"; // for SannfsBFCT etc 47 LIB "dmodapp.lib"; // for initialIdeal etc 48 49 50 /////////////////////////////////////////////////////////////////////////////// 51 // testing for consistency of the library: 46 52 proc testbfctlib () 47 53 { 48 54 // tests all procs for consistency 49 55 "AUXILIARY PROCEDURES:"; 50 example ispositive; 51 example isin; 52 example scalarprod; 56 example allPositive; 57 example scalarProd; 53 58 example vec2poly; 54 59 "MAIN PROCEDURES:"; 55 60 example bfct; 56 example bfct syz;57 example bfct ann;58 example bfct onestep;59 example bfct ideal;60 example p intersect;61 example p intersectsyz;62 example lin reduce;63 example ncsolve;64 } 65 61 example bfctSyz; 62 example bfctAnn; 63 example bfctOneGB; 64 example bfctIdeal; 65 example pIntersect; 66 example pIntersectSyz; 67 example linReduce; 68 example linReduceIdeal; 69 example linSyzSolve; 70 } 66 71 67 72 // auxiliary procedures  … … 71 76 RETURN: a ring, the associated graded ring of the basering w.r.t. u and v 72 77 PURPOSE: compute the associated graded ring of the basering w.r.t. u and v 78 ASSUME: basering is a Weyl algebra 73 79 EXAMPLE: example gradedWeyl; shows examples 74 80 NOTE: u[i] is the weight of x(i), v[i] the weight of D(i). … … 122 128 { 123 129 "EXAMPLE:"; echo = 2; 124 LIB "bfct.lib";125 130 ring @D = 0,(x,y,z,Dx,Dy,Dz),dp; 126 131 def D = Weyl(); … … 132 137 133 138 134 proc ispositive (intvec v)135 "USAGE: ispositive(v); v an intvec136 RETURN: 1 if all components of v are positive, or 0 otherwise139 proc allPositive (intvec v) 140 "USAGE: allPositive(v); v an intvec 141 RETURN: int, 1 if all components of v are positive, or 0 otherwise 137 142 PURPOSE: check whether all components of an intvec are positive 138 EXAMPLE: example ispositive; shows an example143 EXAMPLE: example allPositive; shows an example 139 144 " 140 145 { … … 154 159 "EXAMPLE:"; echo = 2; 155 160 intvec v = 1,2,3; 156 ispositive(v);161 allPositive(v); 157 162 intvec w = 1,2,3; 158 ispositive(w);159 } 160 161 proc isin(list l, i)162 "USAGE: isin(l,i); l a list, i an argument of any type163 RETURN: anint, the position of the first appearance of i in l, or 0 if i is not a member of l163 allPositive(w); 164 } 165 166 static proc findFirst (list l, i) 167 "USAGE: findFirst(l,i); l a list, i an argument of any type 168 RETURN: int, the position of the first appearance of i in l, or 0 if i is not a member of l 164 169 PURPOSE: check whether the second argument is a member of a list 165 EXAMPLE: example isin; shows an example170 EXAMPLE: example findFirst; shows an example 166 171 " 167 172 { … … 177 182 return(0); 178 183 } 184 // isin(list(1, 2, list(1)),2); // seems to be a bit buggy, 185 // isin(list(1, 2, list(1)),3); // but works for the purposes 186 // isin(list(1, 2, list(1)),list(1)); // of this library 187 // isin(l,list(2)); 179 188 example 180 189 { … … 182 191 ring r = 0,x,dp; 183 192 list l = 1,2,3; 184 isin(l,4);185 isin(l,2);186 } 187 188 proc scalar prod (intvec v, intvec w)189 "USAGE: scalar prod(v,w); v,w intvecs190 RETURN: anint, the standard scalar product of v and w191 PURPOSE: compute the scalar product of two intvecs192 NOTE: the arguments must havethe same size193 EXAMPLE: example scalar prod; shows examples193 findFirst(l,4); 194 findFirst(l,2); 195 } 196 197 proc scalarProd (intvec v, intvec w) 198 "USAGE: scalarProd(v,w); v,w intvecs 199 RETURN: int, the standard scalar product of v and w 200 PURPOSE: computes the scalar product of two intvecs 201 ASSUME: the arguments are of the same size 202 EXAMPLE: example scalarProd; shows examples 194 203 " 195 204 { … … 213 222 intvec v = 1,2,3; 214 223 intvec w = 4,5,6; 215 scalarprod(v,w); 216 } 217 224 scalarProd(v,w); 225 } 218 226 219 227 // main procedures  220 228 221 222 proc linreduce(poly f, ideal I, list #) 223 "USAGE: linreduce(f, I [,s,t]); f a poly, I an ideal, s,t optional ints 224 RETURN: a poly obtained by linear reductions of the given poly with the given ideal 225 PURPOSE: reduce a poly only by linear reductions 226 NOTE: If s<>0, a list consisting of the reduced poly and the coefficient vector of the used 227 @* reductions is returned. 228 @* If t<>0, only leading monomials are reduced, otherwise, and by default, all monomials 229 @* are reduced, if possible. 230 EXAMPLE: example linreduce; shows examples 229 proc linReduceIdeal(ideal I, list #) 230 "USAGE: linReduceIdeal(I [,s,t]); I an ideal, s,t optional ints 231 RETURN: ideal/list, linear reductum (over field) of f by elements from I 232 PURPOSE: reduce a poly only by linear reductions (no monomial multiplications) 233 NOTE: If s<>0, a list consisting of the reduced ideal and the coefficient 234 @* vectors of the used reductions given as module is returned. 235 @* Otherwise (and by default), only the reduced ideal is returned. 236 @* If t=0 (and by default) all monomials are reduced (if possible), 237 @* otherwise, only leading monomials are reduced. 238 DISPLAY: If @code{printlevel}>=1, all debug messages will be printed. 239 EXAMPLE: example linReduceIdeal; shows examples 231 240 " 232 241 { 242 // #[1] = remembercoeffs 243 // #[2] = redtail 233 244 int ppl = printlevel  voice + 2; 234 245 int remembercoeffs = 0; // default 235 int red lm= 0; // default246 int redtail = 0; // default 236 247 if (size(#)>0) 237 248 { … … 244 255 if (typeof(#[2])=="int"  typeof(#[2])=="number") 245 256 { 246 redlm = #[2]; 257 redtail = #[2]; 258 } 259 } 260 } 261 int sI = ncols(I); 262 int sZeros = sI  size(I); 263 dbprint(ppl,"ideal contains "+string(sZeros)+" zero(s)"); 264 int i,j; 265 ideal J,lmJ,ordJ; 266 list lJ = sort(I); 267 module M; // for the coefficients 268 if (sZeros > 0) // I contains zeros 269 { 270 if (remembercoeffs <> 0) 271 { 272 j = 0; 273 for (i=1; i<=sI; i++) 274 { 275 if (I[i] == 0) 276 { 277 j++; 278 J[j] = 0; 279 ordJ[j] = 1; 280 M[j] = gen(i); 281 } 282 else 283 { 284 M[i+sZerosj] = gen(lJ[2][ij]+j); 285 } 286 } 287 } 288 else 289 { 290 for (i=1; i<=sZeros; i++) 291 { 292 J[i] = 0; 293 ordJ[i] = 1; 294 } 295 } 296 I = J,lJ[1]; 297 } 298 else 299 { 300 I = lJ[1]; 301 if (remembercoeffs <> 0) 302 { 303 for (i=1; i<=size(lJ[1]); i++) { M[i+sZeros] = gen(lJ[2][i]); } 304 } 305 } 306 dbprint(ppl,"initially sorted ideal:", I); 307 if (remembercoeffs <> 0) { dbprint(ppl," used permutations:", M); } 308 poly lm,c,redpoly,maxlmJ; 309 J[sZeros+1] = I[sZeros+1]; 310 lm = I[sZeros+1]; 311 maxlmJ = leadmonom(lm); 312 lmJ[sZeros+1] = maxlmJ; 313 int ordlm,reduction,maxordJ; 314 maxordJ = ord(lm); 315 ordJ[sZeros+1] = maxordJ; 316 vector v; 317 for (i=sZeros+2; i<=sI; i++) 318 { 319 redpoly = I[i]; 320 lm = leadmonom(redpoly); 321 ordlm = ord(lm); 322 reduction = 1; 323 if (remembercoeffs <> 0) { v = M[i]; } 324 while (reduction == 1) // while there was a reduction 325 { 326 reduction = 0; 327 for (j=sZeros+1; j<i; j++) 328 { 329 if (lm == 0) { break; } 330 if (ordlm > maxordJ) { break; } 331 if (ordlm == ordJ[j]) 332 { 333 if (lm > maxlmJ) { break; } 334 if (lm == lmJ[j]) 335 { 336 dbprint(ppl,"reducing " + string(redpoly)); 337 dbprint(ppl," with " + string(J[j])); 338 c = leadcoef(redpoly)/leadcoef(J[j]); 339 redpoly = redpoly  c*J[j]; 340 dbprint(ppl," to " + string(redpoly)); 341 lm = leadmonom(redpoly); 342 ordlm = ord(lm); 343 if (remembercoeffs <> 0) { M[i] = M[i]  c * M[j]; } 344 reduction = 1; 345 } 346 } 347 } 348 } 349 for (j=sZeros+1; j<i; j++) 350 { 351 if (redpoly < J[j]) { break; } 352 } 353 J = insertGenerator(J,redpoly,j); 354 lmJ = insertGenerator(lmJ,lm,j); 355 ordJ = insertGenerator(ordJ,poly(ordlm),j); 356 if (remembercoeffs <> 0) 357 { 358 v = M[i]; 359 M = deleteGenerator(M,i); 360 M = insertGenerator(M,v,j); 361 } 362 maxordJ = ord(J[i]); 363 maxlmJ = lmJ[i]; 364 } 365 if (redtail <> 0) 366 { 367 dbprint(ppl,"finished reducing leading monomials:",J); 368 if (remembercoeffs <> 0) { dbprint(ppl,"used reductions:",M); } 369 int k; 370 for (i=sZeros+1; i<=sI; i++) 371 { 372 lm = lmJ[i]; 373 for (j=i+1; j<=sI; j++) 374 { 375 for (k=2; k<=size(J[j]); k++) // run over all terms in J[j] 376 { 377 if (ordJ[i] == ord(J[j][k])) 378 { 379 if (lm == normalize(J[j][k])) 380 { 381 c = leadcoef(J[j][k])/leadcoef(J[i]); 382 dbprint(ppl,"reducing " + string(J[j])); 383 dbprint(ppl," with " + string(J[i])); 384 J[j] = J[j]  c*J[i]; 385 dbprint(ppl," to " + string(J[j])); 386 if (remembercoeffs <> 0) { M[j] = M[j]  c * M[i]; } 387 } 388 } 389 } 390 } 391 } 392 } 393 if (remembercoeffs <> 0) { return(list(J,M)); } 394 else { return(J); } 395 } 396 example 397 { 398 "EXAMPLE:"; echo = 2; 399 ring r = 0,(x,y),dp; 400 ideal I = 3,x+9,y4,y4+7x+2; 401 linReduceIdeal(I); 402 linReduceIdeal(I,0,1); 403 list l = linReduceIdeal(I,1,1); l; 404 module M = I; 405 l[1]  ideal(M*l[2]); 406 } 407 408 proc linReduce(poly f, ideal I, list #) 409 "USAGE: linReduce(f, I [,s,t]); f a poly, I an ideal, s,t optional ints 410 RETURN: poly/list, linear reductum (over field) of f by elements from I 411 PURPOSE: reduce a poly only by linear reductions (no monomial multiplications) 412 NOTE: If s<>0, a list consisting of the reduced poly and the coefficient 413 @* vector of the used reductions is returned, otherwise (and by default) 414 @* only reduced poly is returned. 415 @* If t=0 (and by default) all monomials are reduced (if possible), 416 @* otherwise, only leading monomials are reduced. 417 DISPLAY: If @code{printlevel}>=1, all debug messages will be printed. 418 EXAMPLE: example linReduce; shows examples 419 " 420 { 421 int ppl = printlevel  voice + 2; 422 int remembercoeffs = 0; // default 423 int redtail = 0; // default 424 int prepareideal = 1; // default 425 if (size(#)>0) 426 { 427 if (typeof(#[1])=="int"  typeof(#[1])=="number") 428 { 429 remembercoeffs = #[1]; 430 } 431 if (size(#)>1) 432 { 433 if (typeof(#[2])=="int"  typeof(#[2])=="number") 434 { 435 redtail = #[2]; 436 } 437 if (size(#)>2) 438 { 439 if (typeof(#[3])=="int"  typeof(#[3])=="number") 440 { 441 prepareideal = #[3]; 442 } 247 443 } 248 444 } … … 250 446 int i,j,k; 251 447 int sI = ncols(I); 252 ideal lmI,lcI; 448 // prereduce I: 449 module M; 450 if (prepareideal <> 0) 451 { 452 if (remembercoeffs <> 0) 453 { 454 list sortedI = linReduceIdeal(I,1,redtail); 455 I = sortedI[1]; 456 M = sortedI[2]; 457 dbprint(ppl,"prepared ideal:",I); 458 dbprint(ppl," with operations:",M); 459 } 460 else 461 { 462 I = linReduceIdeal(I,0,redtail); 463 } 464 } 465 else 466 { 467 if (remembercoeffs <> 0) 468 { 469 for (i=1; i<=sI; i++) 470 { 471 M[i] = gen(i); 472 } 473 } 474 } 475 ideal lmI,lcI,ordI; 253 476 for (i=1; i<=sI; i++) 254 477 { 255 478 lmI[i] = leadmonom(I[i]); 256 479 lcI[i] = leadcoef(I[i]); 480 ordI[i] = ord(lmI[i]); 257 481 } 258 482 vector v; 259 483 poly c; 260 int reduction = 1; 261 if (redlm == 0) 262 { 263 ideal monomf; 264 for (k=1; k<=size(f); k++) 265 { 266 monomf[k] = normalize(f[k]); 267 } 268 while (reduction == 1) // while there was a reduction 269 { 270 reduction = 0; 271 for (i=sI; i>=1; i) 272 { 273 dbprint(ppl,"testing ideal entry:",i); 274 for (j=1; j<=size(f); j++) 484 // === reduce leading monomials === 485 poly lm = leadmonom(f); 486 int ordf = ord(lm); 487 for (i=sI; i>=1; i) // I is sorted: smallest lm's on top 488 { 489 if (lm == 0) 490 { 491 break; 492 } 493 else 494 { 495 if (ordf == ordI[i]) 496 { 497 if (lm == lmI[i]) 275 498 { 276 if (monomf[j] == lmI[i]) 499 c = leadcoef(f)/lcI[i]; 500 f = f  c*I[i]; 501 lm = leadmonom(f); 502 ordf = ord(lm); 503 if (remembercoeffs <> 0) 504 { 505 v = v  c * M[i]; 506 } 507 } 508 } 509 } 510 } 511 // === reduce tails as well === 512 if (redtail <> 0) 513 { 514 dbprint(ppl,"poly after reduction of leading monomials:",f); 515 for (i=1; i<=sI; i++) 516 { 517 dbprint(ppl,"testing ideal entry:",i); 518 for (j=1; j<=size(f); j++) 519 { 520 if (ord(f[j]) == ordI[i]) 521 { 522 if (normalize(f[j]) == lmI[i]) 277 523 { 278 524 c = leadcoef(f[j])/lcI[i]; 279 525 f = f  c*I[i]; 280 526 dbprint(ppl,"reducing poly to ",f); 281 monomf = 0;282 for (k=1; k<=size(f); k++)283 {284 monomf[k] = normalize(f[k]);285 }286 reduction = 1;287 527 if (remembercoeffs <> 0) 288 528 { 289 v = v  c * gen(i);529 v = v  c * M[i]; 290 530 } 291 break;292 }293 }294 if (reduction == 1)295 {296 break;297 }298 }299 }300 }301 else // reduce only leading monomials302 {303 poly lm = leadmonom(f);304 while (reduction == 1) // while there was a reduction305 {306 reduction = 0;307 for (i=sI;i>=1;i)308 {309 if (lm <> 0 && lm == lmI[i])310 {311 c = leadcoef(f)/lcI[i];312 f = f  c*I[i];313 lm = leadmonom(f);314 reduction = 1;315 if (remembercoeffs <> 0)316 {317 v = v  c * gen(i);318 531 } 319 532 } … … 338 551 poly f = 5xy+7y+3; 339 552 poly g = 7x+5y+3; 340 linreduce(g,I); 341 linreduce(g,I,0,1); 342 linreduce(f,I,1); 343 } 344 345 proc ncsolve (ideal I, list #) 346 "USAGE: ncsolve(I[,s]); I an ideal, s an optional int 347 RETURN: coefficient vector of a linear combination of 0 in the elements of I 348 PURPOSE: compute a linear dependency between the elements of an ideal if such one exists 553 linReduce(g,I); 554 linReduce(g,I,0,1); 555 linReduce(f,I,1); 556 f = x3 + y2 + x2 + y + x; 557 I = x3  y3, y3  x2, x3  y2, x2  y, y2x; 558 list l = linReduce(f, I, 1); 559 l; 560 module M = I; 561 f  (l[1]  (M*l[2])[1,1]); 562 } 563 564 proc linSyzSolve (ideal I, list #) 565 "USAGE: linSyzSolve(I[,s]); I an ideal, s an optional int 566 RETURN: vector, coefficient vector of a linear combination of 0 in the elements of I 567 PURPOSE: compute a linear dependency between the elements of an ideal 568 @* if such one exists 349 569 NOTE: If s<>0, @code{std} is used for Groebner basis computations, 350 570 @* otherwise, @code{slimgb} is used. 351 571 @* By default, @code{slimgb} is used in char 0 and @code{std} in char >0. 352 @*If printlevel=1, progress debug messages will be printed,572 DISPLAY: If printlevel=1, progress debug messages will be printed, 353 573 @* if printlevel>=2, all the debug messages will be printed. 354 EXAMPLE: example ncsolve; shows examples574 EXAMPLE: example linSyzSolve; shows examples 355 575 " 356 576 { … … 374 594 else 375 595 { 376 // 1. introduce undefined coeffs596 //  1. introduce undefined coeffs  377 597 def save = basering; 378 598 int p = char(save); … … 386 606 ring @A = p,(@a(1..sI)),lp; 387 607 ring @aA = (p,@a(1..sI)), (@z),dp; 608 // TODO: BUG! WHAT IF THE ORIGINAL RING ALREADY HAS SUCH VARIABLES/PARAMETERS!!!? 609 // TODO: YOU CAN OVERCOME THIS BY MEANS SIMILAR TO "chooseSafeVarName" FROM NEW "matrix.lib" 388 610 def @B = save + @aA; 389 611 setring @B; 390 612 ideal I = imap(save,I); 391 // 2. form the linear system for the undef coeffs613 //  2. form the linear system for the undef coeffs  392 614 int i; poly W; ideal QQ; 393 615 for (i=1; i<=sI; i++) … … 403 625 setring @A; 404 626 ideal QQ = imap(@B,QQ); 405 // 3. this QQ is a polynomial ideal, so "solve" the system406 dbprint(ppl, " ncsolve: starting Groebner basis computation with engine:", whichengine);627 //  3. this QQ is a polynomial ideal, so "solve" the system  628 dbprint(ppl, "linSyzSolve: starting Groebner basis computation with engine:", whichengine); 407 629 QQ = engine(QQ,whichengine); 408 630 dbprint(ppl, "QQ after engine:", QQ); 409 631 if (dim(QQ) == 1) 410 632 { 411 dbprint(ppl+1, "no solutions by ncsolve");633 dbprint(ppl+1, "no solutions by linSyzSolve"); 412 634 // output zeroes 413 635 setring save; … … 415 637 return(v); 416 638 } 417 // 4. in order to get the numeric values639 //  4. in order to get the numeric values  418 640 matrix AA = matrix(maxideal(1)); 419 641 module MQQ = std(module(QQ)); … … 439 661 ring r = 0,x,dp; 440 662 ideal I = x,2x; 441 ncsolve(I);663 linSyzSolve(I); 442 664 ideal J = x,x2; 443 ncsolve(J); 444 } 445 446 proc pintersect (poly s, ideal I) 447 "USAGE: pintersect(f, I); f a poly, I an ideal 448 RETURN: coefficient vector of the monic generator of the intersection of the ideal generated by f with I 449 PURPOSE: compute the intersection of an ideal with a principal ideal defined by f 665 linSyzSolve(J); 666 } 667 668 proc pIntersect (poly s, ideal I) 669 "USAGE: pIntersect(f, I); f a poly, I an ideal 670 RETURN: vector, coefficient vector of the monic polynomial 671 PURPOSE: compute the intersection of ideal I with the subalgebra K[f] 672 ASSUME: I is given as Groebner basis. 450 673 NOTE: If the intersection is zero, this proc might not terminate. 451 @* I should be given as standard basis. 452 @* If printlevel=1, progress debug messages will be printed, 674 DISPLAY: If printlevel=1, progress debug messages will be printed, 453 675 @* if printlevel>=2, all the debug messages will be printed. 454 EXAMPLE: example p intersect; shows examples676 EXAMPLE: example pIntersect; shows examples 455 677 " 456 678 { 457 679 // assume I is given in Groebner basis 458 attrib(I,"isSB",1); // set attribute for suppressing NF messages 680 if (attrib(I,"isSB") <> 1) 681 { 682 print("WARNING: The input has no SB attribute!"); 683 print(" Treating it as if it were a Groebner basis and proceeding..."); 684 attrib(I,"isSB",1); // set attribute for suppressing NF messages 685 } 459 686 int ppl = printlevelvoice+2; 460 687 // case 1: I = basering 461 688 if (size(I) == 1) 462 689 { 463 if (simplify(I, 1) == ideal(1))690 if (simplify(I,2)[1] == 1) 464 691 { 465 692 return(gen(2)); // = s … … 512 739 poly toNF,tobracket,newNF,rednewNF,oldNF,secNF; 513 740 i = 1; 514 dbprint(ppl+1,"p intersect starts...");741 dbprint(ppl+1,"pIntersect starts..."); 515 742 while (1) 516 743 { … … 535 762 secNF = newNF; 536 763 } 537 ll = lin reduce(newNF,redNI,1);764 ll = linReduce(newNF,redNI,1); 538 765 rednewNF = ll[1]; 539 766 l[i+1] = ll[2]; … … 578 805 v = imap(@R,v); 579 806 kill @R; 580 dbprint(ppl+1,"p intersect finished");807 dbprint(ppl+1,"pIntersect finished"); 581 808 return(v); 582 809 } … … 586 813 ring r = 0,(x,y),dp; 587 814 poly f = x^2+y^3+x*y^2; 588 def D = initial malgrange(f);815 def D = initialMalgrange(f); 589 816 setring D; 590 817 inF; 591 pintersect(t*Dt,inF); 592 } 593 594 proc pintersectsyz (poly s, ideal II, list #) 595 "USAGE: pintersectsyz(f, I [,p,s,t]); f a poly, I an ideal, p, t optial ints, p a prime number 596 RETURN: coefficient vector of the monic generator of the intersection of the ideal generated by f with I 597 PURPOSE: compute the intersection of an ideal with a principal ideal defined by f 598 NOTE: If the intersection is zero, this proc might not terminate. 599 @* I should be given as standard basis. 600 @* If p>0 is given, this proc computes the generator of the intersection in char p first and 601 @* then only searches for a generator of the obtained degree in the basering. 602 @* Otherwise, it searched for all degrees. 603 @* This is done by computing syzygies. 818 pIntersect(t*Dt,inF); 819 } 820 821 proc pIntersectSyz (poly s, ideal II, list #) 822 "USAGE: pIntersectSyz(f, I [,p,s,t]); f a poly, I an ideal, p, t optial ints, p a prime number 823 RETURN: vector, coefficient vector of the monic polynomial 824 PURPOSE: compute the intersection of an ideal I with the subalgebra K[f] 825 ASSUME: I is given as Groebner basis. 826 NOTE: If the intersection is zero, this procedure might not terminate. 827 @* If p>0 is given, this proc computes the generator of the intersection in char p first 828 @* and then only searches for a generator of the obtained degree in the basering. 829 @* Otherwise, it searched for all degrees by computing syzygies. 604 830 @* If s<>0, @code{std} is used for Groebner basis computations in char 0, 605 831 @* otherwise, and by default, @code{slimgb} is used. 606 832 @* If t<>0 and by default, @code{std} is used for Groebner basis computations in char >0, 607 833 @* otherwise, @code{slimgb} is used. 608 @*If printlevel=1, progress debug messages will be printed,834 DISPLAY: If printlevel=1, progress debug messages will be printed, 609 835 @* if printlevel>=2, all the debug messages will be printed. 610 EXAMPLE: example p intersectsyz; shows examples836 EXAMPLE: example pIntersectSyz; shows examples 611 837 " 612 838 { 613 839 // assume I is given in Groebner basis 614 840 ideal I = II; 615 attrib(I,"isSB",1); // set attribute for suppressing NF messages 841 if (attrib(I,"isSB") <> 1) 842 { 843 print("WARNING: The input has no SB attribute!"); 844 print(" Treating it as if it were a Groebner basis and proceeding..."); 845 attrib(I,"isSB",1); // set attribute for suppressing NF messages 846 } 616 847 int ppl = printlevelvoice+2; 617 848 int whichengine = 0; // default … … 668 899 } 669 900 i = 1; 670 dbprint(ppl+1,"p intersectsyz starts...");901 dbprint(ppl+1,"pIntersectSyz starts..."); 671 902 dbprint(ppl+1,"with ideal I=", I); 672 903 while (1) … … 688 919 } 689 920 // look for a solution 690 dbprint(ppl," ncsolve starts with: "+string(matrix(NI)));921 dbprint(ppl,"linSyzSolve starts with: "+string(matrix(NI))); 691 922 if (solveincharp<>0) // modular method 692 923 { 693 924 setring Rp; 694 925 NI[i+1] = phi(newNF); 695 v = ncsolve(NI,modengine);926 v = linSyzSolve(NI,modengine); 696 927 if (v!=0) // there is a modular solution 697 928 { 698 929 dbprint(ppl,"got solution in char ",solveincharp," of degree " ,i); 699 930 setring save; 700 v = ncsolve(NI,whichengine);931 v = linSyzSolve(NI,whichengine); 701 932 if (v==0) 702 933 { … … 712 943 else // nonmodular method 713 944 { 714 v = ncsolve(NI,whichengine);945 v = linSyzSolve(NI,whichengine); 715 946 } 716 947 matrix MM[1][nrows(v)] = matrix(v); 717 dbprint(ppl," ncsolve ready with: "+string(MM));948 dbprint(ppl,"linSyzSolve ready with: "+string(MM)); 718 949 kill MM; 719 // " ncsolve ready with"; print(v);950 // "linSyzSolve ready with"; print(v); 720 951 if (v!=0) 721 952 { … … 729 960 if (p!=0) 730 961 { 731 dbprint(ppl," ncsolve: bad solution!");962 dbprint(ppl,"linSyzSolve: bad solution!"); 732 963 } 733 964 else 734 965 { 735 dbprint(ppl," ncsolve: got solution!");966 dbprint(ppl,"linSyzSolve: got solution!"); 736 967 // "got solution!"; 737 968 break; … … 741 972 i++; 742 973 } 743 dbprint(ppl+1,"p intersectsyz finished");974 dbprint(ppl+1,"pIntersectSyz finished"); 744 975 return(v); 745 976 } … … 749 980 ring r = 0,(x,y),dp; 750 981 poly f = x^2+y^3+x*y^2; 751 def D = initial malgrange(f);982 def D = initialMalgrange(f); 752 983 setring D; 753 984 inF; 754 985 poly s = t*Dt; 755 p intersectsyz(s,inF);986 pIntersectSyz(s,inF); 756 987 int p = prime(20000); 757 p intersectsyz(s,inF,p,0,0);988 pIntersectSyz(s,inF,p,0,0); 758 989 } 759 990 760 991 proc vec2poly (list #) 761 992 "USAGE: vec2poly(v [,i]); v a vector or an intvec, i an optional int 762 RETURN: a poly with coefficient vector v 763 PURPOSE: convert a coefficient vector to a poly 764 NOTE: If i>0 is given, the returned poly is an element of K[var(i)], 765 @* otherwise, and by default, @code{i=1} is used. 766 @* The first entry of v is the coefficient of 1. 993 RETURN: poly, an univariate poly in ith variable with coefficients given by v 994 PURPOSE: constructs an univariate poly in K[var(i)] with given coefficients, 995 @* such that the coefficient at var(i)^{j1} is v[j]. 996 NOTE: The optional argument i must be positive, by default i is 1. 767 997 EXAMPLE: example vec2poly; shows examples 768 998 " … … 858 1088 } 859 1089 860 static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodp intersect, int pintersectchar, int modengine, intvec u0)1090 static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodpIntersect, int pIntersectchar, int modengine, intvec u0) 861 1091 { 862 1092 int ppl = printlevel  voice +2; … … 866 1096 if (inorann == 0) // bfct using initial ideal 867 1097 { 868 def D = initial malgrange(f,whichengine,methodord,1,u0);1098 def D = initialMalgrange(f,whichengine,methodord,1,u0); 869 1099 setring D; 870 1100 ideal J = inF; … … 881 1111 vector b; 882 1112 // try it modular 883 if (methodp intersect <> 0) // pintersectsyz884 { 885 if (p intersectchar == 0) // pintersectsyz::modular1113 if (methodpIntersect <> 0) // pIntersectSyz 1114 { 1115 if (pIntersectchar == 0) // pIntersectSyz::modular 886 1116 { 887 1117 int lb = 30000; … … 893 1123 dbprint(ppl,"number of run in the loop: "+string(i)); 894 1124 int q = prime(random(lb,ub)); 895 if ( isin(usedprimes,q)==0) // if q was not already used1125 if (findFirst(usedprimes,q)==0) // if q was not already used 896 1126 { 897 1127 usedprimes = usedprimes,q; 898 1128 dbprint(ppl,"used prime is: "+string(q)); 899 b = p intersectsyz(s,J,q,whichengine,modengine);1129 b = pIntersectSyz(s,J,q,whichengine,modengine); 900 1130 } 901 1131 i++; 902 1132 } 903 1133 } 904 else // p intersectsyz::nonmodular905 { 906 b = p intersectsyz(s,J,0,whichengine);907 } 908 } 909 else // p intersect: linreduce910 { 911 b = p intersect(s,J);1134 else // pIntersectSyz::nonmodular 1135 { 1136 b = pIntersectSyz(s,J,0,whichengine); 1137 } 1138 } 1139 else // pIntersect: linReduce 1140 { 1141 b = pIntersect(s,J); 912 1142 } 913 1143 setring save; … … 926 1156 proc bfct (poly f, list #) 927 1157 "USAGE: bfct(f [,s,t,v]); f a poly, s,t optional ints, v an optional intvec 928 RETURN: list of roots of the BernsteinSato polynomial bs(f) and their multiplicies 929 PURPOSE: compute the global BernsteinSato polynomial for a hypersurface, defined by f, according to the algorithm by Masayuki Noro 930 NOTE: In this proc, the initial Malgrange ideal is computed. 931 @* Further, a system of linear equations is solved by linear reductions. 932 @* If s<>0, @code{std} is used for Groebner basis computations, 1158 RETURN: list of ideal and intvec 1159 PURPOSE: computes the roots of the BernsteinSato polynomial b(s) 1160 @* for the hypersurface defined by f. 1161 ASSUME: The basering is a commutative polynomial ring in char 0. 1162 BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm 1163 @* by Masayuki Noro and then a system of linear equations is solved by linear reductions. 1164 NOTE: In the output list, the ideal contains all the roots 1165 @* and the intvec their multiplicities. 1166 @* If s<>0, @code{std} is used for GB computations, 933 1167 @* otherwise, and by default, @code{slimgb} is used. 934 1168 @* If t<>0, a matrix ordering is used for Groebner basis computations, 935 1169 @* otherwise, and by default, a block ordering is used. 936 @* If v is a positive weight vector, v is used for homogenization computations, 1170 @* If v is a positive weight vector, v is used for homogenization computations, 937 1171 @* otherwise and by default, no weights are used. 938 @*If printlevel=1, progress debug messages will be printed,1172 DISPLAY: If printlevel=1, progress debug messages will be printed, 939 1173 @* if printlevel>=2, all the debug messages will be printed. 940 1174 EXAMPLE: example bfct; shows examples … … 965 1199 if (size(#)>2) 966 1200 { 967 if (typeof(#[3])=="intvec" && size(#[3])==n && ispositive(#[3])==1)1201 if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1) 968 1202 { 969 1203 u0 = #[3]; … … 985 1219 } 986 1220 987 proc bfctsyz (poly f, list #) 988 "USAGE: bfctsyz(f [,r,s,t,u,v]); f a poly, r,s,t,u optional ints, v an optional intvec 989 RETURN: list of roots of the BernsteinSato polynomial bs(f) and its multiplicies 990 PURPOSE: compute the global BernsteinSato polynomial for a hypersurface, defined by f, according to the algorithm by Masayuki Noro 991 NOTE: In this proc, the initial Malgrange ideal is computed. 992 @* Further, a system of linear equations is solved by computing syzygies. 993 @* If r<>0, @code{std} is used for Groebner basis computations in characteristic 0, 1221 proc bfctSyz (poly f, list #) 1222 "USAGE: bfctSyz(f [,r,s,t,u,v]); f a poly, r,s,t,u optional ints, v an optional intvec 1223 RETURN: list of ideal and intvec 1224 PURPOSE: computes the roots of the BernsteinSato polynomial b(s) 1225 @* for the hypersurface defined by f 1226 ASSUME: The basering is a commutative polynomial ring in char 0. 1227 BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm 1228 @* by Masayuki Noro and then a system of linear equations is solved by computing syzygies. 1229 NOTE: In the output list, the ideal contains all the roots 1230 @* and the intvec their multiplicities. 1231 @* If r<>0, @code{std} is used for GB computations in characteristic 0, 994 1232 @* otherwise, and by default, @code{slimgb} is used. 995 @* If s<>0, a matrix ordering is used for G roebner basis computations,996 @* otherwise,and by default, a block ordering is used.997 @* If t<>0, the computation of the intersection is solely performed over charasteristic 0,998 @* otherwise and by default, a modular method is used.999 @* If u<>0 and by default, @code{std} is used for G roebner basis computations in characteristic >0,1000 @* otherwise, @code{slimgb} is used.1001 @* If v is a positive weight vector, v is used for homogenization computations,1002 @* otherwise and by default, no weights are used.1003 @*If printlevel=1, progress debug messages will be printed,1233 @* If s<>0, a matrix ordering is used for GB computations, otherwise, 1234 @* and by default, a block ordering is used. 1235 @* If t<>0, the computation of the intersection is solely performed over 1236 @* charasteristic 0, otherwise and by default, a modular method is used. 1237 @* If u<>0 and by default, @code{std} is used for GB computations in 1238 @* characteristic >0, otherwise, @code{slimgb} is used. 1239 @* If v is a positive weight vector, v is used for homogenization 1240 @* computations, otherwise and by default, no weights are used. 1241 DISPLAY: If printlevel=1, progress debug messages will be printed, 1004 1242 @* if printlevel>=2, all the debug messages will be printed. 1005 1243 EXAMPLE: example bfct; shows examples … … 1018 1256 int whichengine = 0; // default 1019 1257 int methodord = 0; // default 1020 int p intersectchar = 0; // default1258 int pIntersectchar = 0; // default 1021 1259 int modengine = 1; // default 1022 1260 intvec u0 = 0; // default … … 1037 1275 if (typeof(#[3])=="int"  typeof(#[3])=="number") 1038 1276 { 1039 p intersectchar = int(#[3]);1277 pIntersectchar = int(#[3]); 1040 1278 } 1041 1279 if (size(#)>3) … … 1047 1285 if (size(#)>4) 1048 1286 { 1049 if (typeof(#[5])=="intvec" && size(#[5])==n && ispositive(#[5])==1)1287 if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1) 1050 1288 { 1051 1289 u0 = #[5]; … … 1056 1294 } 1057 1295 } 1058 list b = bfctengine(f,0,whichengine,methodord,1,p intersectchar,modengine,u0);1296 list b = bfctengine(f,0,whichengine,methodord,1,pIntersectchar,modengine,u0); 1059 1297 return(b); 1060 1298 } … … 1064 1302 ring r = 0,(x,y),dp; 1065 1303 poly f = x^2+y^3+x*y^2; 1066 bfct syz(f);1304 bfctSyz(f); 1067 1305 intvec v = 3,2; 1068 bfctsyz(f,0,1,1,0,v); 1069 } 1070 1071 proc bfctideal (ideal I, intvec w, list #) 1072 "USAGE: bfctideal(I,w[,s,t]); I an ideal, w an intvec, s,t optional ints 1073 RETURN: list of roots and their multiplicies of the global bfunction of I w.r.t. the weight vector (w,w) 1074 PURPOSE: compute the global bfunction of an ideal according to the algorithm by M. Noro 1075 NOTE: Assume, I is an ideal in the nth Weyl algebra where the sequence of the 1076 @* variables is x(1),...,x(n),D(1),...,D(n). 1077 @* If s<>0, @code{std} is used for Groebner basis computations in characteristic 0, 1306 bfctSyz(f,0,1,1,0,v); 1307 } 1308 1309 proc bfctIdeal (ideal I, intvec w, list #) 1310 "USAGE: bfctIdeal(I,w[,s,t]); I an ideal, w an intvec, s,t optional ints 1311 RETURN: list of ideal and intvec 1312 PURPOSE: computes the roots of the global bfunction of I wrt the weight vector (w,w). 1313 ASSUME: The basering is the nth Weyl algebra in char 0, where the sequence of 1314 @* the variables is x(1),...,x(n),D(1),...,D(n). 1315 BACKGROUND: In this proc, the initial ideal of I is computed according to the algorithm by 1316 @* Masayuki Noro and then a system of linear equations is solved by linear reductions. 1317 NOTE: In the output list, the ideal contains all the roots 1318 @* and the intvec their multiplicities. 1319 @* If s<>0, @code{std} is used for GB computations in characteristic 0, 1078 1320 @* otherwise, and by default, @code{slimgb} is used. 1079 @* If t<>0, a matrix ordering is used for G roebner basis computations,1080 @* otherwise,and by default, a block ordering is used.1081 @*If printlevel=1, progress debug messages will be printed,1321 @* If t<>0, a matrix ordering is used for GB computations, otherwise, 1322 @* and by default, a block ordering is used. 1323 DISPLAY: If printlevel=1, progress debug messages will be printed, 1082 1324 @* if printlevel>=2, all the debug messages will be printed. 1083 1325 EXAMPLE: example bfctideal; shows examples … … 1104 1346 } 1105 1347 } 1106 ideal J = initial ideal(I,w,w,whichengine,methodord);1348 ideal J = initialIdeal(I,w,w,whichengine,methodord); 1107 1349 poly s; 1108 1350 for (i=1; i<=n; i++) … … 1110 1352 s = s + w[i]*var(i)*var(n+i); 1111 1353 } 1112 vector b = p intersect(s,J);1354 vector b = pIntersect(s,J); 1113 1355 list l = listofroots(b,0); 1114 1356 return(l); … … 1120 1362 def D = Weyl(); 1121 1363 setring D; 1122 ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6;1364 ideal I = std(3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6); I; 1123 1365 intvec w1 = 1,1; 1124 1366 intvec w2 = 1,2; 1125 1367 intvec w3 = 2,3; 1126 bfctideal(I,w1); 1127 bfctideal(I,w2,1); 1128 bfctideal(I,w3,0,1); 1129 } 1130 1131 1132 proc bfctonestep (poly f,list #) 1133 "USAGE: bfctonestep(f [,s,t]); f a poly, s,t optional ints 1134 RETURN: list of roots of the BernsteinSato polynomial bs(f) and its multiplicies 1135 PURPOSE: compute the global BernsteinSato polynomial for a hypersurface, defined by f, using only one Groebner basis computation 1136 NOTE: If s<>0, @code{std} is used for the Groebner basis computation, otherwise, 1368 bfctIdeal(I,w1); 1369 bfctIdeal(I,w2,1); 1370 bfctIdeal(I,w3,0,1); 1371 } 1372 1373 proc bfctOneGB (poly f,list #) 1374 "USAGE: bfctOneGB(f [,s,t]); f a poly, s,t optional ints 1375 RETURN: list of ideal and intvec 1376 PURPOSE: computes the roots of the BernsteinSato polynomial b(s) for the 1377 @* hypersurface defined by f, using only one GB computation 1378 ASSUME: The basering is a commutative polynomial ring in char 0. 1379 BACKGROUND: In this proc, the initial Malgrange ideal is computed based on the 1380 @* algorithm by Masayuki Noro and combined with an elimination ordering. 1381 NOTE: In the output list, the ideal contains all the roots 1382 @* and the intvec their multiplicities. 1383 @* If s<>0, @code{std} is used for the GB computation, otherwise, 1137 1384 @* and by default, @code{slimgb} is used. 1138 @* If t<>0, a matrix ordering is used for G roebner basiscomputations,1385 @* If t<>0, a matrix ordering is used for GB computations, 1139 1386 @* otherwise, and by default, a block ordering is used. 1140 @*If printlevel=1, progress debug messages will be printed,1387 DISPLAY: If printlevel=1, progress debug messages will be printed, 1141 1388 @* if printlevel>=2, all the debug messages will be printed. 1142 EXAMPLE: example bfct onestep; shows examples1389 EXAMPLE: example bfctOneGB; shows examples 1143 1390 " 1144 1391 { … … 1146 1393 def save = basering; 1147 1394 int n = nvars(save); 1395 int noofvars = 2*n+4; 1148 1396 int i; 1149 1397 int whichengine = 0; // default … … 1163 1411 } 1164 1412 } 1165 def DDh = initialidealengine("bfctonestep", whichengine, methodord, f); 1413 intvec uv; 1414 uv[n+3] = 1; 1415 ring r = 0,(x(1..n)),dp; 1416 poly f = fetch(save,f); 1417 uv[1] = 1; uv[noofvars] = 0; 1418 // for the ordering 1419 intvec @a; @a = 1:noofvars; @a[2] = 2; 1420 intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0; 1421 if (methodord == 0) // default: block ordering 1422 { 1423 ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(noofvars1),lp(1)); 1424 } 1425 else // M() ordering 1426 { 1427 intmat @Ord[noofvars][noofvars]; 1428 @Ord[1,1..noofvars] = uv; 1429 @Ord[2,1..noofvars] = 1:(noofvars1); 1430 for (i=1; i<=noofvars2; i++) 1431 { 1432 @Ord[2+i,noofvars  i] = 1; 1433 } 1434 dbprint(ppl,"weights for ordering:",transpose(@a)); 1435 dbprint(ppl,"the ordering matrix:",@Ord); 1436 ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord)); 1437 } 1438 dbprint(ppl,"the ring Dh:",Dh); 1439 // noncommutative relations 1440 matrix @relD[noofvars][noofvars]; 1441 @relD[1,2] = t*h^2;// s*t = t*s+t*h^2 1442 @relD[2,n+3] = Dt*h^2;// Dt*s = s*Dt+h^2 1443 @relD[1,n+3] = h^2; 1444 for (i=1; i<=n; i++) 1445 { 1446 @relD[i+2,n+3+i] = h^2; 1447 } 1448 dbprint(ppl,"nc relations:",@relD); 1449 def DDh = nc_algebra(1,@relD); 1166 1450 setring DDh; 1167 dbprint(ppl, "the initial ideal:", string(matrix(inF))); 1451 dbprint(ppl,"computing in ring",DDh); 1452 ideal I; 1453 poly f = imap(r,f); 1454 kill r; 1455 f = homog(f,h); 1456 I = t  f, t*Dt  s; // defining the Malgrange ideal 1457 for (i=1; i<=n; i++) 1458 { 1459 I = I, D(i)+diff(f,x(i))*Dt; 1460 } 1461 dbprint(ppl, "starting Groebner basis computation with engine:", whichengine); 1462 I = engine(I, whichengine); 1463 dbprint(ppl, "finished Groebner basis computation"); 1464 dbprint(ppl, "I before dehomogenization is" ,I); 1465 I = subst(I,h,1); // dehomogenization 1466 dbprint(ppl, "I after dehomogenization is" ,I); 1467 I = inForm(I,uv); // we are only interested in the initial form w.r.t. uv 1468 dbprint(ppl, "the initial ideal:", string(matrix(I))); 1168 1469 intvec tonselect = 1; 1169 for (i=3; i<=2*n+4; i++) 1170 { 1171 tonselect = tonselect,i; 1172 } 1173 inF = nselect(inF,tonselect); 1174 dbprint(ppl, "generators containing only s:", string(matrix(inF))); 1175 inF = engine(inF, whichengine); // is now a principal ideal; 1470 for (i=3; i<=2*n+4; i++) { tonselect = tonselect,i; } 1471 I = nselect(I,tonselect); 1472 dbprint(ppl, "generators containing only s:", string(matrix(I))); 1473 I = engine(I, whichengine); // is now a principal ideal; 1176 1474 setring save; 1177 1475 ideal J; J[2] = var(1); 1178 1476 map @m = DDh,J; 1179 ideal inF = @m(inF);1180 poly p = inF[1];1477 ideal I = @m(I); 1478 poly p = I[1]; 1181 1479 list l = listofroots(p,1); 1182 1480 return(l); … … 1187 1485 ring r = 0,(x,y),dp; 1188 1486 poly f = x^2+y^3+x*y^2; 1189 bfctonestep(f); 1190 bfctonestep(f,1,1); 1191 } 1192 1193 proc bfctann (poly f, list #) 1194 "USAGE: bfctann(f [,r]); f a poly, r an optional int 1195 RETURN: list of roots of the BernsteinSato polynomial bs(f) and their multiplicies 1196 PURPOSE: compute the global BernsteinSato polynomial for a hypersurface, defined by f 1197 NOTE: In this proc, ann(f^s) is computed. 1198 @* If r<>0, @code{std} is used for Groebner basis computations, 1487 bfctOneGB(f); 1488 bfctOneGB(f,1,1); 1489 } 1490 1491 proc bfctAnn (poly f, list #) 1492 "USAGE: bfctAnn(f [,r]); f a poly, r an optional int 1493 RETURN: list of ideal and intvec 1494 PURPOSE: computes the roots of the BernsteinSato polynomial b(s) 1495 @* for the hypersurface defined by f 1496 ASSUME: The basering is a commutative polynomial ring in char 0. 1497 BACKGROUND: In this proc, ann(f^s) is computed and then a system of linear 1498 @* equations is solved by linear reductions. 1499 NOTE: In the output list, the ideal contains all the roots 1500 @* and the intvec their multiplicities. 1501 @* If r<>0, @code{std} is used for GB computations, 1199 1502 @* otherwise, and by default, @code{slimgb} is used. 1200 @*If printlevel=1, progress debug messages will be printed,1503 DISPLAY: If printlevel=1, progress debug messages will be printed, 1201 1504 @* if printlevel>=2, all the debug messages will be printed. 1202 1505 EXAMPLE: example bfctann; shows examples … … 1221 1524 ring r = 0,(x,y),dp; 1222 1525 poly f = x^2+y^3+x*y^2; 1223 bfctann(f); 1224 } 1225 1226 static proc hardexamples () 1526 bfctAnn(f); 1527 } 1528 1529 /* 1530 //static proc hardexamples () 1227 1531 { 1228 1532 // some hard examples … … 1246 1550 bfct(xyzreiffen45); 1247 1551 } 1248 1552 */ 
Singular/LIB/dmod.lib
r6045c3 r565e86 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: dmod.lib,v 1.3 2 20081209 16:50:21 levandov Exp $";2 version="$Id: dmod.lib,v 1.33 20081223 21:39:31 levandov Exp $"; 3 3 category="Noncommutative"; 4 4 info=" … … 27 27 @* OT) the classical Ann F^s algorithm from Oaku and Takayama (J. Pure 28 28 Applied Math., 1999), 29 @* LOT) Levandovskyy's modification of the OakuTakayama algorithm ( unpublished)29 @* LOT) Levandovskyy's modification of the OakuTakayama algorithm (ISSAC 2007) 30 30 @* BM) the Ann F^s algorithm by Briancon and Maisonobe (Remarques sur 31 31 l'ideal de Bernstein associe a des polynomes, preprint, 2002) 32 32 33 33 GUIDE: 34 @*  Ann F^s = I = I(F^s) = LD in D(R)[s] can be computed by SannfsBM SannfsOT, SannfsLOT34 @*  Ann F^s = I = I(F^s) = LD in D(R)[s] can be computed by SannfsBM, SannfsOT, SannfsLOT 35 35 @*  Ann^(1) F^s in D(R)[s] can be computed by Sannfslog 36 36 @*  global Bernstein polynomial bs resp. BS in K[s] can be computed by bernsteinBM 37 @*  Ann F^s0 = I(F^s0) = LD0 in D(R) can be computed by annfs0, annfsBM, annfsOT, annfsLOT 37 @*  Ann F^s0 = I(F^s0) = LD0 in D(R) can be computed by annfs0, annfsBM, annfsOT, annfsLOT, annfs2 38 38 @*  all the relevant data (LD, LD0, bs, PS) are computed by operatorBM 39 39 … … 49 49 annfsBMI(F[,eng]); compute Ann F^s and Bernstein ideal for a poly F=f1*..*fP (multivariate algorithm of BrianconMaisonobe) 50 50 checkRoot(F,a[,S,eng]); check if a given rational is a root of the global Bernstein polynomial of F and compute its multiplicity 51 52 SECONDARY PROCEDURES FOR DMODULES:53 54 51 annfsBM(F[,eng]); compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of BrianconMaisonobe) 55 52 annfsLOT(F[,eng]); compute Ann F^s0 in D and Bernstein poly for a poly F (Levandovskyy modification of the OakuTakayama algorithm) 56 53 annfsOT(F[,eng]); compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of OakuTakayama) 57 54 SannfsBM(F[,eng]); compute Ann F^s in D[s] for a poly F (algorithm of BrianconMaisonobe) 58 SannfsBFCT(F[,eng]); compute Ann F^s in D[s] for a poly F (algorithm of BrianconMaisonobe )55 SannfsBFCT(F[,eng]); compute Ann F^s in D[s] for a poly F (algorithm of BrianconMaisonobe, other output ordering) 59 56 SannfsLOT(F[,eng]); compute Ann F^s in D[s] for a poly F (Levandovskyy modification of the OakuTakayama algorithm) 60 57 SannfsOT(F[,eng]); compute Ann F^s in D[s] for a poly F (algorithm of OakuTakayama) 61 annfs0(I,F [,eng]); compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s] 58 annfs0(I,F [,eng]); compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s] 59 annfs2(I,F [,eng]); compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s] by using a trick of Noro 60 annfsRB(I,F [,eng]); compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s] by reduceB strategy of Macaulay 62 61 checkRoot1(I,F,a[,eng]); check whether a rational is a root of the global Bernstein polynomial of F from the known Ann F^s in D[s] 63 62 checkRoot2(I,F,a[,eng]); check whether a rational is a root of the global Bernstein polynomial of F and compute its multiplicity from the known Ann F^s in D[s] 64 checkFactor(I,F,q s[,eng]); check whether a polynomial qsin K[s] is a factor of the global Bernstein polynomial of F from the known Ann F^s in D[s]63 checkFactor(I,F,q[,eng]); check whether a polynomial q in K[s] is a factor of the global Bernstein polynomial of F from the known Ann F^s in D[s] 65 64 66 65 AUXILIARY PROCEDURES: … … 72 71 minIntRoot(P,fact); minimal integer root among the roots of a maximal ideal P 73 72 varnum(s); the number of the variable with the name s 73 74 74 75 75 SEE ALSO: gmssing_lib, bfct_lib, dmodapp_lib … … 92 92 "MAIN PROCEDURES:"; 93 93 example annfs; 94 example annfs0;95 94 example Sannfs; 96 95 example Sannfslog; … … 101 100 example annfsBMI; 102 101 example checkRoot; 102 example annfs0; 103 example annfs2; 104 example annfsRB; 105 example annfs2; 103 106 "SECONDARY DMOD PROCEDURES:"; 104 107 example annfsBM; … … 108 111 example SannfsLOT; 109 112 example SannfsOT; 113 example SannfsBFCT; 110 114 example checkRoot1; 111 115 example checkRoot2; … … 130 134 @*  the ideal LD (which is a Groebner basis) is the needed Dmodule structure, 131 135 @*  the list BS is the list of roots and multiplicities of a Bernstein polynomial of f. 132 @*If @code{printlevel}=1, progress debug messages will be printed,133 @* if @code{printlevel}>=2, all the debug messages will be printed.136 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 137 @* if @code{printlevel}>=2, all the debug messages will be printed. 134 138 EXAMPLE: example annfs; shows examples 135 139 " … … 252 256 @* In the output ring: 253 257 @*  the ideal LD is the needed Dmodule structure. 254 @*If @code{printlevel}=1, progress debug messages will be printed,255 @* if @code{printlevel}>=2, all the debug messages will be printed.258 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 259 @* if @code{printlevel}>=2, all the debug messages will be printed. 256 260 EXAMPLE: example Sannfs; shows examples 257 261 " … … 365 369 @* If eng <>0, @code{std} is used for Groebner basis computations, 366 370 @* otherwise, and by default @code{slimgb} is used. 367 @* If printlevel=1, progress debug messages will be printed,368 @* if printlevel>=2, all the debug messages will be printed.371 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 372 @* if @code{printlevel}>=2, all the debug messages will be printed. 369 373 EXAMPLE: example Sannfslog; shows examples 370 374 " … … 506 510 // alternative code for SannfsBM, renamed from annfsBM to ALTannfsBM 507 511 // is superfluos  will not be included in the official documentation 508 proc ALTannfsBM (poly F, list #)512 static proc ALTannfsBM (poly F, list #) 509 513 "USAGE: ALTannfsBM(f [,eng]); f a poly, eng an optional int 510 514 RETURN: ring … … 514 518 @* If eng <>0, @code{std} is used for Groebner basis computations, 515 519 @* otherwise, and by default @code{slimgb} is used. 516 @* If printlevel=1, progress debug messages will be printed,517 @* if printlevel>=2, all the debug messages will be printed.520 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 521 @* if @code{printlevel}>=2, all the debug messages will be printed. 518 522 EXAMPLE: example ALTannfsBM; shows examples 519 523 " … … 703 707 NOTE: If eng <>0, @code{std} is used for Groebner basis computations, 704 708 @* otherwise, and by default @code{slimgb} is used. 705 @* If printlevel=1, progress debug messages will be printed,706 @* if printlevel>=2, all the debug messages will be printed.709 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 710 @* if @code{printlevel}>=2, all the debug messages will be printed. 707 711 EXAMPLE: example bernsteinBM; shows examples 708 712 " … … 933 937 @* If eng <>0, @code{std} is used for Groebner basis computations, 934 938 @* otherwise, and by default @code{slimgb} is used. 935 @* If printlevel=1, progress debug messages will be printed,936 @* if printlevel>=2, all the debug messages will be printed.939 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 940 @* if @code{printlevel}>=2, all the debug messages will be printed. 937 941 EXAMPLE: example annfsBM; shows examples 938 942 " … … 1217 1221 1218 1222 1219 // try to replaces with s1 => data is shorter1223 // replacing s with s1 => data is shorter 1220 1224 // analogue of annfs0 1221 1225 proc annfs2(ideal I, poly F, list #) … … 1223 1227 RETURN: ring 1224 1228 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based on the 1225 output of procedures SannfsBM, SannfsOT or SannfsLOT1229 output of Sannfslike procedure 1226 1230 NOTE: activate this ring with the @code{setring} command. In this ring, 1227 1231 @*  the ideal LD (which is a Groebner basis) is the annihilator of f^s, … … 1229 1233 @* If eng <>0, @code{std} is used for Groebner basis computations, 1230 1234 @* otherwise and by default @code{slimgb} is used. 1231 @* Uses the shorter form of expressions in the variable s (the idea of Noro).1232 @* If printlevel=1, progress debug messages will be printed,1233 @* if printlevel>=2, all the debug messages will be printed.1235 @* annfs2 uses the shorter form of expressions in the variable s (the idea of Noro). 1236 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 1237 @* if @code{printlevel}>=2, all the debug messages will be printed. 1234 1238 EXAMPLE: example annfs2; shows examples 1235 1239 " … … 1390 1394 @* If eng <>0, @code{std} is used for Groebner basis computations, 1391 1395 @* otherwise and by default @code{slimgb} is used. 1392 @* Uses the shorter form of expressions in the variable s (the idea of Noro).1393 @* If printlevel=1, progress debug messages will be printed,1394 @* if printlevel>=2, all the debug messages will be printed.1396 @* This procedure follows the 'reduceB' strategy, used in Macaulay2. 1397 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 1398 @* if @code{printlevel}>=2, all the debug messages will be printed. 1395 1399 EXAMPLE: example annfsRB; shows examples 1396 1400 " … … 1561 1565 poly F = x^3+y^3+z^3; 1562 1566 printlevel = 0; 1563 def A = SannfsBM(F); 1564 setring A; 1565 LD; 1567 def A = SannfsBM(F); setring A; 1568 LD; // sparametric ahhinilator 1566 1569 poly F = imap(r,F); 1567 def B = annfsRB(LD,F); 1568 setring B; 1570 def B = annfsRB(LD,F); setring B; 1569 1571 LD; 1570 1572 BS; … … 1584 1586 @* If eng <>0, @code{std} is used for Groebner basis computations, 1585 1587 @* otherwise and by default @code{slimgb} is used. 1586 @* If printlevel=1, progress debug messages will be printed,1587 @* if printlevel>=2, all the debug messages will be printed.1588 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 1589 @* if @code{printlevel}>=2, all the debug messages will be printed. 1588 1590 EXAMPLE: example operatorBM; shows examples 1589 1591 " … … 1889 1891 // for this, do special modulo with emphasis on the 2comp 1890 1892 // of the transp matrix, there must appear 1 in the GB 1893 // then use lift to get the expression... 1891 1894 // need: (c,<) ordering for such comp's 1892 1895 1893 // proc operatorModulo(poly F, ideal I, poly b) 1894 // "USAGE: operatorModulo(f,I,b); f a poly, I an ideal, b a poly 1895 // RETURN: poly 1896 // PURPOSE: compute the Boperator from the poly F, ideal I = Ann f^s and BernsteinSato 1897 // polynomial b using modulo 1898 // NOTE: The computations take place in the ring, similar to the one returned by Sannfs procedure. 1899 // @* If printlevel=1, progress debug messages will be printed, 1900 // @* if printlevel>=2, all the debug messages will be printed. 1901 // EXAMPLE: example operatorModulo; shows examples 1902 // " 1903 // { 1904 // // with hom_kernel; 1905 // matrix AA[1][2] = F,b; 1906 // // matrix M = matrix(subst(I,s,s+1)*freemodule(2)); // ann f^{s+1} 1907 // matrix N = matrix(I); // ann f^s 1908 // // matrix K = hom_kernel(AA,M,N); 1909 // matrix K = modulo(AA,N); 1910 // K = transpose(K); 1911 // print(K); 1912 // ideal J; 1913 // int i; 1914 // poly t; number n; 1915 // for(i=1; i<=nrows(K); i++) 1916 // { 1917 // if (K[i,2]!=0) 1918 // { 1919 // if ( leadmonom(K[i,2]) == 1) 1920 // { 1921 // t = K[i,1]; 1922 // n = leadcoef(K[i,2]); 1923 // t = t/n; 1924 // // J = J, K[i][2]; 1925 // break; 1926 // } 1927 // } 1928 // } 1929 // ideal J = groebner(subst(I,s,s+1)); // for NF 1930 // t = NF(t,J); 1931 // "candidate:"; t; 1932 // J = subst(J,s,s1); 1933 // // test: 1934 // if ( NF(t*Fb,J) !=0) 1935 // { 1936 // "Problem: PS does not work on F"; 1937 // } 1938 // return(t); 1939 // } 1940 // example 1941 // { 1942 // "EXAMPLE:"; echo = 2; 1943 // // ring r = 0,(x,y,z),Dp; 1944 // // poly F = x^3+y^3+z^3; 1945 // LIB "dmod.lib"; option(prot); option(mem); 1946 // ring r = 0,(x,y),Dp; 1947 // // poly F = x^3+y^3+x*y^2; 1948 // poly F = x^4 + y^4 + x*y^4; 1949 // def A = Sannfs(F); // here we get LD = ann f^s 1950 // setring A; 1951 // poly F = imap(r,F); 1952 // def B = annfs0(LD,F); // to obtain BS polynomial 1953 // list BS = imap(B,BS); 1954 // poly b = fl2poly(BS,"s"); 1955 // LD = groebner(LD); 1956 // poly PS = operatorModulo(F,LD,b); 1957 // PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD 1958 // reduce(PS*Fbs,LD); // check the property of PS 1959 // } 1960 1961 proc fl2poly(list L, string s) 1962 "USAGE: fl2poly(L,s); L a list, s a string 1896 /* 1897 proc operatorModulo(poly F, ideal I, poly b) 1898 "USAGE: operatorModulo(f,I,b); f a poly, I an ideal, b a poly 1963 1899 RETURN: poly 1964 PURPOSE: reconstruct a monic polynomial in one variable from its factorization 1965 ASSUME: s is a string with the name of some variable and L is supposed to consist of two entries: 1966 @* L[1] of the type ideal with the roots of a polynomial 1967 @* L[2] of the type intvec with the multiplicities of corr. roots 1968 EXAMPLE: example fl2poly; shows examples 1900 PURPOSE: compute the Boperator from the poly F, ideal I = Ann f^s and BernsteinSato 1901 polynomial b using modulo 1902 NOTE: The computations take place in the ring, similar to the one returned by Sannfs procedure. 1903 @* If printlevel=1, progress debug messages will be printed, 1904 @* if printlevel>=2, all the debug messages will be printed. 1905 EXAMPLE: example operatorModulo; shows examples 1969 1906 " 1970 1907 { 1971 if (varnum(s)==0) 1972 { 1973 ERROR("no such variable found"); return(0); 1974 } 1975 poly x = var(varnum(s)); 1976 poly P = 1; 1977 int sl = size(L[1]); 1978 ideal RR = L[1]; 1979 intvec IV = L[2]; 1980 for(int i=1; i<= sl; i++) 1981 { 1982 P = P*((xRR[i])^IV[i]); 1983 } 1984 return(P); 1908 // with hom_kernel; 1909 matrix AA[1][2] = F,b; 1910 // matrix M = matrix(subst(I,s,s+1)*freemodule(2)); // ann f^{s+1} 1911 matrix N = matrix(I); // ann f^s 1912 // matrix K = hom_kernel(AA,M,N); 1913 matrix K = modulo(AA,N); 1914 K = transpose(K); 1915 print(K); 1916 ideal J; 1917 int i; 1918 poly t; number n; 1919 for(i=1; i<=nrows(K); i++) 1920 { 1921 if (K[i,2]!=0) 1922 { 1923 if ( leadmonom(K[i,2]) == 1) 1924 { 1925 t = K[i,1]; 1926 n = leadcoef(K[i,2]); 1927 t = t/n; 1928 // J = J, K[i][2]; 1929 break; 1930 } 1931 } 1932 } 1933 ideal J = groebner(subst(I,s,s+1)); // for NF 1934 t = NF(t,J); 1935 "candidate:"; t; 1936 J = subst(J,s,s1); 1937 // test: 1938 if ( NF(t*Fb,J) !=0) 1939 { 1940 "Problem: PS does not work on F"; 1941 } 1942 return(t); 1985 1943 } 1986 1944 example 1987 1945 { 1988 1946 "EXAMPLE:"; echo = 2; 1989 ring r = 0,(x,y,z,s),Dp; 1990 ideal I = 1,4/3,5/3,2; 1991 intvec mI = 2,1,1,1; 1992 list BS = I,mI; 1993 poly p = fl2poly(BS,"s"); 1994 p; 1995 factorize(p,2); 1996 } 1947 // ring r = 0,(x,y,z),Dp; 1948 // poly F = x^3+y^3+z^3; 1949 LIB "dmod.lib"; option(prot); option(mem); 1950 ring r = 0,(x,y),Dp; 1951 // poly F = x^3+y^3+x*y^2; 1952 poly F = x^4 + y^4 + x*y^4; 1953 def A = Sannfs(F); // here we get LD = ann f^s 1954 setring A; 1955 poly F = imap(r,F); 1956 def B = annfs0(LD,F); // to obtain BS polynomial 1957 list BS = imap(B,BS); 1958 poly b = fl2poly(BS,"s"); 1959 LD = groebner(LD); 1960 poly PS = operatorModulo(F,LD,b); 1961 PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD 1962 reduce(PS*Fbs,LD); // check the property of PS 1963 } 1964 */ 1965 1966 1997 1967 1998 1968 proc annfsParamBM (poly F, list #) … … 2005 1975 @* If eng <>0, @code{std} is used for Groebner basis computations, 2006 1976 @* otherwise, and by default @code{slimgb} is used. 2007 @* If printlevel=1, progress debug messages will be printed,2008 @* if printlevel>=2, all the debug messages will be printed.1977 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, 1978 @* if @code{printlevel}>=2, all the debug messages will be printed. 2009 1979 EXAMPLE: example annfsParamBM; shows examples 2010 1980 " … … 3445 3415 @* If eng <>0, @code{std} is used for Groebner basis computations, 3446 3416 @* otherwise, and by default @code{slimgb} is used. 3447 @*If printlevel=1, progress debug messages will be printed,3417 DISPLAY: If printlevel=1, progress debug messages will be printed, 3448 3418 @* if printlevel>=2, all the debug messages will be printed. 3449 3419 EXAMPLE: example SannfsBFCT; shows examples … … 3637 3607 poly F = x^3+y^3+z^3*w; 3638 3608 printlevel = 0; 3639 def A = SannfsBFCT(F); 3640 setring A; 3609 def A = SannfsBFCT(F); setring A; 3641 3610 intvec v = 1,2,3,4,5,6,7,8; 3642 3611 // are there polynomials, depending on s only? 3643 3612 nselect(LD,v); 3644 3613 // a fancier example: 3645 def R = reiffen(4,5); 3646 setring R; 3614 def R = reiffen(4,5); setring R; 3647 3615 v = 1,2,3,4; 3616 RC; // the Reiffen curve in 4,5 3648 3617 def B = SannfsBFCT(RC); 3649 3618 setring B; … … 3885 3854 } 3886 3855 3887 3856 /* 3888 3857 proc SannfsLOTold(poly F, list #) 3889 3858 "USAGE: SannfsLOT(f [,eng]); f a poly, eng an optional int … … 4109 4078 } 4110 4079 4080 */ 4111 4081 4112 4082 proc annfsLOT(poly F, list #) … … 4158 4128 RETURN: ring 4159 4129 PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based on the 4160 output of procedures SannfsBM, SannfsOT or SannfsLOT4130 output of Sannfslike procedure 4161 4131 NOTE: activate this ring with the @code{setring} command. In this ring, 4162 4132 @*  the ideal LD (which is a Groebner basis) is the annihilator of f^s, … … 4299 4269 poly F = x^3+y^3+z^3; 4300 4270 printlevel = 0; 4301 def A = SannfsBM(F); 4271 def A = SannfsBM(F); setring A; 4302 4272 // alternatively, one can use SannfsOT or SannfsLOT 4303 setring A;4304 4273 LD; 4305 4274 poly F = imap(r,F); 4306 def B = annfs0(LD,F); 4307 setring B; 4275 def B = annfs0(LD,F); setring B; 4308 4276 LD; 4309 4277 BS; … … 4656 4624 @* We compute the real annihilator for any rational value of n (both generic and exceptional). 4657 4625 @* The implementation goes along the lines of SaitoSturmfelsTakayama, Alg. 5.3.15 4658 @*If printlevel=1, progress debug messages will be printed,4626 DISPLAY: If printlevel=1, progress debug messages will be printed, 4659 4627 @* if printlevel>=2, all the debug messages will be printed. 4660 4628 EXAMPLE: example annfspecial; shows examples … … 4724 4692 } 4725 4693 4726 static proc new_ex_annfspecial() 4694 /* 4695 //static proc new_ex_annfspecial() 4727 4696 { 4728 4697 //another example for annfspecial: x3+y3+z3 … … 4740 4709 annfspecial(LD,F,1,1); // exceptional root 4741 4710 } 4711 */ 4742 4712 4743 4713 proc minIntRoot(ideal P, int fact) … … 4817 4787 "USAGE: isHolonomic(M); M an ideal/module/matrix 4818 4788 RETURN: int, 1 if M is holonomic and 0 otherwise 4789 ASSUME: basering is a Weyl algebra 4819 4790 PURPOSE: check the modules for the property of holonomy 4820 4791 NOTE: M is holonomic, if 2*dim(M) = dim(R), where R is a … … 4956 4927 @* 'alg1' (default)  for the algorithm 1 of J. MartinMorales (unpublished) 4957 4928 @* 'alg2'  for the algorithm 2 of J. MartinMorales (unpublished) 4958 @* The output int is:4959 @*  if the algorithm 1 is chosen: 1 if alpha is a root of the global Bernstein polynomial and 0 otherwise4960 @*  if the algorithm 2 is chosen:the multiplicity of alpha as a root of the global Bernstein polynomial of f.4961 @* 4929 @* The output of type int is: 4930 @* 'alg1': 1 if alpha is a root of the global Bernstein polynomial and 0 otherwise 4931 @* 'alg2': the multiplicity of alpha as a root of the global Bernstein polynomial of f. 4932 @* (If alpha is not a root, the output is 0) 4962 4933 @* If eng <>0, @code{std} is used for Groebner basis computations, 4963 4934 @* otherwise (and by default) @code{slimgb} is used. 4964 @*If printlevel=1, progress debug messages will be printed,4965 @* if printlevel>=2, all the debug messages will be printed.4935 DISPLAY: If printlevel=1, progress debug messages will be printed, 4936 @* if printlevel>=2, all the debug messages will be printed. 4966 4937 EXAMPLE: example checkRoot; shows examples 4967 4938 " … … 5068 5039 NOTE: If eng <>0, @code{std} is used for Groebner basis computations, 5069 5040 @* otherwise (and by default) @code{slimgb} is used. 5070 @*If printlevel=1, progress debug messages will be printed,5071 @* if printlevel>=2, all the debug messages will be printed.5041 DISPLAY: If printlevel=1, progress debug messages will be printed, 5042 @* if printlevel>=2, all the debug messages will be printed. 5072 5043 EXAMPLE: example checkRoot1; shows examples 5073 5044 " … … 5192 5163 NOTE: If eng <>0, @code{std} is used for Groebner basis computations, 5193 5164 @* otherwise (and by default) @code{slimgb} is used. 5194 @*If printlevel=1, progress debug messages will be printed,5195 @* if printlevel>=2, all the debug messages will be printed.5165 DISPLAY: If printlevel=1, progress debug messages will be printed, 5166 @* if printlevel>=2, all the debug messages will be printed. 5196 5167 EXAMPLE: example checkRoot2; shows examples 5197 5168 " … … 5324 5295 NOTE: If eng <>0, @code{std} is used for Groebner basis computations, 5325 5296 @* otherwise (and by default) @code{slimgb} is used. 5326 @*If printlevel=1, progress debug messages will be printed,5327 @* if printlevel>=2, all the debug messages will be printed.5297 DISPLAY: If printlevel=1, progress debug messages will be printed, 5298 @* if printlevel>=2, all the debug messages will be printed. 5328 5299 EXAMPLE: example checkFactor; shows examples 5329 5300 " … … 5461 5432 } 5462 5433 5463 5464 static proc exCheckGenericity() 5434 /// ****** EXAMPLES ************ 5435 5436 /* 5437 5438 //static proc exCheckGenericity() 5465 5439 { 5466 5440 LIB "control.lib"; … … 5484 5458 } 5485 5459 5486 static proc exOT_17()5460 //static proc exOT_17() 5487 5461 { 5488 5462 // OakuTakayama, p.208 … … 5499 5473 } 5500 5474 5501 static proc exOT_16()5475 //static proc exOT_16() 5502 5476 { 5503 5477 // OakuTakayama, p.208 … … 5513 5487 } 5514 5488 5515 static proc ex_bcheck()5489 //static proc ex_bcheck() 5516 5490 { 5517 5491 ring R = 0,(x,y),dp; … … 5526 5500 } 5527 5501 5528 static proc ex_bcheck2()5502 //static proc ex_bcheck2() 5529 5503 { 5530 5504 ring R = 0,(x,y),dp; … … 5536 5510 } 5537 5511 5538 static proc ex_BMI()5512 //static proc ex_BMI() 5539 5513 { 5540 5514 // a hard example … … 5549 5523 } 5550 5524 5551 static proc ex2_BMI()5525 //static proc ex2_BMI() 5552 5526 { 5553 5527 // this example was believed to be intractable in 2005 by GagoVargas, Castro and Ucha … … 5563 5537 } 5564 5538 5565 static proc ex_operatorBM()5539 //static proc ex_operatorBM() 5566 5540 { 5567 5541 ring r = 0,(x,y,z,w),Dp; … … 5578 5552 reduce(PS*Fbs,LD); // check the property of PS 5579 5553 } 5554 5555 */ 
Singular/LIB/dmodapp.lib
r6045c3 r565e86 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: dmodapp.lib,v 1.1 2 20081209 17:52:35levandov Exp $";2 version="$Id: dmodapp.lib,v 1.13 20081223 21:39:31 levandov Exp $"; 3 3 category="Noncommutative"; 4 4 info=" 5 5 LIBRARY: dmodapp.lib Applications of algebraic Dmodules 6 AUTHORS: Viktor Levandovskyy, levandov@math.rwthaachen.de6 AUTHORS: Viktor Levandovskyy, levandov@math.rwthaachen.de 7 7 @* Daniel Andres, daniel.andres@math.rwthaachen.de 8 8 … … 18 18 rational function G/F from Quot(R). These can be computed via annPoly resp. annRat. 19 19 20 @*  initial form and initial ideals in Weyl algebras with respect to a given weight vector can be computed with the help of inForm, initial malgrange, initialideal.20 @*  initial form and initial ideals in Weyl algebras with respect to a given weight vector can be computed with the help of inForm, initialMalgrange, initialIdeal. 21 21 22 22 @*  appelF1, appelF2 and appelF4 return ideals in parametric Weyl algebra, which … … 32 32 DLoc0(I, F); presentation of the localization of D/I w.r.t. f^s, based on the procedure SDLoc 33 33 34 initial malgrange(f[,s,t,u,v]); Groebner basis of the initial Malgrange ideal for a given poly35 initial ideal(I,u,v[,s,t]); initial ideal of a given ideal w.r.t. given weights34 initialMalgrange(f[,s,t,u,v]); Groebner basis of the initial Malgrange ideal for a given poly 35 initialIdeal(I,u,v[,s,t]); initial ideal of a given ideal w.r.t. given weights 36 36 inForm(f,w); initial form of a poly/ideal w.r.t. a given weight 37 37 isFsat(I, F); check whether the ideal I is Fsaturated … … 44 44 engine(I,i); computes a Groebner basis with the algorithm given by i 45 45 poly2list(f); decompose a poly to a list of terms and exponents 46 fl2poly(L,s); reconstruct a monic univariate polynomial from its factorization 47 insertGenerator(id,p[,k]); insert an element into an ideal/module 48 deleteGenerator(id,k); delete the kth element from an ideal/module 49 46 50 47 51 SEE ALSO: dmod_lib, gmssing_lib, bfct_lib … … 58 62 // charInfo(); ??? 59 63 64 65 /////////////////////////////////////////////////////////////////////////////// 66 // testing for consistency of the library: 60 67 proc testdmodapp() 61 68 { 62 example initial ideal;63 example initial malgrange;69 example initialIdeal; 70 example initialMalgrange; 64 71 example DLoc; 65 72 example DLoc0; … … 73 80 example appelF4; 74 81 example poly2list; 75 } 76 77 78 proc initialideal (ideal I, intvec u, intvec v, list #) 79 "USAGE: initialideal(I,u,v [,s,t]); I an ideal, u,v intvecs, s,t optional ints 80 RETURN: an ideal, a Broebner basis of the initial ideal of the input ideal w.r.t. the weights u and v 81 PURPOSE: compute the initial ideal 82 NOTE: Assume, I is an ideal in the nth Weyl algebra where the sequence of the 83 @* indeterminates is x(1),...,x(n),D(1),...,D(n). 84 @* Further assume that u is the weight for the x(i) and v the weight for the D(i). 82 example fl2poly; 83 example insertGenerator; 84 example deleteGenerator; 85 } 86 87 88 proc initialIdeal (ideal I, intvec u, intvec v, list #) 89 "USAGE: initialIdeal(I,u,v [,s,t]); I ideal, u,v intvecs, s,t optional ints 90 RETURN: ideal, GB of initial ideal of the input ideal wrt the weights u and v 91 PURPOSE: computes the initial ideal 92 NOTE: Assume, I is an ideal in the nth Weyl algebra D, where the sequence 93 @* of the indeterminates is x(1),...,x(n),D(1),...,D(n). 94 @* Further assume that u is the weight for the x(i) and v the weight for 95 @* the D(i). 85 96 @* Note that the returned ideal is not a Dideal but an ideal in the associated 86 97 @* graded ring while the Groebner basis is a subset of D. … … 89 100 @* If t<>0, a matrix ordering is used for Groebner basis computations, 90 101 @* otherwise, and by default, a block ordering is used. 91 @*If printlevel=1, progress debug messages will be printed,102 DISPLAY: If printlevel=1, progress debug messages will be printed, 92 103 @* if printlevel>=2, all the debug messages will be printed. 93 EXAMPLE: example initial ideal; shows examples104 EXAMPLE: example initialIdeal; shows examples 94 105 " 95 106 { … … 113 124 } 114 125 } 115 def D = initial idealengine("initialideal", whichengine, methodord, I, u, v);116 ideal inF = fetch(D,inF); 126 def D = initialIdealEngine("initialIdeal", whichengine, methodord, I, u, v); 127 ideal inF = fetch(D,inF); attrib(inF,"isSB",1); 117 128 return(inF); 118 129 } … … 125 136 intvec u = 1; intvec v = 2; 126 137 ideal I = x^2*Dx^2,x*Dx^4; 127 ideal J = initialideal(I,u,v); J; 128 } 129 130 proc initialmalgrange (poly f,list #) 131 "USAGE: initialmalgrange(f, [,s,t,u,v]); f a poly, s,t,u optional ints, v an optional intvec 132 RETURN: a ring, the Weyl algebra induced by the basering, extended in the indeterminates t and Dt 133 PURPOSE: compute the initial Malgrange ideal of a given poly w.r.t. the weight vector 134 @* (1,0...,0,1,0,...,0) such that the weight of t is 1 and the weight of Dt is 1. 138 ideal J = initialIdeal(I,u,v); J; 139 } 140 141 proc initialMalgrange (poly f,list #) 142 "USAGE: initialMalgrange(f, [,s,t,u,v]); f poly, s,t,u optional ints, v optional intvec 143 RETURN: ring, the Weyl algebra induced by basering, extended with vars t and Dt, 144 @* containing the ideal \"inF\", being the initial ideal of the Malgrange 145 @* ideal of f. 146 PURPOSE: computes the initial Malgrange ideal of a given poly wrt the weight 147 @* vector (1,0...,0,1,0,...,0) such that the weight of t is 1 and the 148 @* weight of Dt is 1. 135 149 NOTE: Activate the output ring with the @code{setring} command. 136 150 @* Varnames of the basering should not include t and Dt. 137 @* In the ouput ring, inF is the initial Malgrange ideal. 138 @* If s<>0, @code{std} is used for Groebner basis computations, otherwise, 139 @* and by default, @code{slimgb} is used. 151 @* If s<>0, @code{std} is used for Groebner basis computations, 152 @* otherwise, and by default, @code{slimgb} is used. 140 153 @* If t<>0, a matrix ordering is used for Groebner basis computations, 141 154 @* otherwise, and by default, a block ordering is used. 142 155 @* If u<>0, the order of the variables is reversed. 143 @* If v is a positive weight vector, v is used for homogenization computations,144 @* otherwise and by default, no weights are used.145 @*If printlevel=1, progress debug messages will be printed,156 @* If v is a positive weight vector, v is used for homogenization 157 @* computations, otherwise and by default, no weights are used. 158 DISPLAY: If printlevel=1, progress debug messages will be printed, 146 159 @* if printlevel>=2, all the debug messages will be printed. 147 EXAMPLE: example initial malgrange; shows examples160 EXAMPLE: example initialMalgrange; shows examples 148 161 " 149 162 { … … 175 188 if (size(#)>3) 176 189 { 177 if (typeof(#[4])=="intvec" && size(#[4])==n && ispositive(#[4])==1)190 if (typeof(#[4])=="intvec" && size(#[4])==n && allPositive(#[4])==1) 178 191 { 179 192 u0 = #[4]; … … 187 200 u0 = 1:n; 188 201 } 189 def D = initial idealengine("initialmalgrange",whichengine, methodord, f, 0, 0, u0, reversevars);190 setring D;202 def D = initialIdealEngine("initialMalgrange",whichengine, methodord, f, 0, 0, u0, reversevars); 203 setring save; 191 204 return(D); 192 205 } … … 196 209 ring r = 0,(x,y),dp; 197 210 poly f = x^2+y^3+x*y^2; 198 def D = initial malgrange(f);211 def D = initialMalgrange(f); 199 212 setring D; 200 213 inF; 201 214 setring r; 202 215 intvec v = 3,2; 203 def D2 = initial malgrange(f,1,0,1,v);216 def D2 = initialMalgrange(f,1,0,1,v); 204 217 setring D2; 205 218 inF; 206 219 } 207 220 208 proc initialidealengine(string calledfrom, int whichengine, int blockord, list #)221 static proc initialIdealEngine(string calledfrom, int whichengine, int blockord, list #) 209 222 { 210 223 //#[1] = f or I … … 218 231 n = nvars(save); 219 232 intvec uv; 220 if (calledfrom == "initial ideal")233 if (calledfrom == "initialIdeal") 221 234 { 222 235 ideal I = #[1]; … … 227 240 noofvars = 2*n+1; 228 241 } 229 else 242 else // initialMalgrange 230 243 { 231 244 poly f = #[1]; 232 if (calledfrom == "initialmalgrange") 233 { 234 uv[n+2] = 1; 235 noofvars = 2*n+3; 236 intvec u0 = #[4]; 237 int reversevars = #[5]; 238 ring r = 0,(x(1..n)),wp(u0); 239 } 240 else // bfctonestep 241 { 242 uv[n+3] = 1; 243 noofvars = 2*n+4; 244 ring r = 0,(x(1..n)),dp; 245 } 245 uv[n+2] = 1; 246 noofvars = 2*n+3; 247 intvec u0 = #[4]; 248 int reversevars = #[5]; 249 ring r = 0,(x(1..n)),wp(u0); 246 250 poly f = fetch(save,f); 247 251 uv[1] = 1; uv[noofvars] = 0; … … 249 253 // for the ordering 250 254 intvec @a; 251 if (calledfrom == "initial malgrange")255 if (calledfrom == "initialMalgrange") 252 256 { 253 257 int d = deg(f); … … 264 268 @a = weighttx,weightD,1; 265 269 } 266 else 270 else // initialIdeal 267 271 { 268 272 @a = 1:noofvars; 269 if (calledfrom == "bfctonestep")270 {271 @a[2] = 2;272 intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0;273 }274 273 } 275 274 if (blockord == 0) // default: blockordering 276 275 { 277 if (calledfrom == "initial malgrange")276 if (calledfrom == "initialMalgrange") 278 277 { 279 278 ring Dh = 0,(t,x(n..1),Dt,D(n..1),h),(a(@a),a(uv),dp(noofvars1),lp(1)); 280 279 } 281 else 282 { 283 if (calledfrom == "initialideal") 284 { 285 ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),dp(noofvars1),lp(1)); 286 } 287 else // bfctonestep 288 { 289 ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(noofvars1),lp(1)); 290 } 280 else // initialIdeal 281 { 282 ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),dp(noofvars1),lp(1)); 291 283 } 292 284 } … … 302 294 dbprint(ppl,"weights for ordering:",transpose(@a)); 303 295 dbprint(ppl,"the ordering matrix:",@Ord); 304 if (calledfrom == "initial malgrange")296 if (calledfrom == "initialMalgrange") 305 297 { 306 298 ring Dh = 0,(t,x(n..1),Dt,D(n..1),h),(a(@a),M(@Ord)); 307 299 } 308 else 309 { 310 if (calledfrom == "initialideal") 311 { 312 ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),M(@Ord)); 313 } 314 else // bfctonestep 315 { 316 ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord)); 317 } 300 else // initialIdeal 301 { 302 ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),M(@Ord)); 318 303 } 319 304 } … … 321 306 // noncommutative relations 322 307 matrix @relD[noofvars][noofvars]; 323 if (calledfrom == "initialmalgrange") 324 { 325 for (i=1; i<=n+1; i++) 326 { 327 @relD[i,n+1+i] = h^(d+1); 328 } 329 } 330 else 331 { 332 if (calledfrom == "initialideal") 333 { 334 for (i=1; i<=n; i++) 335 { 336 @relD[i,n+i] = h^2; 337 } 338 } 339 else // bfctonestep 340 { 341 @relD[1,2] = t*h^2;// s*t = t*s+t*h^2 342 @relD[2,n+3] = Dt*h^2;// Dt*s = s*Dt+h^2 343 @relD[1,n+3] = h^2; 344 for (i=1; i<=n; i++) 345 { 346 @relD[i+2,n+3+i] = h^2; 347 } 308 if (calledfrom == "initialMalgrange") 309 { 310 for (i=1; i<=n+1; i++) { @relD[i,n+1+i] = h^(d+1); } 311 } 312 else // initialIdeal 313 { 314 for (i=1; i<=n; i++) 315 { 316 @relD[i,n+i] = h^2; 348 317 } 349 318 } … … 353 322 dbprint(ppl,"computing in ring",DDh); 354 323 ideal I; 355 if (calledfrom == "initial ideal")324 if (calledfrom == "initialIdeal") 356 325 { 357 326 I = fetch(save,I); … … 367 336 { 368 337 I = I, D(i)+diff(f,x(i))*Dt; 369 }370 if (calledfrom == "bfctonestep")371 {372 I = I,t*Dts;373 338 } 374 339 } … … 380 345 dbprint(ppl, "I after dehomogenization is" ,I); 381 346 I = inForm(I,uv); // we are only interested in the initial form w.r.t. uv 382 if (calledfrom == "initial malgrange")347 if (calledfrom == "initialMalgrange") 383 348 { 384 349 // keep the names of the variables of the basering … … 425 390 } 426 391 } 427 else 428 { 429 if (calledfrom == "initialideal") 430 { 431 ring @D = 0,(x(1..n),D(1..n)),dp; 432 def D = Weyl(@D); 433 setring D; 434 ideal I = imap(DDh,I); 435 kill Dh,DDh; 436 } 437 } 438 if (calledfrom <> "bfctonestep") 439 { 440 dbprint(ppl, "starting cosmetic Groebner basis computation with engine:", whichengine); 441 I = engine(I, whichengine); 442 dbprint(ppl,"finished cosmetic Groebner basis computation:"); 443 dbprint(ppl,"the initial ideal is:", I); 444 } 445 ideal inF = I; 392 else // initialIdeal 393 { 394 ring @D = 0,(x(1..n),D(1..n)),dp; 395 def D = Weyl(@D); 396 setring D; 397 ideal I = imap(DDh,I); 398 kill Dh,DDh; 399 } 400 dbprint(ppl, "starting cosmetic Groebner basis computation with engine:", whichengine); 401 I = engine(I, whichengine); 402 dbprint(ppl,"finished cosmetic Groebner basis computation:"); 403 dbprint(ppl,"the initial ideal is:", I); 404 ideal inF = I; attrib(inF,"isSB",1); 446 405 export(inF); 447 406 return(basering); 448 407 } 449 408 450 proc inForm (list #) 451 "USAGE: inForm(I,w); I an ideal, w an intvec 452 RETURN: the initial form of I w.r.t. the weight vector w 453 PURPOSE: compute the initial form of an ideal w.r.t a given weight vector 454 NOTE: the size of the weight vector must be equal to the number of variables of the basering 409 proc inForm (ideal I, intvec w) 410 "USAGE: inForm(I,w); I ideal, w intvec 411 RETURN: the initial form of I wrt the weight vector w 412 PURPOSE: computes the initial form of an ideal wrt a given weight vector 413 NOTE: the size of the weight vector must be equal to the number of variables 414 @* of the basering. 455 415 EXAMPLE: example inForm; shows examples 456 416 " 457 417 { 458 if (size(#)<2)459 {460 ERROR("inForm needs two arguments of type poly/ideal,intvec");461 }462 if (typeof(#[1])=="poly"  typeof(#[1])=="ideal")463 {464 ideal I = #[1];465 }466 else467 {468 ERROR("first argument must be of type poly or ideal");469 }470 if (typeof(#[2])=="intvec")471 {472 def w = #[2];473 }474 else475 {476 ERROR("second argument must be of type intvec");477 }478 418 if (size(w) != nvars(basering)) 479 419 { 480 420 ERROR("weight vector has wrong dimension"); 481 421 } 422 if (I == 0) 423 { 424 return(I); 425 } 482 426 int j,i,s,m; 483 427 list l; … … 487 431 { 488 432 l = poly2list(I[j]); 489 m = scalar prod(w,l[1][1]);490 g = 0;433 m = scalarProd(w,l[1][1]); 434 g = l[1][2]; 491 435 for (i=2; i<=size(l); i++) 492 436 { 493 s = scalarprod(w,l[i][1]); 494 m = Max(list(s,m)); 495 } 496 for (i=1; i<=size(l); i++) 497 { 498 if (scalarprod(w,l[i][1]) == m) 437 s = scalarProd(w,l[i][1]); 438 if (s == m) 499 439 { 500 g = g +l[i][2];440 g = g + l[i][2]; 501 441 } 442 else 443 { 444 if (s > m) 445 { 446 m = s; 447 g = l[i][2]; 448 } 449 } 502 450 } 503 451 J[j] = g; 504 452 } 505 if (typeof(#[1])=="poly") 506 { 507 return(J[1]); // if the input was a poly, return a poly 508 } 509 else 510 { 511 return(J); // otherwise, return an ideal 512 } 453 return(J); 513 454 } 514 455 example … … 527 468 inForm(I,w2); 528 469 inForm(I,w3); 529 inForm(F,w1); 530 } 531 532 static proc charVariety(ideal I) 470 } 471 472 /* 473 474 proc charVariety(ideal I) 533 475 "USAGE: charVariety(I); I an ideal 534 476 RETURN: ring … … 615 557 } 616 558 559 /* 560 617 561 // TODO 618 static proc charInfo(ideal I) 562 563 /* 564 proc charInfo(ideal I) 619 565 "USAGE: charInfo(I); I an ideal 620 566 RETURN: ring … … 638 584 // run primdec 639 585 } 586 */ 640 587 641 588 … … 666 613 "EXAMPLE:"; echo = 2; 667 614 ring r = 0,x,dp; 668 def A = appelF1(); //(1,2,3,4);615 def A = appelF1(); 669 616 setring A; 670 617 IAppel1; … … 696 643 "EXAMPLE:"; echo = 2; 697 644 ring r = 0,x,dp; 698 def A = appelF2(); //(1,2,3,4);645 def A = appelF2(); 699 646 setring A; 700 647 IAppel2; … … 726 673 "EXAMPLE:"; echo = 2; 727 674 ring r = 0,x,dp; 728 def A = appelF4(); //(1,2,3,4);675 def A = appelF4(); 729 676 setring A; 730 677 IAppel4; … … 740 687 list l; 741 688 int i = 1; 742 if (f ==0) // just for the zero polynomial689 if (f == 0) // just for the zero polynomial 743 690 { 744 691 l[1] = list(leadexp(f), lead(f)); 745 692 } 746 while (f!=0) 693 else { l[size(f)] = list(); } // memory preallocation 694 while (f != 0) 747 695 { 748 696 l[i] = list(leadexp(f), lead(f)); 749 f = f lead(f);697 f = f  lead(f); 750 698 i++; 751 699 } … … 803 751 "USAGE: DLoc(I, F); I an ideal, F a poly 804 752 RETURN: nothing (exports objects instead) 805 ASSUME: the basering is a Weyl algebra and I is Fsaturated753 ASSUME: the basering is a Weyl algebra 806 754 PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s 807 755 NOTE: In the basering, the following objects are exported: … … 860 808 { 861 809 /* assume: to be run in the output ring of SDLoc */ 862 /* todo: add F, eliminate vars*Dvars, factorize BS */810 /* doing: add F, eliminate vars*Dvars, factorize BS */ 863 811 /* analogue to annfs0 */ 864 812 def @R2 = basering; … … 1292 1240 NOTE: activate the output ring with the @code{setring} command. 1293 1241 @* In the output ring: 1294 @*  the ideal RLD (which is given in a Groebner basis) is the annihilator.1242 @*  the ideal LD (which is given in a Groebner basis) is the annihilator. 1295 1243 @* If @code{printlevel}=1, progress debug messages will be printed, 1296 1244 @* if @code{printlevel}>=2, all the debug messages will be printed. … … 1369 1317 kill @R1; kill @R2; 1370 1318 dbprint(ppl,"// 4[annRat] running modulo"); 1371 ideal RLD = modulo(G,LL);1319 ideal LD = modulo(G,LL); 1372 1320 dbprint(ppl,"// 4[annRat] running GB on the final result"); 1373 RLD = engine(RLD,0);1374 export RLD;1321 LD = engine(LD,0); 1322 export LD; 1375 1323 return(@R3); 1376 1324 } … … 1382 1330 def B = annRat(g,f); 1383 1331 setring B; 1384 RLD;1332 LD; 1385 1333 // Now, compare with the output of Macaulay2: 1386 1334 ideal tst = 3*x*Dx + 2*y*Dy + 1, y^3*Dy^2  x^2*Dy^2 + 6*y^2*Dy + 6*y, 1387 1335 9*y^2*Dx^2*Dy  4*y*Dy^3 + 27*y*Dx^2 + 2*Dy^2, 9*y^3*Dx^2  4*y^2*Dy^2 + 10*y*Dy 10; 1388 1336 option(redSB); option(redTail); 1389 RLD = groebner(RLD);1337 LD = groebner(LD); 1390 1338 tst = groebner(tst); 1391 print(matrix(NF( RLD,tst))); print(matrix(NF(tst,RLD)));1339 print(matrix(NF(LD,tst))); print(matrix(NF(tst,LD))); 1392 1340 // So, these two answers are the same 1393 1341 } 1394 1342 1395 static proc ex_annRat() 1396 { 1397 // more complicated example 1343 /* 1344 //static proc ex_annRat() 1345 { 1346 // more complicated example for annRat 1398 1347 ring r = 0,(x,y,z),dp; 1399 1348 poly f = x3+y3+z3; // mir = 2 … … 1402 1351 setring A; 1403 1352 } 1353 */ 1404 1354 1405 1355 proc annPoly(poly f) … … 1409 1359 NOTE: activate the output ring with the @code{setring} command. 1410 1360 @* In the output ring: 1411 @*  the ideal RLD (which is given in a Groebner basis) is the annihilator.1361 @*  the ideal LD (which is given in a Groebner basis) is the annihilator. 1412 1362 @* If @code{printlevel}=1, progress debug messages will be printed, 1413 1363 @* if @code{printlevel}>=2, all the debug messages will be printed. … … 1455 1405 } 1456 1406 1457 static proc exCusp() 1407 /* DIFFERENT EXAMPLES 1408 1409 //static proc exCusp() 1458 1410 { 1459 1411 "EXAMPLE:"; echo = 2; … … 1474 1426 } 1475 1427 1476 static proc exWalther1()1428 //static proc exWalther1() 1477 1429 { 1478 1430 // p.18 Rem 3.10 … … 1494 1446 } 1495 1447 1496 static proc exWalther2()1448 //static proc exWalther2() 1497 1449 { 1498 1450 // p.19 Rem 3.10 cont'd … … 1518 1470 } 1519 1471 1520 static proc exWalther3()1472 //static proc exWalther3() 1521 1473 { 1522 1474 // can check with annFs too :) … … 1550 1502 } 1551 1503 1504 */ 1505 1552 1506 proc engine(ideal I, int i) 1553 1507 "USAGE: engine(I,i); I an ideal, i an int … … 1581 1535 engine(I,1); // uses std 1582 1536 } 1537 1538 proc insertGenerator (list #) 1539 "USAGE: insertGenerator(id,p[,k]); id an ideal/module, p a poly/vector, k an optional int 1540 RETURN: same as id 1541 PURPOSE: insert an element into an ideal or a module 1542 NOTE: If k is given, p is inserted at position k, otherwise (and by default), 1543 @* p is inserted at the beginning. 1544 EXAMPLE: example insertGenerator; shows examples 1545 " 1546 { 1547 if (size(#) < 2) 1548 { 1549 ERROR("insertGenerator has to be called with at least 2 arguments (ideal/module,poly/vector)"); 1550 } 1551 string inp1 = typeof(#[1]); 1552 if (inp1 == "ideal"  inp1 == "module") 1553 { 1554 if (inp1 == "ideal") { ideal id = #[1]; } 1555 else { module id = #[1]; } 1556 } 1557 else { ERROR("first argument has to be of type ideal or module"); } 1558 string inp2 = typeof(#[2]); 1559 if (inp2 == "poly"  inp2 == "vector") 1560 { 1561 if (inp2 =="poly") { poly f = #[2]; } 1562 else 1563 { 1564 if (inp1 == "ideal") 1565 { 1566 ERROR("second argument has to be a poly if first argument is an ideal"); 1567 } 1568 else { vector f = #[2]; } 1569 } 1570 } 1571 else { ERROR("second argument has to be of type poly/vector"); } 1572 int n = ncols(id); 1573 int k = 1; // default 1574 if (size(#)>=3) 1575 { 1576 if (typeof(#[3]) == "int") 1577 { 1578 k = #[3]; 1579 if (k<=0) 1580 { 1581 ERROR("third argument has to be positive"); 1582 } 1583 } 1584 else { ERROR("third argument has to be of type int"); } 1585 } 1586 execute(inp1 +" J;"); 1587 if (k == 1) { J = f,id; } 1588 else 1589 { 1590 if (k>n) 1591 { 1592 J = id; 1593 J[k] = f; 1594 } 1595 else // 1<k<=n 1596 { 1597 J[1..k1] = id[1..k1]; 1598 J[k] = f; 1599 J[k+1..n+1] = id[k..n]; 1600 } 1601 } 1602 return(J); 1603 } 1604 example 1605 { 1606 "EXAMPLE:"; echo = 2; 1607 ring r = 0,(x,y,z),dp; 1608 ideal I = x^2,z^4; 1609 insertGenerator(I,y^3); 1610 insertGenerator(I,y^3,2); 1611 module M = I; 1612 insertGenerator(M,[x^3,y^2,z],2); 1613 } 1614 1615 proc deleteGenerator (list #) 1616 "USAGE: deleteGenerator(id,k); id an ideal/module, k an int 1617 RETURN: same as id 1618 PURPOSE: deletes the kth element from an ideal or a module 1619 EXAMPLE: example insertGenerator; shows examples 1620 " 1621 { 1622 if (size(#) < 2) 1623 { 1624 ERROR("deleteGenerator has to be called with 2 arguments (ideal/module,int)"); 1625 } 1626 string inp1 = typeof(#[1]); 1627 if (inp1 == "ideal"  inp1 == "module") 1628 { 1629 if (inp1 == "ideal") { ideal id = #[1]; } 1630 else { module id = #[1]; } 1631 } 1632 else { ERROR("first argument has to be of type ideal or module"); } 1633 string inp2 = typeof(#[2]); 1634 if (inp2 == "int"  inp2 == "number") { int k = int(#[2]); } 1635 else { ERROR("second argument has to be of type int"); } 1636 int n = ncols(id); 1637 if (n == 1) { ERROR(inp1+" must have more than one generator"); } 1638 if (k<=0  k>n) { ERROR("second argument has to be in the range 1,...,"+string(n)); } 1639 execute(inp1 +" J;"); 1640 if (k == 1) { J = id[2..n]; } 1641 else 1642 { 1643 if (k == n) { J = id[1..n1]; } 1644 else 1645 { 1646 J[1..k1] = id[1..k1]; 1647 J[k..n1] = id[k+1..n]; 1648 } 1649 } 1650 return(J); 1651 } 1652 example 1653 { 1654 "EXAMPLE:"; echo = 2; 1655 ring r = 0,(x,y,z),dp; 1656 ideal I = x^2,y^3,z^4; 1657 deleteGenerator(I,2); 1658 module M = [x,y,z],[x2,y2,z2],[x3,y3,z3];; 1659 deleteGenerator(M,2); 1660 } 1661 1662 proc fl2poly(list L, string s) 1663 "USAGE: fl2poly(L,s); L a list, s a string 1664 RETURN: poly 1665 PURPOSE: reconstruct a monic polynomial in one variable from its factorization 1666 ASSUME: s is a string with the name of some variable and L is supposed to consist of two entries: 1667 @* L[1] of the type ideal with the roots of a polynomial 1668 @* L[2] of the type intvec with the multiplicities of corr. roots 1669 EXAMPLE: example fl2poly; shows examples 1670 " 1671 { 1672 if (varnum(s)==0) 1673 { 1674 ERROR("no such variable found"); return(0); 1675 } 1676 poly x = var(varnum(s)); 1677 poly P = 1; 1678 int sl = size(L[1]); 1679 ideal RR = L[1]; 1680 intvec IV = L[2]; 1681 for(int i=1; i<= sl; i++) 1682 { 1683 P = P*((xRR[i])^IV[i]); 1684 } 1685 return(P); 1686 } 1687 example 1688 { 1689 "EXAMPLE:"; echo = 2; 1690 ring r = 0,(x,y,z,s),Dp; 1691 ideal I = 1,4/3,5/3,2; 1692 intvec mI = 2,1,1,1; 1693 list BS = I,mI; 1694 poly p = fl2poly(BS,"s"); 1695 p; 1696 factorize(p,2); 1697 }
Note: See TracChangeset
for help on using the changeset viewer.