////////////////////////////////////////////////////////////////////////////// version="$Id: bfun.lib,v 1.6 2009-03-10 16:27:55 Singular Exp $"; category="Noncommutative"; info=" LIBRARY: bfun.lib Algorithms for b-functions and Bernstein-Sato polynomial AUTHORS: Daniel Andres, daniel.andres@math.rwth-aachen.de @* Viktor Levandovskyy, levandov@math.rwth-aachen.de THEORY: Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R, @* one is interested in the global b-function (also known as Bernstein-Sato @* polynomial) b(s) in K[s], defined to be the monic polynomial of minimal @* degree, satisfying a functional identity L * F^{s+1} = b(s) F^s, @* for some operator L in D[s] (* stands for the action of differential operator) @* By D one denotes the n-th Weyl algebra @* K. @* One is interested in the following data: @* - Bernstein-Sato polynomial b(s) in K[s], @* - the list of its roots, which are known to be rational @* - the multiplicities of the roots. @* @* There is a general definition of a b-function of a holonomic ideal [SST] @* (that is, an ideal I in a Weyl algebra D, such that D/I is holonomic module) @* with respect to the given weight vector w: For a poly p in D, its initial @* form w.r.t. (-w,w) is defined as the sum of all terms of p which have @* maximal weighted total degree where the weight of x_i is -w_i and the weight @* of d_i is w_i. Let J be the initial ideal of I w.r.t. (-w,w), i.e. the @* K-vector space generated by all initial forms w.r.t (-w,w) of elements of I. @* Put s = w_1 x_1 d_1 + ... + w_n x_n d_n. Then the monic generator b_w(s) of @* the intersection of J with the PID K[s] is called the b-function of I w.r.t. w. @* Unlike Bernstein-Sato polynomial, general b-function with respect to @* arbitrary weights need not have rational roots at all. However, b-function @* of a holonomic ideal is known to be non-zero as well. @* @* References: [SST] Saito, Sturmfels, Takayama: Groebner Deformations of @* Hypergeometric Differential Equations (2000), @* Noro: An Efficient Modular Algorithm for Computing the Global b-function, @* (2002). MAIN PROCEDURES: bfct(f[,s,t,v]); compute the BS polynomial of f with linear reductions bfctSyz(f[,r,s,t,u,v]); compute the BS polynomial of f with syzygy-solver bfctAnn(f[,s]); compute the BS polynomial of f via Ann f^s + f bfctOneGB(f[,s,t]); compute the BS polynomial of f by just one GB computation bfctIdeal(I,w[,s,t]); compute the b-function of ideal wrt weight pIntersect(f,I[,s]); intersection of ideal with subalgebra K[f] for a poly f pIntersectSyz(f,I[,p,s,t]); intersection of ideal with subalgebra K[f] with syz-solver linReduce(f,I[,s]); reduce a poly by linear reductions wrt ideal linReduceIdeal(I,[s,t]); interreduce generators of ideal by linear reductions linSyzSolve(I[,s]); compute a linear dependency of elements of ideal I AUXILIARY PROCEDURES: allPositive(v); checks whether all entries of an intvec are positive scalarProd(v,w); computes the standard scalar product of two intvecs vec2poly(v[,i]); constructs an univariate poly with given coefficients SEE ALSO: dmod_lib, dmodapp_lib, gmssing_lib "; LIB "qhmoduli.lib"; // for Max LIB "dmod.lib"; // for SannfsBFCT etc LIB "dmodapp.lib"; // for initialIdealW etc LIB "nctools.lib"; // for isWeyl etc /////////////////////////////////////////////////////////////////////////////// // testing for consistency of the library: proc testbfunlib () { // tests all procs for consistency "AUXILIARY PROCEDURES:"; example allPositive; example scalarProd; example vec2poly; "MAIN PROCEDURES:"; example bfct; example bfctSyz; example bfctAnn; example bfctOneGB; example bfctIdeal; example pIntersect; example pIntersectSyz; example linReduce; example linReduceIdeal; example linSyzSolve; } //--------------- auxiliary procedures ---------------------------------------- static proc gradedWeyl (intvec u,intvec v) "USAGE: gradedWeyl(u,v); u,v intvecs RETURN: a ring, the associated graded ring of the basering w.r.t. u and v PURPOSE: compute the associated graded ring of the basering w.r.t. u and v ASSUME: basering is a Weyl algebra EXAMPLE: example gradedWeyl; shows examples NOTE: u[i] is the weight of x(i), v[i] the weight of D(i). @* u+v has to be a non-negative intvec. " { // todo check assumption int i; def save = basering; int n = nvars(save)/2; if (nrows(u)<>n || nrows(v)<>n) { ERROR("weight vectors have wrong dimension"); } intvec uv,gr; uv = u+v; for (i=1; i<=n; i++) { if (uv[i]>=0) { if (uv[i]==0) { gr[i] = 0; } else { gr[i] = 1; } } else { ERROR("the sum of the weight vectors has to be a non-negative intvec"); } } list l = ringlist(save); list l2 = l[2]; matrix l6 = l[6]; for (i=1; i<=n; i++) { if (gr[i] == 1) { l2[n+i] = "xi("+string(i)+")"; l6[i,n+i] = 0; } } l[2] = l2; l[6] = l6; def G = ring(l); return(G); } example { "EXAMPLE:"; echo = 2; ring @D = 0,(x,y,z,Dx,Dy,Dz),dp; def D = Weyl(); setring D; intvec u = -1,-1,1; intvec v = 2,1,1; def G = gradedWeyl(u,v); setring G; G; } static proc safeVarName (string s) { string S = "," + charstr(basering) + "," + varstr(basering) + ","; s = "," + s + ","; while (find(S,s) <> 0) { s[1] = "@"; s = "," + s; } s = s[2..size(s)-1]; return(s) } proc allPositive (intvec v) "USAGE: allPositive(v); v an intvec RETURN: int, 1 if all components of v are positive, or 0 otherwise PURPOSE: check whether all components of an intvec are positive EXAMPLE: example allPositive; shows an example " { int i; for (i=1; i<=size(v); i++) { if (v[i]<=0) { return(0); break; } } return(1); } example { "EXAMPLE:"; echo = 2; intvec v = 1,2,3; allPositive(v); intvec w = 1,-2,3; allPositive(w); } static proc findFirst (list l, i) "USAGE: findFirst(l,i); l a list, i an argument of any type RETURN: int, the position of the first appearance of i in l, @* or 0 if i is not a member of l PURPOSE: check whether the second argument is a member of a list EXAMPLE: example findFirst; shows an example " { int j; for (j=1; j<=size(l); j++) { if (l[j]==i) { return(j); break; } } return(0); } // findFirst(list(1, 2, list(1)),2); // seems to be a bit buggy, // findFirst(list(1, 2, list(1)),3); // but works for the purposes // findFirst(list(1, 2, list(1)),list(1)); // of this library // findFirst(l,list(2)); example { "EXAMPLE:"; echo = 2; ring r = 0,x,dp; list l = 1,2,3; findFirst(l,4); findFirst(l,2); } proc scalarProd (intvec v, intvec w) "USAGE: scalarProd(v,w); v,w intvecs RETURN: int, the standard scalar product of v and w PURPOSE: computes the scalar product of two intvecs ASSUME: the arguments are of the same size EXAMPLE: example scalarProd; shows examples " { int i; int sp; if (size(v)!=size(w)) { ERROR("non-matching dimensions"); } else { for (i=1; i<=size(v);i++) { sp = sp + v[i]*w[i]; } } return(sp); } example { "EXAMPLE:"; echo = 2; intvec v = 1,2,3; intvec w = 4,5,6; scalarProd(v,w); } //-------------- main procedures ------------------------------------------------------- proc linReduceIdeal(ideal I, list #) "USAGE: linReduceIdeal(I [,s,t,u]); I an ideal, s,t,u optional ints RETURN: ideal or list, linear reductum (over field) of I by its elements PURPOSE: reduces a list of polys only by linear reductions (no monomial @* multiplications) NOTE: If s<>0, a list consisting of the reduced ideal and the coefficient @* vectors of the used reductions given as module is returned. @* Otherwise (and by default), only the reduced ideal is returned. @* If t<>0 (and by default) all monomials are reduced (if possible), @* otherwise, only leading monomials are reduced. @* If u<>0 (and by default), the ideal is first sorted in increasing order. @* If u is set to 0 and the given ideal is not sorted in the way described, @* the result might not be as expected. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example linReduceIdeal; shows examples " { // #[1] = remembercoeffs // #[2] = redtail // #[3] = sortideal int ppl = printlevel - voice + 2; int remembercoeffs = 0; // default int redtail = 1; // default int sortideal = 1; // default if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { remembercoeffs = #[1]; } if (size(#)>1) { if (typeof(#[2])=="int" || typeof(#[2])=="number") { redtail = #[2]; } if (size(#)>2) { if (typeof(#[3])=="int" || typeof(#[3])=="number") { sortideal = #[3]; } } } } int sI = ncols(I); int sZeros = sI - size(I); int i,j; ideal J,lmJ,ordJ; list lJ = sort(I); module M; // for the coefficients // step 1: prepare, e.g. sort I if (sortideal <> 0) { if (sZeros > 0) // I contains zeros { if (remembercoeffs <> 0) { j = 0; for (i=1; i<=sI; i++) { if (I[i] == 0) { j++; J[j] = 0; ordJ[j] = -1; M[j] = gen(i); } else { M[i+sZeros-j] = gen(lJ[2][i-j]+j); } } } else { for (i=1; i<=sZeros; i++) { J[i] = 0; ordJ[i] = -1; } } I = J,lJ[1]; } else // I contains no zeros { I = lJ[1]; if (remembercoeffs <> 0) { for (i=1; i<=size(lJ[1]); i++) { M[i+sZeros] = gen(lJ[2][i]); } } } } else // assume I is already sorted { if (remembercoeffs <> 0) { for (i=1; i<=ncols(I); i++) { M[i] = gen(i); } } } dbprint(ppl-1,"// initially sorted ideal:", I); if (remembercoeffs <> 0) { dbprint(ppl-1,"// used permutations:", M); } // step 2: reduce leading monomials by linear reductions poly lm,c,redpoly,maxlmJ; J[sZeros+1] = I[sZeros+1]; lm = I[sZeros+1]; maxlmJ = leadmonom(lm); lmJ[sZeros+1] = maxlmJ; int ordlm,reduction; ordJ[sZeros+1] = ord(lm); vector v; int noRedPast; for (i=sZeros+2; i<=sI; i++) { redpoly = I[i]; lm = leadmonom(redpoly); ordlm = ord(lm); if (remembercoeffs <> 0) { v = M[i]; } reduction = 1; while (reduction == 1) // while there was a reduction { noRedPast = i; reduction = 0; for (j=sZeros+1; j maxlmJ) { break; } //lm is bigger than maximal monomial to reduce with if (ordlm == ordJ[j]) { if (lm == lmJ[j]) { dbprint(ppl-1,"// reducing " + string(redpoly)); dbprint(ppl-1,"// with " + string(J[j])); c = leadcoef(redpoly)/leadcoef(J[j]); redpoly = redpoly - c*J[j]; dbprint(ppl-1,"// to " + string(redpoly)); lm = leadmonom(redpoly); ordlm = ord(lm); if (remembercoeffs <> 0) { M[i] = M[i] - c * M[j]; } noRedPast = j; reduction = 1; } } } } for (j=sZeros+1; j 0) { v = M[i]; M = deleteGenerator(M,i); M = insertGenerator(M,v,j); } maxlmJ = lmJ[i]; } // step 3: reduce tails by linear reductions as well if (redtail <> 0) { dbprint(ppl,"// finished reducing leading monomials"); dbprint(ppl-1,string(J)); if (remembercoeffs <> 0) { dbprint(ppl-1,"// used reductions:" + string(M)); } int k; for (i=sZeros+1; i<=sI; i++) { lm = lmJ[i]; for (j=i+1; j<=sI; j++) { for (k=2; k<=size(J[j]); k++) // run over all terms in J[j] { if (ordJ[i] == ord(J[j][k])) { if (lm == normalize(J[j][k])) { c = leadcoef(J[j][k])/leadcoef(J[i]); dbprint(ppl-1,"// reducing " + string(J[j])); dbprint(ppl-1,"// with " + string(J[i])); J[j] = J[j] - c*J[i]; dbprint(ppl-1,"// to " + string(J[j])); if (remembercoeffs <> 0) { M[j] = M[j] - c * M[i]; } } } } } } } if (remembercoeffs <> 0) { return(list(J,M)); } else { return(J); } } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; ideal I = 3,x+9,y4+5x,2y4+7x+2; linReduceIdeal(I); // reduces tails linReduceIdeal(I,0,0); // no reductions of tails list l = linReduceIdeal(I,1); // reduces tails and shows reductions used l; module M = I; l[1] - ideal(M*l[2]); } proc linReduce(poly f, ideal I, list #) "USAGE: linReduce(f, I [,s,t,u]); f a poly, I an ideal, s,t,u optional ints RETURN: poly or list, linear reductum (over field) of f by elements from I PURPOSE: reduce a poly only by linear reductions (no monomial multiplications) NOTE: If s<>0, a list consisting of the reduced poly and the coefficient @* vector of the used reductions is returned, otherwise (and by default) @* only reduced poly is returned. @* If t<>0 (and by default) all monomials are reduced (if possible), @* otherwise, only leading monomials are reduced. @* If u<>0 (and by default), the ideal is linearly pre-reduced, i.e. @* instead of the given ideal, the output of @code{linReduceIdeal} is used. @* If u is set to 0 and the given ideal does not equal the output of @* @code{linReduceIdeal}, the result might not be as expected. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example linReduce; shows examples " { int ppl = printlevel - voice + 2; int remembercoeffs = 0; // default int redtail = 1; // default int prepareideal = 1; // default if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { remembercoeffs = #[1]; } if (size(#)>1) { if (typeof(#[2])=="int" || typeof(#[2])=="number") { redtail = #[2]; } if (size(#)>2) { if (typeof(#[3])=="int" || typeof(#[3])=="number") { prepareideal = #[3]; } } } } int i,j,k; int sI = ncols(I); // pre-reduce I: module M; if (prepareideal <> 0) { if (remembercoeffs <> 0) { list sortedI = linReduceIdeal(I,1,redtail); I = sortedI[1]; M = sortedI[2]; dbprint(ppl-1,"// prepared ideal:" +string(I)); dbprint(ppl-1,"// with operations:" +string(M)); } else { I = linReduceIdeal(I,0,redtail); } } else { if (remembercoeffs <> 0) { for (i=1; i<=sI; i++) { M[i] = gen(i); } } } ideal lmI,lcI,ordI; for (i=1; i<=sI; i++) { lmI[i] = leadmonom(I[i]); lcI[i] = leadcoef(I[i]); ordI[i] = ord(lmI[i]); } vector v; poly c; // === reduce leading monomials === poly lm = leadmonom(f); int ordf = ord(lm); for (i=sI; i>=1; i--) // I is sorted: smallest lm's on top { if (lm == 0) { break; } else { if (ordf == ordI[i]) { if (lm == lmI[i]) { c = leadcoef(f)/lcI[i]; f = f - c*I[i]; lm = leadmonom(f); ordf = ord(lm); if (remembercoeffs <> 0) { v = v - c * M[i]; } } } } } // === reduce tails as well === if (redtail <> 0) { dbprint(ppl,"// finished reducing leading monomials"); dbprint(ppl-1,string(f)); if (remembercoeffs <> 0) { dbprint(ppl-1,"// used reductions:" + string(v)); } for (i=1; i<=sI; i++) { dbprint(ppl,"// testing ideal entry "+string(i)); for (j=1; j<=size(f); j++) { if (ord(f[j]) == ordI[i]) { if (normalize(f[j]) == lmI[i]) { c = leadcoef(f[j])/lcI[i]; f = f - c*I[i]; dbprint(ppl-1,"// reducing with " + string(I[i])); dbprint(ppl-1,"// to " + string(f)); if (remembercoeffs <> 0) { v = v - c * M[i]; } } } } } } if (remembercoeffs <> 0) { list l = f,v; return(l); } else { return(f); } } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; ideal I = 1,y,xy; poly f = 5xy+7y+3; poly g = 7x+5y+3; linReduce(g,I); // reduces tails linReduce(g,I,0,0); // no reductions of tails linReduce(f,I,1); // reduces tails and shows reductions used f = x3+y2+x2+y+x; I = x3-y3, y3-x2,x3-y2,x2-y,y2-x; list l = linReduce(f,I,1); l; module M = I; f - (l[1]-(M*l[2])[1,1]); } proc linSyzSolve (ideal I, list #) "USAGE: linSyzSolve(I[,s]); I an ideal, s an optional int RETURN: vector, coefficient vector of linear combination of 0 in elements of I PURPOSE: compute a linear dependency between the elements of an ideal @* if such one exists NOTE: If s<>0, @code{std} is used for Groebner basis computations, @* otherwise, @code{slimgb} is used. @* By default, @code{slimgb} is used in char 0 and @code{std} in char >0. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example linSyzSolve; shows examples " { int whichengine = 0; // default int enginespecified = 0; // default if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { whichengine = int( #[1]); enginespecified = 1; } } int ppl = printlevel - voice +2; int sI = ncols(I); // check if we are done if (I[sI]==0) { vector v = gen(sI); } else { // ------- 1. introduce undefined coeffs ------------------ def save = basering; list RL = ringlist(save); int nv = nvars(save); int np = npars(save); int p = char(save); string cs = "(" + charstr(save) + ")"; if (enginespecified == 0) { if (p <> 0) { whichengine = 1; } } int i; list Lvar; for (i=1; i<=sI; i++) { Lvar[i] = safeVarName("@a(" + string(i) + ")"); } list L@A = RL[1..4]; L@A[2] = Lvar; L@A[3] = list(list("lp",intvec(1:sI)),list("C",intvec(0))); def @A = ring(L@A); L@A[2] = list(safeVarName("@z")); L@A[3][1] = list("dp",intvec(1)); if (size(L@A[1])>1) { L@A[1][2] = L@A[1][2] + Lvar; L@A[1][3][size(L@A[1][3])+1] = list("lp",intvec(1:sI)); } else { L@A[1] = list(p,Lvar,list(list("lp",intvec(1:sI))),ideal(0)); } def @aA = ring(L@A); def @B = save + @aA; setring @B; ideal I = imap(save,I); // ------- 2. form the linear system for the undef coeffs --- poly W; ideal QQ; for (i=1; i<=sI; i++) { W = W + par(np+i)*I[i]; } while (W!=0) { QQ = QQ,leadcoef(W); W = W - lead(W); } // QQ consists of polynomial expressions in @a(i) of type number setring @A; ideal QQ = imap(@B,QQ); // ------- 3. this QQ is a polynomial ideal, so "solve" the system ----- dbprint(ppl, "// linSyzSolve: starting Groebner basis computation with engine:", whichengine); QQ = engine(QQ,whichengine); dbprint(ppl, "// QQ after engine:", QQ); if (dim(QQ) == -1) { dbprint(ppl+1, "// no solutions by linSyzSolve"); // output zeroes setring save; kill @A,@aA,@B; return(v); } // ------- 4. in order to get the numeric values ------- matrix AA = matrix(maxideal(1)); module MQQ = std(module(QQ)); AA = NF(AA,MQQ); // todo: we still receive NF warnings dbprint(ppl, "// AA after NF:",AA); // "AA after NF:"; print(AA); for(i=1; i<=sI; i++) { AA = subst(AA,var(i),1); } dbprint(ppl, "// AA after subst:",AA); // "AA after subst: "; print(AA); vector v = (module(transpose(AA)))[1]; setring save; vector v = imap(@A,v); kill @A,@aA,@B; } return(v); } example { "EXAMPLE:"; echo = 2; ring r = 0,x,dp; ideal I = x,2x; linSyzSolve(I); ideal J = x,x2; linSyzSolve(J); } proc pIntersect (poly s, ideal I, list #) "USAGE: pIntersect(f, I [,s]); f a poly, I an ideal, s an optional int RETURN: vector, coefficient vector of the monic polynomial PURPOSE: compute the intersection of ideal I with the subalgebra K[f] ASSUME: I is given as Groebner basis, basering is not a qring. NOTE: If the intersection is zero, this proc might not terminate. @* If s>0 is given, it is searched for the generator of the intersection @* only up to degree s. Otherwise (and by default), no bound is assumed. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example pIntersect; shows examples " { def save = basering; int n = nvars(save); list RL = ringlist(save); int i,j,k; if (RL[4]<>0) { ERROR ("basering must not be a qring"); } // assume I is given in Groebner basis if (attrib(I,"isSB") <> 1) { print("// WARNING: The input has no SB attribute!"); print("// Treating it as if it were a Groebner basis and proceeding..."); attrib(I,"isSB",1); // set attribute for suppressing NF messages } int bound = 0; // default if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { bound = #[1]; } } int ppl = printlevel-voice+2; // ---case 1: I = basering--- if (size(I) == 1) { if (simplify(I,2)[1] == 1) { return(gen(1)); // = 1 } } // ---case 2: intersection is zero--- intvec degs = leadexp(s); intvec possdegbounds; list degI; i = 1; while (i <= ncols(I)) { if (i == ncols(I)+1) { break; } degI[i] = leadexp(I[i]); for (j=1; j<=n; j++) { if (degs[j] == 0) { if (degI[i][j] <> 0) { break; } } if (j == n) { k++; possdegbounds[k] = Max(degI[i]); } } i++; } int degbound = Min(possdegbounds); if (bound>0 && bound0 && i>bound) { return(vector(0)); } dbprint(ppl,"// Testing degree "+string(i)); if (i>1) { oldNF = newNF; tobracket = s^(i-1) - oldNF; if (tobracket==0) // todo bug in bracket? { toNF = 0; } else { toNF = bracket(tobracket,secNF); } newNF = NF(toNF+oldNF*secNF,I); // = NF(s^i,I) } else { newNF = NF(s,I); secNF = newNF; } ll = linReduce(newNF,redNI,1); rednewNF = ll[1]; l[i+1] = ll[2]; dbprint(ppl-1,"// newNF is: "+string(newNF)); dbprint(ppl-1,"// rednewNF is: "+string(rednewNF)); if (rednewNF != 0) // no linear dependency { redNI[i+1] = rednewNF; i++; } else // there is a linear dependency, hence we are done { dbprint(ppl,"// degree of the generator of the intersection is: "+string(i)); break; } } dbprint(ppl-1,"// used linear reductions:", l); // we obtain the coefficients of the generator by the used reductions: list Lvar; for (j=1; j<=i+1; j++) { Lvar[j] = safeVarName("a("+string(j)+")"); } list Lord = list("dp",intvec(1:(i+1))),list("C",intvec(0)); list L@R = RL[1..4]; L@R[2] = Lvar; L@R[3] = Lord; def @R = ring(L@R); setring @R; list l = imap(save,l); ideal C; for (j=1;j<=i+1;j++) { C[j] = 0; for (k=1;k<=j;k++) { C[j] = C[j]+l[j][k]*var(k); } } for (j=i;j>=1;j--) { C[i+1]= subst(C[i+1],var(j),var(j)+C[j]); } matrix m = coeffs(C[i+1],maxideal(1)); vector v = gen(i+1); for (j=1;j<=i+1;j++) { v = v + m[j,1]*gen(j); } setring save; v = imap(@R,v); kill @R; return(v); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; poly f = x^2+y^3+x*y^2; def D = initialMalgrange(f); setring D; inF; pIntersect(t*Dt,inF); pIntersect(t*Dt,inF,1); } proc pIntersectSyz (poly s, ideal I, list #) "USAGE: pIntersectSyz(f, I [,p,s,t]); f poly, I ideal, p,t optial ints, p prime RETURN: vector, coefficient vector of the monic polynomial PURPOSE: compute the intersection of an ideal I with the subalgebra K[f] ASSUME: I is given as Groebner basis. NOTE: If the intersection is zero, this procedure might not terminate. @* If p>0 is given, this proc computes the generator of the intersection in @* char p first and then only searches for a generator of the obtained @* degree in the basering. Otherwise, it searches for all degrees by @* computing syzygies. @* If s<>0, @code{std} is used for Groebner basis computations in char 0, @* otherwise, and by default, @code{slimgb} is used. @* If t<>0 and by default, @code{std} is used for Groebner basis @* computations in char >0, otherwise, @code{slimgb} is used. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example pIntersectSyz; shows examples " { // assume I is given in Groebner basis if (attrib(I,"isSB") <> 1) { print("// WARNING: The input has no SB attribute!"); print("// Treating it as if it were a Groebner basis and proceeding..."); attrib(I,"isSB",1); // set attribute for suppressing NF messages } int ppl = printlevel-voice+2; int whichengine = 0; // default int modengine = 1; // default int solveincharp = 0; // default def save = basering; int n = nvars(save); if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { solveincharp = int(#[1]); } if (size(#)>1) { if (typeof(#[2])=="int" || typeof(#[2])=="number") { whichengine = int(#[2]); } if (size(#)>2) { if (typeof(#[3])=="int" || typeof(#[3])=="number") { modengine = int(#[3]); } } } } int i,j; vector v; poly tobracket,toNF,newNF,p; ideal NI = 1; newNF = NF(s,I); NI[2] = newNF; list RL = ringlist(save); if (solveincharp) { int psolveincharp = prime(solveincharp); if (solveincharp <> psolveincharp) { print("// " + string(solveincharp) + " is invalid characteristic of ground field."); print("// " + string(psolveincharp) + " is used."); solveincharp = psolveincharp; kill psolveincharp; } list RLp = RL[1..4]; if (typeof(RL[1]) == "int") { RLp[1] = solveincharp; } else { RLp[1][1] = solveincharp; } // parameters def @Rp = ring(RLp); setring @Rp; number c; setring save; int shortSave = short; // workaround for maps Q(a_i) -> Z/p(a_i) short = 0; string str; int badprime; i=1; while (badprime == 0 && i<=size(s)) // detect bad primes { str = string(denominator(leadcoef(s[i]))); str = "c = " + str + ";"; setring @Rp; execute(str); if (c == 0) { badprime = 1; } setring save; i++; } str = "poly s = " + string(s) + ";"; if (size(RL) > 4) // basering is NC-algebra { string RL5 = "@C = " + string(RL[5]) + ";"; string RL6 = "@D = " + string(RL[6]) + ";"; setring @Rp; matrix @C[n][n]; matrix @D[n][n]; execute(RL5); execute(RL6); def Rp = nc_algebra(@C,@D); } else { def Rp = @Rp; } setring Rp; kill @Rp; dbprint(ppl-1,"// solving in ring ", Rp); execute(str); vector v; number c; ideal NI = 1; setring save; i=1; while (badprime == 0 && i<=size(I)) // detect bad primes { str = string(leadcoef(cleardenom(I[i]))); str = "c = " + str + ";"; setring Rp; execute(str); if (c == 0) { badprime = 1; } setring save; i++; } if (badprime == 1) { print("// WARNING: bad prime"); short = shortSave; return(vector(0)); } } i = 1; dbprint(ppl,"// pIntersectSyz starts..."); dbprint(ppl-1,"// with ideal I=", I); while (1) { dbprint(ppl,"// testing degree: "+string(i)); if (i>1) { tobracket = s^(i-1)-NI[i]; if (tobracket!=0) { toNF = bracket(tobracket,NI[2]) + NI[i]*NI[2]; } else { toNF = NI[i]*NI[2]; } newNF = NF(toNF,I); NI[i+1] = newNF; } // look for a solution dbprint(ppl-1,"// linSyzSolve starts with: "+string(matrix(NI))); if (solveincharp) // modular method { for (j=1; j<=size(newNF); j++) { str = string(denominator(leadcoef(s[i]))); str = "c = " + str + ";"; setring Rp; execute(str); if (c == 0) { print("// WARNING: bad prime"); setring save; short = shortSave; return(vector(0)); } setring save; } str = "NI[" + string(i) + "+1] = " + string(newNF) + ";"; setring Rp; execute(str); // NI[i+1] = [newNF]_{solveincharp} v = linSyzSolve(NI,modengine); if (v!=0) // there is a modular solution { dbprint(ppl,"// got solution in char "+string(solveincharp)+" of degree "+string(i)); setring save; v = linSyzSolve(NI,whichengine); if (v==0) { break; } } else // no modular solution { setring save; v = 0; } } else // non-modular method { v = linSyzSolve(NI,whichengine); } matrix MM[1][nrows(v)] = matrix(v); dbprint(ppl-1,"// linSyzSolve ready with: "+string(MM)); kill MM; // "linSyzSolve ready with"; print(v); if (v!=0) { // a solution: //check for the reality of the solution p = 0; for (j=1; j<=i+1; j++) { p = p + v[j]*NI[j]; } if (p!=0) { dbprint(ppl,"// linSyzSolve: bad solution!"); } else { dbprint(ppl,"// linSyzSolve: got solution!"); // "got solution!"; break; } } // no solution: i++; } dbprint(ppl,"// pIntersectSyz finished"); if (solveincharp) { short = shortSave; } return(v); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; poly f = x^2+y^3+x*y^2; def D = initialMalgrange(f); setring D; inF; poly s = t*Dt; pIntersectSyz(s,inF); int p = prime(20000); pIntersectSyz(s,inF,p,0,0); } proc vec2poly (list #) "USAGE: vec2poly(v [,i]); v a vector or an intvec, i an optional int RETURN: poly, an univariate poly in i-th variable with coefficients given by v PURPOSE: constructs an univariate poly in K[var(i)] with given coefficients, @* such that the coefficient at var(i)^{j-1} is v[j]. NOTE: The optional argument i must be positive, by default i is 1. EXAMPLE: example vec2poly; shows examples " { def save = basering; int i,ringvar; ringvar = 1; // default if (size(#) > 0) { if (typeof(#[1])=="vector" || typeof(#[1])=="intvec") { def v = #[1]; } else { ERROR("wrong input: expected vector/intvec expression"); } if (size(#) > 1) { if (typeof(#[2])=="int" || typeof(#[2])=="number") { ringvar = int(#[2]); } } } if (ringvar > nvars(save)) { ERROR("var out of range"); } poly p; for (i=1; i<=nrows(v); i++) { p = p + v[i]*(var(ringvar))^(i-1); } return(p); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; vector v = gen(1) + 3*gen(3) + 22/9*gen(4); intvec iv = 3,2,1; vec2poly(v,2); vec2poly(iv); } static proc listofroots (list #) { def save = basering; int n = nvars(save); int i; poly p; if (typeof(#[1])=="vector") { vector b = #[1]; for (i=1; i<=nrows(b); i++) { p = p + b[i]*(var(1))^(i-1); } } else { p = #[1]; } int substitution = int(#[2]); string s = safeVarName("s"); list RL = ringlist(save); RL = RL[1..4]; RL[2] = list(s); RL[3] = list(list("dp",intvec(1)),list("C",0)); def S = ring(RL); setring S; ideal J; for (i=1; i<=n; i++) { J[i] = var(1); } map @m = save,J; poly p = @m(p); if (substitution == 1) { p = subst(p,var(1),-var(1)-1); } // the rest of this proc is nicked from bernsteinBM from dmod.lib list P = factorize(p);//with constants and multiplicities ideal bs; intvec m; //the BS poly is monic, so we are not interested in constants for (i=2; i<= size(P[1]); i++) //we delete P[1][1] and P[2][1] { bs[i-1] = P[1][i]; m[i-1] = P[2][i]; } bs = normalize(bs); bs = -subst(bs,var(1),0); setring save; ideal bs = imap(S,bs); kill S; list BS = bs,m; return(BS); } static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodpIntersect, int pIntersectchar, int modengine, intvec u0) { int ppl = printlevel - voice +2; int i; def save = basering; int n = nvars(save); if (isCommutative() == 0) { ERROR("basering must be commutative"); } if (char(save) <> 0) { ERROR("characteristic of basering has to be 0"); } list L = ringlist(save); int qr; if (L[4] <> 0) // qring { print("// basering is qring:"); print("// discarding the quotient and proceeding..."); L[4] = ideal(0); qr = 1; def save2 = ring(L); setring save2; poly f = imap(save,f); } if (size(variables(f)) == 0) // f is constant { return(list(ideal(0),intvec(0))); } if (inorann == 0) // bfct using initial ideal { def D = initialMalgrange(f,whichengine,methodord,u0); setring D; ideal J = inF; kill inF; poly s = t*Dt; } else // bfct using Ann(f^s) { def D = SannfsBFCT(f,whichengine); setring D; ideal J = LD; kill LD; poly s = var(2*n+1); } vector b; // try it modular if (methodpIntersect <> 0) // pIntersectSyz { if (pIntersectchar == 0) // pIntersectSyz::modular { list L = ringlist(D); int lb = 10000; int ub = 536870909; // bounds for random primes list usedprimes; int sJ = size(J); int sLJq; ideal LJ; for (i=1; i<=sJ; i++) { LJ[i] = leadcoef(cleardenom(J[i])); } int short_save = short; // workaround for map Q(a_i) -> Z/q(a_i) short = 0; string strLJq = "ideal LJq = " + string(LJ) + ";"; int nD = nvars(D); string L5 = "matrix @C[nD][nD]; @C = " + string(L[5]) + ";"; string L6 = "matrix @D[nD][nD]; @D = " + string(L[6]) + ";"; L = L[1..4]; i = 1; while (b == 0) { dbprint(ppl,"// number of run in the loop: "+string(i)); int q = prime(random(lb,ub)); if (findFirst(usedprimes,q)==0) // if q has not been used already { usedprimes = usedprimes,q; dbprint(ppl,"// using prime: "+string(q)); if (typeof(L[1]) == "int") { L[1] = q; } else { L[1][1] = q; } // parameters def @Rq = ring(L); setring @Rq; execute(L5); execute(L6); def Rq = nc_algebra(@C,@D); // def Rq = nc_algebra(1,@D); setring Rq; kill @Rq; execute(strLJq); sLJq = size(LJq); setring D; kill Rq; if (sLJq <> sJ ) // detect unlucky prime { dbprint(ppl,"// " +string(q) + " is unlucky"); b = 0; } else { b = pIntersectSyz(s,J,q,whichengine,modengine); } } i++; } short = short_save; } else // pIntersectSyz::non-modular { b = pIntersectSyz(s,J,0,whichengine); } } else // pIntersect: linReduce { b = pIntersect(s,J); } if (inorann == 0) { list l = listofroots(b,1); } else { list l = listofroots(b,0); } setring save; list l = imap(D,l); return(l); } proc bfct (poly f, list #) "USAGE: bfct(f [,s,t,v]); f a poly, s,t optional ints, v an optional intvec RETURN: list of ideal and intvec PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) @* for the hypersurface defined by f. ASSUME: The basering is commutative and of characteristic 0. BACKGROUND: In this proc, the initial Malgrange ideal is computed according to @* the algorithm by Masayuki Noro and then a system of linear equations is @* solved by linear reductions. NOTE: In the output list, the ideal contains all the roots @* and the intvec their multiplicities. @* If s<>0, @code{std} is used for GB computations, @* otherwise, and by default, @code{slimgb} is used. @* If t<>0, a matrix ordering is used for Groebner basis computations, @* otherwise, and by default, a block ordering is used. @* If v is a positive weight vector, v is used for homogenization @* computations, otherwise and by default, no weights are used. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example bfct; shows examples " { int ppl = printlevel - voice +2; int i; int n = nvars(basering); // in # we have two switches: // one for the engine used for Groebner basis computations, // one for M() ordering or its realization // in # can also be the optional weight vector int whichengine = 0; // default int methodord = 0; // default intvec u0 = 0; // default if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { whichengine = int(#[1]); } if (size(#)>1) { if (typeof(#[2])=="int" || typeof(#[2])=="number") { methodord = int(#[2]); } if (size(#)>2) { if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1) { u0 = #[3]; } } } } list b = bfctengine(f,0,whichengine,methodord,0,0,0,u0); return(b); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; poly f = x^2+y^3+x*y^2; bfct(f); intvec v = 3,2; bfct(f,1,0,v); } proc bfctSyz (poly f, list #) "USAGE: bfctSyz(f [,r,s,t,u,v]); f poly, r,s,t,u optional ints, v opt. intvec RETURN: list of ideal and intvec PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) @* for the hypersurface defined by f ASSUME: The basering is commutative and of characteristic 0. BACKGROUND: In this proc, the initial Malgrange ideal is computed according to @* the algorithm by Masayuki Noro and then a system of linear equations is @* solved by computing syzygies. NOTE: In the output list, the ideal contains all the roots and the intvec @* their multiplicities. @* If r<>0, @code{std} is used for GB computations in characteristic 0, @* otherwise, and by default, @code{slimgb} is used. @* If s<>0, a matrix ordering is used for GB computations, otherwise, @* and by default, a block ordering is used. @* If t<>0, the computation of the intersection is solely performed over @* charasteristic 0, otherwise and by default, a modular method is used. @* If u<>0 and by default, @code{std} is used for GB computations in @* characteristic >0, otherwise, @code{slimgb} is used. @* If v is a positive weight vector, v is used for homogenization @* computations, otherwise and by default, no weights are used. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example bfctSyz; shows examples " { int ppl = printlevel - voice +2; int i; // in # we have four switches: // one for the engine used for Groebner basis computations in char 0, // one for M() ordering or its realization // one for a modular method when computing the intersection // and one for the engine used for Groebner basis computations in char >0 // in # can also be the optional weight vector int n = nvars(basering); int whichengine = 0; // default int methodord = 0; // default int pIntersectchar = 0; // default int modengine = 1; // default intvec u0 = 0; // default if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { whichengine = int(#[1]); } if (size(#)>1) { if (typeof(#[2])=="int" || typeof(#[2])=="number") { methodord = int(#[2]); } if (size(#)>2) { if (typeof(#[3])=="int" || typeof(#[3])=="number") { pIntersectchar = int(#[3]); } if (size(#)>3) { if (typeof(#[4])=="int" || typeof(#[4])=="number") { modengine = int(#[4]); } if (size(#)>4) { if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1) { u0 = #[5]; } } } } } } list b = bfctengine(f,0,whichengine,methodord,1,pIntersectchar,modengine,u0); return(b); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; poly f = x^2+y^3+x*y^2; bfctSyz(f); intvec v = 3,2; bfctSyz(f,0,1,1,0,v); } proc bfctIdeal (ideal I, intvec w, list #) "USAGE: bfctIdeal(I,w[,s,t]); I an ideal, w an intvec, s,t optional ints RETURN: list of ideal and intvec PURPOSE: computes the roots of the global b-function of I wrt the weight (-w,w). ASSUME: The basering is the n-th Weyl algebra in characteristic 0 and for all @* 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the @* sequence of variables is given by x(1),...,x(n),D(1),...,D(n), @* where D(i) is the differential operator belonging to x(i). @* Further we assume that I is holonomic. BACKGROUND: In this proc, the initial ideal of I is computed according to the @* algorithm by Masayuki Noro and then a system of linear equations is @* solved by linear reductions. NOTE: In the output list, say L, @* - L[1] of type ideal contains all the rational roots of a b-function, @* - L[2] of type intvec contains the multiplicities of above roots, @* - optional L[3] of type string is the part of b-function without rational roots. @* Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and @* L[3] is 1 (for nonzero constant) or 0 (for zero b-function). @* If s<>0, @code{std} is used for GB computations in characteristic 0, @* otherwise, and by default, @code{slimgb} is used. @* If t<>0, a matrix ordering is used for GB computations, otherwise, @* and by default, a block ordering is used. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example bfctIdeal; shows examples " { int ppl = printlevel - voice +2; int i; def save = basering; int n = nvars(save)/2; int whichengine = 0; // default int methodord = 0; // default if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { whichengine = int(#[1]); } if (size(#)>1) { if (typeof(#[2])=="int" || typeof(#[2])=="number") { methodord = int(#[2]); } } } if (isWeyl()==0) { ERROR("basering is not a Weyl algebra"); } for (i=1; i<=n; i++) { if (bracket(var(i+n),var(i))<>1) { ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i))); } } int isH = isHolonomic(I); if (isH<>1) { print("WARNING: given ideal is not holonomic"); print("... setting bound for degree of b-function to 10 and proceeding"); isH = 10; } else { isH = 0; } // no degree bound for pIntersect ideal J = initialIdealW(I,-w,w,whichengine,methodord); poly s; for (i=1; i<=n; i++) { s = s + w[i]*var(i)*var(n+i); } vector b = pIntersect(s,J,isH); list RL = ringlist(save); RL = RL[1..4]; RL[2] = list(safeVarName("s")); RL[3] = list(list("dp",intvec(1)),list("C",intvec(0))); def @S = ring(RL); setring @S; vector b = imap(save,b); poly bs = vec2poly(b); list l = bFactor(bs); setring save; list l = imap(@S,l); return(l); } example { "EXAMPLE:"; echo = 2; ring @D = 0,(x,y,Dx,Dy),dp; def D = Weyl(); setring D; ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6; I = std(I); intvec w1 = 0,1; intvec w2 = 2,3; bfctIdeal(I,w1); bfctIdeal(I,w2,0,1); ideal J = I[size(I)]; // J is not holonomic by construction bfctIdeal(J,w1); // b-function of D/J w.r.t. w1 is non-zero bfctIdeal(J,w2); // b-function of D/J w.r.t. w2 is zero } proc bfctOneGB (poly f,list #) "USAGE: bfctOneGB(f [,s,t]); f a poly, s,t optional ints RETURN: list of ideal and intvec PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the @* hypersurface defined by f, using only one GB computation ASSUME: The basering is commutative and of characteristic 0. BACKGROUND: In this proc, the initial Malgrange ideal is computed based on the @* algorithm by Masayuki Noro and combined with an elimination ordering. NOTE: In the output list, the ideal contains all the roots and the intvec @* their multiplicities. @* If s<>0, @code{std} is used for the GB computation, otherwise, @* and by default, @code{slimgb} is used. @* If t<>0, a matrix ordering is used for GB computations, @* otherwise, and by default, a block ordering is used. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example bfctOneGB; shows examples " { int ppl = printlevel - voice +2; if (!isCommutative()) { ERROR("Basering must be commutative"); } def save = basering; int n = nvars(save); if (char(save) <> 0) { ERROR("characteristic of basering has to be 0"); } list L = ringlist(save); int qr; if (L[4] <> 0) // qring? { print("// basering is qring:"); print("// discarding the quotient and proceeding..."); L[4] = ideal(0); qr = 1; def save2 = ring(L); setring save2; poly f = imap(save,f); } int N = 2*n+4; int i; int whichengine = 0; // default int methodord = 0; // default if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { whichengine = int(#[1]); } if (size(#)>1) { if (typeof(#[2])=="int" || typeof(#[2])=="number") { methodord = int(#[2]); } } } // creating the homogenized extended Weyl algebra // create names for vars list Lvar; Lvar[1] = safeVarName("t"); Lvar[2] = safeVarName("s"); Lvar[n+3] = safeVarName("D"+Lvar[1]); Lvar[N] = safeVarName("h"); for (i=1; i<=n; i++) { Lvar[i+2] = string(var(i)); Lvar[i+n+3] = safeVarName("D" + string(var(i))); } // create ordering intvec uv = -1; uv[n+3] = 1; uv[N] = 0; intvec @a = 1:N; @a[2] = 2; intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0; list Lord; Lord[1] = list("a",@a); Lord[2] = list("a",@a2); if (methodord == 0) // default: block ordering { //ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(N-1),lp(1)); Lord[3] = list("a",uv); Lord[4] = list("dp",intvec(1:(N-1))); Lord[5] = list("lp",intvec(1)); Lord[6] = list("C",intvec(0)); } else // M() ordering { intmat @Ord[N][N]; @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1); for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; } dbprint(ppl,"// weights for ordering: "+string(transpose(@a))); dbprint(ppl,"// the ordering matrix:",@Ord); //ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord)); Lord[3] = list("M",intvec(@Ord)); Lord[4] = list("C",intvec(0)); } // create commutative ring list L@Dh = ringlist(basering); L@Dh = L@Dh[1..4]; // if basering is commutative nc_algebra L@Dh[2] = Lvar; L@Dh[3] = Lord; def @Dh = ring(L@Dh); setring @Dh; dbprint(ppl,"// the ring @Dh:",@Dh); // create non-commutative relations matrix @relD[N][N]; @relD[1,2] = var(1)*var(N)^2; // s*t = t*s + t*h^2 @relD[2,n+3] = var(n+3)*var(N)^2; // Dt*s = s*Dt+Dt*h^2 @relD[1,n+3] = var(N)^2; for (i=1; i<=n; i++) { @relD[i+2,n+3+i] = var(N)^2; } dbprint(ppl,"// nc relations:",@relD); def Dh = nc_algebra(1,@relD); setring Dh; kill @Dh; dbprint(ppl,"// computing in ring",Dh); poly f = imap(save,f); f = homog(f,h); // create the Malgrange ideal ideal I = var(1) - f, var(1)*var(n+3) - var(2); for (i=1; i<=n; i++) { I[3+i] = var(i+n+3)+diff(f,var(i+2))*var(n+3); } dbprint(ppl-1, "// the Malgrange ideal: " +string(I)); // the hard part: Groebner basis computation dbprint(ppl, "// starting Groebner basis computation with engine: "+string(whichengine)); I = engine(I, whichengine); dbprint(ppl, "// finished Groebner basis computation"); I = subst(I,h,1); // dehomogenization dbprint(ppl-1,string(I)); // 3.3 the initial form I = inForm(I,uv); dbprint(ppl, "// the initial ideal:", string(matrix(I))); // read off the solution intvec tonselect = 1; for (i=3; i<=2*n+4; i++) { tonselect = tonselect,i; } I = nselect(I,tonselect); dbprint(ppl, "// generators containing only s:", string(matrix(I))); I = engine(I, whichengine); // is now a principal ideal; if (qr == 1) { setring save2; } else { setring save; } ideal J; J[2] = var(1); map @m = Dh,J; ideal I = @m(I); poly p = I[1]; list l = listofroots(p,1); if (qr == 1) { setring save; list l = imap(save2,l); } return(l); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; poly f = x^2+y^3+x*y^2; bfctOneGB(f); bfctOneGB(f,1,1); } proc bfctAnn (poly f, list #) "USAGE: bfctAnn(f [,r]); f a poly, r an optional int RETURN: list of ideal and intvec PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the @* hypersurface defined by f. ASSUME: The basering is commutative and of characteristic 0. BACKGROUND: In this proc, ann(f^s) is computed and then a system of linear @* equations is solved by linear reductions. NOTE: In the output list, the ideal contains all the roots and the intvec @* their multiplicities. @* If r<>0, @code{std} is used for GB computations, @* otherwise, and by default, @code{slimgb} is used. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example bfctAnn; shows examples " { def save = basering; int ppl = printlevel - voice + 2; int whichengine = 0; // default if (size(#)>0) { if (typeof(#[1])=="int" || typeof(#[1])=="number") { whichengine = int(#[1]); } } list b = bfctengine(f,1,whichengine,0,0,0,0,0); return(b); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; poly f = x^2+y^3+x*y^2; bfctAnn(f); } /* //static proc hardexamples () { // some hard examples ring r1 = 0,(x,y,z,w),dp; // ab34 poly ab34 = (z3+w4)*(3z2x+4w3y); bfct(ab34); // ha3 poly ha3 = xyzw*(x+y)*(x+z)*(x+w)*(y+z+w); bfct(ha3); // ha4 poly ha4 = xyzw*(x+y)*(x+z)*(x+w)*(y+z)*(y+w); bfct(ha4); // chal4: reiffen(4,5)*reiffen(5,4) ring r2 = 0,(x,y),dp; poly chal4 = (x4+xy4+y5)*(x5+x4y+y4); bfct(chal4); // (xy+z)*reiffen(4,5) ring r3 = 0,(x,y,z),dp; poly xyzreiffen45 = (xy+z)*(y4+yz4+z5); bfct(xyzreiffen45); // sparse ideal as suggested by Alex; gives 1 as std ideal I1 = 28191*y^2+14628*x*Dy, 24865*x^2+24072*x*Dx+17756*Dy^2; std(I1); } */