- Timestamp:
- Mar 10, 2009, 5:12:52 PM (15 years ago)
- Branches:
- (u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '4bd32dfef92ec9f5ed8dceee82d14318ae147107')
- Children:
- af9f1bc5804a0d34be405baa2e1604f18c96b4d4
- Parents:
- 094f80e33c853c56c092642444a1dbea668c09f8
- Location:
- Singular/LIB
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/bfun.lib
r094f80 reda414 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: bfun.lib,v 1. 4 2009-03-09 18:34:51levandov Exp $";2 version="$Id: bfun.lib,v 1.5 2009-03-10 16:12:52 levandov Exp $"; 3 3 category="Noncommutative"; 4 4 info=" … … 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 b-Function (also known as Bernstein-Sato 11 @* polynomial) b(s) in K[s], defined to be the monic polynomial, satisfying 12 @* a functional identity L * F^{s+1} = b(s) F^s, 13 @* for some operator L in D[s]. Here, D stands for an n-th Weyl algebra 10 @* one is interested in the global b-function (also known as Bernstein-Sato 11 @* polynomial) b(s) in K[s], defined to be the monic polynomial of minimal 12 @* degree, satisfying a functional identity L * F^{s+1} = b(s) F^s, 13 @* for some operator L in D[s] (* stands for the action of differential operator) 14 @* By D one denotes the n-th Weyl algebra 14 15 @* K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>. 15 16 @* One is interested in the following data: … … 18 19 @* - the multiplicities of the roots. 19 20 @* 20 @* There is a general definition of a b-function of a holonomic ideal [SST] 21 @* (that is, an ideal I in a Weyl algebra D, such that D/I is holonomic module) 22 @* with respect to the given weight vector. B-function b_w(s) is in K[s]. 23 @* Unlike Bernstein-Sato polynomial, general B-function with respect to 24 @* arbitrary weights need not have rational roots at all. 21 @* There is a general definition of a b-function of a holonomic ideal [SST] 22 @* (that is, an ideal I in a Weyl algebra D, such that D/I is holonomic module) 23 @* with respect to the given weight vector w: For a poly p in D, its initial 24 @* form w.r.t. (-w,w) is defined as the sum of all terms of p which have 25 @* maximal weighted total degree where the weight of x_i is -w_i and the weight 26 @* of d_i is w_i. Let J be the initial ideal of I w.r.t. (-w,w), i.e. the 27 @* K-vector space generated by all initial forms w.r.t (-w,w) of elements of I. 28 @* Put s = w_1 x_1 d_1 + ... + w_n x_n d_n. Then the monic generator b_w(s) of 29 @* the intersection of J with the PID K[s] is called the b-function of I w.r.t. w. 30 @* Unlike Bernstein-Sato polynomial, general b-function with respect to 31 @* arbitrary weights need not have rational roots at all. However, b-function 32 @* of a holonomic ideal is known to be non-zero as well. 25 33 @* 26 34 @* References: [SST] Saito, Sturmfels, Takayama: Groebner Deformations of … … 32 40 MAIN PROCEDURES: 33 41 34 bfct(f[,s,t,v]); compute s the Bernstein-Sato polynomial of poly f35 bfctSyz(f[,r,s,t,u,v]); computes the Bernstein-Sato polynomial of poly f36 bfctAnn(f[,s]); computes the Bernstein-Sato polynomial of polyf37 bfctOneGB(f[,s,t]); computes the Bernstein-Sato polynomial of poly f38 bfctIdeal(I,w[,s,t]); computes the globalb-function of ideal wrt weight39 pIntersect(f,I[,s]); intersection of ideal with subalgebra K[f] forpoly f40 pIntersectSyz(f,I[,p,s,t]); intersection of ideal with subalgebra K[f] for poly f41 linReduce(f,I[,s]); reduce s a poly by linear reductions only42 linReduceIdeal(I,[s,t]); reduces generators of ideal by linear reductions only43 linSyzSolve(I[,s]); compute sa linear dependency of elements of ideal I42 bfct(f[,s,t,v]); compute the BS polynomial of f with linear reductions 43 bfctSyz(f[,r,s,t,u,v]); compute the BS polynomial of f with syzygy-solver 44 bfctAnn(f[,s]); compute the BS polynomial of f via Ann f^s + f 45 bfctOneGB(f[,s,t]); compute the BS polynomial of f by just one GB computation 46 bfctIdeal(I,w[,s,t]); compute the b-function of ideal wrt weight 47 pIntersect(f,I[,s]); intersection of ideal with subalgebra K[f] for a poly f 48 pIntersectSyz(f,I[,p,s,t]); intersection of ideal with subalgebra K[f] with syz-solver 49 linReduce(f,I[,s]); reduce a poly by linear reductions wrt ideal 50 linReduceIdeal(I,[s,t]); interreduce generators of ideal by linear reductions 51 linSyzSolve(I[,s]); compute a linear dependency of elements of ideal I 44 52 45 53 AUXILIARY PROCEDURES: … … 911 919 @* If p>0 is given, this proc computes the generator of the intersection in 912 920 @* char p first and then only searches for a generator of the obtained 913 @* degree in the basering. Otherwise, it searched for all degrees by 914 @* computing syzygies. Note that the basering must be free of parameters to 915 @* use this modular method. 921 @* degree in the basering. Otherwise, it searches for all degrees by 922 @* computing syzygies. 916 923 @* If s<>0, @code{std} is used for Groebner basis computations in char 0, 917 924 @* otherwise, and by default, @code{slimgb} is used. … … 935 942 int solveincharp = 0; // default 936 943 def save = basering; 944 int n = nvars(save); 937 945 if (size(#)>0) 938 946 { … … 962 970 newNF = NF(s,I); 963 971 NI[2] = newNF; 964 list l = ringlist(save); 965 if (typeof(l[1]) == "int") { l[1] = solveincharp; } 966 else // parameters in basering: problems with mappings from Q(a_i)->Z/p(a_i) 967 { 968 if (solveincharp <> 0) 969 { 970 print("// basering must be free of parameters to use modular method"); 971 print("// solving in characteristic of basering only..."); 972 solveincharp = 0; 973 } 974 } 975 if (solveincharp<>0) 976 { 977 if (size(l)>4) 978 { 979 matrix l5 = l[5]; 980 matrix l6 = l[6]; 981 } 982 def @Rp = ring(l); 972 list RL = ringlist(save); 973 if (solveincharp) 974 { 975 int psolveincharp = prime(solveincharp); 976 if (solveincharp <> psolveincharp) 977 { 978 print("// " + string(solveincharp) + " is invalid characteristic of ground field."); 979 print("// " + string(psolveincharp) + " is used."); 980 solveincharp = psolveincharp; 981 kill psolveincharp; 982 } 983 list RLp = RL[1..4]; 984 if (typeof(RL[1]) == "int") { RLp[1] = solveincharp; } 985 else { RLp[1][1] = solveincharp; } // parameters 986 def @Rp = ring(RLp); 983 987 setring @Rp; 984 list l = ringlist(@Rp);988 number c; 985 989 setring save; 986 if (size(l)>4) 987 { 990 int shortSave = short; // workaround for maps Q(a_i) -> Z/p(a_i) 991 short = 0; 992 string str; 993 int badprime; 994 i=1; 995 while (badprime == 0 && i<=size(s)) // detect bad primes 996 { 997 str = string(denominator(leadcoef(s[i]))); 998 str = "c = " + str + ";"; 988 999 setring @Rp; 989 l[5] = fetch(save,l5); 990 l[6] = fetch(save,l6); 991 } 992 setring @Rp; 993 def Rp = ring(l); 1000 execute(str); 1001 if (c == 0) { badprime = 1; } 1002 setring save; 1003 i++; 1004 } 1005 str = "poly s = " + string(s) + ";"; 1006 if (size(RL) > 4) // basering is NC-algebra 1007 { 1008 string RL5 = "@C = " + string(RL[5]) + ";"; 1009 string RL6 = "@D = " + string(RL[6]) + ";"; 1010 setring @Rp; 1011 matrix @C[n][n]; matrix @D[n][n]; 1012 execute(RL5); execute(RL6); 1013 def Rp = nc_algebra(@C,@D); 1014 } 1015 else { def Rp = @Rp; } 994 1016 setring Rp; 995 1017 kill @Rp; 996 1018 dbprint(ppl-1,"// solving in ring ", Rp); 1019 execute(str); 997 1020 vector v; 998 map phi = save,maxideal(1); 999 poly s = phi(s); 1021 number c; 1000 1022 ideal NI = 1; 1001 1023 setring save; 1024 i=1; 1025 while (badprime == 0 && i<=size(I)) // detect bad primes 1026 { 1027 str = string(leadcoef(cleardenom(I[i]))); 1028 str = "c = " + str + ";"; 1029 setring Rp; 1030 execute(str); 1031 if (c == 0) { badprime = 1; } 1032 setring save; 1033 i++; 1034 } 1035 if (badprime == 1) 1036 { 1037 print("// WARNING: bad prime"); 1038 short = shortSave; 1039 return(vector(0)); 1040 } 1002 1041 } 1003 1042 i = 1; … … 1023 1062 // look for a solution 1024 1063 dbprint(ppl-1,"// linSyzSolve starts with: "+string(matrix(NI))); 1025 if (solveincharp<>0) // modular method 1026 { 1064 if (solveincharp) // modular method 1065 { 1066 for (j=1; j<=size(newNF); j++) 1067 { 1068 str = string(denominator(leadcoef(s[i]))); 1069 str = "c = " + str + ";"; 1070 setring Rp; 1071 execute(str); 1072 if (c == 0) 1073 { 1074 print("// WARNING: bad prime"); 1075 setring save; 1076 short = shortSave; 1077 return(vector(0)); 1078 } 1079 setring save; 1080 } 1081 str = "NI[" + string(i) + "+1] = " + string(newNF) + ";"; 1027 1082 setring Rp; 1028 NI[i+1] = phi(newNF);1083 execute(str); // NI[i+1] = [newNF]_{solveincharp} 1029 1084 v = linSyzSolve(NI,modengine); 1030 1085 if (v!=0) // there is a modular solution … … 1076 1131 } 1077 1132 dbprint(ppl,"// pIntersectSyz finished"); 1133 if (solveincharp) { short = shortSave; } 1078 1134 return(v); 1079 1135 } … … 1240 1296 { 1241 1297 list L = ringlist(D); 1242 if (typeof(L[1])=="int") 1243 { 1244 int lb = 10000; 1245 int ub = 536870909; 1246 i = 1; 1247 list usedprimes; 1248 matrix L5 = L[5]; matrix L6 = L[6]; 1249 int sJ = size(J); int sJq; 1250 while (b == 0) 1298 int lb = 10000; int ub = 536870909; // bounds for random primes 1299 list usedprimes; 1300 int sJ = size(J); 1301 int sLJq; 1302 ideal LJ; 1303 for (i=1; i<=sJ; i++) 1304 { 1305 LJ[i] = leadcoef(cleardenom(J[i])); 1306 } 1307 int short_save = short; // workaround for map Q(a_i) -> Z/q(a_i) 1308 short = 0; 1309 string strLJq = "ideal LJq = " + string(LJ) + ";"; 1310 int nD = nvars(D); 1311 string L5 = "matrix @C[nD][nD]; @C = " + string(L[5]) + ";"; 1312 string L6 = "matrix @D[nD][nD]; @D = " + string(L[6]) + ";"; 1313 L = L[1..4]; 1314 i = 1; 1315 while (b == 0) 1316 { 1317 dbprint(ppl,"// number of run in the loop: "+string(i)); 1318 int q = prime(random(lb,ub)); 1319 if (findFirst(usedprimes,q)==0) // if q has not been used already 1251 1320 { 1252 dbprint(ppl,"// number of run in the loop: "+string(i)); 1253 int q = prime(random(lb,ub)); 1254 if (findFirst(usedprimes,q)==0) // if q was not already used 1321 usedprimes = usedprimes,q; 1322 dbprint(ppl,"// using prime: "+string(q)); 1323 if (typeof(L[1]) == "int") { L[1] = q; } 1324 else { L[1][1] = q; } // parameters 1325 def @Rq = ring(L); setring @Rq; 1326 execute(L5); execute(L6); 1327 def Rq = nc_algebra(@C,@D); // def Rq = nc_algebra(1,@D); 1328 setring Rq; kill @Rq; 1329 execute(strLJq); 1330 sLJq = size(LJq); 1331 setring D; kill Rq; 1332 if (sLJq <> sJ ) // detect unlucky prime 1255 1333 { 1256 usedprimes = usedprimes,q; 1257 dbprint(ppl,"// using prime: "+string(q)); 1258 L[1] = q; 1259 def @Rq = ring(L); setring @Rq; 1260 list lq = ringlist(@Rq); 1261 lq[5] = fetch(D,L5); lq[6] = fetch(D,L6); 1262 def Rq = ring(lq); 1263 setring Rq; kill @Rq; 1264 ideal Jq = fetch(D,J); 1265 sJq = size(lead(Jq)); 1266 setring D; kill Rq; 1267 if (sJq <> sJ) 1268 { 1269 dbprint(ppl,"// " +string(q) + " is unlucky"); 1270 b = 0; 1271 } 1272 else 1273 { 1274 b = pIntersectSyz(s,J,q,whichengine,modengine); 1275 } 1334 dbprint(ppl,"// " +string(q) + " is unlucky"); 1335 b = 0; 1276 1336 } 1277 i++; 1337 else 1338 { 1339 b = pIntersectSyz(s,J,q,whichengine,modengine); 1340 } 1278 1341 } 1279 } 1280 else //parameters: problem when mapping Q(a_i) -> Z/p(a_i) 1281 { 1282 print("// basering must be free of parameters to use modular method"); 1283 print("// solving only in char 0..."); 1284 b = pIntersectSyz(s,J,0,whichengine); 1285 } 1286 } 1287 else // pIntersectSyz::non-modular else // parameters in basering: problems with mappings 1342 i++; 1343 } 1344 short = short_save; 1345 } 1346 else // pIntersectSyz::non-modular 1288 1347 { 1289 1348 b = pIntersectSyz(s,J,0,whichengine); … … 1384 1443 @* If t<>0, the computation of the intersection is solely performed over 1385 1444 @* charasteristic 0, otherwise and by default, a modular method is used. 1386 @* Note that in this case, the basering must be free of parameters.1387 1445 @* If u<>0 and by default, @code{std} is used for GB computations in 1388 1446 @* characteristic >0, otherwise, @code{slimgb} is used. … … 1460 1518 RETURN: list of ideal and intvec 1461 1519 PURPOSE: computes the roots of the global b-function of I wrt the weight (-w,w). 1462 ASSUME: The basering is the n-th Weyl algebra in characteristic 0 and 1520 ASSUME: The basering is the n-th Weyl algebra in characteristic 0 and for all 1463 1521 @* 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the 1464 1522 @* sequence of variables is given by x(1),...,x(n),D(1),...,D(n), 1465 1523 @* where D(i) is the differential operator belonging to x(i). 1466 @* Further assume that I is holonomic.1524 @* Further we assume that I is holonomic. 1467 1525 BACKGROUND: In this proc, the initial ideal of I is computed according to the 1468 1526 @* algorithm by Masayuki Noro and then a system of linear equations is 1469 1527 @* solved by linear reductions. 1470 NOTE: In the output list, the ideal contains all the roots and the intvec 1471 @* their multiplicities. 1528 NOTE: In the output list, say L, 1529 @* - L[1] of type ideal contains all the rational roots of a b-function, 1530 @* - L[2] of type intvec contains the multiplicities of above roots, 1531 @* - optional L[3] of type string is the part of b-function without rational roots. 1532 @* Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and 1533 @* L[3] is 1 (for nonzero constant) or 0 (for zero b-function). 1472 1534 @* If s<>0, @code{std} is used for GB computations in characteristic 0, 1473 1535 @* otherwise, and by default, @code{slimgb} is used. … … 1510 1572 if (isH<>1) 1511 1573 { 1512 print(" //given ideal is not holonomic");1513 print(" //setting bound for degree of b-function to 10 and proceeding");1574 print("WARNING: given ideal is not holonomic"); 1575 print("... setting bound for degree of b-function to 10 and proceeding"); 1514 1576 isH = 10; 1515 1577 } … … 1540 1602 setring D; 1541 1603 ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6; I = std(I); 1542 intvec w1 = 1,1; 1543 intvec w2 = 1,2; 1544 intvec w3 = 2,3; 1604 intvec w1 = 0,1; 1605 intvec w2 = 2,3; 1545 1606 bfctIdeal(I,w1); 1546 bfctIdeal(I,w2,1); 1547 bfctIdeal(I,w3,0,1); 1607 bfctIdeal(I,w2,0,1); 1608 ideal J = I[size(I)]; // J is not holonomic by construction 1609 bfctIdeal(J,w1); // b-function of D/J w.r.t. w1 is non-zero 1610 bfctIdeal(J,w2); // b-function of D/J w.r.t. w2 is zero 1548 1611 } 1549 1612 … … 1768 1831 poly xyzreiffen45 = (xy+z)*(y4+yz4+z5); 1769 1832 bfct(xyzreiffen45); 1833 1834 // sparse ideal as suggested by Alex; gives 1 as std 1835 ideal I1 = 28191*y^2+14628*x*Dy, 24865*x^2+24072*x*Dx+17756*Dy^2; 1836 std(I1); 1837 1838 1770 1839 } 1771 1840 */ -
Singular/LIB/dmodapp.lib
r094f80 reda414 1 1 ////////////////////////////////////////////////////////////////////////////// 2 version="$Id: dmodapp.lib,v 1.2 2 2009-03-09 18:34:52 levandov Exp $";2 version="$Id: dmodapp.lib,v 1.23 2009-03-10 16:12:52 levandov Exp $"; 3 3 category="Noncommutative"; 4 4 info=" … … 93 93 example insertGenerator; 94 94 example deleteGenerator; 95 example bFactor; 95 96 } 96 97
Note: See TracChangeset
for help on using the changeset viewer.