////////////////////////////////////////////////////////////////////////////// version="$Id: dmod.lib,v 1.40 2009-04-14 12:00:14 Singular Exp $"; category="Noncommutative"; info=" LIBRARY: dmod.lib Algorithms for algebraic D-modules AUTHORS: Viktor Levandovskyy, levandov@math.rwth-aachen.de @* Jorge Martin Morales, jorge@unizar.es THEORY: Let K be a field of characteristic 0. Given a polynomial ring @* R = K[x_1,...,x_n] and a polynomial F in R, @* one is interested in the R[1/F]-module of rank one, generated by @* the symbol F^s for a symbolic discrete variable s. @* In fact, the module R[1/F]*F^s has a structure of a D(R)[s]-module, where D(R) @* is an n-th Weyl algebra K and @* D(R)[s] = D(R) tensored with K[s] over K. @* Constructively, one needs to find a left ideal I = I(F^s) in D(R), such @* that K[x_1,...,x_n,1/F]*F^s is isomorphic to D(R)/I as a D(R)-module. @* We often write just D for D(R) and D[s] for D(R)[s]. @* One is interested in the following data: @* - Ann F^s = I = I(F^s) in D(R)[s], denoted by LD in the output @* - global Bernstein polynomial in K[s], denoted by bs, @* - its minimal integer root s0, the list of all roots of bs, which are known @* to be rational, with their multiplicities, which is denoted by BS @* - Ann F^s0 = I(F^s0) in D(R), denoted by LD0 in the output @* (LD0 is a holonomic ideal in D(R)) @* - Ann^(1) F^s in D(R)[s], denoted by LD1 (logarithmic derivations) @* - an operator in D(R)[s], denoted by PS, such that the functional equality @* PS*F^(s+1) = bs*F^s holds in K[x_1,...,x_n,1/F]*F^s. REFERENCES: @* We provide the following implementations of algorithms: @* (OT) the classical Ann F^s algorithm from Oaku and Takayama (Journal of @* Pure and Applied Math., 1999), @* (LOT) Levandovskyy's modification of the Oaku-Takayama algorithm (ISSAC 2007) @* (BM) the Ann F^s algorithm by Briancon and Maisonobe (Remarques sur @* l'ideal de Bernstein associe a des polynomes, preprint, 2002) @* (LM08) V. Levandovskyy and J. Martin-Morales, ISSAC 2008 @* (C) Countinho, A Primer of Algebraic D-Modules, @* (SST) Saito, Sturmfels, Takayama 'Groebner Deformations of Hypergeometric @* Differential Equations', Springer, 2000 GUIDE: @* - Ann F^s = I(F^s) = LD in D(R)[s] can be computed by Sannfs [BM, OT, LOT] @* - Ann^(1) F^s in D(R)[s] can be computed by Sannfslog @* - global Bernstein polynomial bs in K[s] can be computed by bernsteinBM @* - Ann F^s0 = I(F^s0) = LD0 in D(R) can be computed by annfs0, annfs, annfsBM, @* annfsOT, annfsLOT, annfs2, annfsRB etc. @* - all the relevant data to F^s (LD, LD0, bs, PS) are computed by operatorBM @* @* - annihilator of F^{s1} for a number s1 is computed with annfspecial @* - annihilator of F_1^s_1 * ... * F_p^s_p is computed with annfsBMI @* - computing the multiplicity of a rational number r in the Bernstein poly @* of a given ideal goes with checkRoot @* - check, whether a given univariate polynomial divides the Bernstein poly @* goes with checkFactor MAIN PROCEDURES: annfs(F[,S,eng]); compute Ann F^s0 in D and Bernstein poly for a poly F annfspecial(I, F, m, n); compute Ann F^n from Ann F^s for a poly F and a number n Sannfs(F[,S,eng]); compute Ann F^s in D[s] for a poly F Sannfslog(F[,eng]); compute Ann^(1) F^s in D[s] for a poly F bernsteinBM(F[,eng]); compute global Bernstein poly for a poly F (algorithm of Briancon-Maisonobe) operatorBM(F[,eng]); compute Ann F^s, Ann F^s0, BS and PS for a poly F (algorithm of Briancon-Maisonobe) annfsParamBM(F[,eng]); compute the generic Ann F^s (algorithm by Briancon and Maisonobe) and exceptional parametric constellations for a poly F with parametric coefficients annfsBMI(F[,eng]); compute Ann F^s and Bernstein ideal for a poly F=f1*..*fP (multivariate algorithm of Briancon-Maisonobe) checkRoot(F,a[,S,eng]); check if a given rational is a root of the global Bernstein polynomial of F and compute its multiplicity annfsBM(F[,eng]); compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of Briancon-Maisonobe) annfsLOT(F[,eng]); compute Ann F^s0 in D and Bernstein poly for a poly F (Levandovskyy modification of the Oaku-Takayama algorithm) annfsOT(F[,eng]); compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of Oaku-Takayama) SannfsBM(F[,eng]); compute Ann F^s in D[s] for a poly F (algorithm of Briancon-Maisonobe) SannfsBFCT(F[,eng]); compute Ann F^s in D[s] for a poly F (algorithm of Briancon-Maisonobe, other output ordering) SannfsLOT(F[,eng]); compute Ann F^s in D[s] for a poly F (Levandovskyy modification of the Oaku-Takayama algorithm) SannfsOT(F[,eng]); compute Ann F^s in D[s] for a poly F (algorithm of Oaku-Takayama) annfs0(I,F [,eng]); compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s] 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 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 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] 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] 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] AUXILIARY PROCEDURES: arrange(p); create a poly, describing a full hyperplane arrangement reiffen(p,q); create a poly, describing a Reiffen curve isHolonomic(M); check whether a module is holonomic convloc(L); replace global orderings with local in the ringlist L minIntRoot(P,fact); minimal integer root among the roots of a maximal ideal P varNum(s); the number of the variable with the name s isRational(n); check whether n is a rational number SEE ALSO: gmssing_lib, bfun_lib, dmodapp_lib "; LIB "matrix.lib"; // for submat LIB "nctools.lib"; LIB "elim.lib"; // for nselect LIB "qhmoduli.lib"; // for Max LIB "gkdim.lib"; LIB "gmssing.lib"; LIB "control.lib"; // for genericity LIB "dmodapp.lib"; // for e.g. engine LIB "poly.lib"; proc testdmodlib() { /* tests all procs for consistency */ /* adding the new proc, add it here */ "MAIN PROCEDURES:"; example annfs; example Sannfs; example Sannfslog; example bernsteinBM; example operatorBM; example annfspecial; example annfsParamBM; example annfsBMI; example checkRoot; example annfs0; example annfs2; example annfsRB; example annfs2; "SECONDARY D-MOD PROCEDURES:"; example annfsBM; example annfsLOT; example annfsOT; example SannfsBM; example SannfsLOT; example SannfsOT; example SannfsBFCT; example checkRoot1; example checkRoot2; example checkFactor; } // new top-level procedures proc annfs(poly F, list #) "USAGE: annfs(f [,S,eng]); f a poly, S a string, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s with the algorithm @* given in S and with the Groebner basis engine given in ''eng'' NOTE: activate the output ring with the @code{setring} command. @* String S; S can be one of the following: @* 'bm' (default) - for the algorithm of Briancon and Maisonobe, @* 'ot' - for the algorithm of Oaku and Takayama, @* 'lot' - for the Levandovskyy's modification of the algorithm of OT. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise and by default @code{slimgb} is used. @* In the output ring: @* - the ideal LD (which is a Groebner basis) is the needed D-module structure, @* - the list BS contains roots and multiplicities of a BS-polynomial of f. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example annfs; shows examples " { int eng = 0; int chs = 0; // choice string algo = "bm"; string st; if ( size(#)>0 ) { if ( typeof(#[1]) == "string" ) { st = string(#[1]); if ( (st == "BM") || (st == "Bm") || (st == "bM") ||(st == "bm")) { algo = "bm"; chs = 1; } if ( (st == "OT") || (st == "Ot") || (st == "oT") || (st == "ot")) { algo = "ot"; chs = 1; } if ( (st == "LOT") || (st == "lOT") || (st == "Lot") || (st == "lot")) { algo = "lot"; chs = 1; } if (chs != 1) { // incorrect value of S print("Incorrect algorithm given, proceed with the default BM"); algo = "bm"; } // second arg if (size(#)>1) { // exists 2nd arg if ( typeof(#[2]) == "int" ) { // the case: given alg, given engine eng = int(#[2]); } else { eng = 0; } } else { // no second arg eng = 0; } } else { if ( typeof(#[1]) == "int" ) { // the case: default alg, engine eng = int(#[1]); // algo = "bm"; //is already set } else { // incorr. 1st arg algo = "bm"; } } } // size(#)=0, i.e. there is no algorithm and no engine given // eng = 0; algo = "bm"; //are already set // int ppl = printlevel-voice+2; int old_printlevel = printlevel; printlevel=printlevel+1; def save = basering; if ( algo=="ot") { def @A = annfsOT(F,eng); } else { if ( algo=="lot") { def @A = annfsLOT(F,eng); } else { // bm = default def @A = annfsBM(F,eng); } } printlevel = old_printlevel; return(@A); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = z*x^2+y^3; def A = annfs(F); // here, the default BM algorithm will be used setring A; // the Weyl algebra in (x,y,z,Dx,Dy,Dz) LD; //the annihilator of F^{-1} over A BS; // roots with multiplicities of BS polynomial } proc Sannfs(poly F, list #) "USAGE: Sannfs(f [,S,eng]); f a poly, S a string, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[f^s] with the algorithm @* given in S and with the Groebner basis engine given in eng NOTE: activate the output ring with the @code{setring} command. @* The value of a string S can be @* 'bm' (default) - for the algorithm of Briancon and Maisonobe, @* 'lot' - for the Levandovskyy's modification of the algorithm of OT, @* 'ot' - for the algorithm of Oaku and Takayama. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. @* In the output ring: @* - the ideal LD is the needed D-module structure. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example Sannfs; shows examples " { int eng = 0; int chs = 0; // choice string algo = "bm"; string st; if ( size(#)>0 ) { if ( typeof(#[1]) == "string" ) { st = string(#[1]); if ( (st == "BM") || (st == "Bm") || (st == "bM") ||(st == "bm")) { algo = "bm"; chs = 1; } if ( (st == "OT") || (st == "Ot") || (st == "oT") || (st == "ot")) { algo = "ot"; chs = 1; } if ( (st == "LOT") || (st == "lOT") || (st == "Lot") || (st == "lot")) { algo = "lot"; chs = 1; } if (chs != 1) { // incorrect value of S print("Incorrect algorithm given, proceed with the default BM"); algo = "bm"; } // second arg if (size(#)>1) { // exists 2nd arg if ( typeof(#[2]) == "int" ) { // the case: given alg, given engine eng = int(#[2]); } else { eng = 0; } } else { // no second arg eng = 0; } } else { if ( typeof(#[1]) == "int" ) { // the case: default alg, engine eng = int(#[1]); // algo = "bm"; //is already set } else { // incorr. 1st arg algo = "bm"; } } } // size(#)=0, i.e. there is no algorithm and no engine given // eng = 0; algo = "bm"; //are already set // int ppl = printlevel-voice+2; printlevel=printlevel+1; def save = basering; if ( algo=="ot") { def @A = SannfsOT(F,eng); } else { if ( algo=="lot") { def @A = SannfsLOT(F,eng); } else { // bm = default def @A = SannfsBM(F,eng); } } printlevel=printlevel-1; return(@A); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^3+y^3+z^3; printlevel = 0; def A = Sannfs(F); // here, the default BM algorithm will be used setring A; LD; } proc Sannfslog (poly F, list #) "USAGE: Sannfslog(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s NOTE: activate the output ring with the @code{setring} command. @* In the output ring D[s], the ideal LD1 is generated by the elements @* in Ann F^s in D[s], coming from logarithmic derivations. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example Sannfslog; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } int ppl = printlevel-voice+2; def save = basering; int N = nvars(basering); int Nnew = 2*N+1; int i; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether vars have admissible names list Name = RL[2]; for (i=1; i<=N; i++) { if (Name[i] == "s") { ERROR("Variable names should not include s"); } } // the ideal I ideal I = -F, jacob(F); dbprint(ppl,"// -1-1- starting the computation of syz(-F,_Dx(F))"); dbprint(ppl-1, I); matrix M = syz(I); M = transpose(M); // it is more usefull working with columns dbprint(ppl,"// -1-2- the module syz(-F,_Dx(F)) has been computed"); dbprint(ppl-1, M); // ------------ the ring @R ------------ // _x, _Dx, s; elim.ord for _x,_Dx. // now, create the names for new vars list DName; for (i=1; i<=N; i++) { DName[i] = "D"+Name[i]; // concat } tmp[1] = "s"; list NName; for (i=1; i<=N; i++) { NName[2*i-1] = Name[i]; NName[2*i] = DName[i]; //NName[2*i-1] = DName[i]; //NName[2*i] = Name[i]; } NName[Nnew] = tmp[1]; L[2] = NName; tmp = 0; // block ord (a(1,1),a(0,0,1,1),...,dp); //list("a",intvec(1,1)), list("a",intvec(0,0,1,1)), ... tmp[1] = "a"; // string for (i=1; i<=N; i++) { iv[2*i-1] = 1; iv[2*i] = 1; tmp[2] = iv; iv = 0; // intvec Lord[i] = tmp; } //list("dp",intvec(1,1,1,1,1,...)) s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); kill s; tmp[1] = "dp"; // string tmp[2] = iv; // intvec Lord[N+1] = tmp; //list("C",intvec(0)) tmp[1] = "C"; // string iv = 0; tmp[2] = iv; // intvec Lord[N+2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[2*i-1,2*i]=1; //@D[2*i-1,2*i]=-1; } def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -2-1- the ring @R(_x,_Dx,s) is ready"); dbprint(ppl-1, @R); matrix M = imap(save,M); // now, create the vector [-s,_Dx] vector v = [-s]; // now s is a variable for (i=1; i<=N; i++) { v = v + var(2*i)*gen(i+1); //v = v + var(2*i-1)*gen(i+1); } ideal J = ideal(M*v); // make leadcoeffs positive for (i=1; i<= ncols(J); i++) { if ( leadcoef(J[i])<0 ) { J[i] = -J[i]; } } ideal LD1 = J; kill J; export LD1; return(@R); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),Dp; poly F = x^4+y^5+x*y^4; printlevel = 0; def A = Sannfslog(F); setring A; LD1; } // alternative code for SannfsBM, renamed from annfsBM to ALTannfsBM // is superfluos - will not be included in the official documentation static proc ALTannfsBM (poly F, list #) "USAGE: ALTannfsBM(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the annihilator ideal of f^s in D[s], where D is the Weyl @* algebra, according to the algorithm by Briancon and Maisonobe NOTE: activate the output ring with the @code{setring} command. In this ring, @* - the ideal LD is the annihilator of f^s. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example ALTannfsBM; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*N+2; int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; //char L[4] = RL[4]; //char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "t"; RName[2] = "s"; for (i=1; i<=N; i++) { for(j=1; j<=size(RName); j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include t,s"); } } } // now, create the names for new vars list DName; for (i=1; i<=N; i++) { DName[i] = "D"+Name[i]; //concat } tmp[1] = "t"; tmp[2] = "s"; list NName = tmp + Name + DName; L[2] = NName; // Name, Dname will be used further kill NName; // block ord (lp(2),dp); tmp[1] = "lp"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill s; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2]=t; for(i=1; i<=N; i++) { @D[2+i,N+2+i]=1; } // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = t*F+s; poly p; for(i=1; i<=N; i++) { p = t; //t p = diff(F,var(2+i))*p; I = I, var(N+2+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of t in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1); kill I,J; dbprint(ppl,"// -1-3- t is eliminated"); dbprint(ppl-1, K); //K is without t // create Dn[s], where Dn is the ordinary Weyl Algebra, and put the result into it, // thus creating the ring @R2 // keep: N, i,j,s, tmp, RL setring save; Nnew = 2*N+1; // list RL = ringlist(save); //is defined earlier kill Lord, tmp, iv; L = 0; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2] // DName is defined earlier tmp[1] = "s"; list NName = Name + DName + tmp; L[2] = NName; // dp ordering; string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)] = ";"; execute(s); kill s; tmp = 0; tmp[1] = "dp"; //string tmp[2] = iv; //intvec Lord[1] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list // Add: Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal K = imap(@R,K); option(redSB); //dbprint(ppl,"// -2-2- the final cosmetic std"); //K = engine(K,eng); //std does the job too // total cleanup kill @R; ideal LD = K; export LD; return(@R2); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z,w),Dp; poly F = x^3+y^3+z^2*w; printlevel = 0; def A = ALTannfsBM(F); setring A; LD; } proc bernsteinBM(poly F, list #) "USAGE: bernsteinBM(f [,eng]); f a poly, eng an optional int RETURN: list (of roots of the Bernstein polynomial b and their multiplicies) PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, @* defined by f, according to the algorithm by Briancon and Maisonobe NOTE: If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example bernsteinBM; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*N+2; int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; //char L[4] = RL[4]; //char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "t"; RName[2] = "s"; for (i=1; i<=N; i++) { for(j=1; j<=size(RName); j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include t,s"); } } } // now, create the names for new vars list DName; for (i=1; i<=N; i++) { DName[i] = "D"+Name[i]; //concat } tmp[1] = "t"; tmp[2] = "s"; list NName = tmp + Name + DName; L[2] = NName; // Name, Dname will be used further kill NName; // block ord (lp(2),dp); tmp[1] = "lp"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill s; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2]=t; for(i=1; i<=N; i++) { @D[2+i,N+2+i]=1; } // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = t*F+s; poly p; for(i=1; i<=N; i++) { p = t; //t p = diff(F,var(2+i))*p; I = I, var(N+2+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of t in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1); kill I,J; dbprint(ppl,"// -1-3- t is eliminated"); dbprint(ppl-1, K); //K is without t // ----------- the ring @R2 ------------ // _x, _Dx,s; elim.ord for _x,_Dx. // keep: N, i,j,s, tmp, RL setring save; Nnew = 2*N+1; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var tmp[1] = "s"; // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (dp(N),dp); string s = "iv="; for (i=1; i<=Nnew-1; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); tmp[1] = "dp"; //string tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; //string s[size(s)] = ","; s = s+"1;"; execute(s); kill s; kill NName; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); MM = 0,s,MM; map R01 = @R, MM; ideal K = R01(K); kill @R, R01; poly F = imap(save,F); K = K,F; dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2"); dbprint(ppl-1, K); ideal M = engine(K,eng); ideal K2 = nselect(M,1..Nnew-1); kill K,M; dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2"); dbprint(ppl-1, K2); // the ring @R3 and the search for minimal negative int s ring @R3 = 0,s,dp; dbprint(ppl,"// -3-1- the ring @R3(s) is ready"); ideal K3 = imap(@R2,K2); kill @R2; poly p = K3[1]; dbprint(ppl,"// -3-2- factorization"); list P = factorize(p); //with constants and multiplicities ideal bs; intvec m; //the Bernstein polynomial 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]; } // "--------- b-function factorizes into ---------"; P; //int sP = minIntRoot(bs,1); //dbprint(ppl,"// -3-3- minimal integer root found"); //dbprint(ppl-1, sP); // convert factors to a list of their roots and multiplicities bs = normalize(bs); bs = -subst(bs,s,0); setring save; ideal bs = imap(@R3,bs); kill @R3; list BS = bs,m; return(BS); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z,w),Dp; poly F = x^3+y^3+z^2*w; printlevel = 0; bernsteinBM(F); } // some changes proc annfsBM (poly F, list #) "USAGE: annfsBM(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s, according @* to the algorithm by Briancon and Maisonobe NOTE: activate the output ring with the @code{setring} command. In this ring, @* - the ideal LD (which is a Groebner basis) is the needed D-module structure, @* which is obtained by substituting the minimal integer root of a Bernstein @* polynomial into the s-parametric ideal; @* - the list BS is the list of roots and multiplicities of a Bernstein polynomial of f. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example annfsBM; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*N+2; int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; //char L[4] = RL[4]; //char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "t"; RName[2] = "s"; for (i=1; i<=N; i++) { for(j=1; j<=size(RName); j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include t,s"); } } } // now, create the names for new vars list DName; for (i=1; i<=N; i++) { DName[i] = "D"+Name[i]; //concat } tmp[1] = "t"; tmp[2] = "s"; list NName = tmp + Name + DName; L[2] = NName; // Name, Dname will be used further kill NName; // block ord (lp(2),dp); tmp[1] = "lp"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill s; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2]=t; for(i=1; i<=N; i++) { @D[2+i,N+2+i]=1; } // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = t*F+s; poly p; for(i=1; i<=N; i++) { p = t; //t p = diff(F,var(2+i))*p; I = I, var(N+2+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of t in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1); kill I,J; dbprint(ppl,"// -1-3- t is eliminated"); dbprint(ppl-1, K); //K is without t setring save; // ----------- the ring @R2 ------------ // _x, _Dx,s; elim.ord for _x,_Dx. // keep: N, i,j,s, tmp, RL Nnew = 2*N+1; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var tmp[1] = "s"; // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (dp(N),dp); string s = "iv="; for (i=1; i<=Nnew-1; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); tmp[1] = "dp"; //string tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; //string s[size(s)] = ","; s = s+"1;"; execute(s); kill s; kill NName; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); MM = 0,s,MM; map R01 = @R, MM; ideal K = R01(K); poly F = imap(save,F); K = K,F; dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2"); dbprint(ppl-1, K); ideal M = engine(K,eng); ideal K2 = nselect(M,1..Nnew-1); kill K,M; dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2"); dbprint(ppl-1, K2); // the ring @R3 and the search for minimal negative int s ring @R3 = 0,s,dp; dbprint(ppl,"// -3-1- the ring @R3(s) is ready"); ideal K3 = imap(@R2,K2); poly p = K3[1]; dbprint(ppl,"// -3-2- factorization"); list P = factorize(p); //with constants and multiplicities ideal bs; intvec m; //the Bernstein polynomial is monic, so we are not interested in constants for (i=2; i<= size(P[1]); i++) //we ignore P[1][1] (constant) and P[2][1] (its mult.) { bs[i-1] = P[1][i]; m[i-1] = P[2][i]; } // "--------- b-function factorizes into ---------"; P; int sP = minIntRoot(bs,1); dbprint(ppl,"// -3-3- minimal integer root found"); dbprint(ppl-1, sP); // convert factors to a list of their roots bs = normalize(bs); bs = -subst(bs,s,0); list BS = bs,m; //TODO: sort BS! // --------- substitute s found in the ideal --------- // --------- going back to @R and substitute --------- setring @R; ideal K2 = subst(K,s,sP); kill K; // create the ordinary Weyl algebra and put the result into it, // thus creating the ring @R4 // keep: N, i,j,s, tmp, RL setring save; Nnew = 2*N; // list RL = ringlist(save); //is defined earlier kill Lord, tmp, iv; L = 0; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2]M // DName is defined earlier list NName = Name + DName; L[2] = NName; // dp ordering; string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)] = ";"; execute(s); kill s; tmp = 0; tmp[1] = "dp"; //string tmp[2] = iv; //intvec Lord[1] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list // Add: Plural part def @R4@ = ring(L); setring @R4@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R4 = nc_algebra(1,@D); setring @R4; kill @R4@; dbprint(ppl,"// -4-1- the ring @R4(_x,_Dx) is ready"); dbprint(ppl-1, @R4); ideal K4 = imap(@R,K2); option(redSB); dbprint(ppl,"// -4-2- the final cosmetic std"); K4 = engine(K4,eng); //std does the job too // total cleanup kill @R; kill @R2; list BS = imap(@R3,BS); export BS; kill @R3; ideal LD = K4; export LD; return(@R4); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = z*x^2+y^3; printlevel = 0; def A = annfsBM(F); setring A; LD; BS; } // replacing s with -s-1 => data is shorter // analogue of annfs0 proc annfs2(ideal I, poly F, list #) "USAGE: annfs2(I, F [,eng]); I an ideal, F a poly, eng an optional int RETURN: ring PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, @* based on the output of Sannfs-like procedure @* annfs2 uses shorter expressions in the variable s (the idea of Noro). NOTE: activate the output ring with the @code{setring} command. In this ring, @* - the ideal LD (which is a Groebner basis) is the annihilator of f^s, @* - the list BS contains the roots with multiplicities of the BS polynomial. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise and by default @code{slimgb} is used. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example annfs2; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } def @R2 = basering; // we're in D_n[s], where the elim ord for s is set ideal J = NF(I,std(F)); // make leadcoeffs positive int i; J = subst(J,s,-s-1); for (i=1; i<= ncols(J); i++) { if (leadcoef(J[i]) <0 ) { J[i] = -J[i]; } } J = J,F; ideal M = engine(J,eng); int Nnew = nvars(@R2); ideal K2 = nselect(M,1..Nnew-1); int ppl = printlevel-voice+2; dbprint(ppl,"// -1-1- _x,_Dx are eliminated in basering"); dbprint(ppl-1, K2); // the ring @R3 and the search for minimal negative int s ring @R3 = 0,s,dp; dbprint(ppl,"// -2-1- the ring @R3 i.e. K[s] is ready"); ideal K3 = imap(@R2,K2); poly p = K3[1]; dbprint(ppl,"// -2-2- factorization"); // ideal P = factorize(p,1); //without constants and multiplicities // "--------- b-function factorizes into ---------"; P; // convert factors to the list of their roots with mults // assume all factors are linear // ideal BS = normalize(P); // BS = subst(BS,s,0); // BS = -BS; list P = factorize(p); //with constants and multiplicities ideal bs; intvec m; //the Bernstein polynomial 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]; bs[i-1] = subst(bs[i-1],s,-s-1); m[i-1] = P[2][i]; } int sP = minIntRoot(bs,1); bs = normalize(bs); bs = -subst(bs,s,0); dbprint(ppl,"// -2-3- minimal integer root found"); dbprint(ppl-1, sP); //TODO: sort BS! // --------- substitute s found in the ideal --------- // --------- going back to @R and substitute --------- setring @R2; K2 = subst(I,s,sP); // create the ordinary Weyl algebra and put the result into it, // thus creating the ring @R5 // keep: N, i,j,s, tmp, RL Nnew = Nnew - 1; // former 2*N; // list RL = ringlist(save); // is defined earlier // kill Lord, tmp, iv; list L = 0; list Lord, tmp; intvec iv; list RL = ringlist(basering); L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2]M // DName is defined earlier list NName; // = RL[2]; // skip the last var 's' for (i=1; i<=Nnew; i++) { NName[i] = RL[2][i]; } L[2] = NName; // dp ordering; string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)] = ";"; execute(s); tmp = 0; tmp[1] = "dp"; // string tmp[2] = iv; // intvec Lord[1] = tmp; kill s; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list // Add: Plural part def @R4@ = ring(L); setring @R4@; int N = Nnew/2; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R4 = nc_algebra(1,@D); setring @R4; kill @R4@; dbprint(ppl,"// -3-1- the ring @R4 is ready"); dbprint(ppl-1, @R4); ideal K4 = imap(@R2,K2); option(redSB); dbprint(ppl,"// -3-2- the final cosmetic std"); K4 = engine(K4,eng); // std does the job too // total cleanup ideal bs = imap(@R3,bs); kill @R3; list BS = bs,m; export BS; ideal LD = K4; export LD; return(@R4); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^3+y^3+z^3; printlevel = 0; def A = SannfsBM(F); setring A; LD; poly F = imap(r,F); def B = annfs2(LD,F); setring B; LD; BS; } // try to replace s with -s-1 => data is shorter as in annfs2 // and use Macaulay's reduceB strategy, that is add // not F but ; the resulting B-function // has to be multiplied with (s+1) at the very end proc annfsRB(ideal I, poly F, list #) "USAGE: annfsRB(I, F [,eng]); I an ideal, F a poly, eng an optional int RETURN: ring PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, @* based on the output of Sannfs like procedure NOTE: activate the output ring with the @code{setring} command. In this ring, @* - the ideal LD (which is a Groebner basis) is the annihilator of f^s, @* - the list BS contains the roots with multiplicities of a Bernstein polynomial of f. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise and by default @code{slimgb} is used. @* This procedure follows the 'reduceB' strategy, used in Macaulay2. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example annfsRB; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } def @R2 = basering; int ppl = printlevel-voice+2; // we're in D_n[s], where the elim ord for s is set // switch to comm. ring in X's and compute the GB of Tjurina ideal dbprint(ppl,"// -1-0- creating K[x] and Tjurina ideal"); list nL = ringlist(@R2); list temp,t2; temp[1] = nL[1]; temp[4] = nL[4]; int @n = int((nvars(@R2)-1)/2); // # of x's int i; for (i=1; i<=@n; i++) { t2[i] = nL[2][i]; } temp[2] = t2; t2 = 0; t2[1] = nL[3][1]; // more weights than vars? t2[2] = nL[3][3]; temp[3] = t2; def @R22 = ring(temp); setring @R22; poly F = imap(@R2,F); ideal J = F; for (i=1; i<=@n; i++) { J = J, diff(F,var(i)); } J = engine(J,eng); dbprint(ppl,"// -1-1- finished computing the GB of Tjurina ideal"); dbprint(ppl-1, J); setring @R2; ideal JF = imap(@R22,J); kill @R22; attrib(JF,"isSB",1); // embedded comm ring is used ideal J = NF(I,JF); dbprint(ppl,"// -1-2- finished computing the NF of I wrt Tjurina ideal"); dbprint(ppl-1, J2); // make leadcoeffs positive J = subst(J,s,-s-1); for (i=1; i<= ncols(J); i++) { if (leadcoef(J[i]) <0 ) { J[i] = -J[i]; } } J = J,JF; ideal M = engine(J,eng); int Nnew = nvars(@R2); ideal K2 = nselect(M,1..Nnew-1); dbprint(ppl,"// -2-0- _x,_Dx are eliminated in basering"); dbprint(ppl-1, K2); // the ring @R3 and the search for minimal negative int s ring @R3 = 0,s,dp; dbprint(ppl,"// -2-1- the ring @R3 i.e. K[s] is ready"); ideal K3 = imap(@R2,K2); poly p = K3[1]; p = s*p; // mult with the lost (s+1) factor dbprint(ppl,"// -2-2- factorization"); // ideal P = factorize(p,1); //without constants and multiplicities // "--------- b-function factorizes into ---------"; P; // convert factors to the list of their roots with mults // assume all factors are linear // ideal BS = normalize(P); // BS = subst(BS,s,0); // BS = -BS; list P = factorize(p); //with constants and multiplicities ideal bs; intvec m; //the Bernstein polynomial 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]; bs[i-1] = subst(bs[i-1],s,-s-1); m[i-1] = P[2][i]; } int sP = minIntRoot(bs,1); bs = normalize(bs); bs = -subst(bs,s,0); dbprint(ppl,"// -2-3- minimal integer root found"); dbprint(ppl-1, sP); //TODO: sort BS! // --------- substitute s found in the ideal --------- // --------- going back to @R and substitute --------- setring @R2; K2 = subst(I,s,sP); // create the ordinary Weyl algebra and put the result into it, // thus creating the ring @R5 // keep: N, i,j,s, tmp, RL Nnew = Nnew - 1; // former 2*N; // list RL = ringlist(save); // is defined earlier // kill Lord, tmp, iv; list L = 0; list Lord, tmp; intvec iv; list RL = ringlist(basering); L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2]M // DName is defined earlier list NName; // = RL[2]; // skip the last var 's' for (i=1; i<=Nnew; i++) { NName[i] = RL[2][i]; } L[2] = NName; // dp ordering; string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)] = ";"; execute(s); tmp = 0; tmp[1] = "dp"; // string tmp[2] = iv; // intvec Lord[1] = tmp; kill s; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list // Add: Plural part def @R4@ = ring(L); setring @R4@; int N = Nnew/2; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R4 = nc_algebra(1,@D); setring @R4; kill @R4@; dbprint(ppl,"// -3-1- the ring @R4 is ready"); dbprint(ppl-1, @R4); ideal K4 = imap(@R2,K2); option(redSB); dbprint(ppl,"// -3-2- the final cosmetic std"); K4 = engine(K4,eng); // std does the job too // total cleanup ideal bs = imap(@R3,bs); kill @R3; list BS = bs,m; export BS; ideal LD = K4; export LD; return(@R4); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^3+y^3+z^3; printlevel = 0; def A = SannfsBM(F); setring A; LD; // s-parametric ahhinilator poly F = imap(r,F); def B = annfsRB(LD,F); setring B; LD; BS; } proc operatorBM(poly F, list #) "USAGE: operatorBM(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the B-operator and other relevant data for Ann F^s, @* using e.g. algorithm by Briancon and Maisonobe for Ann F^s and BS. NOTE: activate the output ring with the @code{setring} command. In the output ring D[s] @* - the polynomial F is the same as the input, @* - the ideal LD is the annihilator of f^s in Dn[s], @* - the ideal LD0 is the needed D-mod structure, where LD0 = LD|s=s0, @* - the polynomial bs is the global Bernstein polynomial of f in the variable s, @* - the list BS contains all the roots with multiplicities of the global Bernstein polynomial of f, @* - the polynomial PS is an operator in Dn[s] such that PS*f^(s+1) = bs*f^s. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise and by default @code{slimgb} is used. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example operatorBM; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*N+2; int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; //char L[4] = RL[4]; //char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "t"; RName[2] = "s"; for (i=1; i<=N; i++) { for(j=1; j<=size(RName); j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include t,s"); } } } // now, create the names for new vars list DName; for (i=1; i<=N; i++) { DName[i] = "D"+Name[i]; //concat } tmp[1] = "t"; tmp[2] = "s"; list NName = tmp + Name + DName; L[2] = NName; // Name, Dname will be used further kill NName; // block ord (lp(2),dp); tmp[1] = "lp"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill s; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2]=t; for(i=1; i<=N; i++) { @D[2+i,N+2+i]=1; } // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = t*F+s; poly p; for(i=1; i<=N; i++) { p = t; //t p = diff(F,var(2+i))*p; I = I, var(N+2+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of t in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1); kill I,J; dbprint(ppl,"// -1-3- t is eliminated"); dbprint(ppl-1, K); //K is without t setring save; // ----------- the ring @R2 ------------ // _x, _Dx,s; elim.ord for _x,_Dx. // keep: N, i,j,s, tmp, RL Nnew = 2*N+1; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var tmp[1] = "s"; // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (dp(N),dp); string s = "iv="; for (i=1; i<=Nnew-1; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); tmp[1] = "dp"; //string tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; //string s[size(s)] = ","; s = s+"1;"; execute(s); kill s; kill NName; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); MM = 0,s,MM; map R01 = @R, MM; ideal K = R01(K); poly F = imap(save,F); K = K,F; dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2"); dbprint(ppl-1, K); ideal M = engine(K,eng); ideal K2 = nselect(M,1..Nnew-1); kill K,M; dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2"); dbprint(ppl-1, K2); // the ring @R3 and the search for minimal negative int s ring @R3 = 0,s,dp; dbprint(ppl,"// -3-1- the ring @R3(s) is ready"); ideal K3 = imap(@R2,K2); kill @R2; poly p = K3[1]; dbprint(ppl,"// -3-2- factorization"); list P = factorize(p); //with constants and multiplicities ideal bs; intvec m; //the Bernstein polynomial 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]; } // "--------- b-function factorizes into ---------"; P; int sP = minIntRoot(bs,1); dbprint(ppl,"// -3-3- minimal integer root found"); dbprint(ppl-1, sP); // convert factors to a list of their roots with multiplicities bs = normalize(bs); bs = -subst(bs,s,0); list BS = bs,m; //TODO: sort BS! // --------- substitute s found in the ideal --------- // --------- going back to @R and substitute --------- setring @R; ideal K2 = subst(K,s,sP); // create Dn[s], where Dn is the ordinary Weyl algebra, and put the result into it, // thus creating the ring @R4 // keep: N, i,j,s, tmp, RL setring save; Nnew = 2*N+1; // list RL = ringlist(save); //is defined earlier kill Lord, tmp, iv; L = 0; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2] // DName is defined earlier tmp[1] = "s"; list NName = Name + DName + tmp; L[2] = NName; // dp ordering; string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)] = ";"; execute(s); kill s; tmp = 0; tmp[1] = "dp"; //string tmp[2] = iv; //intvec Lord[1] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list // Add: Plural part def @R4@ = ring(L); setring @R4@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R4 = nc_algebra(1,@D); setring @R4; kill @R4@; dbprint(ppl,"// -4-1- the ring @R4(_x,_Dx,s) is ready"); dbprint(ppl-1, @R4); ideal LD0 = imap(@R,K2); ideal LD = imap(@R,K); kill @R; poly bs = imap(@R3,p); list BS = imap(@R3,BS); kill @R3; bs = normalize(bs); poly F = imap(save,F); dbprint(ppl,"// -4-2- starting the computation of PS via lift"); //better liftstd, I didn't knot it works also for Plural, liftslimgb? // liftstd may give extra coeffs in the resulting ideal matrix T = lift(F+LD,bs); poly PS = T[1,1]; dbprint(ppl,"// -4-3- an operator PS found, PS*f^(s+1) = b(s)*f^s"); dbprint(ppl-1,PS); option(redSB); dbprint(ppl,"// -4-4- the final cosmetic std"); LD0 = engine(LD0,eng); //std does the job too LD = engine(LD,eng); export F,LD,LD0,bs,BS,PS; return(@R4); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^3+y^3+z^3; printlevel = 0; def A = operatorBM(F); setring A; F; // the original polynomial itself LD; // generic annihilator LD0; // annihilator bs; // normalized Bernstein poly BS; // roots and multiplicities of the Bernstein poly PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD reduce(PS*F-bs,LD); // check the property of PS } // more interesting: // ring r = 0,(x,y,z,w),Dp; // poly F = x^3+y^3+z^2*w; // TODO: 1 has to appear in the 1nd column of transp. matrix // this does not happen automatically // for this, do special modulo with emphasis on the 1comp (c,<) // of the transp matrix, there must appear 1 in the GB // then use lift to get the expression... // need: (c,<) ordering for such comp's /* proc operatorModulo(poly F, ideal I, poly b) "USAGE: operatorModulo(f,I,b); f a poly, I an ideal, b a poly RETURN: poly PURPOSE: compute the B-operator from the poly f, ideal I = Ann f^s and Bernstein-Sato polynomial b using modulo i.e. kernel of module homomorphism NOTE: The computations take place in the ring, similar to the one returned by Sannfs procedure. @* If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example operatorModulo; shows examples " { // with hom_kernel; matrix AA[1][2] = F,-b; // matrix M = matrix(subst(I,s,s+1)*freemodule(2)); // ann f^{s+1} matrix N = matrix(I); // ann f^s // matrix K = hom_kernel(AA,M,N); matrix K = modulo(AA,N); K = transpose(K); print(K); ideal J; int i; poly t; number n; for(i=1; i<=nrows(K); i++) { if (K[i,2]!=0) { if ( leadmonom(K[i,2]) == 1) { t = K[i,1]; n = leadcoef(K[i,2]); t = t/n; // J = J, K[i][2]; break; } } } ideal J = groebner(subst(I,s,s+1)); // for NF t = NF(t,J); "candidate:"; t; J = subst(J,s,s-1); // test: if ( NF(t*F-b,J) !=0) { "Problem: PS does not work on F"; } return(t); } example { "EXAMPLE:"; echo = 2; // ring r = 0,(x,y,z),Dp; // poly F = x^3+y^3+z^3; LIB "dmod.lib"; option(prot); option(mem); ring r = 0,(x,y),Dp; // poly F = x^3+y^3+x*y^2; poly F = x^4 + y^4 + x*y^4; def A = Sannfs(F); // here we get LD = ann f^s setring A; poly F = imap(r,F); def B = annfs0(LD,F); // to obtain BS polynomial list BS = imap(B,BS); poly b = fl2poly(BS,"s"); LD = groebner(LD); poly PS = operatorModulo(F,LD,b); PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD reduce(PS*F-bs,LD); // check the property of PS } */ proc annfsParamBM (poly F, list #) "USAGE: annfsParamBM(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the generic Ann F^s and exceptional parametric constellations @* of a polynomial with parametric coefficients with the BM algorithm NOTE: activate the output ring with the @code{setring} command. In this ring, @* - the ideal LD is the D-module structure oa Ann F^s @* - the ideal Param contains special parameters as entries @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example annfsParamBM; shows examples " { //PURPOSE: compute the list of all possible Bernstein-Sato polynomials for a polynomial with parametric coefficients, according to the algorithm by Briancon and Maisonobe // @* - the list BS is the list of roots and multiplicities of a Bernstein polynomial of f. // ***** not implented yet **** int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*N+2; int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; //char L[4] = RL[4]; //char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "t"; RName[2] = "s"; for (i=1; i<=N; i++) { for(j=1; j<=size(RName); j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include t,s"); } } } // now, create the names for new vars list DName; for (i=1; i<=N; i++) { DName[i] = "D"+Name[i]; //concat } tmp[1] = "t"; tmp[2] = "s"; list NName = tmp + Name + DName; L[2] = NName; // Name, Dname will be used further kill NName; // block ord (lp(2),dp); tmp[1] = "lp"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill s; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2]=t; for(i=1; i<=N; i++) { @D[2+i,N+2+i]=1; } // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = t*F+s; poly p; for(i=1; i<=N; i++) { p = t; //t p = diff(F,var(2+i))*p; I = I, var(N+2+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of t in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1); dbprint(ppl,"// -1-3- t is eliminated"); dbprint(ppl-1, K); //K is without t // ----- looking for special parameters ----- dbprint(ppl,"// -2-1- starting the computation of the transformation matrix (via lift)"); J = normalize(J); matrix T = lift(I,J); //try also with liftstd kill I,J; dbprint(ppl,"// -2-2- the transformation matrix has been computed"); dbprint(ppl-1, T); //T is the transformation matrix dbprint(ppl,"// -2-3- genericity does the job"); list lParam = genericity(T); int ip = size(lParam); int cip; string sParam; if (sParam[1]=="-") { sParam=""; } //genericity returns "-" // if no parameters exist in a basering for (cip=1; cip <= ip; cip++) { sParam = sParam + "," +lParam[cip]; } if (size(sParam) >=2) { sParam = sParam[2..size(sParam)]; // removes the 1st colon } export sParam; kill T; dbprint(ppl,"// -2-4- the special parameters has been computed"); dbprint(ppl, sParam); // create Dn[s], where Dn is the ordinary Weyl Algebra, and put the result into it, // thus creating the ring @R2 // keep: N, i,j,s, tmp, RL setring save; Nnew = 2*N+1; // list RL = ringlist(save); //is defined earlier kill Lord, tmp, iv; L = 0; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2]M // DName is defined earlier tmp[1] = "s"; list NName = Name + DName + tmp; L[2] = NName; // dp ordering; string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)] = ";"; execute(s); kill s; tmp = 0; tmp[1] = "dp"; //string tmp[2] = iv; //intvec Lord[1] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list // Add: Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -3-1- the ring @R2(_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal K = imap(@R,K); kill @R; option(redSB); dbprint(ppl,"// -3-2- the final cosmetic std"); K = engine(K,eng); //std does the job too ideal LD = K; export LD; if (sParam[1] == ",") { sParam = sParam[2..size(sParam)]; } // || ((sParam[1] == " ") && (sParam[2] == ","))) execute("ideal Param ="+sParam+";"); export Param; kill sParam; return(@R2); } example { "EXAMPLE:"; echo = 2; ring r = (0,a,b),(x,y),Dp; poly F = x^2 - (y-a)*(y-b); printlevel = 0; def A = annfsParamBM(F); setring A; LD; Param; setring r; poly G = x2-(y-a)^2; // try the exceptional value b=a of parameters def B = annfsParamBM(G); setring B; LD; Param; } // *** the following example is nice, but too complicated for the documentation *** // ring r = (0,a),(x,y,z),Dp; // poly F = x^4+y^4+z^2+a*x*y*z; // printlevel = 2; //0 // def A = annfsParamBM(F); // setring A; // LD; // Param; proc annfsBMI(ideal F, list #) "USAGE: annfsBMI(F [,eng]); F an ideal, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s where @* f = F[1]*..*F[P], according to the algorithm by Briancon and Maisonobe. NOTE: activate the output ring with the @code{setring} command. In this ring, @* - the ideal LD is the needed D-mod structure, @* - the list BS is the Bernstein ideal of a polynomial f = F[1]*..*F[P]. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example annfsBMI; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int P = size(F); //if F has some generators which are zero, int P = ncols(I); int Nnew = 2*N+2*P; int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; //char L[4] = RL[4]; //char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; for (j=1; j<=P; j++) { RName[j] = "t("+string(j)+")"; RName[j+P] = "s("+string(j)+")"; } for(i=1; i<=N; i++) { for(j=1; j<=size(RName); j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include t(i),s(i)"); } } } // now, create the names for new vars list DName; for(i=1; i<=N; i++) { DName[i] = "D"+Name[i]; //concat } list NName = RName + Name + DName; L[2] = NName; // Name, Dname will be used further kill NName; // block ord (lp(P),dp); tmp[1] = "lp"; //string s = "iv="; for (i=1; i<=2*P; i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; //string s = "iv="; for (i=1; i<=Nnew; i++) //actually i<=2*N { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill s; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; for (i=1; i<=P; i++) { @D[i,i+P] = t(i); } for(i=1; i<=N; i++) { @D[2*P+i,2*P+N+i] = 1; } // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(_t,_s,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I ideal F = imap(save,F); ideal I = t(1)*F[1]+s(1); for (j=2; j<=P; j++) { I = I, t(j)*F[j]+s(j); } poly p,q; for (i=1; i<=N; i++) { p=0; for (j=1; j<=P; j++) { q = t(j); q = diff(F[j],var(2*P+i))*q; p = p + q; } I = I, var(2*P+N+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of "+string(t(1..P))+" in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1..P); kill I,J; dbprint(ppl,"// -1-3- all t(i) are eliminated"); dbprint(ppl-1, K); //K is without t(i) // ----------- the ring @R2 ------------ // _x, _Dx,s; elim.ord for _x,_Dx. // keep: N, i,j,s, tmp, RL setring save; Nnew = 2*N+P; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; //char L[4] = RL[4]; //char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var for (j=1; j<=P; j++) { tmp[j] = "s("+string(j)+")"; } // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (dp(N),dp); string s = "iv="; for (i=1; i<=Nnew-P; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); tmp[1] = "dp"; //string tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; //string s[size(s)] = ","; for (j=1; j<=P; j++) { s = s+"1,"; } s[size(s)]=";"; execute(s); kill s; kill NName; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,_s) is ready"); dbprint(ppl-1, @R2); // ideal MM = maxideal(1); // MM = 0,s,MM; // map R01 = @R, MM; // ideal K = R01(K); ideal F = imap(save,F); // maybe ideal F = R01(I); ? ideal K = imap(@R,K); // maybe ideal K = R01(I); ? poly f=1; for (j=1; j<=P; j++) { f = f*F[j]; } K = K,f; // to compute B (Bernstein-Sato ideal) //j=2; // for example //K = K,F[j]; // to compute Bj (see "On the computation of Bernstein-Sato ideals"; Castro, Ucha) //K = K,F; // to compute Bsigma (see "On the computation of Bernstein-Sato ideals"; Castro, Ucha) dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2"); dbprint(ppl-1, K); ideal M = engine(K,eng); ideal K2 = nselect(M,1..Nnew-P); kill K,M; dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2"); dbprint(ppl-1, K2); // the ring @R3 and factorize ring @R3 = 0,s(1..P),dp; dbprint(ppl,"// -3-1- the ring @R3(_s) is ready"); ideal K3 = imap(@R2,K2); if (size(K3)==1) { poly p = K3[1]; dbprint(ppl,"// -3-2- factorization"); // Warning: now P is an integer list Q = factorize(p); //with constants and multiplicities ideal bs; intvec m; for (i=2; i<=size(Q[1]); i++) //we delete Q[1][1] and Q[2][1] { bs[i-1] = Q[1][i]; m[i-1] = Q[2][i]; } // "--------- Q-ideal factorizes into ---------"; list(bs,m); list BS = bs,m; } else { // conjecture: the Bernstein ideal is principal dbprint(ppl,"// -3-2- the Bernstein ideal is not principal"); ideal BS = K3; } // create the ring @R4(_x,_Dx,_s) and put the result into it, // _x, _Dx,s; ord "dp". // keep: N, i,j,s, tmp, RL setring save; Nnew = 2*N+P; // list RL = ringlist(save); //is defined earlier kill Lord, tmp, iv; L = 0; list Lord, tmp; intvec iv; L[1] = RL[1]; //char L[4] = RL[4]; //char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var for (j=1; j<=P; j++) { tmp[j] = "s("+string(j)+")"; } // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // dp ordering; string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); kill s; kill NName; tmp[1] = "dp"; //string tmp[2] = iv; //intvec Lord[1] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R4@ = ring(L); setring @R4@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R4 = nc_algebra(1,@D); setring @R4; kill @R4@; dbprint(ppl,"// -4-1- the ring @R4i(_x,_Dx,_s) is ready"); dbprint(ppl-1, @R4); ideal K4 = imap(@R,K); option(redSB); dbprint(ppl,"// -4-2- the final cosmetic std"); K4 = engine(K4,eng); //std does the job too // total cleanup kill @R; kill @R2; def BS = imap(@R3,BS); export BS; kill @R3; ideal LD = K4; export LD; return(@R4); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),Dp; ideal F = x,y,x+y; printlevel = 0; def A = annfsBMI(F); setring A; LD; BS; } proc annfsOT(poly F, list #) "USAGE: annfsOT(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s, @* according to the algorithm by Oaku and Takayama NOTE: activate the output ring with the @code{setring} command. In this ring, @* - the ideal LD (which is a Groebner basis) is the needed D-module structure, @* which is obtained by substituting the minimal integer root of a Bernstein @* polynomial into the s-parametric ideal; @* - the list BS contains roots with multiplicities of a Bernstein polynomial of f. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example annfsOT; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*(N+2); int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "u"; RName[2] = "v"; RName[3] = "t"; RName[4] = "Dt"; for(i=1;i<=N;i++) { for(j=1; j<=size(RName);j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include u,v,t,Dt"); } } } // now, create the names for new vars tmp[1] = "u"; tmp[2] = "v"; list UName = tmp; list DName; for(i=1;i<=N;i++) { DName[i] = "D"+Name[i]; // concat } tmp = 0; tmp[1] = "t"; tmp[2] = "Dt"; list NName = UName + tmp + Name + DName; L[2] = NName; tmp = 0; // Name, Dname will be used further kill UName; kill NName; // block ord (a(1,1),dp); tmp[1] = "a"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for(i=1;i<=Nnew;i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[3,4]=1; for(i=1; i<=N; i++) { @D[4+i,N+4+i]=1; } // @D[N+3,2*(N+2)]=1; old t,Dt stuff // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(u,v,t,Dt,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = u*F-t,u*v-1; poly p; for(i=1; i<=N; i++) { p = u*Dt; // u*Dt p = diff(F,var(4+i))*p; I = I, var(N+4+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of u,v in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1..2); dbprint(ppl,"// -1-3- u,v are eliminated"); dbprint(ppl-1, K); // K is without u,v setring save; // ------------ new ring @R2 ------------------ // without u,v and with the elim.ord for t,Dt // tensored with the K[s] // keep: N, i,j,s, tmp, RL Nnew = 2*N+2+1; // list RL = ringlist(save); // is defined earlier L = 0; // kill L; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2]; list RName; RName[1] = "t"; RName[2] = "Dt"; // now, create the names for new var (here, s only) tmp[1] = "s"; // DName is defined earlier list NName = RName + Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (a(1,1),dp); tmp[1] = "a"; iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with a(1,1,1,1)... tmp[1] = "dp"; s = "iv="; for(i=1; i<= Nnew; i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill NName; tmp[2] = iv; Lord[2] = tmp; // extra block for s // tmp[1] = "dp"; iv = 1; // s[size(s)]= ","; s = s + "1,1,1;"; execute(s); tmp[2] = iv; // Lord[3] = tmp; kill s; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now, add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; @D[1,2] = 1; for(i=1; i<=N; i++) { @D[2+i,2+N+i] = 1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(t,Dt,_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); MM = 0,0,MM; map R01 = @R, MM; ideal K = R01(K); // ideal K = imap(@R,K); // names of vars are important! poly G = t*Dt+s+1; // s is a variable here K = NF(K,std(G)),G; // -------- the ideal K_(@R2) is ready ---------- dbprint(ppl,"// -2-2- starting the elimination of t,Dt in @R2"); dbprint(ppl-1, K); ideal M = engine(K,eng); ideal K2 = nselect(M,1..2); dbprint(ppl,"// -2-3- t,Dt are eliminated"); dbprint(ppl-1, K2); // dbprint(ppl-1+1," -2-4- std of K2"); // option(redSB); option(redTail); K2 = std(K2); // K2; // without t,Dt, and with s // -------- the ring @R3 ---------- // _x, _Dx, s; elim.ord for _x,_Dx. // keep: N, i,j,s, tmp, RL setring save; Nnew = 2*N+1; // list RL = ringlist(save); // is defined earlier // kill L; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2]; // now, create the names for new var (here, s only) tmp[1] = "s"; // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (a(1,1...),dp); string s = "iv="; for(i=1; i<=Nnew-1; i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); tmp[1] = "a"; // string tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s[size(s)]=","; s= s+"1;"; execute(s); kill s; kill NName; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R3@ = ring(L); setring @R3@; matrix @D[Nnew][Nnew]; for(i=1; i<=N; i++) { @D[i,N+i]=1; } def @R3 = nc_algebra(1,@D); setring @R3; kill @R3@; dbprint(ppl,"// -3-1- the ring @R3(_x,_Dx,s) is ready"); dbprint(ppl-1, @R3); ideal MM = maxideal(1); MM = 0,0,MM; map R12 = @R2, MM; ideal K = R12(K2); poly F = imap(save,F); K = K,F; dbprint(ppl,"// -3-2- starting the elimination of _x,_Dx in @R3"); dbprint(ppl-1, K); ideal M = engine(K,eng); ideal K3 = nselect(M,1..Nnew-1); dbprint(ppl,"// -3-3- _x,_Dx are eliminated in @R3"); dbprint(ppl-1, K3); // the ring @R4 and the search for minimal negative int s ring @R4 = 0,(s),dp; dbprint(ppl,"// -4-1- the ring @R4 is ready"); ideal K4 = imap(@R3,K3); poly p = K4[1]; dbprint(ppl,"// -4-2- factorization"); //// ideal P = factorize(p,1); // without constants and multiplicities list P = factorize(p); // with constants and multiplicities ideal bs; intvec m; // the Bernstein polynomial 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]; } // "------ b-function factorizes into ----------"; P; //// int sP = minIntRoot(P, 1); int sP = minIntRoot(bs,1); dbprint(ppl,"// -4-3- minimal integer root found"); dbprint(ppl-1, sP); // convert factors to a list of their roots // assume all factors are linear //// ideal BS = normalize(P); //// BS = subst(BS,s,0); //// BS = -BS; bs = normalize(bs); bs = subst(bs,s,0); bs = -bs; list BS = bs,m; // TODO: sort BS! // ------ substitute s found in the ideal ------ // ------- going back to @R2 and substitute -------- setring @R2; ideal K3 = subst(K2,s,sP); // create the ordinary Weyl algebra and put the result into it, // thus creating the ring @R5 // keep: N, i,j,s, tmp, RL setring save; Nnew = 2*N; // list RL = ringlist(save); // is defined earlier kill Lord, tmp, iv; L = 0; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2]; // DName is defined earlier list NName = Name + DName; L[2] = NName; // dp ordering; string s = "iv="; for(i=1;i<=Nnew;i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); tmp = 0; tmp[1] = "dp"; // string tmp[2] = iv; //intvec Lord[1] = tmp; kill s; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list // Add: Plural part def @R5@ = ring(L); setring @R5@; matrix @D[Nnew][Nnew]; for(i=1; i<=N; i++) { @D[i,N+i]=1; } def @R5 = nc_algebra(1,@D); setring @R5; kill @R5@; dbprint(ppl,"// -5-1- the ring @R5 is ready"); dbprint(ppl-1, @R5); ideal K5 = imap(@R2,K3); option(redSB); dbprint(ppl,"// -5-2- the final cosmetic std"); K5 = engine(K5,eng); // std does the job too // total cleanup kill @R; kill @R2; kill @R3; //// ideal BS = imap(@R4,BS); list BS = imap(@R4,BS); export BS; ideal LD = K5; kill @R4; export LD; return(@R5); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^2+y^3+z^5; printlevel = 0; def A = annfsOT(F); setring A; LD; BS; } proc SannfsOT(poly F, list #) "USAGE: SannfsOT(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the @* 1st step of the algorithm by Oaku and Takayama in the ring D[s] NOTE: activate the output ring with the @code{setring} command. @* In the output ring D[s], the ideal LD (which is NOT a Groebner basis) @* is the needed D-module structure. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example SannfsOT; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*(N+2); int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "u"; RName[2] = "v"; RName[3] = "t"; RName[4] = "Dt"; for(i=1;i<=N;i++) { for(j=1; j<=size(RName);j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include u,v,t,Dt"); } } } // now, create the names for new vars tmp[1] = "u"; tmp[2] = "v"; list UName = tmp; list DName; for(i=1;i<=N;i++) { DName[i] = "D"+Name[i]; // concat } tmp = 0; tmp[1] = "t"; tmp[2] = "Dt"; list NName = UName + tmp + Name + DName; L[2] = NName; tmp = 0; // Name, Dname will be used further kill UName; kill NName; // block ord (a(1,1),dp); tmp[1] = "a"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for(i=1;i<=Nnew;i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[3,4]=1; for(i=1; i<=N; i++) { @D[4+i,N+4+i]=1; } // @D[N+3,2*(N+2)]=1; old t,Dt stuff // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(u,v,t,Dt,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = u*F-t,u*v-1; poly p; for(i=1; i<=N; i++) { p = u*Dt; // u*Dt p = diff(F,var(4+i))*p; I = I, var(N+4+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of u,v in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1..2); dbprint(ppl,"// -1-3- u,v are eliminated"); dbprint(ppl-1, K); // K is without u,v setring save; // ------------ new ring @R2 ------------------ // without u,v and with the elim.ord for t,Dt // tensored with the K[s] // keep: N, i,j,s, tmp, RL Nnew = 2*N+2+1; // list RL = ringlist(save); // is defined earlier L = 0; // kill L; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2]; list RName; RName[1] = "t"; RName[2] = "Dt"; // now, create the names for new var (here, s only) tmp[1] = "s"; // DName is defined earlier list NName = RName + Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (a(1,1),dp); tmp[1] = "a"; iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with a(1,1,1,1)... tmp[1] = "dp"; s = "iv="; for(i=1; i<= Nnew; i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill NName; tmp[2] = iv; Lord[2] = tmp; // extra block for s // tmp[1] = "dp"; iv = 1; // s[size(s)]= ","; s = s + "1,1,1;"; execute(s); tmp[2] = iv; // Lord[3] = tmp; kill s; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now, add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; @D[1,2] = 1; for(i=1; i<=N; i++) { @D[2+i,2+N+i] = 1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(t,Dt,_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); MM = 0,0,MM; map R01 = @R, MM; ideal K = R01(K); // ideal K = imap(@R,K); // names of vars are important! poly G = t*Dt+s+1; // s is a variable here K = NF(K,std(G)),G; // -------- the ideal K_(@R2) is ready ---------- dbprint(ppl,"// -2-2- starting the elimination of t,Dt in @R2"); dbprint(ppl-1, K); ideal M = engine(K,eng); ideal K2 = nselect(M,1..2); dbprint(ppl,"// -2-3- t,Dt are eliminated"); dbprint(ppl-1, K2); K2 = engine(K2,eng); setring save; // ----------- the ring @R3 ------------ // _x, _Dx,s; elim.ord for _x,_Dx. // keep: N, i,j,s, tmp, RL Nnew = 2*N+1; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var tmp[1] = "s"; // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (dp(N),dp); string s = "iv="; for (i=1; i<=Nnew-1; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); tmp[1] = "dp"; // string tmp[2] = iv; // intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s[size(s)] = ","; s = s+"1;"; execute(s); kill s; kill NName; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2]=iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R3@ = ring(L); setring @R3@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R3 = nc_algebra(1,@D); setring @R3; kill @R3@; dbprint(ppl,"// -3-1- the ring @R3(_x,_Dx,s) is ready"); dbprint(ppl-1, @R3); ideal MM = maxideal(1); MM = 0,s,MM; map R01 = @R2, MM; ideal K2 = R01(K2); // total cleanup ideal LD = K2; // make leadcoeffs positive for (i=1; i<= ncols(LD); i++) { if (leadcoef(LD[i]) <0 ) { LD[i] = -LD[i]; } } export LD; kill @R; kill @R2; return(@R3); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^3+y^3+z^3; printlevel = 0; def A = SannfsOT(F); setring A; LD; } proc SannfsBM(poly F, list #) "USAGE: SannfsBM(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the @* 1st step of the algorithm by Briancon and Maisonobe in the ring D[s]. NOTE: activate the output ring with the @code{setring} command. @* In the output ring D[s], the ideal LD (which is NOT a Groebner basis) is @* the needed D-module structure. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example SannfsBM; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*N+2; int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "t"; RName[2] = "s"; for(i=1;i<=N;i++) { for(j=1; j<=size(RName);j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include t,s"); } } } // now, create the names for new vars list DName; for(i=1;i<=N;i++) { DName[i] = "D"+Name[i]; // concat } tmp[1] = "t"; tmp[2] = "s"; list NName = tmp + Name + DName; L[2] = NName; // Name, Dname will be used further kill NName; // block ord (lp(2),dp); tmp[1] = "lp"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for(i=1;i<=Nnew;i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill s; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2]=t; for(i=1; i<=N; i++) { @D[2+i,N+2+i]=1; } // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = t*F+s; poly p; for(i=1; i<=N; i++) { p = t; // t p = diff(F,var(2+i))*p; I = I, var(N+2+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of t in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1); dbprint(ppl,"// -1-3- t is eliminated"); dbprint(ppl-1, K); // K is without t K = engine(K,eng); // std does the job too // now, we must change the ordering // and create a ring without t, Dt // setring S; // ----------- the ring @R3 ------------ // _x, _Dx,s; elim.ord for _x,_Dx. // keep: N, i,j,s, tmp, RL Nnew = 2*N+1; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; list L=imap(save,L); list RL=imap(save,RL); L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var tmp[1] = "s"; // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (dp(N),dp); string s = "iv="; for (i=1; i<=Nnew-1; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); tmp[1] = "dp"; // string tmp[2] = iv; // intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s[size(s)] = ","; s = s+"1;"; execute(s); kill s; kill NName; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2]=iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); MM = 0,s,MM; map R01 = @R, MM; ideal K = R01(K); // total cleanup ideal LD = K; // make leadcoeffs positive for (i=1; i<= ncols(LD); i++) { if (leadcoef(LD[i]) <0 ) { LD[i] = -LD[i]; } } export LD; kill @R; return(@R2); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^3+y^3+z^3; printlevel = 0; def A = SannfsBM(F); setring A; LD; } proc SannfsBFCT(poly F, list #) "USAGE: SannfsBFCT(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute Ann f^s and Groebner basis of Ann f^s+f in D[s] NOTE: activate the output ring with the @code{setring} command. @* This procedure, unlike SannfsBM, returns a ring with the degrevlex @* ordering in all variables. @* In the ring D[s], the ideal LD is the ideal needed (which is returned as a Groebner basis). @* If eng <>0, @code{std} is used for Groebner basis 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 SannfsBFCT; shows examples " { // DEBUG INFO: ordering on the output ring = dp, engine on the sum of ideals is used // hence it's the situation for slimgb int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*N+2; int i,j; list RL = ringlist(basering); list L, Lord; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "@t"; RName[2] = "@s"; list PName; int BadName; if (size(RL[1])>1) { PName = RL[1][2]; } // PName; for(i=1;i<=2;i++) { for(j=1; j<=N; j++) { if (Name[j] == RName[i]) { BadName = 1; break; } } for(j=1; j<=size(PName); j++) { if (PName[j] == RName[i]) { BadName = 1; break; } } } if (BadName == 1) { ERROR("Variable/parameter names should not include @t,@s"); } kill PName,BadName; // now, create the names for new vars list DName; for(i=1;i<=N;i++) { DName[i] = "D"+Name[i]; // concat } list NName = RName + DName + Name ; L[2] = NName; // Name, Dname will be used further kill NName; // block ord (lp(2),dp); Lord[1] = list("lp",intvec(1,1)); // continue with dp 1,1,1,1... Lord[2] = list("dp",intvec(1:Nnew)); Lord[3] = list("C",intvec(0)); L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2] = var(1); for(i=1; i<=N; i++) { @D[2+i,N+2+i]=-1; } // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(@t,@s,_Dx,_x) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = var(1)*F+var(2); poly p; for(i=1; i<=N; i++) { p = var(1); // t p = diff(F,var(N+2+i))*p; I = I, var(2+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of @t in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1); dbprint(ppl,"// -1-3- @t is eliminated"); dbprint(ppl-1, K); // K is without @t K = engine(K,eng); // std does the job too // now, we must change the ordering // and create a ring without @t // setring S; // ----------- the ring @R3 ------------ // _Dx,_x,s; +fast ord ! // keep: N, i,j,RL Nnew = 2*N+1; kill Lord, RName; list Lord; list L=imap(save,L); list RL=imap(save,RL); L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars hava admissible names -> done earlier // DName is defined earlier list NName = DName + Name + list(string(var(2))); L[2] = NName; kill NName; // just dp // block ord (dp(N),dp); Lord[1] = list("dp",intvec(1:Nnew)); Lord[2]= list("C",intvec(0)); L[3] = Lord; // we are done with the list. Now add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=-1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_Dx,_x,@s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); MM = 0,var(Nnew),MM; map R01 = @R, MM; ideal K = R01(K); // total cleanup poly F = imap(save, F); ideal LD = K,F; dbprint(ppl,"// -2-2- start GB computations for Ann f^s + f"); dbprint(ppl-1, LD); LD = engine(LD,eng); dbprint(ppl,"// -2-3- finished GB computations for Ann f^s + f"); dbprint(ppl-1, LD); // make leadcoeffs positive for (i=1; i<= ncols(LD); i++) { if (leadcoef(LD[i]) <0 ) { LD[i] = -LD[i]; } } export LD; kill @R; return(@R2); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z,w),Dp; poly F = x^3+y^3+z^3*w; printlevel = 0; def A = SannfsBFCT(F); setring A; intvec v = 1,2,3,4,5,6,7,8; // are there polynomials, depending on @s only? nselect(LD,v); // a fancier example: def R = reiffen(4,5); setring R; v = 1,2,3,4; RC; // the Reiffen curve in 4,5 def B = SannfsBFCT(RC); setring B; // Are there polynomials, depending on @s only? nselect(LD,v); // It is not the case. Are there leading terms in @s only? nselect(lead(LD),v); } proc SannfsBFCTstd(poly F, list #) "USAGE: SannfsBFCTstd(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute Ann f^s and Groebner basis of Ann f^s+f in D[s] NOTE: activate the output ring with the @code{setring} command. @* This procedure, unlike SannfsBM, returns a ring with the degrevlex @* ordering in all variables. @* In the ring D[s], the ideal LD (which IS a Groebner basis) is the needed ideal. @* In this procedure @code{std} is used for Groebner basis computations. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example SannfsBFCTstd; shows examples " { // DEBUG INFO: ordering on the output ring = dp, // use std(K,F); for reusing the std property of K int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); int Nnew = 2*N+2; int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; RName[1] = "@t"; RName[2] = "@s"; for(i=1;i<=N;i++) { for(j=1; j<=size(RName);j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include @t,@s"); } } } // now, create the names for new vars list DName; for(i=1;i<=N;i++) { DName[i] = "D"+Name[i]; // concat } tmp[1] = "t"; tmp[2] = "s"; list NName = tmp + DName + Name ; L[2] = NName; // Name, Dname will be used further kill NName; // block ord (lp(2),dp); tmp[1] = "lp"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for(i=1;i<=Nnew;i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); kill s; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2]=t; for(i=1; i<=N; i++) { @D[2+i,N+2+i]=-1; } // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(t,s,_Dx,_x) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); ideal I = t*F+s; poly p; for(i=1; i<=N; i++) { p = t; // t p = diff(F,var(N+2+i))*p; I = I, var(2+i) + p; } // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of t in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1); dbprint(ppl,"// -1-3- t is eliminated"); dbprint(ppl-1, K); // K is without t K = engine(K,eng); // std does the job too // now, we must change the ordering // and create a ring without t // setring S; // ----------- the ring @R3 ------------ // _Dx,_x,s; +fast ord ! // keep: N, i,j,s, tmp, RL Nnew = 2*N+1; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; list L=imap(save,L); list RL=imap(save,RL); L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var tmp[1] = "s"; // DName is defined earlier list NName = DName + Name + tmp; L[2] = NName; tmp = 0; // just dp // block ord (dp(N),dp); string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); tmp[1] = "dp"; // string tmp[2] = iv; // intvec Lord[1] = tmp; kill s; kill NName; tmp[1] = "C"; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=-1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_Dx,_x,s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); MM = 0,s,MM; map R01 = @R, MM; ideal K = R01(K); // total cleanup poly F = imap(save, F); // ideal LD = K,F; dbprint(ppl,"// -2-2- start GB computations for Ann f^s + f"); // dbprint(ppl-1, LD); ideal LD = std(K,F); // LD = engine(LD,eng); dbprint(ppl,"// -2-3- finished GB computations for Ann f^s + f"); dbprint(ppl-1, LD); // make leadcoeffs positive for (i=1; i<= ncols(LD); i++) { if (leadcoef(LD[i]) <0 ) { LD[i] = -LD[i]; } } export LD; kill @R; return(@R2); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z,w),Dp; poly F = x^3+y^3+z^3*w; printlevel = 0; def A = SannfsBFCT(F); setring A; intvec v = 1,2,3,4,5,6,7,8; // are there polynomials, depending on s only? nselect(LD,v); // a fancier example: def R = reiffen(4,5); setring R; v = 1,2,3,4; RC; // the Reiffen curve in 4,5 def B = SannfsBFCT(RC); setring B; // Are there polynomials, depending on s only? nselect(LD,v); // It is not the case. Are there leading monomials in s only? nselect(lead(LD),v); } // use a finer ordering proc SannfsLOT(poly F, list #) "USAGE: SannfsLOT(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the @* Levandovskyy's modification of the algorithm by Oaku and Takayama in D[s] NOTE: activate the output ring with the @code{setring} command. @* In the ring D[s], the ideal LD (which is NOT a Groebner basis) is @* the needed D-module structure. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example SannfsLOT; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); // int Nnew = 2*(N+2); int Nnew = 2*(N+1)+1; //removed u,v; added s int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; // RName[1] = "u"; // RName[2] = "v"; RName[1] = "t"; RName[2] = "Dt"; for(i=1;i<=N;i++) { for(j=1; j<=size(RName);j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include t,Dt"); } } } // now, create the names for new vars // tmp[1] = "u"; // tmp[2] = "v"; // list UName = tmp; list DName; for(i=1;i<=N;i++) { DName[i] = "D"+Name[i]; // concat } tmp = 0; tmp[1] = "t"; tmp[2] = "Dt"; list SName ; SName[1] = "s"; // list NName = tmp + Name + DName + SName; list NName = tmp + SName + Name + DName; L[2] = NName; tmp = 0; // Name, Dname will be used further // kill UName; kill NName; // block ord (a(1,1),dp); tmp[1] = "a"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with a(0,0,1) tmp[1] = "a"; // string iv = 0,0,1; tmp[2] = iv; //intvec Lord[2] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for(i=1;i<=Nnew;i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); tmp[2] = iv; Lord[3] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[4] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2]=1; for(i=1; i<=N; i++) { @D[3+i,N+3+i]=1; } // ADD [s,t]=-t, [s,Dt]=Dt @D[1,3] = -var(1); @D[2,3] = var(2); // @D[N+3,2*(N+2)]=1; old t,Dt stuff // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(t,Dt,s,_x,_Dx) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); // ideal I = u*F-t,u*v-1; ideal I = F-t; poly p; for(i=1; i<=N; i++) { // p = u*Dt; // u*Dt p = Dt; p = diff(F,var(3+i))*p; I = I, var(N+3+i) + p; } // I = I, var(1)*var(2) + var(Nnew) +1; // reduce it with t-f!!! // t*Dt + s +1 reduced with t-f gives f*Dt + s I = I, F*var(2) + var(3); // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of t,Dt in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1..2); dbprint(ppl,"// -1-3- t,Dt are eliminated"); dbprint(ppl-1, K); // K is without t, Dt K = engine(K,eng); // std does the job too // now, we must change the ordering // and create a ring without t, Dt setring save; // ----------- the ring @R3 ------------ // _x, _Dx,s; elim.ord for _x,_Dx. // keep: N, i,j,s, tmp, RL Nnew = 2*N+1; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var tmp[1] = "s"; // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (dp(N),dp); // string s is already defined s = "iv="; for (i=1; i<=Nnew-1; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); tmp[1] = "dp"; // string tmp[2] = iv; // intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s[size(s)] = ","; s = s+"1;"; execute(s); kill s; kill NName; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2]=iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); // MM = 0,s,MM; MM = 0,0,s,MM[1..size(MM)-1]; map R01 = @R, MM; ideal K = R01(K); // total cleanup ideal LD = K; // make leadcoeffs positive for (i=1; i<= ncols(LD); i++) { if (leadcoef(LD[i]) <0 ) { LD[i] = -LD[i]; } } export LD; kill @R; return(@R2); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^3+y^3+z^3; printlevel = 0; def A = SannfsLOT(F); setring A; LD; } /* proc SannfsLOTold(poly F, list #) "USAGE: SannfsLOT(f [,eng]); f a poly, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the Levandovskyy's modification of the algorithm by Oaku and Takayama in the ring D[s], where D is the Weyl algebra NOTE: activate the output ring with the @code{setring} command. @* In the ring D[s], the ideal LD (which is NOT a Groebner basis) is the needed D-module structure. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise, and by default @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example SannfsLOT; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } // returns a list with a ring and an ideal LD in it int ppl = printlevel-voice+2; // printf("plevel :%s, voice: %s",printlevel,voice); def save = basering; int N = nvars(basering); // int Nnew = 2*(N+2); int Nnew = 2*(N+1)+1; //removed u,v; added s int i,j; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether vars have admissible names list Name = RL[2]; list RName; // RName[1] = "u"; // RName[2] = "v"; RName[1] = "t"; RName[2] = "Dt"; for(i=1;i<=N;i++) { for(j=1; j<=size(RName);j++) { if (Name[i] == RName[j]) { ERROR("Variable names should not include t,Dt"); } } } // now, create the names for new vars // tmp[1] = "u"; // tmp[2] = "v"; // list UName = tmp; list DName; for(i=1;i<=N;i++) { DName[i] = "D"+Name[i]; // concat } tmp = 0; tmp[1] = "t"; tmp[2] = "Dt"; list SName ; SName[1] = "s"; // list NName = UName + tmp + Name + DName; list NName = tmp + Name + DName + SName; L[2] = NName; tmp = 0; // Name, Dname will be used further // kill UName; kill NName; // block ord (a(1,1),dp); tmp[1] = "a"; // string iv = 1,1; tmp[2] = iv; //intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s = "iv="; for(i=1;i<=Nnew;i++) { s = s+"1,"; } s[size(s)]= ";"; execute(s); tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; @D[1,2]=1; for(i=1; i<=N; i++) { @D[2+i,N+2+i]=1; } // ADD [s,t]=-t, [s,Dt]=Dt @D[1,Nnew] = -var(1); @D[2,Nnew] = var(2); // @D[N+3,2*(N+2)]=1; old t,Dt stuff // L[5] = matrix(UpOneMatrix(Nnew)); // L[6] = @D; def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(t,Dt,_x,_Dx,s) is ready"); dbprint(ppl-1, @R); // create the ideal I poly F = imap(save,F); // ideal I = u*F-t,u*v-1; ideal I = F-t; poly p; for(i=1; i<=N; i++) { // p = u*Dt; // u*Dt p = Dt; p = diff(F,var(2+i))*p; I = I, var(N+2+i) + p; } // I = I, var(1)*var(2) + var(Nnew) +1; // reduce it with t-f!!! // t*Dt + s +1 reduced with t-f gives f*Dt + s I = I, F*var(2) + var(Nnew); // -------- the ideal I is ready ---------- dbprint(ppl,"// -1-2- starting the elimination of t,Dt in @R"); dbprint(ppl-1, I); ideal J = engine(I,eng); ideal K = nselect(J,1..2); dbprint(ppl,"// -1-3- t,Dt are eliminated"); dbprint(ppl-1, K); // K is without t, Dt K = engine(K,eng); // std does the job too // now, we must change the ordering // and create a ring without t, Dt setring save; // ----------- the ring @R3 ------------ // _x, _Dx,s; elim.ord for _x,_Dx. // keep: N, i,j,s, tmp, RL Nnew = 2*N+1; kill Lord, tmp, iv, RName; list Lord, tmp; intvec iv; L[1] = RL[1]; L[4] = RL[4]; // char, minpoly // check whether vars hava admissible names -> done earlier // now, create the names for new var tmp[1] = "s"; // DName is defined earlier list NName = Name + DName + tmp; L[2] = NName; tmp = 0; // block ord (dp(N),dp); // string s is already defined s = "iv="; for (i=1; i<=Nnew-1; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); tmp[1] = "dp"; // string tmp[2] = iv; // intvec Lord[1] = tmp; // continue with dp 1,1,1,1... tmp[1] = "dp"; // string s[size(s)] = ","; s = s+"1;"; execute(s); kill s; kill NName; tmp[2] = iv; Lord[2] = tmp; tmp[1] = "C"; iv = 0; tmp[2]=iv; Lord[3] = tmp; tmp = 0; L[3] = Lord; // we are done with the list. Now add a Plural part def @R2@ = ring(L); setring @R2@; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R2 = nc_algebra(1,@D); setring @R2; kill @R2@; dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); ideal MM = maxideal(1); MM = 0,s,MM; map R01 = @R, MM; ideal K = R01(K); // total cleanup ideal LD = K; // make leadcoeffs positive for (i=1; i<= ncols(LD); i++) { if (leadcoef(LD[i]) <0 ) { LD[i] = -LD[i]; } } export LD; kill @R; return(@R2); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^3+y^3+z^3; printlevel = 0; def A = SannfsLOTold(F); setring A; LD; } */ proc annfsLOT(poly F, list #) "USAGE: annfsLOT(F [,eng]); F a poly, eng an optional int RETURN: ring PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to @* the Levandovskyy's modification of the algorithm by Oaku and Takayama NOTE: activate the output ring with the @code{setring} command. In this ring, @* - the ideal LD (which is a Groebner basis) is the needed D-module structure, @* which is obtained by substituting the minimal integer root of a Bernstein @* polynomial into the s-parametric ideal; @* - the list BS contains the roots with multiplicities of BS polynomial of f. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise and by default @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example annfsLOT; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } printlevel=printlevel+1; def save = basering; def @A = SannfsLOT(F,eng); setring @A; poly F = imap(save,F); def B = annfs0(LD,F,eng); return(B); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = z*x^2+y^3; printlevel = 0; def A = annfsLOT(F); setring A; LD; BS; } proc annfs0(ideal I, poly F, list #) "USAGE: annfs0(I, F [,eng]); I an ideal, F a poly, eng an optional int RETURN: ring PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based @* on the output of Sannfs-like procedure NOTE: activate the output ring with the @code{setring} command. In this ring, @* - the ideal LD (which is a Groebner basis) is the annihilator of f^s, @* - the list BS contains the roots with multiplicities of BS polynomial of f. @* If eng <>0, @code{std} is used for Groebner basis computations, @* otherwise and by default @code{slimgb} is used. @* If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example annfs0; shows examples " { int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } def @R2 = basering; // we're in D_n[s], where the elim ord for s is set ideal J = NF(I,std(F)); // make leadcoeffs positive int i; for (i=1; i<= ncols(J); i++) { if (leadcoef(J[i]) <0 ) { J[i] = -J[i]; } } J = J,F; ideal M = engine(J,eng); int Nnew = nvars(@R2); ideal K2 = nselect(M,1..Nnew-1); int ppl = printlevel-voice+2; dbprint(ppl,"// -1-1- _x,_Dx are eliminated in basering"); dbprint(ppl-1, K2); // the ring @R3 and the search for minimal negative int s ring @R3 = 0,s,dp; dbprint(ppl,"// -2-1- the ring @R3 i.e. K[s] is ready"); ideal K3 = imap(@R2,K2); poly p = K3[1]; dbprint(ppl,"// -2-2- factorization"); // ideal P = factorize(p,1); //without constants and multiplicities // "--------- b-function factorizes into ---------"; P; // convert factors to the list of their roots with mults // assume all factors are linear // ideal BS = normalize(P); // BS = subst(BS,s,0); // BS = -BS; list P = factorize(p); //with constants and multiplicities ideal bs; intvec m; //the Bernstein polynomial 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]; } int sP = minIntRoot(bs,1); bs = normalize(bs); bs = -subst(bs,s,0); dbprint(ppl,"// -2-3- minimal integer root found"); dbprint(ppl-1, sP); //TODO: sort BS! // --------- substitute s found in the ideal --------- // --------- going back to @R and substitute --------- setring @R2; K2 = subst(I,s,sP); // create the ordinary Weyl algebra and put the result into it, // thus creating the ring @R5 // keep: N, i,j,s, tmp, RL Nnew = Nnew - 1; // former 2*N; // list RL = ringlist(save); // is defined earlier // kill Lord, tmp, iv; list L = 0; list Lord, tmp; intvec iv; list RL = ringlist(basering); L[1] = RL[1]; L[4] = RL[4]; //char, minpoly // check whether vars have admissible names -> done earlier // list Name = RL[2]M // DName is defined earlier list NName; // = RL[2]; // skip the last var 's' for (i=1; i<=Nnew; i++) { NName[i] = RL[2][i]; } L[2] = NName; // dp ordering; string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)] = ";"; execute(s); tmp = 0; tmp[1] = "dp"; // string tmp[2] = iv; // intvec Lord[1] = tmp; kill s; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list // Add: Plural part def @R4@ = ring(L); setring @R4@; int N = Nnew/2; matrix @D[Nnew][Nnew]; for (i=1; i<=N; i++) { @D[i,N+i]=1; } def @R4 = nc_algebra(1,@D); setring @R4; kill @R4@; dbprint(ppl,"// -3-1- the ring @R4 is ready"); dbprint(ppl-1, @R4); ideal K4 = imap(@R2,K2); option(redSB); dbprint(ppl,"// -3-2- the final cosmetic std"); K4 = engine(K4,eng); // std does the job too // total cleanup ideal bs = imap(@R3,bs); kill @R3; list BS = bs,m; export BS; ideal LD = K4; export LD; return(@R4); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x^3+y^3+z^3; printlevel = 0; def A = SannfsBM(F); setring A; // alternatively, one can use SannfsOT or SannfsLOT LD; poly F = imap(r,F); def B = annfs0(LD,F); setring B; LD; BS; } // proc annfsgms(poly F, list #) // "USAGE: annfsgms(f [,eng]); f a poly, eng an optional int // ASSUME: f has an isolated critical point at 0 // RETURN: ring // PURPOSE: compute the D-module structure of basering[1/f]*f^s // NOTE: activate the output ring with the @code{setring} command. In this ring, // @* - the ideal LD is the needed D-mod structure, // @* - the ideal BS is the list of roots of a Bernstein polynomial of f. // @* If eng <>0, @code{std} is used for Groebner basis computations, // @* otherwise (and by default) @code{slimgb} is used. // @* If printlevel=1, progress debug messages will be printed, // @* if printlevel>=2, all the debug messages will be printed. // EXAMPLE: example annfsgms; shows examples // " // { // LIB "gmssing.lib"; // int eng = 0; // if ( size(#)>0 ) // { // if ( typeof(#[1]) == "int" ) // { // eng = int(#[1]); // } // } // int ppl = printlevel-voice+2; // // returns a ring with the ideal LD in it // def save = basering; // // compute the Bernstein poly from gmssing.lib // list RL = ringlist(basering); // // in the descr. of the ordering, replace "p" by "s" // list NL = convloc(RL); // // create a ring with the ordering, converted to local // def @LR = ring(NL); // setring @LR; // poly F = imap(save, F); // ideal B = bernstein(F)[1]; // // since B may not contain (s+1) [following gmssing.lib] // // add it! // B = B,-1; // B = simplify(B,2+4); // erase zero and repeated entries // // find the minimal integer value // int S = minIntRoot(B,0); // dbprint(ppl,"// -0- minimal integer root found"); // dbprint(ppl-1,S); // setring save; // int N = nvars(basering); // int Nnew = 2*(N+2); // int i,j; // string s; // // list RL = ringlist(basering); // list L, Lord; // list tmp; // intvec iv; // L[1] = RL[1]; // char // L[4] = RL[4]; // char, minpoly // // check whether vars have admissible names // list Name = RL[2]; // list RName; // RName[1] = "u"; // RName[2] = "v"; // RName[3] = "t"; // RName[4] = "Dt"; // for(i=1;i<=N;i++) // { // for(j=1; j<=size(RName);j++) // { // if (Name[i] == RName[j]) // { // ERROR("Variable names should not include u,v,t,Dt"); // } // } // } // // now, create the names for new vars // // tmp[1] = "u"; tmp[2] = "v"; tmp[3] = "t"; tmp[4] = "Dt"; // list UName = RName; // list DName; // for(i=1;i<=N;i++) // { // DName[i] = "D"+Name[i]; // concat // } // list NName = UName + Name + DName; // L[2] = NName; // tmp = 0; // // Name, Dname will be used further // kill UName; // kill NName; // // block ord (a(1,1),dp); // tmp[1] = "a"; // string // iv = 1,1; // tmp[2] = iv; //intvec // Lord[1] = tmp; // // continue with dp 1,1,1,1... // tmp[1] = "dp"; // string // s = "iv="; // for(i=1; i<=Nnew; i++) // need really all vars! // { // s = s+"1,"; // } // s[size(s)]= ";"; // execute(s); // tmp[2] = iv; // Lord[2] = tmp; // tmp[1] = "C"; // iv = 0; // tmp[2] = iv; // Lord[3] = tmp; // tmp = 0; // L[3] = Lord; // // we are done with the list // def @R = ring(L); // setring @R; // matrix @D[Nnew][Nnew]; // @D[3,4] = 1; // t,Dt // for(i=1; i<=N; i++) // { // @D[4+i,4+N+i]=1; // } // // L[5] = matrix(UpOneMatrix(Nnew)); // // L[6] = @D; // nc_algebra(1,@D); // dbprint(ppl,"// -1-1- the ring @R is ready"); // dbprint(ppl-1,@R); // // create the ideal // poly F = imap(save,F); // ideal I = u*F-t,u*v-1; // poly p; // for(i=1; i<=N; i++) // { // p = u*Dt; // u*Dt // p = diff(F,var(4+i))*p; // I = I, var(N+4+i) + p; // Dx, Dy // } // // add the relations between t,Dt and s // // I = I, t*Dt+1+S; // // -------- the ideal I is ready ---------- // dbprint(ppl,"// -1-2- starting the elimination of u,v in @R"); // ideal J = engine(I,eng); // ideal K = nselect(J,1..2); // dbprint(ppl,"// -1-3- u,v are eliminated in @R"); // dbprint(ppl-1,K); // without u,v: not yet our answer // //----- create a ring with elim.ord for t,Dt ------- // setring save; // // ------------ new ring @R2 ------------------ // // without u,v and with the elim.ord for t,Dt // // keep: N, i,j,s, tmp, RL // Nnew = 2*N+2; // // list RL = ringlist(save); // is defined earlier // kill Lord,tmp,iv, RName; // L = 0; // list Lord, tmp; // intvec iv; // L[1] = RL[1]; // char // L[4] = RL[4]; // char, minpoly // // check whether vars have admissible names -> done earlier // // list Name = RL[2]; // list RName; // RName[1] = "t"; // RName[2] = "Dt"; // // DName is defined earlier // list NName = RName + Name + DName; // L[2] = NName; // tmp = 0; // // block ord (a(1,1),dp); // tmp[1] = "a"; // string // iv = 1,1; // tmp[2] = iv; //intvec // Lord[1] = tmp; // // continue with dp 1,1,1,1... // tmp[1] = "dp"; // string // s = "iv="; // for(i=1;i<=Nnew;i++) // { // s = s+"1,"; // } // s[size(s)]= ";"; // execute(s); // kill s; // kill NName; // tmp[2] = iv; // Lord[2] = tmp; // tmp[1] = "C"; // iv = 0; // tmp[2] = iv; // Lord[3] = tmp; // tmp = 0; // L[3] = Lord; // // we are done with the list // // Add: Plural part // def @R2 = ring(L); // setring @R2; // matrix @D[Nnew][Nnew]; // @D[1,2]=1; // for(i=1; i<=N; i++) // { // @D[2+i,2+N+i]=1; // } // nc_algebra(1,@D); // dbprint(ppl,"// -2-1- the ring @R2 is ready"); // dbprint(ppl-1,@R2); // ideal MM = maxideal(1); // MM = 0,0,MM; // map R01 = @R, MM; // ideal K2 = R01(K); // // add the relations between t,Dt and s // // K2 = K2, t*Dt+1+S; // poly G = t*Dt+S+1; // K2 = NF(K2,std(G)),G; // dbprint(ppl,"// -2-2- starting elimination for t,Dt in @R2"); // ideal J = engine(K2,eng); // ideal K = nselect(J,1..2); // dbprint(ppl,"// -2-3- t,Dt are eliminated"); // dbprint(ppl-1,K); // // "------- produce a final result --------"; // // ----- create the ordinary Weyl algebra and put // // ----- the result into it // // ------ create the ring @R5 // // keep: N, i,j,s, tmp, RL // setring save; // Nnew = 2*N; // // list RL = ringlist(save); // is defined earlier // kill Lord, tmp, iv; // // kill L; // L = 0; // list Lord, tmp; // intvec iv; // L[1] = RL[1]; // char // L[4] = RL[4]; // char, minpoly // // check whether vars have admissible names -> done earlier // // list Name = RL[2]; // // DName is defined earlier // list NName = Name + DName; // L[2] = NName; // // dp ordering; // string s = "iv="; // for(i=1;i<=2*N;i++) // { // s = s+"1,"; // } // s[size(s)]= ";"; // execute(s); // tmp = 0; // tmp[1] = "dp"; // string // tmp[2] = iv; //intvec // Lord[1] = tmp; // kill s; // tmp[1] = "C"; // iv = 0; // tmp[2] = iv; // Lord[2] = tmp; // tmp = 0; // L[3] = Lord; // // we are done with the list // // Add: Plural part // def @R5 = ring(L); // setring @R5; // matrix @D[Nnew][Nnew]; // for(i=1; i<=N; i++) // { // @D[i,N+i]=1; // } // nc_algebra(1,@D); // dbprint(ppl,"// -3-1- the ring @R5 is ready"); // dbprint(ppl-1,@R5); // ideal K5 = imap(@R2,K); // option(redSB); // dbprint(ppl,"// -3-2- the final cosmetic std"); // K5 = engine(K5,eng); // std does the job too // // total cleanup // kill @R; // kill @R2; // ideal LD = K5; // ideal BS = imap(@LR,B); // kill @LR; // export BS; // export LD; // return(@R5); // } // example // { // "EXAMPLE:"; echo = 2; // ring r = 0,(x,y,z),Dp; // poly F = x^2+y^3+z^5; // def A = annfsgms(F); // setring A; // LD; // print(matrix(BS)); // } proc convloc(list @NL) "USAGE: convloc(L); L a list RETURN: list PURPOSE: convert a ringlist L into another ringlist, @* where all the 'p' orderings are replaced with the 's' orderings, e.g. @code{dp} by @code{ds}. ASSUME: L is a result of a ringlist command EXAMPLE: example convloc; shows examples " { list NL = @NL; // given a ringlist, returns a new ringlist, // where all the p-orderings are replaced with s-ord's int i,j,k,found; int nblocks = size(NL[3]); for(i=1; i<=nblocks; i++) { for(j=1; j<=size(NL[3][i]); j++) { if (typeof(NL[3][i][j]) == "string" ) { found = 0; for (k=1; k<=size(NL[3][i][j]); k++) { if (NL[3][i][j][k]=="p") { NL[3][i][j][k]="s"; found = 1; // printf("replaced at %s,%s,%s",i,j,k); } } } } } return(NL); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),(Dp(2),dp(1)); list L = ringlist(r); list N = convloc(L); def rs = ring(N); setring rs; rs; } proc annfspecial(ideal I, poly F, int mir, number n) "USAGE: annfspecial(I,F,mir,n); I an ideal, F a poly, int mir, number n RETURN: ideal PURPOSE: compute the annihilator ideal of F^n in the Weyl Algebra @* for the given rational number n ASSUME: The basering is D[s] and contains 's' explicitly as a variable, @* the ideal I is the Ann F^s in D[s] (obtained with e.g. SannfsBM), @* the integer 'mir' is the minimal integer root of the BS polynomial of F, @* and the number n is rational. NOTE: We compute the real annihilator for any rational value of n (both @* generic and exceptional). The implementation goes along the lines of @* the Algorithm 5.3.15 from Saito-Sturmfels-Takayama. DISPLAY: If printlevel=1, progress debug messages will be printed, @* if printlevel>=2, all the debug messages will be printed. EXAMPLE: example annfspecial; shows examples " { if (!isRational(n)) { "ERROR: n must be a rational number!"; } int ppl = printlevel-voice+2; // int mir = minIntRoot(L[1],0); int ns = varNum("s"); if (!ns) { ERROR("Variable s expected in the ideal AnnFs"); } int d; ideal P = subst(I,var(ns),n); if (denominator(n) == 1) { // n is fraction-free d = int(numerator(n)); if ( (!d) && (n!=0)) { ERROR("no parametric special values are allowed"); } d = d - mir; if (d>0) { dbprint(ppl,"// -1-1- starting syzygy computations"); matrix J[1][1] = F^d; dbprint(ppl-1,"// -1-1-1- of the polynomial ideal"); dbprint(ppl-1,J); matrix K[1][size(I)] = subst(I,var(ns),mir); dbprint(ppl-1,"// -1-1-2- modulo ideal:"); dbprint(ppl-1, K); module M = modulo(J,K); dbprint(ppl-1,"// -1-1-3- getting the result:"); dbprint(ppl-1, M); P = P, ideal(M); dbprint(ppl,"// -1-2- finished syzygy computations"); } else { dbprint(ppl,"// -1-1- d<=0, no syzygy computations needed"); dbprint(ppl-1,"// -1-2- for d ="); dbprint(ppl-1, d); } } // also the else case: d<=0 or n is not an integer dbprint(ppl,"// -2-1- starting final Groebner basis"); P = groebner(P); dbprint(ppl,"// -2-2- finished final Groebner basis"); return(P); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),dp; poly F = x3-y2; def B = annfs(F); setring B; minIntRoot(BS[1],0); // So, the minimal integer root is -1 setring r; def A = SannfsBM(F); setring A; poly F = x3-y2; annfspecial(LD,F,-1,3/4); // generic root annfspecial(LD,F,-1,-2); // integer but still generic root annfspecial(LD,F,-1,1); // exceptional root } /* //static proc new_ex_annfspecial() { //another example for annfspecial: x3+y3+z3 ring r = 0,(x,y,z),dp; poly F = x3+y3+z3; def B = annfs(F); setring B; minIntRoot(BS[1],0); // So, the minimal integer root is -1 setring r; def A = SannfsBM(F); setring A; poly F = x3+y3+z3; annfspecial(LD,F,-1,3/4); // generic root annfspecial(LD,F,-1,-2); // integer but still generic root annfspecial(LD,F,-1,1); // exceptional root } */ proc minIntRoot(ideal P, int fact) "USAGE: minIntRoot(P, fact); P an ideal, fact an int RETURN: int PURPOSE: minimal integer root of a maximal ideal P NOTE: if fact==1, P is the result of some 'factorize' call, @* else P is treated as the result of bernstein::gmssing.lib @* in both cases without constants and multiplicities EXAMPLE: example minIntRoot; shows examples " { // ideal P = factorize(p,1); // or ideal P = bernstein(F)[1]; intvec vP; number nP; int i; if ( fact ) { // the result comes from "factorize" P = normalize(P); // now leadcoef = 1 // TODO: hunt for units and kill then !!! P = lead(P)-P; // nP = leadcoef(P[i]-lead(P[i])); // for 1 var only, extract the coeff } else { // bernstein deletes -1 usually P = P, -1; } // for both situations: // now we have an ideal of fractions of type "number" int sP = size(P); for(i=1; i<=sP; i++) { nP = leadcoef(P[i]); if ( (nP - int(nP)) == 0 ) { vP = vP,int(nP); } } if ( size(vP)>=2 ) { vP = vP[2..size(vP)]; } sP = -Max(-vP); if (sP == 0) { "Warning: zero root!"; } return(sP); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),ds; poly f1 = x*y*(x+y); ideal I1 = bernstein(f1)[1]; // a local Bernstein poly I1; minIntRoot(I1,0); poly f2 = x2-y3; ideal I2 = bernstein(f2)[1]; I2; minIntRoot(I2,0); // now we illustrate the behaviour of factorize // together with a global ordering ring r2 = 0,x,dp; poly f3 = 9*(x+2/3)*(x+1)*(x+4/3); //global b-poly of f1=x*y*(x+y) ideal I3 = factorize(f3,1); I3; minIntRoot(I3,1); // and a more interesting situation ring s = 0,(x,y,z),ds; poly f = x3 + y3 + z3; ideal I = bernstein(f)[1]; I; minIntRoot(I,0); } proc isHolonomic(def M) "USAGE: isHolonomic(M); M an ideal/module/matrix RETURN: int, 1 if M is holonomic over the base ring, and 0 otherwise ASSUME: basering is a Weyl algebra in characteristic 0 PURPOSE: check whether M is holonomic over the base ring NOTE: M is holonomic if 2*dim(M) = dim(R), where R is the base ring; dim stays for Gelfand-Kirillov dimension EXAMPLE: example isHolonomic; shows examples " { if (dmodappassumeViolation()) { ERROR("Basering is inappropriate: characteristic>0 or qring present"); } if (!isWeyl(basering)) { ERROR("Basering is not a Weyl algebra"); } if ( (typeof(M) != "ideal") && (typeof(M) != "module") && (typeof(M) != "matrix") ) { // print(typeof(M)); ERROR("an argument of type ideal/module/matrix expected"); } if (attrib(M,"isSB")!=1) { M = std(M); } int dimR = gkdim(std(0)); int dimM = gkdim(M); return( (dimR==2*dimM) ); } example { "EXAMPLE:"; echo = 2; ring R = 0,(x,y),dp; poly F = x*y*(x+y); def A = annfsBM(F,0); setring A; LD; isHolonomic(LD); ideal I = std(LD[1]); I; isHolonomic(I); } proc reiffen(int p, int q) "USAGE: reiffen(p, q); int p, int q RETURN: ring PURPOSE: set up the polynomial, describing a Reiffen curve NOTE: activate the output ring with the @code{setring} command and @* find the curve as a polynomial RC. @* A Reiffen curve is defined as RC = x^p + y^q + xy^{q-1}, q >= p+1 >= 5 EXAMPLE: example reiffen; shows examples " { // we allow also other numbers, p \geq 1, q\geq 1 // a Reiffen curve is defined as // F = x^p + y^q +x*y^{q-1}, q \geq p+1 \geq 5 // ASSUME: q >= p+1 >= 5. Otherwise an error message is returned // if ( (p<4) || (q<5) || (q-p<1) ) // { // ERROR("Some of conditions p>=4, q>=5 or q>=p+1 is not satisfied!"); // } if ( (p<1) || (q<1) ) { ERROR("Some of conditions p>=1, q>=1 is not satisfied!"); } ring @r = 0,(x,y),dp; poly RC = y^q +x^p + x*y^(q-1); export RC; return(@r); } example { "EXAMPLE:"; echo = 2; def r = reiffen(4,5); setring r; RC; } proc arrange(int p) "USAGE: arrange(p); int p RETURN: ring PURPOSE: set up the polynomial, describing a hyperplane arrangement NOTE: must be executed in a commutative ring ASSUME: basering is present and it is commutative EXAMPLE: example arrange; shows examples " { int UseBasering = 0 ; if (p<2) { ERROR("p>=2 is required!"); } if ( nameof(basering)!="basering" ) { if (p > nvars(basering)) { ERROR("too big p"); } else { def @r = basering; UseBasering = 1; } } else { // no basering found ERROR("no basering found!"); // ring @r = 0,(x(1..p)),dp; } int i,j,sI; poly q = 1; list ar; ideal tmp; tmp = ideal(var(1)); ar[1] = tmp; for (i = 2; i<=p; i++) { // add i-nary sums to the product ar = indAR(ar,i); } for (i = 1; i<=p; i++) { tmp = ar[i]; sI = size(tmp); for (j = 1; j<=sI; j++) { q = q*tmp[j]; } } if (UseBasering) { return(q); } // poly AR = q; export AR; // return(@r); } example { "EXAMPLE:"; echo = 2; ring X = 0,(x,y,z,t),dp; poly q = arrange(3); factorize(q,1); } proc checkRoot(poly F, number a, list #) "USAGE: checkRoot(f,alpha [,S,eng]); poly f, number alpha, string S, int eng RETURN: int ASSUME: Basering is a commutative ring, alpha is a rational number. PURPOSE: check whether a rational number alpha is a root of the global @* Bernstein-Sato polynomial of f and compute its multiplicity, @* with the algorithm given by S and with the Groebner basis engine given by eng. NOTE: The annihilator of f^s in D[s] is needed and hence it is computed with the @* algorithm by Briancon and Maisonobe. The value of a string S can be @* 'alg1' (default) - for the algorithm 1 of [LM08] @* 'alg2' - for the algorithm 2 of [LM08] @* Depending on the value of S, the output of type int is: @* 'alg1': 1 only if -alpha is a root of the global Bernstein-Sato polynomial @* 'alg2': the multiplicity of -alpha as a root of the global Bernstein-Sato @* polynomial of f. If -alpha is not a root, the output is 0. @* If eng <>0, @code{std} is used for Groebner basis 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 checkRoot; shows examples " { int eng = 0; int chs = 0; // choice string algo = "alg1"; string st; if ( size(#)>0 ) { if ( typeof(#[1]) == "string" ) { st = string(#[1]); if ( (st == "alg1") || (st == "ALG1") || (st == "Alg1") ||(st == "aLG1")) { algo = "alg1"; chs = 1; } if ( (st == "alg2") || (st == "ALG2") || (st == "Alg2") || (st == "aLG2")) { algo = "alg2"; chs = 1; } if (chs != 1) { // incorrect value of S print("Incorrect algorithm given, proceed with the default alg1"); algo = "alg1"; } // second arg if (size(#)>1) { // exists 2nd arg if ( typeof(#[2]) == "int" ) { // the case: given alg, given engine eng = int(#[2]); } else { eng = 0; } } else { // no second arg eng = 0; } } else { if ( typeof(#[1]) == "int" ) { // the case: default alg, engine eng = int(#[1]); // algo = "alg1"; //is already set } else { // incorr. 1st arg algo = "alg1"; } } } // size(#)=0, i.e. there is no algorithm and no engine given // eng = 0; algo = "alg1"; //are already set // int ppl = printlevel-voice+2; printlevel=printlevel+1; def save = basering; def @A = SannfsBM(F); setring @A; poly F = imap(save,F); number a = imap(save,a); if ( algo=="alg1") { int output = checkRoot1(LD,F,a,eng); } else { if ( algo=="alg2") { int output = checkRoot2(LD,F,a,eng); } } printlevel=printlevel-1; return(output); } example { "EXAMPLE:"; echo = 2; printlevel=0; ring r = 0,(x,y),Dp; poly F = x^4+y^5+x*y^4; checkRoot(F,11/20); // -11/20 is a root of bf poly G = x*y; checkRoot(G,1,"alg2"); // -1 is a root of bg with multiplicity 2 } proc checkRoot1(ideal I, poly F, number a, list #) "USAGE: checkRoot1(I,f,alpha [,eng]); ideal I, poly f, number alpha, int eng ASSUME: Basering is D[s], I is the annihilator of f^s in D[s], @* that is basering and I are the output of Sannfs-like procedure, @* f is a polynomial in K[x] and alpha is a rational number. RETURN: int, 1 if -alpha is a root of the Bernstein-Sato polynomial of f PURPOSE: check, whether alpha is a root of the global Bernstein-Sato polynomial of f NOTE: If eng <>0, @code{std} is used for Groebner basis 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 checkRoot1; shows examples " { // to check: alpha is rational (has char 0 check inside) if (!isRational(a)) { "ERROR: alpha must be a rational number!"; } // no qring if ( size(ideal(basering)) >0 ) { "ERROR: no qring is allowed"; } int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } int ppl = printlevel-voice+2; dbprint(ppl,"// -0-1- starting the procedure checkRoot1"); def save = basering; int N = nvars(basering); int Nnew = N-1; int n = Nnew / 2; int i; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether basering is D[s]=K(_x,_Dx,s) list Name = RL[2]; // for (i=1; i<=n; i++) // { // if ( bracket(var(i+n),var(i))!=1 ) // { // ERROR("basering should be D[s]=K(_x,_Dx,s)"); // } // } if ( Name[N]!="s" ) { ERROR("the last variable of basering should be s"); } // now, create the new vars list NName; for (i=1; i<=Nnew; i++) { NName[i] = Name[i]; } L[2] = NName; kill Name,NName; // block ord (dp); tmp[1] = "dp"; // string s = "iv="; for (i=1; i<=Nnew; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); kill s; tmp[2] = iv; Lord[1] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[Nnew][Nnew]; for (i=1; i<=n; i++) { @D[i,i+n]=1; } def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(_x,_Dx) is ready"); dbprint(ppl-1, S); // create the ideal K = ann_D[s](f^s)_{s=-alpha} + < f > setring save; ideal K = subst(I,s,-a); dbprint(ppl,"// -1-2- the variable s has been substituted by "+string(-a)); dbprint(ppl-1, K); K = NF(K,std(F)); // make leadcoeffs positive for (i=1; i<=ncols(K); i++) { if ( leadcoef(K[i])<0 ) { K[i] = -K[i]; } } K = K,F; // ------------ the ideal K is ready ------------ setring @R; ideal K = imap(save,K); dbprint(ppl,"// -1-3- starting the computation of a Groebner basis of K in @R"); dbprint(ppl-1, K); ideal G = engine(K,eng); dbprint(ppl,"// -1-4- the Groebner basis has been computed"); dbprint(ppl-1, G); return(G[1]!=1); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),Dp; poly F = x^4+y^5+x*y^4; printlevel = 0; def A = Sannfs(F); setring A; poly F = imap(r,F); checkRoot1(LD,F,31/20); // -31/20 is not a root of bs checkRoot1(LD,F,11/20); // -11/20 is a root of bs } proc checkRoot2 (ideal I, poly F, number a, list #) "USAGE: checkRoot2(I,f,a [,eng]); I an ideal, f a poly, alpha a number, eng an optional int ASSUME: I is the annihilator of f^s in D[s], basering is D[s], @* that is basering and I are the output os Sannfs-like procedure, @* f is a polynomial in K[_x] and alpha is a rational number. RETURN: int, the multiplicity of -alpha as a root of the BS polynomial of f. PURPOSE: check whether a rational number alpha is a root of the global Bernstein- @* Sato polynomial of f and compute its multiplicity from the known Ann F^s in D[s] NOTE: If -alpha is not a root, the output is 0. @* If eng <>0, @code{std} is used for Groebner basis 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 checkRoot2; shows examples " { // to check: alpha is rational (has char 0 check inside) if (!isRational(a)) { "ERROR: alpha must be a rational number!"; } // no qring if ( size(ideal(basering)) >0 ) { "ERROR: no qring is allowed"; } int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } int ppl = printlevel-voice+2; dbprint(ppl,"// -0-1- starting the procedure checkRoot2"); def save = basering; int N = nvars(basering); int n = (N-1) / 2; int i; string s; list RL = ringlist(basering); list L, Lord; list tmp; intvec iv; L[1] = RL[1]; // char L[4] = RL[4]; // char, minpoly // check whether basering is D[s]=K(_x,_Dx,s) list Name = RL[2]; for (i=1; i<=n; i++) { if ( bracket(var(i+n),var(i))!=1 ) { ERROR("basering should be D[s]=K(_x,_Dx,s)"); } } if ( Name[N]!="s" ) { ERROR("the last variable of basering should be s"); } // now, create the new vars L[2] = Name; kill Name; // block ord (dp); tmp[1] = "dp"; // string s = "iv="; for (i=1; i<=N; i++) { s = s+"1,"; } s[size(s)]=";"; execute(s); kill s; tmp[2] = iv; Lord[1] = tmp; tmp[1] = "C"; iv = 0; tmp[2] = iv; Lord[2] = tmp; tmp = 0; L[3] = Lord; // we are done with the list def @R@ = ring(L); setring @R@; matrix @D[N][N]; for (i=1; i<=n; i++) { @D[i,i+n]=1; } def @R = nc_algebra(1,@D); setring @R; kill @R@; dbprint(ppl,"// -1-1- the ring @R(_x,_Dx,s) is ready"); dbprint(ppl-1, @R); // now, continue with the algorithm ideal I = imap(save,I); poly F = imap(save,F); number a = imap(save,a); ideal II = NF(I,std(F)); // make leadcoeffs positive for (i=1; i<=ncols(II); i++) { if ( leadcoef(II[i])<0 ) { II[i] = -II[i]; } } ideal J,G; int m; // the output (multiplicity) dbprint(ppl,"// -2- starting the bucle"); for (i=0; i<=n; i++) // the multiplicity has to be <= n { // create the ideal Ji = ann_D[s](f^s) + < f, (s+alpha)^{i+1} > // (s+alpha)^i in Ji <==> -alpha is a root with multiplicity >= i J = II,F,(s+a)^(i+1); // ------------ the ideal Ji is ready ----------- dbprint(ppl,"// -2-"+string(i+1)+"-1- starting the computation of a Groebner basis of J"+string(i)+" in @R"); dbprint(ppl-1, J); G = engine(J,eng); dbprint(ppl,"// -2-"+string(i+1)+"-2- the Groebner basis has been computed"); dbprint(ppl-1, G); if ( NF((s+a)^i,G)==0 ) { dbprint(ppl,"// -2-"+string(i+1)+"-3- the number "+string(-a)+" has not multiplicity "+string(i+1)); m = i; break; } dbprint(ppl,"// -2-"+string(i+1)+"-3- the number "+string(-a)+" has multiplicity at least "+string(i+1)); } dbprint(ppl,"// -3- the bucle has finished"); return(m); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z),Dp; poly F = x*y*z; printlevel = 0; def A = Sannfs(F); setring A; poly F = imap(r,F); checkRoot2(LD,F,1); // -1 is a root of bs with multiplicity 3 checkRoot2(LD,F,1/3); // -1/3 is not a root } proc checkFactor(ideal I, poly F, poly q, list #) "USAGE: checkFactor(I,f,qs [,eng]); I an ideal, f a poly, qs a poly, eng an optional int ASSUME: checkFactor is called from the basering, created by Sannfs-like proc, @* that is, from the Weyl algebra in x1,..,xN,d1,..,dN tensored with K[s]. @* The ideal I is the annihilator of f^s in D[s], that is the ideal, computed @* by Sannfs-like procedure (usually called LD there). @* Moreover, f is a polynomial in K[x1,..,xN] and qs is a polynomial in K[s]. RETURN: int, 1 if qs is a factor of the global Bernstein polynomial of f and 0 otherwise PURPOSE: check whether a univariate polynomial qs is a factor of the @* Bernstein-Sato polynomial of f without explicit knowledge of the latter. NOTE: If eng <>0, @code{std} is used for Groebner basis 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 checkFactor; shows examples " { // ASSUME too complicated, cannot check it. int eng = 0; if ( size(#)>0 ) { if ( typeof(#[1]) == "int" ) { eng = int(#[1]); } } int ppl = printlevel-voice+2; def @R2 = basering; int N = nvars(@R2); int i; // we're in D_n[s], where the elim ord for s is set dbprint(ppl,"// -0-1- starting the procedure checkFactor"); dbprint(ppl,"// -1-1- the ring @R2(_x,_Dx,s) is ready"); dbprint(ppl-1, @R2); // create the ideal J = ann_D[s](f^s) + < f,q > ideal J = NF(I,std(F)); // make leadcoeffs positive for (i=1; i<=ncols(J); i++) { if ( leadcoef(J[i])<0 ) { J[i] = -J[i]; } } J = J,F,q; // ------------ the ideal J is ready ----------- dbprint(ppl,"// -1-2- starting the elimination of _x,_Dx in @R2"); dbprint(ppl-1, J); ideal G = engine(J,eng); ideal K = nselect(G,1..N-1); kill J,G; dbprint(ppl,"// -1-3- _x,_Dx are eliminated"); dbprint(ppl-1, K); //q is a factor of bs iff K = < q > //K = normalize(K); //q = normalize(q); //return( (K[1]==q) ); return( NF(K[1],std(q))==0 ); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y),Dp; poly F = x^4+y^5+x*y^4; printlevel = 0; def A = Sannfs(F); setring A; poly F = imap(r,F); checkFactor(LD,F,20*s+31); // -31/20 is not a root of bs checkFactor(LD,F,20*s+11); // -11/20 is a root of bs checkFactor(LD,F,(20*s+11)^2); // the multiplicity of -11/20 is 1 } proc varNum(string s) "USAGE: varNum(s); string s RETURN: int PURPOSE: returns the number of the variable with the name s @* among the variables of basering or 0 if there is no such variable EXAMPLE: example varNum; shows examples " { int i; for (i=1; i<= nvars(basering); i++) { if ( string(var(i)) == s ) { return(i); } } return(0); } example { "EXAMPLE:"; echo = 2; ring X = 0,(x,y1,t,z(0),z,tTa),dp; varNum("z"); varNum("t"); varNum("xyz"); } static proc indAR(list L, int n) "USAGE: indAR(L,n); list L, int n RETURN: list PURPOSE: computes arrangement inductively, using L and @* var(n) as the next variable ASSUME: L has a structure of an arrangement EXAMPLE: example indAR; shows examples " { if ( (n<2) || (n>nvars(basering)) ) { ERROR("incorrect n"); } int sl = size(L); list K; ideal tmp; poly @t = L[sl][1] + var(n); //1 elt K[sl+1] = ideal(@t); tmp = L[1]+var(n); K[1] = tmp; tmp = 0; int i,j,sI; ideal I; for(i=sl; i>=2; i--) { I = L[i-1]; sI = size(I); for(j=1; j<=sI; j++) { I[j] = I[j] + var(n); } tmp = L[i],I; K[i] = tmp; I = 0; tmp = 0; } kill I; kill tmp; return(K); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,y,z,t,v),dp; list L; L[1] = ideal(x); list K = indAR(L,2); K; list M = indAR(K,3); M; M = indAR(M,4); M; } proc isRational(number n) "USAGE: isRational(n); n number RETURN: int PURPOSE: determine whether n is a rational number, @* that is it does not contain parameters. ASSUME: ground field is of characteristic 0 EXAMPLE: example indAR; shows examples " { if (char(basering) != 0) { ERROR("The ground field must be of characteristic 0!"); } number dn = denominator(n); number nn = numerator(n); return( ((int(dn)==dn) && (int(nn)==nn)) ); } example { "EXAMPLE:"; echo = 2; ring r = (0,a),(x,y),dp; number n1 = 11/73; isRational(n1); number n2 = (11*a+3)/72; isRational(n2); } /// ****** EXAMPLES ************ /* //static proc exCheckGenericity() { LIB "control.lib"; ring r = (0,a,b,c),x,dp; poly p = (x-a)*(x-b)*(x-c); def A = annfsBM(p); setring A; ideal J = slimgb(LD); matrix T = lift(LD,J); T = normalize(T); genericity(T); // Ann =x^3*Dx+3*x^2*t*Dt+(-a-b-c)*x^2*Dx+(-2*a-2*b-2*c)*x*t*Dt+3*x^2+(a*b+a*c+b*c)*x*Dx+(a*b+a*c+b*c)*t*Dt+(-2*a-2*b-2*c)*x+(-a*b*c)*Dx+(a*b+a*c+b*c) // genericity: g = a2-ab-ac+b2-bc+c2 =0 // g = (a -(b+c)/2)^2 + (3/4)*(b-c)^2; // g ==0 <=> a=b=c // indeed, Ann = (x-a)^2*(x*Dx+3*t*Dt+(-a)*Dx+3) // -------------------------------------------- // BUT a direct computation shows // when a=b=c, // Ann = x*Dx+3*t*Dt+(-a)*Dx+3 } //static proc exOT_17() { // Oaku-Takayama, p.208 ring R = 0,(x,y),dp; poly F = x^3-y^2; // x^2+x*y+y^2; option(prot); option(mem); // option(redSB); def A = annfsOT(F,0); setring A; LD; gkdim(LD); // a holonomic check // poly F = x^3-y^2; // = x^7 - y^5; // x^3-y^4; // x^5 - y^4; } //static proc exOT_16() { // Oaku-Takayama, p.208 ring R = 0,(x),dp; poly F = x*(1-x); option(prot); option(mem); // option(redSB); def A = annfsOT(F,0); setring A; LD; gkdim(LD); // a holonomic check } //static proc ex_bcheck() { ring R = 0,(x,y),dp; poly F = x*y*(x+y); option(prot); option(mem); int eng = 0; // option(redSB); def A = annfsOT(F,eng); setring A; LD; } //static proc ex_bcheck2() { ring R = 0,(x,y),dp; poly F = x*y*(x+y); int eng = 0; def A = annfsBM(F,eng); setring A; LD; } //static proc ex_BMI() { // a hard example ring r = 0,(x,y),Dp; poly F1 = (x2-y3)*(x3-y2); poly F2 = (x2-y3)*(xy4+y5+x4); ideal F = F1,F2; def A = annfsBMI(F); setring A; LD; BS; } //static proc ex2_BMI() { // this example was believed to be intractable in 2005 by Gago-Vargas, Castro and Ucha ring r = 0,(x,y),Dp; option(prot); option(mem); ideal F = x2+y3,x3+y2; printlevel = 2; def A = annfsBMI(F); setring A; LD; BS; } //static proc ex_operatorBM() { ring r = 0,(x,y,z,w),Dp; poly F = x^3+y^3+z^2*w; printlevel = 0; def A = operatorBM(F); setring A; F; // the original polynomial itself LD; // generic annihilator LD0; // annihilator bs; // normalized Bernstein poly BS; // root and multiplicities of the Bernstein poly PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD reduce(PS*F-bs,LD); // check the property of PS } */