Changeset eaf8f8 in git

Timestamp:

Mar 9, 2009, 7:34:52 PM (15 years ago)

Author:

Viktor Levandovskyy <levandov@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

b22a72d3cf5f625b056ece9285f0442601763cf0

Parents:

e6d283f58657a6d67f7a592a1e62a50f9464d936

Message:

*levandov: further bugfixes, finer work with varnames, doc enhancements, assume procs etc


git-svn-id: file:///usr/local/Singular/svn/trunk@11529 2c84dea3-7e68-4137-9b89-c4e89433aadc

Location:

Singular/LIB

Files:

• bfun.lib (modified) (20 diffs)
• dmod.lib (modified) (42 diffs)
• dmodapp.lib (modified) (5 diffs)

Singular/LIB/bfun.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: bfun.lib,v 1.3 2009-03-06 20:32:28 levandov Exp $";
+version="$Id: bfun.lib,v 1.4 2009-03-09 18:34:51 levandov Exp $";
 category="Noncommutative";
 info="
@@ -56 +56 @@
 LIB "dmod.lib";     // for SannfsBFCT etc
 LIB "dmodapp.lib";  // for initialIdealW etc
-LIB "nctools.lib"; // for isWeyl etc
+LIB "nctools.lib";  // for isWeyl etc
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -149 +149 @@
 static proc safeVarName (string s)
 {
-  string cv = "," + charstr(basering) + "," + varstr(basering) + ",";
+  string S = "," + charstr(basering) + "," + varstr(basering) + ",";
   s = "," + s + ",";
-  while (find(cv,s) <> 0)
+  while (find(S,s) <> 0)
   {
     s[1] = "@";
@@ -903 +903 @@
 }
 
-proc pIntersectSyz (poly s, ideal II, list #)
+proc pIntersectSyz (poly s, ideal I, list #)
 "USAGE:  pIntersectSyz(f, I [,p,s,t]);  f poly, I ideal, p,t optial ints, p prime
 RETURN:  vector, coefficient vector of the monic polynomial
 PURPOSE: compute the intersection of an ideal I with the subalgebra K[f]
-ASSUME:  I is given as Groebner basis, basering is free of parameters.
+ASSUME:  I is given as Groebner basis.
 NOTE:    If the intersection is zero, this procedure might not terminate.
 @*       If p>0 is given, this proc computes the generator of the intersection in
 @*       char p first and then only searches for a generator of the obtained 
 @*       degree in the basering. Otherwise, it searched for all degrees by 
-@*       computing syzygies.
+@*       computing syzygies. Note that the basering must be free of parameters to
+@*       use this modular method.
 @*       If s<>0, @code{std} is used for Groebner basis computations in char 0,
 @*       otherwise, and by default, @code{slimgb} is used.
@@ -923 +924 @@
 {
   // assume I is given in Groebner basis
-  ideal I = II;
   if (attrib(I,"isSB") <> 1)
   {
@@ -962 +962 @@
   newNF = NF(s,I);
   NI[2] = newNF;
+  list l = ringlist(save);
+  if (typeof(l[1]) == "int") { l[1] = solveincharp; }
+  else // parameters in basering: problems with mappings from Q(a_i)->Z/p(a_i)
+  {
+    if (solveincharp <> 0)
+    {
+      print("// basering must be free of parameters to use modular method");
+      print("// solving in characteristic of basering only...");
+      solveincharp = 0;
+    }
+  }
   if (solveincharp<>0)
   {
-    list l = ringlist(save);
-    l[1] = solveincharp;
-    matrix l5 = l[5];
-    matrix l6 = l[6];
+    if (size(l)>4)
+    {
+      matrix l5 = l[5];
+      matrix l6 = l[6];
+    }
     def @Rp = ring(l);
     setring @Rp;
     list l = ringlist(@Rp);
-    l[5] = fetch(save,l5);
-    l[6] = fetch(save,l6);
+    setring save;
+    if (size(l)>4)
+    {
+      setring @Rp;
+      l[5] = fetch(save,l5);
+      l[6] = fetch(save,l6);
+    }
+    setring @Rp;
     def Rp = ring(l);
     setring Rp;
     kill @Rp;
-    dbprint(ppl+1,"// solving in ring ", Rp);
+    dbprint(ppl-1,"// solving in ring ", Rp);
     vector v;
     map phi = save,maxideal(1);
@@ -984 +1002 @@
   }
   i = 1;
-  dbprint(ppl+1,"// pIntersectSyz starts...");
-  dbprint(ppl+1,"// with ideal I=", I);
+  dbprint(ppl,"// pIntersectSyz starts...");
+  dbprint(ppl-1,"// with ideal I=", I);
   while (1)
   {
-    dbprint(ppl,"// i:"+string(i));
+    dbprint(ppl,"// testing degree: "+string(i));
     if (i>1)
     {
@@ -1004 +1022 @@
     }
     // look for a solution
-    dbprint(ppl,"// linSyzSolve starts with: "+string(matrix(NI)));
+    dbprint(ppl-1,"// linSyzSolve starts with: "+string(matrix(NI)));
     if (solveincharp<>0) // modular method
     {
@@ -1031 +1049 @@
     }
     matrix MM[1][nrows(v)] = matrix(v);
-    dbprint(ppl,"// linSyzSolve ready  with: "+string(MM));
+    dbprint(ppl-1,"// linSyzSolve ready  with: "+string(MM));
     kill MM;
     //  "linSyzSolve ready with"; print(v);
@@ -1057 +1075 @@
     i++;
   }
-  dbprint(ppl+1,"// pIntersectSyz finished");
+  dbprint(ppl,"// pIntersectSyz finished");
   return(v);
 }
@@ -1190 +1208 @@
     print("// basering is qring:");
     print("// discarding the quotient and proceeding...");
-    L[4] = 0;
+    L[4] = ideal(0);
     qr = 1;
-    def save2 = ring(L);
-    setring save2;
+    def save2 = ring(L); setring save2;
     poly f = imap(save,f);
   }
@@ -1214 +1231 @@
     ideal J = LD;
     kill LD;
+    poly s = var(2*n+1);
   }
   vector b;
@@ -1221 +1239 @@
     if (pIntersectchar == 0) // pIntersectSyz::modular
     {
-      int lb = 30000;
-      int ub = 536870909;
-      i = 1;
-      list usedprimes;
-      while (b == 0)
-      {
-        dbprint(ppl,"// number of run in the loop: "+string(i));
-        int q = prime(random(lb,ub));
-        if (findFirst(usedprimes,q)==0) // if q was not already used
-        {
-          usedprimes = usedprimes,q;
-          dbprint(ppl,"// used prime is: "+string(q));
-          b = pIntersectSyz(s,J,q,whichengine,modengine);
+      list L = ringlist(D);
+      if (typeof(L[1])=="int")
+      {  
+        int lb = 10000;
+        int ub = 536870909;
+        i = 1;
+        list usedprimes;
+        matrix L5 = L[5]; matrix L6 = L[6];
+        int sJ = size(J); int sJq;
+        while (b == 0)
+        {
+          dbprint(ppl,"// number of run in the loop: "+string(i));
+          int q = prime(random(lb,ub));
+          if (findFirst(usedprimes,q)==0) // if q was not already used
+          {
+            usedprimes = usedprimes,q;
+            dbprint(ppl,"// using prime: "+string(q));
+            L[1] = q;
+            def @Rq = ring(L); setring @Rq;
+            list lq = ringlist(@Rq);
+            lq[5] = fetch(D,L5); lq[6] = fetch(D,L6);
+            def Rq = ring(lq);
+            setring Rq; kill @Rq;
+            ideal Jq = fetch(D,J);
+            sJq = size(lead(Jq));
+            setring D; kill Rq;
+            if (sJq <> sJ)
+            {
+              dbprint(ppl,"// " +string(q) + " is unlucky");
+              b = 0;
+            }
+            else
+            {
+              b = pIntersectSyz(s,J,q,whichengine,modengine);
+            }
+          }
+          i++;
         }
-        i++;
-      }
-    }
-    else // pIntersectSyz::non-modular 
+      }
+      else //parameters: problem when mapping Q(a_i) -> Z/p(a_i)
+      {
+        print("// basering must be free of parameters to use modular method");
+        print("// solving only in char 0...");
+        b = pIntersectSyz(s,J,0,whichengine);
+      }
+    }  
+    else // pIntersectSyz::non-modular else // parameters in basering: problems with mappings
     {
       b = pIntersectSyz(s,J,0,whichengine);
@@ -1247 +1294 @@
     b = pIntersect(s,J);
   }
-  if (qr == 1) { setring save2; }
-  else { setring save; }
-  vector b = imap(D,b);
   if (inorann == 0) { list l = listofroots(b,1); }
   else              { list l = listofroots(b,0); }
-  if (qr == 1) 
-  {
-    setring save;
-    list l = imap(save2,l);
-  }
+  setring save;
+  list l = imap(D,l);
   return(l);
 }
@@ -1343 +1384 @@
 @*       If t<>0, the computation of the intersection is solely performed over
 @*       charasteristic 0, otherwise and by default, a modular method is used.
+@*       Note that in this case, the basering must be free of parameters.
 @*       If u<>0 and by default, @code{std} is used for GB computations in
 @*       characteristic >0, otherwise, @code{slimgb} is used.
@@ -1469 +1511 @@
   {
     print("// given ideal is not holonomic");
-    print("// setting bound for degree of b-fubction to 10 and proceeding");
+    print("// setting bound for degree of b-function to 10 and proceeding");
     isH = 10;
   }
@@ -1526 +1568 @@
 {
   int ppl = printlevel - voice +2;
-  if (isCommutative() == 0) { ERROR("basering must be commutative"); }
+  if (!isCommutative()) { ERROR("Basering must be commutative"); }
   def save = basering;
   int n = nvars(save);
@@ -1539 +1581 @@
     print("// basering is qring:");
     print("// discarding the quotient and proceeding...");
-    L[4] = 0;
+    L[4] = ideal(0);
     qr = 1;
     def save2 = ring(L);
@@ -1545 +1587 @@
     poly f = imap(save,f);
   }
-  int noofvars = 2*n+4;
+  int N = 2*n+4;
   int i;
   int whichengine = 0; // default
@@ -1563 +1605 @@
     }
   }
-  intvec uv;
-  uv[n+3] = 1;
-  ring r = 0,(x(1..n)),dp;
-  poly f = fetch(save,f);
-  uv[1] = -1; uv[noofvars] = 0;
-  //  for the ordering
-  intvec @a; @a = 1:noofvars; @a[2] = 2;
+  // creating the homogenized extended Weyl algebra
+  // create names for vars
+  list Lvar;
+  Lvar[1]   = safeVarName("t"); 
+  Lvar[2]   = safeVarName("s");
+  Lvar[n+3] = safeVarName("D"+Lvar[1]);
+  Lvar[N]   = safeVarName("h");
+  for (i=1; i<=n; i++)
+  {
+    Lvar[i+2]   = string(var(i));
+    Lvar[i+n+3] = safeVarName("D" + string(var(i)));
+  }
+  // create ordering
+  intvec uv = -1; uv[n+3] = 1; uv[N] = 0;
+  intvec @a = 1:N; @a[2] = 2;
   intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0;
+  list Lord; 
+  Lord[1] = list("a",@a); Lord[2] = list("a",@a2);
   if (methodord == 0) // default: block ordering
   {
-    ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(noofvars-1),lp(1));
+    //ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(N-1),lp(1));
+    Lord[3] = list("a",uv);
+    Lord[4] = list("dp",intvec(1:(N-1)));
+    Lord[5] = list("lp",intvec(1));
+    Lord[6] = list("C",intvec(0));
   }
   else // M() ordering
   {
-    intmat @Ord[noofvars][noofvars];
-    @Ord[1,1..noofvars] = uv;
-    @Ord[2,1..noofvars] = 1:(noofvars-1);
-    for (i=1; i<=noofvars-2; i++)
-    {
-      @Ord[2+i,noofvars - i] = -1;
-    }
+    intmat @Ord[N][N];
+    @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1);
+    for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; }
     dbprint(ppl,"// weights for ordering: "+string(transpose(@a)));
     dbprint(ppl,"// the ordering matrix:",@Ord);
-    ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord));
-  }
-  dbprint(ppl,"// the ring Dh:",Dh);
-  // non-commutative relations
-  matrix @relD[noofvars][noofvars];
-  @relD[1,2] = t*h^2;// s*t = t*s+t*h^2
-  @relD[2,n+3] = Dt*h^2;// Dt*s = s*Dt+h^2
-  @relD[1,n+3] = h^2;
+    //ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord));
+    Lord[3] = list("M",intvec(@Ord));
+    Lord[4] = list("C",intvec(0));
+  }
+  // create commutative ring
+  list L@Dh = ringlist(basering);
+  L@Dh = L@Dh[1..4]; // if basering is commutative nc_algebra
+  L@Dh[2] = Lvar; L@Dh[3] = Lord;
+  def @Dh = ring(L@Dh); setring @Dh;
+  dbprint(ppl,"// the ring @Dh:",@Dh);
+  // create non-commutative relations
+  matrix @relD[N][N];
+  @relD[1,2] = var(1)*var(N)^2;     //  s*t = t*s + t*h^2
+  @relD[2,n+3] = var(n+3)*var(N)^2; // Dt*s = s*Dt+Dt*h^2
+  @relD[1,n+3] = var(N)^2;
   for (i=1; i<=n; i++)
   {
-    @relD[i+2,n+3+i] = h^2;
+    @relD[i+2,n+3+i] = var(N)^2;
   }
   dbprint(ppl,"// nc relations:",@relD);
   def Dh = nc_algebra(1,@relD);
-  setring Dh;
-  dbprint(ppl,"// computing in ring",DDh);
-  ideal I;
-  poly f = imap(r,f);
-  kill r;
+  setring Dh; kill @Dh;
+  dbprint(ppl,"// computing in ring",Dh);
+  poly f = imap(save,f);
   f = homog(f,h);
-  I = t - f, t*Dt - s;  // defining the Malgrange ideal
+  // create the Malgrange ideal
+  ideal I = var(1) - f, var(1)*var(n+3) - var(2);
   for (i=1; i<=n; i++)
   {
-    I = I, D(i)+diff(f,x(i))*Dt;
-  }
+    I[3+i] = var(i+n+3)+diff(f,var(i+2))*var(n+3);
+  }
+  dbprint(ppl-1, "// the Malgrange ideal: " +string(I));
+  // the hard part: Groebner basis computation
   dbprint(ppl, "// starting Groebner basis computation with engine: "+string(whichengine));
   I = engine(I, whichengine);
   dbprint(ppl, "// finished Groebner basis computation");
-  dbprint(ppl, "// I before dehomogenization: "+string(I));
   I = subst(I,h,1); // dehomogenization
-  dbprint(ppl, "// I after dehomogenization: " +string(I));
-  I = inForm(I,uv); // we are only interested in the initial form w.r.t. uv
+  dbprint(ppl-1,string(I));
+  // 3.3 the initial form
+  I = inForm(I,uv);
   dbprint(ppl, "// the initial ideal:", string(matrix(I)));
+  // read off the solution
   intvec tonselect = 1;
   for (i=3; i<=2*n+4; i++) { tonselect = tonselect,i; }

Singular/LIB/dmod.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: dmod.lib,v 1.36 2009-03-06 21:05:54 levandov Exp $";
+version="$Id: dmod.lib,v 1.37 2009-03-09 18:34:51 levandov Exp $";
 category="Noncommutative";
 info="
@@ -28 +28 @@
 @*     PS*F^(s+1) = bs*F^s holds in K[x_1,...,x_n,1/F]*F^s.
 
-@* We provide the following implementations:
-@* OT) the classical Ann F^s algorithm from Oaku and Takayama (Journal of 
+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
+@* (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:
@@ -78 +84 @@
 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
+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
@@ -97 +103 @@
 LIB "control.lib";  // for genericity
 LIB "dmodapp.lib"; // for e.g. engine
+LIB "poly.lib";
 
 proc testdmodlib()
@@ -136 +143 @@
 "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.
-@*       The value of a string S can be
-@*       '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 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 (which is a Groebner basis) is the needed D-module structure,
-@*       - the list  BS is the list of roots and multiplicities of a Bernstein polynomial of f.
+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.
+@*    The value of a string S can be
+@*    '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.
@@ -259 +267 @@
 "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
+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 Oaku and Takayama,
+@*       '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,
@@ -377 +386 @@
 "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, where D is the Weyl algebra
+PURPOSE: compute the D-module structure of basering[1/f]*f^s
 NOTE:    activate this ring with the @code{setring} command.
-@*       In the ring D[s], the ideal LD1 is generated by the elements in Ann F^s in D[s] coming from logarithmic derivations.
+@*   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.
@@ -526 +536 @@
 "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 this ring with the @code{setring} command. In this 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 this 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,
@@ -716 +727 @@
 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 its multiplicies
-PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, according to the algorithm by Briancon and Maisonobe
+RETURN:  list (of roots of the Bernstein polynomial b and its 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.
@@ -942 +954 @@
 RETURN:  ring
 PURPOSE: compute the D-module structure of basering[1/f]*f^s, according
-to the algorithm by Briancon and Maisonobe
+@* to the algorithm by Briancon and Maisonobe
 NOTE:    activate this ring with the @code{setring} command. In this ring,
 @*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
@@ -1239 +1251 @@
 "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
-NOTE:    activate this 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.
-@*       annfs2 uses the shorter form of expressions in the variable s (the idea of Noro).
+PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, 
+@*           based on the output of Sannfs-like procedure
+NOTE: Activate this 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.
+@* annfs2 uses shorter expressions in the variable s (the idea of Noro).
 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
 @*       if @code{printlevel}>=2, all the debug messages will be printed.
@@ -1400 +1412 @@
 "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
+PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, 
+@* based on the output of Sannfs like procedure
 NOTE:    activate this ring with the @code{setring} command. In this ring,
 @*       - the ideal LD (which is a Groebner basis) is the annihilator of f^s,
@@ -1589 +1601 @@
 "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, according to the algorithm by Briancon and Maisonobe
+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 this ring with the @code{setring} command. In this ring D[s]
 @*       - the polynomial F is the same as the input,
@@ -1982 +1995 @@
 "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, according to the algorithm by Briancon and Maisonobe
+PURPOSE: compute the generic Ann F^s and exceptional parametric constellations
+@* of a polynomial with parametric coefficients with the BM algorithm
 NOTE:    activate this ring with the @code{setring} command. In this ring,
 @*       - the ideal LD is the D-module structure oa Ann F^s
@@ -2224 +2238 @@
 "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.
+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 this ring with the @code{setring} command. In this ring,
 @*       - the ideal LD is the needed D-mod structure,
@@ -2553 +2567 @@
 "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
+PURPOSE: compute the D-module structure of basering[1/f]*f^s, 
+@* according to the algorithm by Oaku and Takayama
 NOTE:    activate this ring with the @code{setring} command. In this ring,
 @*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
@@ -2935 +2949 @@
 "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], where D is the Weyl algebra
+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 this 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.
+@*  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
 "
@@ -3217 +3233 @@
 "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], where D is the Weyl algebra
-NOTE:    activate this 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.
+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 this 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
 "
@@ -3422 +3440 @@
 "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]
+PURPOSE: compute Ann f^s and Groebner basis of Ann f^s+f in D[s]
 NOTE:    activate this 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.
-@*       If eng <>0, @code{std} is used for Groebner basis computations,
-@*       otherwise, and by default @code{slimgb} is used.
+@*  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.
+@*  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 monomials 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 this 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 SannfsBFCT; shows examples
+EXAMPLE: example SannfsBFCTstd; 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
+  // DEBUG INFO: ordering on the output ring = dp, 
+  // use std(K,F); for reusing the std property of K
 
   int eng = 0;
@@ -3599 +3823 @@
   // total cleanup
   poly F = imap(save, F);
-  ideal LD = K,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-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);
@@ -3633 +3858 @@
   def B = SannfsBFCT(RC);
   setring B;
-  // are there polynomials, depending on s only?
+  // 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);
-}
-
-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 this 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?
+  // It is not the case. Are there leading monomials in s only?
   nselect(lead(LD),v);
 }
@@ -3864 +3869 @@
 "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
+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 this 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.
+@*    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
 "
@@ -4319 +4326 @@
 "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
+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 this 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 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.
+@*  - 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
 "
@@ -4363 +4371 @@
 "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
+PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based 
+@* on the output of Sannfs-like procedure
 NOTE:    activate this 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.
-@*       If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
+@* - 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
 "
@@ -4807 +4815 @@
 RETURN:  list
 PURPOSE: convert a ringlist L into another ringlist,
-where all the 'p' orderings are replaced with the 's' orderings.
+@* where all the 'p' orderings are replaced with the 's' orderings.
 ASSUME: L is a result of a ringlist command
 EXAMPLE: example convloc; shows examples
@@ -4852 +4860 @@
 "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 a rational number n
-ASSUME:  the basering contains 's' as a variable
-NOTE:    We assume that the basering is D[s],
-@*          ideal I is the Ann F^s in D[s] (obtained with e.g. SannfsBM, SannfsLOT, SannfsOT)
-@*          integer 'mir' is the minimal integer root of the Bernstein polynomial of F
-@*          and the number n is rational.
-@*       We compute the real annihilator for any rational value of n (both generic and exceptional).
-@*       The implementation goes along the lines of Saito-Sturmfels-Takayama, Alg. 5.3.15
+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.
@@ -4865 +4874 @@
 "
 {
+
+  if (!isRational(n))
+  {
+    "ERROR: n must be a rational number!";
+  }
   int ppl = printlevel-voice+2;
   //  int mir =  minIntRoot(L[1],0);
@@ -5023 +5037 @@
 "USAGE:  isHolonomic(M); M an ideal/module/matrix
 RETURN:  int, 1 if M is holonomic and 0 otherwise
-ASSUME: basering is a Weyl algebra
+ASSUME: basering is a Weyl algebra in characteristic 0
 PURPOSE: check the modules for the property of holonomy
 NOTE:    M is holonomic, if 2*dim(M) = dim(R), where R is a
@@ -5030 +5044 @@
 "
 {
+  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") )
   {
@@ -5061 +5084 @@
 RETURN:  ring
 PURPOSE: set up the polynomial, describing a Reiffen curve
-NOTE:    activate this ring with the @code{setring} command and find the
-         curve as a polynomial RC
-@*   a Reiffen curve is defined as F = x^p + y^q + xy^{q-1}, q >= p+1 >= 5
-ASSUME: q >= p+1 >= 5. Otherwise an error message is returned
+NOTE:  activate this ring with the @code{setring} command and 
+@*       find the curve as a polynomial RC.
+@*  A Reiffen curve is defined as F = 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
 
-  if ( (p<4) || (q<5) || (q-p<1) )
-  {
-    ERROR("Some of conditions p>=4, q>=5 or q>=p+1 is not satisfied!");
+// 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;
@@ -5092 +5122 @@
 RETURN:  ring
 PURPOSE: set up the polynomial, describing a hyperplane arrangement
-NOTE:   must be executed in a ring
-ASSUME: basering is present
+NOTE:  must be executed in a commutative ring
+ASSUME: basering is present and it is commutative
 EXAMPLE: example arrange; shows examples
 "
@@ -5155 +5185 @@
 
 proc checkRoot(poly F, number a, list #)
-"USAGE:  checkRoot(f,alpha [,S,eng]);  f a poly, alpha a number, S a string , eng an optional int
+"USAGE:  checkRoot(f,alpha [,S,eng]);  poly f, number alpha, string S, int eng
 RETURN:  int
-PURPOSE: check whether a rational is a root of the global Bernstein polynomial of f (and compute its multiplicity)
-with the algorithm given in S and with the Groebner basis engine given in eng
-NOTE:    The annihilator of f^s in D[s] is needed and it is computed according to the algorithm by Briancon and Maisonobe
-@*       The value of a string S can be
-@*       'alg1' (default) - for the algorithm 1 of J. Martin-Morales (unpublished)
-@*       'alg2' - for the algorithm 2 of J. Martin-Morales (unpublished)
+ASSUME: Basering is 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]
 @*       The output of type int is:
-@*       'alg1': 1 if -alpha is a root of the global Bernstein polynomial and 0 otherwise
-@*       'alg2':  the multiplicity of -alpha as a root of the global Bernstein polynomial of f.
-@*                (If -alpha is not a root, the output is 0)
+@*       '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.
@@ -5196 +5228 @@
      {
        // incorrect value of S
-       print("Incorrect algorithm given, proceed with the default alg1 of J. Martín-Morales");
+       print("Incorrect algorithm given, proceed with the default alg1");
        algo = "alg1";
      }
@@ -5269 +5301 @@
 
 proc checkRoot1(ideal I, poly F, number a, list #)
-"USAGE:  checkRoot1(I,f,alpha [,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], f is a polynomial in K[_x]
-RETURN:  int, 1 if -alpha is a root of the global Bernstein polynomial of f and 0 otherwise
-PURPOSE: check whether a rational is a root of the global Bernstein polynomial of f
-NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
+"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,
@@ -5280 +5314 @@
 "
 {
+  // 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 )
@@ -5393 +5437 @@
 
 proc checkRoot2 (ideal I, poly F, number a, list #)
-"USAGE:  checkRoot2(I,f,alpha [,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], f is a polynomial in K[_x]
-RETURN:  int, the multiplicity of -alpha as a root of the global Bernstein polynomial of f. If -alpha is not a root, the output is 0
-PURPOSE: 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]
-NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
+"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,
@@ -5404 +5452 @@
 "
 {
+
+  
+  // 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 )
@@ -5525 +5586 @@
 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:  I is the output of Sannfs, SannfsBM, SannfsLOT or SannfsOT,
-         f is a polynomial in K[_x], qs is a polynomial in K[s]
+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 rational is a root of the global Bernstein polynomial of f and compute its multiplicity from the known Ann F^s in D[s]
+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.
@@ -5536 +5601 @@
 "
 {
+
+  // ASSUME too complicated, cannot check it.
+
   int eng = 0;
   if ( size(#)>0 )
@@ -5595 +5663 @@
 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
+@*  among the variables of basering or 0 if there is no such variable
 EXAMPLE: example varnum; shows examples
 "
@@ -5621 +5689 @@
 "USAGE:  indAR(L,n);  list L, int n
 RETURN:  list
-PURPOSE: computes arrangement inductively, using L and var(n) as the
-next variable
+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
@@ -5666 +5734 @@
   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);
 }
 

Singular/LIB/dmodapp.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: dmodapp.lib,v 1.21 2009-03-06 20:32:29 levandov Exp $";
+version="$Id: dmodapp.lib,v 1.22 2009-03-09 18:34:52 levandov Exp $";
 category="Noncommutative";
 info="
@@ -1522 +1522 @@
     }
   }
+  
   // 1. create  homogenized Weyl algebra
   // 1.1 create ordering
@@ -1556 +1557 @@
   def Dh = nc_algebra(1,@relD);
   setring Dh; kill @Dh;
-  dbprint(ppl,"// computing in ring",DDh);
+  dbprint(ppl,"// computing in ring",Dh);
   // 2. Compute the initial ideal
   ideal I = imap(save,I);
@@ -1653 +1654 @@
     u0 = 1:n;
   }
-  if (isCommutative() == 0) { ERROR("basering must be commutative"); }
   if (char(save) <> 0)      { ERROR("characteristic of basering has to be 0"); }
   // creating the homogenized extended Weyl algebra
   list RL = ringlist(save);
+  RL = RL[1..4]; // if basering is commutative nc_algebra
   list C0 = list("C",intvec(0));
   // 1. get the weighted degree of f
@@ -1779 +1780 @@
 }
 
-static proc dmodappassumeViolation()
+proc dmodappassumeViolation()
 {
   // returns Boolean : yes/no [for assume violation]
   // char K = 0
   // no qring
-  // input poly/ideal is nonzero ?
-  if ( (char(basering) !=0 ) || (nvars(basering) != gkdim(std(0)) ) )
-  {
+  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
+  {
+    //    "ERROR: no qring is allowed";
     return(1);
   }