Changeset dfb2c6 in git

Timestamp:

Mar 6, 2009, 9:32:29 PM (14 years ago)

Author:

Viktor Levandovskyy <levandov@…>

Branches:

(u'spielwiese', '5d369c3cbad1a1bf2d5c856a48fb8a30b51cec3b')

Children:

fda698647a5a06bd263906206d3ce89a87b42c0d

Parents:

2a5ce368e94b24f2bbac34de08062477cb206040

Message:

*levandov: bfun and dmodapp seriously reworked


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

Location:

Singular/LIB

Files:

• bfun.lib (modified) (64 diffs)
• dmodapp.lib (modified) (34 diffs)

Singular/LIB/bfun.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: bfun.lib,v 1.2 2009-02-23 15:10:18 Singular Exp $";
+version="$Id: bfun.lib,v 1.3 2009-03-06 20:32:28 levandov Exp $";
 category="Noncommutative";
 info="
@@ -9 +9 @@
 THEORY: Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R,
 @*      one is interested in the global b-Function (also known as Bernstein-Sato
-@*      polynomial) b(s) in K[s], defined to be the monic polynomial, satisfying a functional
-@*      identity L * F^{s+1} = b(s) F^s,   for some operator L in D[s]. Here, D stands for an
-@*      n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>.
+@*      polynomial) b(s) in K[s], defined to be the monic polynomial, satisfying
+@*      a functional identity    L * F^{s+1} = b(s) F^s,   
+@*      for some operator L in D[s]. Here, D stands for an n-th Weyl algebra
+@*      K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>.
 @*      One is interested in the following data:
 @*      - Bernstein-Sato polynomial b(s) in K[s],
 @*      - the list of its roots, which are known to be rational
 @*      - the multiplicities of the roots.
-@*      References: Saito, Sturmfels, Takayama: Groebner Deformations of Hypergeometric
-@*      Differential Equations (2000), Noro: An Efficient Modular Algorithm for Computing
-@*      the Global b-function, (2002).
+@*
+@*    There is a general definition of a b-function of a holonomic ideal [SST]
+@*    (that is, an ideal I in a Weyl algebra D, such that D/I is holonomic module)
+@*    with respect to the given weight vector. B-function b_w(s) is in K[s].
+@*    Unlike Bernstein-Sato polynomial, general B-function with respect to
+@*    arbitrary weights need not have rational roots at all.
+@* 
+@*      References: [SST] Saito, Sturmfels, Takayama: Groebner Deformations of 
+@*      Hypergeometric Differential Equations (2000),
+@*      Noro: An Efficient Modular Algorithm for Computing the Global b-function,
+@*      (2002).
 
 
@@ -27 +36 @@
 bfctAnn(f[,s]);             computes the Bernstein-Sato polynomial of poly f
 bfctOneGB(f[,s,t]);         computes the Bernstein-Sato polynomial of poly f
-bfctIdeal(I,w[,s,t]);       computes the global b-function of ideal I w.r.t. the weight w
-pIntersect(f,I);            intersection of the ideal I with the subalgebra K[f] for a poly f
-pIntersectSyz(f,I[,p,s,t]); intersection of the ideal I with the subalgebra K[f] for a poly f
+bfctIdeal(I,w[,s,t]);       computes the global b-function of ideal wrt weight
+pIntersect(f,I[,s]);        intersection of ideal with subalgebra K[f] for poly f
+pIntersectSyz(f,I[,p,s,t]); intersection of ideal with subalgebra K[f] for poly f
 linReduce(f,I[,s]);         reduces a poly by linear reductions only
 linReduceIdeal(I,[s,t]);    reduces generators of ideal by linear reductions only
-linSyzSolve(I[,s]);         computes a linear dependency of the elements of ideal I
+linSyzSolve(I[,s]);         computes a linear dependency of elements of ideal I
 
 AUXILIARY PROCEDURES:
@@ -47 +56 @@
 LIB "dmod.lib";     // for SannfsBFCT etc
 LIB "dmodapp.lib";  // for initialIdealW etc
-
+LIB "nctools.lib"; // for isWeyl etc
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -83 +92 @@
 "
 {
+  // todo check assumption
   int i;
   def save = basering;
@@ -137 +147 @@
 }
 
+static proc safeVarName (string s)
+{
+  string cv = "," + charstr(basering) + "," + varstr(basering) + ",";
+  s = "," + s + ",";
+  while (find(cv,s) <> 0)
+  {
+    s[1] = "@";
+    s = "," + s;
+  }
+  s = s[2..size(s)-1];
+  return(s)
+}
 
 proc allPositive (intvec v)
@@ -167 +189 @@
 static proc findFirst (list l, i)
 "USAGE:  findFirst(l,i);  l a list, i an argument of any type
-RETURN:  int, the position of the first appearance of i in l, or 0 if i is not a member of l
+RETURN:  int, the position of the first appearance of i in l,
+@*       or 0 if i is not a member of l
 PURPOSE: check whether the second argument is a member of a list
 EXAMPLE: example findFirst; shows an example
@@ -230 +253 @@
 "USAGE:  linReduceIdeal(I [,s,t,u]); I an ideal, s,t,u optional ints
 RETURN:  ideal or list, linear reductum (over field) of I by its elements
-PURPOSE: reduce a list of polys only by linear reductions (no monomial multiplications)
-NOTE:    If s<>0, a list consisting of the reduced ideal and the coefficient 
+PURPOSE: reduces a list of polys only by linear reductions (no monomial
+@*        multiplications)
+NOTE:    If s<>0, a list consisting of the reduced ideal and the coefficient
 @*       vectors of the used reductions given as module is returned.
 @*       Otherwise (and by default), only the reduced ideal is returned.
@@ -237 +261 @@
 @*       otherwise, only leading monomials are reduced.
 @*       If u<>0 (and by default), the ideal is first sorted in increasing order.
-@*       If u is set to 0 and the given ideal is not sorted in the way described, 
+@*       If u is set to 0 and the given ideal is not sorted in the way described,
 @*       the result might not be as expected.
 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
@@ -327 +351 @@
     }
   }  
-  dbprint(ppl-1,"initially sorted ideal:", I);
-  if (remembercoeffs <> 0) { dbprint(ppl-1," used permutations:", M); }
+  dbprint(ppl-1,"// initially sorted ideal:", I);
+  if (remembercoeffs <> 0) { dbprint(ppl-1,"// used permutations:", M); }
   // step 2: reduce leading monomials by linear reductions
   poly lm,c,redpoly,maxlmJ;
@@ -358 +382 @@
           if (lm == lmJ[j])
           {
-            dbprint(ppl-1,"reducing " + string(redpoly));
-            dbprint(ppl-1," with " + string(J[j]));
+            dbprint(ppl-1,"// reducing " + string(redpoly));
+            dbprint(ppl-1,"//  with " + string(J[j]));
             c = leadcoef(redpoly)/leadcoef(J[j]);
             redpoly = redpoly - c*J[j];
-            dbprint(ppl-1," to " + string(redpoly));
+            dbprint(ppl-1,"//  to " + string(redpoly));
             lm = leadmonom(redpoly);
             ordlm = ord(lm);
@@ -390 +414 @@
   if (redtail <> 0)
   {
-    dbprint(ppl,"finished reducing leading monomials");
+    dbprint(ppl,"// finished reducing leading monomials");
     dbprint(ppl-1,string(J));
-    if (remembercoeffs <> 0) { dbprint(ppl-1,"used reductions:" + string(M)); }
+    if (remembercoeffs <> 0) { dbprint(ppl-1,"// used reductions:" + string(M)); }
     int k;
     for (i=sZeros+1; i<=sI; i++)
@@ -406 +430 @@
             {
               c = leadcoef(J[j][k])/leadcoef(J[i]);
-              dbprint(ppl-1,"reducing " + string(J[j]));
-              dbprint(ppl-1," with " + string(J[i]));
+              dbprint(ppl-1,"// reducing " + string(J[j]));
+              dbprint(ppl-1,"//  with " + string(J[i]));
               J[j] = J[j] - c*J[i];
-              dbprint(ppl-1," to " + string(J[j]));
+              dbprint(ppl-1,"//  to " + string(J[j]));
               if (remembercoeffs <> 0) { M[j] = M[j] - c * M[i]; }
             }
@@ -442 +466 @@
 @*       If t<>0 (and by default) all monomials are reduced (if possible),
 @*       otherwise, only leading monomials are reduced.
-@*       If u<>0 (and by default), the ideal is linearly pre-reduced, i.e. instead
-@*       of the given ideal, the output of @code{linReduceIdeal} is used. If u is
-@*       set to 0 and the given ideal does not equal the output of 
-@*       @code{linReduceIdeal}, the result might not be as expected. 
+@*       If u<>0 (and by default), the ideal is linearly pre-reduced, i.e. 
+@*       instead of the given ideal, the output of @code{linReduceIdeal} is used.
+@*       If u is set to 0 and the given ideal does not equal the output of
+@*       @code{linReduceIdeal}, the result might not be as expected.
 DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
@@ -487 +511 @@
       I = sortedI[1];
       M = sortedI[2];
-      dbprint(ppl-1,"prepared ideal:",I);
-      dbprint(ppl-1," with operations:",M);
+      dbprint(ppl-1,"// prepared ideal:" +string(I));
+      dbprint(ppl-1,"//  with operations:" +string(M));
     }
     else { I = linReduceIdeal(I,0,redtail); }
@@ -532 +556 @@
   if (redtail <> 0)
   {
-    dbprint(ppl,"finished reducing leading monomials");
+    dbprint(ppl,"// finished reducing leading monomials");
     dbprint(ppl-1,string(f));
-    if (remembercoeffs <> 0) { dbprint(ppl-1,"used reductions:" + string(v)); }
+    if (remembercoeffs <> 0) { dbprint(ppl-1,"// used reductions:" + string(v)); }
     for (i=1; i<=sI; i++)
     {
-      dbprint(ppl,"testing ideal entry:",i);
+      dbprint(ppl,"// testing ideal entry "+string(i));
       for (j=1; j<=size(f); j++)
       {
@@ -546 +570 @@
             c = leadcoef(f[j])/lcI[i];
             f = f - c*I[i];
-            dbprint(ppl-1,"reducing with " + string(I[i]));
-            dbprint(ppl-1," to " + string(f));
+            dbprint(ppl-1,"// reducing with " + string(I[i]));
+            dbprint(ppl-1,"//  to " + string(f));
             if (remembercoeffs <> 0) { v = v - c * M[i]; }
           }
@@ -581 +605 @@
 proc linSyzSolve (ideal I, list #)
 "USAGE:  linSyzSolve(I[,s]);  I an ideal, s an optional int
-RETURN:  vector, coefficient vector of a linear combination of 0 in the elements of I
+RETURN:  vector, coefficient vector of linear combination of 0 in elements of I
 PURPOSE: compute a linear dependency between the elements of an ideal
 @*       if such one exists
@@ -613 +637 @@
     // ------- 1. introduce undefined coeffs ------------------
     def save = basering;
+    list RL = ringlist(save);
+    int nv = nvars(save);
+    int np = npars(save);
     int p = char(save);
+    string cs = "(" + charstr(save) + ")";
     if (enginespecified == 0)
     {
@@ -621 +649 @@
       }
     }
-    ring @A = p,(@a(1..sI)),lp;
-    ring @aA = (p,@a(1..sI)), (@z),dp;
-    // TODO: BUG! WHAT IF THE ORIGINAL RING ALREADY HAS SUCH VARIABLES/PARAMETERS!!!?
-    // TODO: YOU CAN OVERCOME THIS BY MEANS SIMILAR TO "chooseSafeVarName" FROM NEW "matrix.lib"
+    int i;
+    list Lvar;
+    for (i=1; i<=sI; i++)
+    {
+      Lvar[i] = safeVarName("@a(" + string(i) + ")");
+    }
+    list L@A = RL[1..4];
+    L@A[2] = Lvar;
+    L@A[3] = list(list("lp",intvec(1:sI)),list("C",intvec(0)));
+    def @A = ring(L@A);
+    L@A[2] = list(safeVarName("@z"));
+    L@A[3][1] = list("dp",intvec(1));
+    if (size(L@A[1])>1)
+    {
+      L@A[1][2] = L@A[1][2] + Lvar;
+      L@A[1][3][size(L@A[1][3])+1] = list("lp",intvec(1:sI));
+    }
+    else
+    {
+      L@A[1] = list(p,Lvar,list(list("lp",intvec(1:sI))),ideal(0));
+    }
+    def @aA = ring(L@A);
     def @B = save + @aA;
     setring @B;
     ideal I = imap(save,I);
     // ------- 2. form the linear system for the undef coeffs ---
-    int i;   poly W;  ideal QQ;
+     poly W;  ideal QQ;
     for (i=1; i<=sI; i++)
     {
-      W = W + @a(i)*I[i];
+      W = W + par(np+i)*I[i];
     }
     while (W!=0)
@@ -643 +689 @@
     ideal QQ = imap(@B,QQ);
     // ------- 3. this QQ is a polynomial ideal, so "solve" the system -----
-    dbprint(ppl, "linSyzSolve: starting Groebner basis computation with engine:", whichengine);
+    dbprint(ppl, "// linSyzSolve: starting Groebner basis computation with engine:", whichengine);
     QQ = engine(QQ,whichengine);
-    dbprint(ppl, "QQ after engine:", QQ);
+    dbprint(ppl, "// QQ after engine:", QQ);
     if (dim(QQ) == -1)
     {
-      dbprint(ppl+1, "no solutions by linSyzSolve");
+      dbprint(ppl+1, "// no solutions by linSyzSolve");
       // output zeroes
       setring save;
@@ -658 +704 @@
     module MQQ = std(module(QQ));
     AA = NF(AA,MQQ); // todo: we still receive NF warnings
-    dbprint(ppl, "AA after NF:",AA);
+    dbprint(ppl, "// AA after NF:",AA);
     //    "AA after NF:"; print(AA);
     for(i=1; i<=sI; i++)
@@ -664 +710 @@
       AA = subst(AA,var(i),1);
     }
-    dbprint(ppl, "AA after subst:",AA);
+    dbprint(ppl, "// AA after subst:",AA);
     //    "AA after subst: "; print(AA);
     vector v = (module(transpose(AA)))[1];
@@ -687 +733 @@
 RETURN:  vector, coefficient vector of the monic polynomial 
 PURPOSE: compute the intersection of ideal I with the subalgebra K[f]
-ASSUME:  I is given as Groebner basis.
+ASSUME:  I is given as Groebner basis, basering is not a qring.
 NOTE:    If the intersection is zero, this proc might not terminate.
 @*       If s>0 is given, it is searched for the generator of the intersection
@@ -696 +742 @@
 "
 {
+  def save = basering;
+  int n = nvars(save);
+  list RL = ringlist(save);
+  int i,j,k;
+  if (RL[4]<>0)
+  {
+    ERROR ("basering must not be a qring");
+  }
   // assume I is given in Groebner basis
   if (attrib(I,"isSB") <> 1)
   {
-    print("WARNING: The input has no SB attribute!");
-    print(" Treating it as if it were a Groebner basis and proceeding...");
+    print("// WARNING: The input has no SB attribute!");
+    print("// Treating it as if it were a Groebner basis and proceeding...");
     attrib(I,"isSB",1); // set attribute for suppressing NF messages
   }
@@ -711 +765 @@
     }
   }
-  int ppl = printlevel-voice+1;
+  int ppl = printlevel-voice+2;
   // ---case 1: I = basering---
   if (size(I) == 1)
@@ -717 +771 @@
     if (simplify(I,2)[1] == 1)
     {
-      return(gen(2)); // = s
-    }
-  }
-  def save = basering;
-  int n = nvars(save);
-  int i,j,k;
+      return(gen(1)); // = 1
+    }
+  }
   // ---case 2: intersection is zero---
   intvec degs = leadexp(s);
@@ -752 +803 @@
     return(vector(0));
   }  
-  dbprint(ppl,"a lower bound for the degree of the insection is:");
-  dbprint(ppl,degbound);
+  dbprint(ppl,"// lower bound for the degree of the insection is " +string(degbound));
   if (degbound == 0) // lm(s) does not appear in lm(I)
   {
@@ -766 +816 @@
   poly toNF,tobracket,newNF,rednewNF,oldNF,secNF;
   i = 1;
-  dbprint(ppl+1,"pIntersect starts...");
   while (1)
   {
     if (bound>0 && i>bound) { return(vector(0)); }
-    dbprint(ppl,"testing degree: "+string(i));
+    dbprint(ppl,"// Testing degree "+string(i));
     if (i>1)
     {
@@ -793 +842 @@
     rednewNF = ll[1];
     l[i+1] = ll[2];
-    dbprint(ppl,"newNF is:", newNF);
-    dbprint(ppl,"rednewNF is:", rednewNF);
+    dbprint(ppl-1,"// newNF is: "+string(newNF));
+    dbprint(ppl-1,"// rednewNF is: "+string(rednewNF));
     if (rednewNF != 0) // no linear dependency
     {
@@ -802 +851 @@
     else // there is a linear dependency, hence we are done
     {
-      dbprint(ppl+1,"the degree of the generator of the intersection is:", i);
+      dbprint(ppl,"// degree of the generator of the intersection is: "+string(i));
       break;
     }
   }
-  dbprint(ppl,"used linear reductions:", l);
-  // we obtain the coefficients of the generator of the intersection by the used reductions:
-  ring @R = 0,(a(1..i+1)),dp;
-  setring @R;
+  dbprint(ppl-1,"// used linear reductions:", l);
+  // we obtain the coefficients of the generator by the used reductions:
+  list Lvar;
+  for (j=1; j<=i+1; j++)
+  {
+    Lvar[j] = safeVarName("a("+string(j)+")");
+  }
+  list Lord = list("dp",intvec(1:(i+1))),list("C",intvec(0));
+  list L@R = RL[1..4];
+  L@R[2] = Lvar; L@R[3] = Lord;
+  def @R = ring(L@R); setring @R;
   list l = imap(save,l);
   ideal C;
@@ -817 +873 @@
     for (k=1;k<=j;k++)
     {
-      C[j] = C[j]+l[j][k]*a(k);
+      C[j] = C[j]+l[j][k]*var(k);
     }
   }
   for (j=i;j>=1;j--)
   {
-    C[i+1]= subst(C[i+1],a(j),a(j)+C[j]);
+    C[i+1]= subst(C[i+1],var(j),var(j)+C[j]);
   }
   matrix m = coeffs(C[i+1],maxideal(1));
@@ -833 +889 @@
   v = imap(@R,v);
   kill @R;
-  dbprint(ppl+1,"pIntersect finished");
   return(v);
 }
@@ -849 +904 @@
 
 proc pIntersectSyz (poly s, ideal II, list #)
-"USAGE:  pIntersectSyz(f, I [,p,s,t]);  f a poly, I an ideal, p, t optial ints, p a prime number
+"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.
+ASSUME:  I is given as Groebner basis, basering is free of parameters.
 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.
+@*       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.
 @*       If s<>0, @code{std} is used for Groebner basis computations in char 0,
 @*       otherwise, and by default, @code{slimgb} is used.
-@*       If t<>0 and by default, @code{std} is used for Groebner basis computations in char >0,
-@*       otherwise, @code{slimgb} is used.
+@*       If t<>0 and by default, @code{std} is used for Groebner basis 
+@*       computations in char >0, otherwise, @code{slimgb} is used.
 DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
@@ -870 +926 @@
   if (attrib(I,"isSB") <> 1)
   {
-    print("WARNING: The input has no SB attribute!");
-    print(" Treating it as if it were a Groebner basis and proceeding...");
+    print("// WARNING: The input has no SB attribute!");
+    print("//  Treating it as if it were a Groebner basis and proceeding...");
     attrib(I,"isSB",1); // set attribute for suppressing NF messages
   }
@@ -920 +976 @@
     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);
@@ -928 +984 @@
   }
   i = 1;
-  dbprint(ppl+1,"pIntersectSyz starts...");
-  dbprint(ppl+1,"with ideal I=", I);
+  dbprint(ppl+1,"// pIntersectSyz starts...");
+  dbprint(ppl+1,"// with ideal I=", I);
   while (1)
   {
-    dbprint(ppl,"i:"+string(i));
+    dbprint(ppl,"// i:"+string(i));
     if (i>1)
     {
@@ -948 +1004 @@
     }
     // look for a solution
-    dbprint(ppl,"linSyzSolve starts with: "+string(matrix(NI)));
+    dbprint(ppl,"// linSyzSolve starts with: "+string(matrix(NI)));
     if (solveincharp<>0) // modular method
     {
@@ -956 +1012 @@
       if (v!=0) // there is a modular solution
       {
-        dbprint(ppl,"got solution in char ",solveincharp," of degree " ,i);
+        dbprint(ppl,"// got solution in char "+string(solveincharp)+" of degree "+string(i));
         setring save;
         v = linSyzSolve(NI,whichengine);
@@ -975 +1031 @@
     }
     matrix MM[1][nrows(v)] = matrix(v);
-    dbprint(ppl,"linSyzSolve ready  with: "+string(MM));
+    dbprint(ppl,"// linSyzSolve ready  with: "+string(MM));
     kill MM;
     //  "linSyzSolve ready with"; print(v);
@@ -989 +1045 @@
       if (p!=0)
       {
-        dbprint(ppl,"linSyzSolve: bad solution!");
+        dbprint(ppl,"// linSyzSolve: bad solution!");
       }
       else
       {
-        dbprint(ppl,"linSyzSolve: got solution!");
+        dbprint(ppl,"// linSyzSolve: got solution!");
         // "got solution!";
         break;
@@ -1001 +1057 @@
     i++;
   }
-  dbprint(ppl+1,"pIntersectSyz finished");
+  dbprint(ppl+1,"// pIntersectSyz finished");
   return(v);
 }
@@ -1088 +1144 @@
   }
   int substitution = int(#[2]);
-  ring S = 0,s,dp;
+  string s = safeVarName("s");
+  list RL = ringlist(save); RL = RL[1..4];
+  RL[2] = list(s); RL[3] = list(list("dp",intvec(1)),list("C",0));
+  def S = ring(RL); setring S;
   ideal J;
   for (i=1; i<=n; i++)
   {
-    J[i] = s;
+    J[i] = var(1);
   }
   map @m = save,J;
@@ -1098 +1157 @@
   if (substitution == 1)
   {
-    p = subst(p,s,-s-1);
+    p = subst(p,var(1),-var(1)-1);
   }
   // the rest of this proc is nicked from bernsteinBM from dmod.lib
   list P = factorize(p);//with constants and multiplicities
-  ideal bs; intvec m;   //the Bernstein polynomial is monic, so we are not interested in constants
+  ideal bs; intvec m;   //the BS poly is monic, so we are not interested in constants
   for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
   {
@@ -1109 +1168 @@
   }
   bs =  normalize(bs);
-  bs = -subst(bs,s,0);
+  bs = -subst(bs,var(1),0);
   setring save;
   ideal bs = imap(S,bs);
@@ -1123 +1182 @@
   def save = basering;
   int n = nvars(save);
+  if (isCommutative() == 0) { ERROR("basering must be commutative"); }
+  if (char(save) <> 0) { ERROR("characteristic of basering has to be 0"); }
+  list L = ringlist(save);
+  int qr;
+  if (L[4] <> 0) // qring
+  {
+    print("// basering is qring:");
+    print("// discarding the quotient and proceeding...");
+    L[4] = 0;
+    qr = 1;
+    def save2 = ring(L);
+    setring save2;
+    poly f = imap(save,f);
+  }
   if (size(variables(f)) == 0) // f is constant
   {
@@ -1129 +1202 @@
   if (inorann == 0) // bfct using initial ideal
   {
-    def D = initialMalgrange(f,whichengine,methodord,1,u0);
+    def D = initialMalgrange(f,whichengine,methodord,u0);
     setring D;
     ideal J = inF;
@@ -1154 +1227 @@
       while (b == 0)
       {
-        dbprint(ppl,"number of run in the loop: "+string(i));
+        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));
+          dbprint(ppl,"// used prime is: "+string(q));
           b = pIntersectSyz(s,J,q,whichengine,modengine);
         }
@@ -1174 +1247 @@
     b = pIntersect(s,J);
   }
-  setring save;
+  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 (inorann == 0) { list l = listofroots(b,1); }
+  else              { list l = listofroots(b,0); }
+  if (qr == 1) 
+  {
+    setring save;
+    list l = imap(save2,l);
   }
   return(l);
@@ -1192 +1265 @@
 PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
 @*       for the hypersurface defined by f. 
-ASSUME:  The basering is a commutative polynomial ring in char 0.
-BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm 
-@*       by Masayuki Noro and then a system of linear equations is solved by linear reductions.
+ASSUME:  The basering is commutative and of characteristic 0.
+BACKGROUND: In this proc, the initial Malgrange ideal is computed according to
+@*       the algorithm by Masayuki Noro and then a system of linear equations is
+@*       solved by linear reductions.
 NOTE:    In the output list, the ideal contains all the roots 
 @*       and the intvec their multiplicities.
@@ -1201 +1275 @@
 @*       If t<>0, a matrix ordering is used for Groebner basis computations,
 @*       otherwise, and by default, a block ordering is used.
-@*       If v is a positive weight vector, v is used for homogenization computations, 
-@*       otherwise and by default, no weights are used.
+@*       If v is a positive weight vector, v is used for homogenization
+@*       computations, otherwise and by default, no weights are used.
 DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
@@ -1253 +1327 @@
 
 proc bfctSyz (poly f, list #)
-"USAGE:  bfctSyz(f [,r,s,t,u,v]);  f a poly, r,s,t,u optional ints, v an optional intvec
+"USAGE:  bfctSyz(f [,r,s,t,u,v]);  f poly, r,s,t,u optional ints, v opt. intvec
 RETURN:  list of ideal and intvec
-PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) 
+PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
 @*       for the hypersurface defined by f
-ASSUME:  The basering is a commutative polynomial ring in char 0.
-BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm 
-@*       by Masayuki Noro and then a system of linear equations is solved by computing syzygies.
-NOTE:    In the output list, the ideal contains all the roots 
-@*       and the intvec their multiplicities.
-@*       If r<>0, @code{std} is used for GB computations in characteristic 0, 
-@*       otherwise, and by default, @code{slimgb} is used. 
-@*       If s<>0, a matrix ordering is used for GB computations, otherwise, 
+ASSUME:  The basering is commutative and of characteristic 0.
+BACKGROUND: In this proc, the initial Malgrange ideal is computed according to
+@*       the algorithm by Masayuki Noro and then a system of linear equations is
+@*       solved by computing syzygies.
+NOTE:    In the output list, the ideal contains all the roots and the intvec 
+@*       their multiplicities.
+@*       If r<>0, @code{std} is used for GB computations in characteristic 0,
+@*       otherwise, and by default, @code{slimgb} is used.
+@*       If s<>0, a matrix ordering is used for GB computations, otherwise,
 @*       and by default, a block ordering is used.
-@*       If t<>0, the computation of the intersection is solely performed over 
+@*       If t<>0, the computation of the intersection is solely performed over
 @*       charasteristic 0, otherwise and by default, a modular method is used.
-@*       If u<>0 and by default, @code{std} is used for GB computations in 
-@*       characteristic >0, otherwise, @code{slimgb} is used. 
-@*       If v is a positive weight vector, v is used for homogenization 
+@*       If u<>0 and by default, @code{std} is used for GB computations in
+@*       characteristic >0, otherwise, @code{slimgb} is used.
+@*       If v is a positive weight vector, v is used for homogenization
 @*       computations, otherwise and by default, no weights are used.
 DISPLAY: If printlevel=1, progress debug messages will be printed,
@@ -1285 +1360 @@
   // and one for the engine used for Groebner basis computations in char >0
   // in # can also be the optional weight vector
-  def save = basering;
-  int n = nvars(save);
+  int n = nvars(basering);
   int whichengine = 0; // default
   int methodord   = 0; // default
@@ -1343 +1417 @@
 "USAGE:  bfctIdeal(I,w[,s,t]);  I an ideal, w an intvec, s,t optional ints
 RETURN:  list of ideal and intvec
-PURPOSE: computes the roots of the global b-function of I wrt the weight vector (-w,w).
-ASSUME:  The basering is the n-th Weyl algebra in char 0, where the sequence of
-@*       the variables is x(1),...,x(n),D(1),...,D(n).
-BACKGROUND:  In this proc, the initial ideal of I is computed according to the algorithm by
-@*       Masayuki Noro and then a system of linear equations is solved by linear reductions.
-NOTE:    In the output list, the ideal contains all the roots 
-@*       and the intvec their multiplicities.
+PURPOSE: computes the roots of the global b-function of I wrt the weight (-w,w).
+ASSUME:  The basering is the n-th Weyl algebra in characteristic 0 and  for all
+@*       1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the
+@*       sequence of variables is given by x(1),...,x(n),D(1),...,D(n), 
+@*       where D(i) is the differential operator belonging to x(i).
+@*       Further assume that I is holonomic.
+BACKGROUND:  In this proc, the initial ideal of I is computed according to the
+@*       algorithm by Masayuki Noro and then a system of linear equations is
+@*       solved by linear reductions.
+NOTE:    In the output list, the ideal contains all the roots and the intvec
+@*       their multiplicities.
 @*       If s<>0, @code{std} is used for GB computations in characteristic 0,
 @*       otherwise, and by default, @code{slimgb} is used. 
@@ -1379 +1457 @@
     }
   }
+  if (isWeyl()==0) { ERROR("basering is not a Weyl algebra"); }
+  for (i=1; i<=n; i++)
+  {
+    if (bracket(var(i+n),var(i))<>1)
+    {
+      ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i)));
+    }
+  }
+  int isH = isHolonomic(I);
+  if (isH<>1)
+  {
+    print("// given ideal is not holonomic");
+    print("// setting bound for degree of b-fubction to 10 and proceeding");
+    isH = 10;
+  }
+  else { isH = 0; } // no degree bound for pIntersect
   ideal J = initialIdealW(I,-w,w,whichengine,methodord);
   poly s;
@@ -1385 +1479 @@
     s = s + w[i]*var(i)*var(n+i);
   }
-  vector b = pIntersect(s,J);
-  list l = listofroots(b,0);
+  vector b = pIntersect(s,J,isH);
+  list RL = ringlist(save); RL = RL[1..4];
+  RL[2] = list(safeVarName("s"));
+  RL[3] = list(list("dp",intvec(1)),list("C",intvec(0)));
+  def @S = ring(RL); setring @S;
+  vector b = imap(save,b);
+  poly bs = vec2poly(b);
+  list l = bFactor(bs);
+  setring save;
+  list l = imap(@S,l);
   return(l);
 }
@@ -1407 +1509 @@
 "USAGE:  bfctOneGB(f [,s,t]);  f a poly, s,t optional ints
 RETURN:  list of ideal and intvec
-PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the 
+PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the
 @*       hypersurface defined by f, using only one GB computation
-ASSUME:  The basering is a commutative polynomial ring in char 0.
-BACKGROUND: In this proc, the initial Malgrange ideal is computed based on the 
+ASSUME:  The basering is commutative and of characteristic 0.
+BACKGROUND: In this proc, the initial Malgrange ideal is computed based on the
 @*       algorithm by Masayuki Noro and combined with an elimination ordering.
-NOTE:    In the output list, the ideal contains all the roots 
-@*       and the intvec their multiplicities.
+NOTE:    In the output list, the ideal contains all the roots and the intvec
+@*       their multiplicities.
 @*       If s<>0, @code{std} is used for the GB computation, otherwise,
 @*       and by default, @code{slimgb} is used. 
@@ -1424 +1526 @@
 {
   int ppl = printlevel - voice +2;
+  if (isCommutative() == 0) { ERROR("basering must be commutative"); }
   def save = basering;
   int n = nvars(save);
+  if (char(save) <> 0) 
+  {
+    ERROR("characteristic of basering has to be 0");
+  }
+  list L = ringlist(save);
+  int qr;
+  if (L[4] <> 0) // qring?
+  {
+    print("// basering is qring:");
+    print("// discarding the quotient and proceeding...");
+    L[4] = 0;
+    qr = 1;
+    def save2 = ring(L);
+    setring save2;
+    poly f = imap(save,f);
+  }
   int noofvars = 2*n+4;
   int i;
@@ -1454 +1573 @@
   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(noofvars-1),lp(1));
   }
   else // M() ordering
@@ -1465 +1584 @@
       @Ord[2+i,noofvars - i] = -1;
     }
-    dbprint(ppl,"weights for ordering:",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);
+    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];
@@ -1479 +1598 @@
     @relD[i+2,n+3+i] = h^2;
   }
-  dbprint(ppl,"nc relations:",@relD);
-  def DDh = nc_algebra(1,@relD);
-  setring DDh;
-  dbprint(ppl,"computing in ring",DDh);
+  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);
@@ -1492 +1611 @@
     I = I, D(i)+diff(f,x(i))*Dt;
   }
-  dbprint(ppl, "starting Groebner basis computation with engine:", whichengine);
+  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 is" ,I);
+  dbprint(ppl, "// finished Groebner basis computation");
+  dbprint(ppl, "// I before dehomogenization: "+string(I));
   I = subst(I,h,1); // dehomogenization
-  dbprint(ppl, "I after dehomogenization is" ,I);
+  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, "the initial ideal:", string(matrix(I)));
+  dbprint(ppl, "// the initial ideal:", string(matrix(I)));
   intvec tonselect = 1;
   for (i=3; i<=2*n+4; i++) { tonselect = tonselect,i; }
   I = nselect(I,tonselect);
-  dbprint(ppl, "generators containing only s:", string(matrix(I)));
+  dbprint(ppl, "// generators containing only s:", string(matrix(I)));
   I = engine(I, whichengine); // is now a principal ideal;
-  setring save;
+  if (qr == 1) { setring save2; }
+  else { setring save; }
   ideal J; J[2] = var(1);
-  map @m = DDh,J;
+  map @m = Dh,J;
   ideal I = @m(I);
   poly p = I[1];
   list l = listofroots(p,1);
+  if (qr == 1)
+  {
+    setring save;
+    list l = imap(save2,l);
+  }
   return(l);
 }
@@ -1525 +1650 @@
 "USAGE:  bfctAnn(f [,r]);  f a poly, r an optional int
 RETURN:  list of ideal and intvec
-PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) 
-@*       for the hypersurface defined by f
-ASSUME:  The basering is a commutative polynomial ring in char 0.
+PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the
+@*       hypersurface defined by f.
+ASSUME:  The basering is commutative and of characteristic 0.
 BACKGROUND: In this proc, ann(f^s) is computed and then a system of linear
 @*       equations is solved by linear reductions.
-NOTE:    In the output list, the ideal contains all the roots 
-@*       and the intvec their multiplicities.
+NOTE:    In the output list, the ideal contains all the roots and the intvec
+@*       their multiplicities.
 @*       If r<>0, @code{std} is used for GB computations,
 @*       otherwise, and by default, @code{slimgb} is used. 
@@ -1549 +1674 @@
     }
   }
-  list b = bfctengine(f,1,whichengine,0,1,0,0,0);
+  list b = bfctengine(f,1,whichengine,0,0,0,0,0);
   return(b);
 }

Singular/LIB/dmodapp.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: dmodapp.lib,v 1.20 2009-03-01 17:41:31 levandov Exp $";
+version="$Id: dmodapp.lib,v 1.21 2009-03-06 20:32:29 levandov Exp $";
 category="Noncommutative";
 info="
@@ -7 +7 @@
 @*             Daniel Andres, daniel.andres@math.rwth-aachen.de
 
-GUIDE: Let R = K[x1,..xN] and D be the Weyl algebra in variables x1,..xN,d1,..dN.
-In this library there are the following procedures for algebraic D-modules:
+GUIDE: Let K be a field of characteristic 0, R = K[x1,..xN] and 
+@* D be the Weyl algebra in variables x1,..xN,d1,..dN.
+@* In this library there are the following procedures for algebraic D-modules:
 @* - localization of a holonomic module D/I with respect to a mult. closed set
-of all powers of a given polynomial F from R. Our aim is to compute an ideal L
-in D, such that D/L is a presentation of a localized module. Such L always exists, since
-such localizations are known to be holonomic and thus cyclic modules.
-The procedures for the localization are DLoc, SDLoc and DLoc0.
-
-@* - annihilator in Weyl algebra of a given polynomial F from R as well as of a given
-rational function G/F from Quot(R). These can be computed via annPoly resp. annRat.
-
-@* - initial form and initial ideals in Weyl algebras with respect to a given weight vector can be computed with the help of inForm, initialMalgrange, initialIdealW.
-
-@* - appelF1, appelF2 and appelF4 return ideals in parametric Weyl algebra, which
-annihilate corresponding Appel hypergeometric functions.
-
+@* of all powers of a given polynomial F from R. Our aim is to compute an 
+@* ideal L in D, such that D/L is a presentation of a localized module. Such L 
+@* always exists, since such localizations are known to be holonomic and thus 
+@* cyclic modules. The procedures for the localization are DLoc,SDLoc and DLoc0.
+@*
+@* - annihilator in Weyl algebra of a given polynomial F from R as well as 
+@* of a given rational function G/F from Quot(R). These can be computed via 
+@* procedures annPoly resp. annRat.
+@*
+@* - initial form and initial ideals in Weyl algebras with respect to a given 
+@* weight vector can be computed with  inForm, initialMalgrange, initialIdealW.
+@*
+@* - appelF1, appelF2 and appelF4 return ideals in parametric Weyl algebras, 
+@* which annihilate corresponding Appel hypergeometric functions.
+
+REFERENCES: 
+@* (SST) Saito, Sturmfels, Takayama 'Groebner Deformations of Hypergeometric 
+@*         Differential Equations', Springer, 2000
+@* (ONW) Oaku, Takayama, Walther 'A Localization Algorithm for D-modules', 2000
 
 MAIN PROCEDURES:
@@ -29 +36 @@
 annRat(f,g);  annihilator of a rational function f/g in the corr. Weyl algebra
 DLoc(I,F);     presentation of the localization of D/I w.r.t. f^s
-SDLoc(I, F);  a generic presentation of the localization of D/I w.r.t. f^s, for D a Weyl algebra
-DLoc0(I, F);  presentation of the localization of D/I w.r.t. f^s, based on the procedure SDLoc
-
-initialMalgrange(f[,s,t,u,v]); Groebner basis of the initial Malgrange ideal for a given poly
+SDLoc(I, F);  a generic presentation of the localization of D/I w.r.t. f^s
+DLoc0(I, F);  presentation of the localization of D/I w.r.t. f^s, based on SDLoc
+
+initialMalgrange(f[,s,t,v]); Groebner basis of the initial Malgrange ideal for f
 initialIdealW(I,u,v[,s,t]); initial ideal of a given ideal w.r.t. given weights
 inForm(f,w);     initial form of a poly/ideal w.r.t. a given weight
@@ -39 +46 @@
 AUXILIARY PROCEDURES:
 
+bFactor(F);  computes the roots of irreducible factors of an univariate poly
 appelF1();      create an ideal annihilating Appel F1 function
 appelF2();      create an ideal annihilating Appel F2 function
@@ -57 +65 @@
 LIB "dmod.lib"; // loads e.g. nctools.lib
 LIB "bfun.lib"; //formerly LIB "bfct.lib";
+LIB "nctools.lib"; // for isWeyl etc
+LIB "gkdim.lib";
 
 // todo: complete and include into above list
@@ -83 +93 @@
   example insertGenerator;
   example deleteGenerator;
-}
-
-
-proc initialIdealW (ideal I, intvec u, intvec v, list #)
-"USAGE:  initialIdealW(I,u,v [,s,t]);  I ideal, u,v intvecs, s,t optional ints
-RETURN:  ideal, GB of initial ideal of the input ideal wrt the weights u and v
-ASSUME: basering is the n-th Weyl algebra and the sequence of the indeterminates is x(1),...,x(n),D(1),...,D(n).
-PURPOSE: computes the initial ideal with respect to given weights.
-NOTE:    Let I be an ideal in the n-th Weyl algebra D, then
-@*       u and v understood as intvecs of weights for x(i) and D(i) respectively.
-@*       Note that the returned ideal is not an ideal in D, but an ideal in the associated
-@*       graded ring while the Groebner basis is a subset of D.
-@*       If s<>0, @code{std} is used for Groebner basis computations, otherwise,
-@*       and by default, @code{slimgb} is used.
-@*       If t<>0, a matrix ordering is used for Groebner basis computations,
-@*       otherwise, and by default, a block ordering is used.
-DISPLAY: If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example initialIdealW; shows examples
-"
-{
-  int ppl = printlevel - voice +2;
-  int i;
-  def save = basering;
-  int whichengine = 0; // default
-  int methodord   = 0; // default
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      whichengine = int(#[1]);
-    }
-    if (size(#)>1)
-    {
-      if (typeof(#[2])=="int" || typeof(#[2])=="number")
-      {
-        methodord = int(#[2]);
-      }
-    }
-  }
-  def D = initialIdealEngine("initialIdeal", whichengine, methodord, I, u, v);
-  ideal inF = fetch(D,inF); attrib(inF,"isSB",1);
-  return(inF);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring @D = 0,(x,Dx),dp;
-  def D = Weyl();
-  setring D;
-  intvec u = -1; intvec v = 2;
-  ideal I = x^2*Dx^2,x*Dx^4;
-  ideal J = initialIdealW(I,u,v); J;
-}
-
-proc initialMalgrange (poly f,list #)
-"USAGE:  initialMalgrange(f, [,s,t,u,v]);  f poly, s,t,u optional ints, v optional intvec
-RETURN:  ring, the Weyl algebra induced by basering, extended with vars t and Dt,
-@*       containing the ideal \"inF\", being the initial ideal of the Malgrange
-@*       ideal of f.
-PURPOSE: computes the initial Malgrange ideal of a given poly wrt the weight
-@*       vector (-1,0...,0,1,0,...,0) such that the weight of t is -1 and the
-@*       weight of Dt is 1.
-NOTE:    Activate the output ring with the @code{setring} command.
-@*       Varnames of the basering should not include t and Dt.
-@*       If s<>0, @code{std} is used for Groebner basis computations,
-@*       otherwise, and by default, @code{slimgb} is used.
-@*       If t<>0, a matrix ordering is used for Groebner basis computations,
-@*       otherwise, and by default, a block ordering is used.
-@*       If u<>0, the order of the variables is reversed.
-@*       If v is a positive weight vector, v is used for homogenization
-@*       computations, otherwise and by default, no weights are used.
-DISPLAY: If printlevel=1, progress debug messages will be printed,
-@*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example initialMalgrange; shows examples
-"
-{
-  int ppl = printlevel - voice +2;
-  def save = basering;
-  int n = nvars(save);
-  int whichengine = 0; // default
-  int methodord   = 0; // default
-  int reversevars = 0; // default
-  intvec u0 = 0;
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      whichengine = int(#[1]);
-    }
-    if (size(#)>1)
-    {
-      if (typeof(#[2])=="int" || typeof(#[2])=="number")
-      {
-        methodord = int(#[2]);
-      }
-      if (size(#)>2)
-      {
-        if (typeof(#[3])=="int" || typeof(#[3])=="number")
-        {
-          reversevars = int(#[3]);
-        }
-        if (size(#)>3)
-        {
-          if (typeof(#[4])=="intvec" && size(#[4])==n && allPositive(#[4])==1)
-          {
-            u0 = #[4];
-          }
-        }
-      }
-    }
-  }
-  if (u0 == 0)
-  {
-    u0 = 1:n;
-  }
-  def D = initialIdealEngine("initialMalgrange",whichengine, methodord, f, 0, 0, u0, reversevars);
-  setring save;
-  return(D);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(x,y),dp;
-  poly f = x^2+y^3+x*y^2;
-  def D = initialMalgrange(f);
-  setring D;
-  inF;
-  setring r;
-  intvec v = 3,2;
-  def D2 = initialMalgrange(f,1,0,1,v);
-  setring D2;
-  inF;
-}
-
-static proc initialIdealEngine(string calledfrom, int whichengine, int blockord, list #)
-{
-  //#[1] = f or I
-  //#[2] = u
-  //#[3] = v
-  //#[4] = u0
-  //#[5] = reversevars
-  int ppl = printlevel - voice +2;
-  def save = basering;
-  int i,n,noofvars;
-  n = nvars(save);
-  intvec uv;
-  if (calledfrom == "initialIdeal")
-  {
-    ideal I = #[1];
-    intvec u = #[2];
-    intvec v = #[3];
-    uv = u,v,0;
-    n = n/2;
-    noofvars = 2*n+1;
-  }
-  else // initialMalgrange
-  {
-    poly f = #[1];
-    uv[n+2] = 1;
-    noofvars = 2*n+3;
-    intvec u0 = #[4];
-    int reversevars = #[5];
-    ring r = 0,(x(1..n)),wp(u0);
-    poly f = fetch(save,f);
-    uv[1] = -1; uv[noofvars] = 0;
-  }
-  //  for the ordering
-  intvec @a;
-  if (calledfrom == "initialMalgrange")
-  {
-    int d = deg(f);
-    intvec weighttx = d;
-    for (i=1; i<=n; i++)
-    {
-      weighttx[i+1] = u0[n-i+1];
-    }
-    intvec weightD = 1;
-    for (i=1; i<=n; i++)
-    {
-      weightD[i+1] = d+1-u0[n-i+1];
-    }
-    @a = weighttx,weightD,1;
-  }
-  else // initialIdeal
-  {
-    @a = 1:noofvars;
-  }
-  if (blockord == 0) // default: blockordering
-  {
-    if (calledfrom == "initialMalgrange")
-    {
-      ring Dh = 0,(t,x(n..1),Dt,D(n..1),h),(a(@a),a(uv),dp(noofvars-1),lp(1));
-    }
-    else // initialIdeal
-    {
-      ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),dp(noofvars-1),lp(1));
-    }
-  }
-  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;
-    }
-    dbprint(ppl,"weights for ordering:",transpose(@a));
-    dbprint(ppl,"the ordering matrix:",@Ord);
-    if (calledfrom == "initialMalgrange")
-    {
-      ring Dh = 0,(t,x(n..1),Dt,D(n..1),h),(a(@a),M(@Ord));
-    }
-    else // initialIdeal
-    {
-      ring Dh = 0,(x(1..n),D(1..n),h),(a(@a),M(@Ord));
-    }
-  }
-  dbprint(ppl,"the ring Dh:",Dh);
-  // non-commutative relations
-  matrix @relD[noofvars][noofvars];
-  if (calledfrom == "initialMalgrange")
-  {
-    for (i=1; i<=n+1; i++) { @relD[i,n+1+i] = h^(d+1); }
-  }
-  else // initialIdeal
-  {
-    for (i=1; i<=n; i++)
-    {
-      @relD[i,n+i] = h^2;
-    }
-  }
-  dbprint(ppl,"nc relations:",@relD);
-  def DDh = nc_algebra(1,@relD);
-  setring DDh;
-  dbprint(ppl,"computing in ring",DDh);
-  ideal I;
-  if (calledfrom == "initialIdeal")
-  {
-    I = fetch(save,I);
-    I = homog(I,h);
-  }
-  else
-  {
-    poly f = imap(r,f);
-    kill r;
-    f = homog(f,h);
-    I = t-f;  // defining the Malgrange ideal
-    for (i=1; i<=n; i++)
-    {
-      I = I, D(i)+diff(f,x(i))*Dt;
-    }
-  }
-  dbprint(ppl, "starting Groebner basis computation with engine:", whichengine);
-  I = engine(I, whichengine);
-  dbprint(ppl, "finished Groebner basis computation");
-  dbprint(ppl, "I before dehomogenization is" ,I);
-  I = subst(I,h,1); // dehomogenization
-  dbprint(ppl, "I after dehomogenization is" ,I);
-  I = inForm(I,uv); // we are only interested in the initial form w.r.t. uv
-  if (calledfrom == "initialMalgrange")
-  {
-    // keep the names of the variables of the basering
-    setring save;
-    list rl = ringlist(save);
-    list varnames = rl[2];
-    for (i=1; i<=n; i++)
-    {
-      if (varnames[i] == "t")
-      {
-        ERROR("Variable names should not include t");
-      }
-    }
-    list newvarnamesrev = "t";
-    newvarnamesrev[n+2] = "Dt";
-    for (i=1; i<=n; i++)
-    {
-      newvarnamesrev[i+1] = varnames[n+1-i];
-      newvarnamesrev[i+n+2] = "D"+varnames[n+1-i];
-    }
-    rl[2]=newvarnamesrev;
-    def @Drev = ring(rl);
-    setring @Drev;
-    def Drev = Weyl(@Drev);
-    setring Drev;
-    ideal I = fetch(DDh,I);
-    kill Dh, DDh;
-    if (reversevars == 0)
-    {
-      setring save;
-      list newvarnames = "t";
-      newvarnames[n+2] = "Dt";
-      for (i=1; i<=n; i++)
-      {
-        newvarnames[i+1] = varnames[i];
-        newvarnames[i+n+2] = "D"+varnames[i];
-      }
-      rl[2] = newvarnames;
-      def @D = ring(rl);
-      setring @D;
-      def D = Weyl(@D);
-      setring D;
-      ideal I = imap(Drev,I);
-    }
-  }
-  else // initialIdeal
-  {
-    ring @D = 0,(x(1..n),D(1..n)),dp;
-    def D = Weyl(@D);
-    setring D;
-    ideal I = imap(DDh,I);
-    kill Dh,DDh;
-  }
-  dbprint(ppl, "starting cosmetic Groebner basis computation with engine:", whichengine);
-  I = engine(I, whichengine);
-  dbprint(ppl,"finished cosmetic Groebner basis computation:");
-  dbprint(ppl,"the initial ideal is:", I);
-  ideal inF = I; attrib(inF,"isSB",1);
-  export(inF);
-  return(basering);
 }
 
@@ -588 +277 @@
 proc appelF1()
 "USAGE:  appelF1();
-RETURN:  ring
-PURPOSE: define the ideal in a parametric Weyl algebra, which annihilates Appel F1 hypergeometric function
+RETURN:  ring  (and exports an ideal into it)
+PURPOSE: define the ideal in a parametric Weyl algebra, 
+@* which annihilates Appel F1 hypergeometric function
 NOTE: the ideal called  IAppel1 is exported to the output ring
 EXAMPLE: example appelF1; shows examples
@@ -618 +308 @@
 proc appelF2() //(number a,b,c)
 "USAGE:  appelF2();
-RETURN:  ring
-PURPOSE: define the ideal in a parametric Weyl algebra, which annihilates Appel F2 hypergeometric function
+RETURN:  ring (and exports an ideal into it)
+PURPOSE: define the ideal in a parametric Weyl algebra, 
+@* which annihilates Appel F2 hypergeometric function
 NOTE: the ideal called  IAppel2 is exported to the output ring
 EXAMPLE: example appelF2; shows examples
@@ -647 +338 @@
 proc appelF4()
 "USAGE:  appelF4();
-RETURN:  ring
-PURPOSE: define the ideal in a parametric Weyl algebra, which annihilates Appel F4 hypergeometric function
+RETURN:  ring  (and exports an ideal into it)
+PURPOSE: define the ideal in a parametric Weyl algebra, 
+@* which annihilates Appel F4 hypergeometric function
 NOTE: the ideal called  IAppel4 is exported to the output ring
 EXAMPLE: example appelF4; shows examples
@@ -750 +442 @@
 PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s
 NOTE:    In the basering, the following objects are exported:
-@*       - the ideal LD0 (which is a Groebner basis) is the presentation of the localization
-@*       - the ideal BS contains the roots with multiplicities of a Bernstein polynomial of D/I w.r.t f.
-@*       If printlevel=1, progress debug messages will be printed,
+@* the ideal LD0 (in Groebner basis) is the presentation of the localization
+@* the list BS contains roots with multiplicities of Bernstein poly of (D/I)_f.
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
 EXAMPLE: example DLoc; shows examples
@@ -759 +451 @@
   /* runs SDLoc and DLoc0 */
   /* assume: run from Weyl algebra */
+  if (dmodappassumeViolation())
+  {
+    ERROR("Basering is inappropriate: characteristic>0 or qring present");
+  }
+  if (!isWeyl())
+  {
+    ERROR("Basering is not a Weyl algebra");
+  }
+  if (defined(LD0) || defined(BS)) 
+  {
+    ERROR("Reserved names LD0 and/or BS are used. Please rename the objects.");
+  }
   int old_printlevel = printlevel;
   printlevel=printlevel+1;
@@ -782 +486 @@
   "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,Dx,Dy),dp;
-  def R = Weyl();    setring R;
+  def R = Weyl();    setring R; // Weyl algebra in variables x,y,Dx,Dy
   poly F = x2-y3;
   ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F;
-  DLoc(I, x2-y3);
-  LD0;
-  BS;
+  // I is not holonomic, since its dimension is not 4/2=2
+  gkdim(I);
+  DLoc(I, x2-y3); // exports LD0 and BS
+  LD0; // localized module (R/I)_f is isomorphic to R/LD0
+  BS; // description of b-function for localization
 }
 
@@ -793 +499 @@
 "USAGE:  DLoc0(I, F);  I an ideal, F a poly
 RETURN:  ring
-PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s, where D is a Weyl Algebra, based on the output of procedure SDLoc
+PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s, 
+@*           where D is a Weyl Algebra, based on the output of procedure SDLoc
 ASSUME: the basering is similar to the output ring of SDLoc procedure
 NOTE:    activate this ring with the @code{setring} command. In this ring,
-@*       - the ideal LD0 (which is a Groebner basis) is the presentation of the localization
-@*       - the ideal BS contains the roots with multiplicities of a Bernstein polynomial of D/I w.r.t f.
-@*       If printlevel=1, progress debug messages will be printed,
+@* the ideal LD0 (in Groebner basis) is the presentation of the localization
+@* the list BS contains roots and multiplicities of Bernstein poly of (D/I)_f.
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
 EXAMPLE: example DLoc0; shows examples
 "
 {
+  if (dmodappassumeViolation())
+  {
+    ERROR("Basering is inappropriate: characteristic>0 or qring present");
+  }
   /* assume: to be run in the output ring of SDLoc */
   /* doing: add F, eliminate vars*Dvars, factorize BS */
@@ -976 +687 @@
     dbprint(ppl-1, @R4);
     ideal K4 = imap(@R2,K2);
+    intvec vopt = option(get);
     option(redSB);
     dbprint(ppl,"// -3-2- the final cosmetic std");
     K4 = groebner(K4);  // std does the job too
+    option(set,vopt);
     // total cleanup
     setring @R2;
@@ -1014 +727 @@
   "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,Dx,Dy),dp;
-  def R = Weyl();    setring R;
+  def R = Weyl();    setring R; // Weyl algebra in variables x,y,Dx,Dy
   poly F = x2-y3;
   ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F;
-  def W = SDLoc(I,F);  setring W; // creates ideal LD
-  def U = DLoc0(LD, x2-y3);  setring U;
-  LD0;
-  BS;
+  // moreover I is not holonomic, since its dimension is not 2 = 4/2
+  gkdim(I); // 3
+  def W = SDLoc(I,F);  setring W; // creates ideal LD in W = R[s]
+  def U = DLoc0(LD, x2-y3);  setring U; // compute in R 
+  LD0; // Groebner basis of the presentation of localization
+  BS; // description of b-function for localization
 }
 
@@ -1027 +742 @@
 "USAGE:  SDLoc(I, F);  I an ideal, F a poly
 RETURN:  ring
-PURPOSE: compute a generic presentation of the localization of D/I w.r.t. f^s, where D is a Weyl Algebra
-ASSUME: the basering is a Weyl algebra
+PURPOSE: compute a generic presentation of the localization of D/I w.r.t. f^s
+ASSUME: the basering D is a Weyl algebra
 NOTE:    activate this ring with the @code{setring} command. In this ring,
-@*       - the ideal LD (which is a Groebner basis) is the presentation of the localization
-@*       If printlevel=1, progress debug messages will be printed,
+@*  the ideal LD (in Groebner basis) is the presentation of the localization
+DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
 EXAMPLE: example SDLoc; shows examples
@@ -1039 +754 @@
   /* printlevel >=4 gives debug info */
   /* assume: we're in the Weyl algebra D  in x1,x2,...,d1,d2,... */
+
+  if (dmodappassumeViolation())
+  {
+    ERROR("Basering is inappropriate: characteristic>0 or qring present");
+  }
+  if (!isWeyl())
+  {
+    ERROR("Basering is not a Weyl algebra");
+  }
   def save = basering;
   /* 1. create D <t, dt, s > as in LOT */
@@ -1058 +782 @@
   RName[1] = "@t";
   RName[2] = "@Dt";
-  RName[3] = "s";
+  RName[3] = "@s";
   for(i=1;i<=N;i++)
   {
@@ -1065 +789 @@
       if (Name[i] == RName[j])
       {
-        ERROR("Variable names should not include @t,@Dt,s");
+        ERROR("Variable names should not include @t,@Dt,@s");
       }
     }
@@ -1073 +797 @@
   tmp[1] = "@t";
   tmp[2] = "@Dt";
-  list SName ; SName[1] = "s";
+  list SName ; SName[1] = "@s";
   list NName = tmp + Name + SName;
   L[2]   = NName;
@@ -1115 +839 @@
   setring @R;
   kill @R@;
-  dbprint(ppl,"// -1-1- the ring @R(t,Dt,_x,_Dx,s) is ready");
+  dbprint(ppl,"// -1-1- the ring @R(@t,@Dt,_x,_Dx,@s) is ready");
   dbprint(ppl-1, @R);
   poly  F = imap(save,F);
@@ -1132 +856 @@
   I = I, @t - F;
   // t*Dt + s +1 reduced with t-f gives f*Dt + s
-  I = I, F*var(2) + var(Nnew);
+  I = I, F*var(2) + var(Nnew); // @s
   // -------- the ideal I is ready ----------
   dbprint(ppl,"// -1-2- starting the elimination of @t,@Dt in @R");
@@ -1221 +945 @@
   "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y,Dx,Dy),dp;
-  def R = Weyl();
+  def R = Weyl(); // Weyl algebra on the variables x,y,Dx,Dy
   setring R;
   poly F = x2-y3;
-  ideal I = Dx*F, Dy*F;
+  ideal I = Dx*F, Dy*F; 
+  // note, that I is not holonomic, since it's dimension is not 2
+  gkdim(I); // 3, while dim R = 4
   def W = SDLoc(I,F);
-  setring W;
-  LD;
+  setring W; // = R[s], where s is a new variable
+  LD; // Groebner basis of s-parametric presentation
 }
 
@@ -1233 +959 @@
 "USAGE:  annRat(g,f);  f, g polynomials
 RETURN:  ring
-PURPOSE: compute the ideal in Weyl algebra, annihilating the rational function g*f^{-1}
-NOTE:    activate the output ring with the @code{setring} command.
-@*       In the output ring:
-@*       - the ideal LD (which is given in a Groebner basis) is the annihilator.
-@*       If @code{printlevel}=1, progress debug messages will be printed,
-@*       if @code{printlevel}>=2, all the debug messages will be printed.
+PURPOSE: compute the annihilator of the rational function g/f in Weyl algebra
+NOTE: activate the output ring with the @code{setring} command.
+@*      In the output ring, the ideal LD (in Groebner basis) is the annihilator.
+@*      The algorithm uses the computation of ann f^{-1} via D-modules.
+DISPLAY: If printlevel=1, progress debug messages will be printed,
+@*       if printlevel>=2, all the debug messages will be printed.
+SEE ALSO: annPoly
 EXAMPLE: example annRat; shows examples
 "
 {
-  // computes the annihilator of g/f
+
+  if (dmodappassumeViolation())
+  {
+    ERROR("Basering is inappropriate: characteristic>0 or qring present");
+  }
+
+  // assumptions: f is not a constant
+  if (f==0) { ERROR("Denominator cannot be zero"); }
+  if (leadexp(f) == 0)
+  {
+    // f = const, so use annPoly
+    g = g/f;
+    def @R = annPoly(g);
+    return(@R);
+  }
+    // computes the annihilator of g/f
   def save = basering;
   int ppl = printlevel-voice+2;
@@ -1329 +1071 @@
   // Now, compare with the output of Macaulay2:
   ideal tst = 3*x*Dx + 2*y*Dy + 1, y^3*Dy^2 - x^2*Dy^2 + 6*y^2*Dy + 6*y,
-9*y^2*Dx^2*Dy - 4*y*Dy^3 + 27*y*Dx^2 + 2*Dy^2, 9*y^3*Dx^2 - 4*y^2*Dy^2 + 10*y*Dy -10;
+9*y^2*Dx^2*Dy-4*y*Dy^3+27*y*Dx^2+2*Dy^2, 9*y^3*Dx^2-4*y^2*Dy^2+10*y*Dy -10;
  option(redSB); option(redTail);
  LD = groebner(LD);
@@ -1352 +1094 @@
 "USAGE:  annPoly(f);  f a poly
 RETURN:  ring
-PURPOSE: compute the ideal in Weyl algebra, annihilating the polynomial f
-NOTE:    activate the output ring with the @code{setring} command.
-@*       In the output ring:
-@*       - the ideal LD (which is given in a Groebner basis) is the annihilator.
-@*       If @code{printlevel}=1, progress debug messages will be printed,
-@*       if @code{printlevel}>=2, all the debug messages will be printed.
+PURPOSE: compute the complete annihilator ideal of f in Weyl algebra 
+NOTE:  activate the output ring with the @code{setring} command.
+@*   In the output ring, the ideal LD (in Groebner basis) is the annihilator.
+DISPLAY: If printlevel=1, progress debug messages will be printed,
+@*       if printlevel>=2, all the debug messages will be printed.
 SEE ALSO: annRat
 EXAMPLE: example annPoly; shows examples
@@ -1395 +1136 @@
   poly f = x^2*z - y^3;
   def A = annPoly(f);
-  setring A;
-  LD;
-  gkdim(LD); // must be 3 since LD is holonomic
+  setring A; // A is the 3rd Weyl algebra in 6 variables
+  LD; // the Groebner basis of annihilator
+  gkdim(LD); // must be 3 = 6/2, since A/LD is holonomic module
   NF(Dy^4, LD); // must be 0 since Dy^4 clearly annihilates f
 }
@@ -1500 +1241 @@
 */
 
-proc engine(ideal I, int i)
-"USAGE:  engine(I,i);  I an ideal, i an int
-RETURN:  ideal
+proc engine(def I, int i)
+"USAGE:  engine(I,i);  I  ideal/module/matrix, i an int
+RETURN:  the same type as I
 PURPOSE: compute the Groebner basis of I with the algorithm, chosen via i
 NOTE: By default and if i=0, slimgb is used; otherwise std does the job.
@@ -1509 +1250 @@
 {
   /* std - slimgb mix */
-  ideal J;
+  def J;
+  //  ideal J;
   if (i==0)
   {
@@ -1517 +1259 @@
   {
     // without options -> strange! (ringlist?)
+    intvec v = option(get);
     option(redSB);
     option(redTail);
     J = std(I);
+    option(set, v);
   }
   return(J);
@@ -1660 +1404 @@
 RETURN:  poly
 PURPOSE: reconstruct a monic polynomial in one variable from its factorization
-ASSUME:  s is a string with the name of some variable and L is supposed to consist of two entries:
+ASSUME:  s is a string with the name of some variable and 
+@*         L is supposed to consist of two entries:
 @*         L[1] of the type ideal with the roots of a polynomial
 @*         L[2] of the type intvec with the multiplicities of corr. roots
@@ -1668 +1413 @@
   if (varnum(s)==0)
   {
-    ERROR("no such variable found"); return(0);
+    ERROR("no such variable found in the basering"); return(0);
   }
   poly x = var(varnum(s));
@@ -1677 +1422 @@
   for(int i=1; i<= sl; i++)
   {
-    P = P*((x-RR[i])^IV[i]);
+    if (IV[i] > 0)
+    {
+      P = P*((x-RR[i])^IV[i]);
+    }
+    else
+    {
+      printf("Ignored the root with incorrect multiplicity %s",string(IV[i]));
+    }
   }
   return(P);
@@ -1692 +1444 @@
   factorize(p,2);
 }
+
+static proc safeVarName (string s, string cv)
+{
+  string S;
+  if (cv == "v")  { S = "," + "," + varstr(basering)  + ","; }
+  if (cv == "c")  { S = "," + "," + charstr(basering) + ","; }
+  if (cv == "cv") { S = "," + charstr(basering) + "," + varstr(basering) + ","; }
+  s = "," + s + ",";
+  while (find(S,s) <> 0)
+  {
+    s[1] = "@";
+    s = "," + s;
+  }
+  s = s[2..size(s)-1];
+  return(s)
+}
+
+proc initialIdealW (ideal I, intvec u, intvec v, list #)
+"USAGE:  initialIdealW(I,u,v [,s,t]);  I ideal, u,v intvecs, s,t optional ints
+RETURN:  ideal, GB of initial ideal of the input ideal wrt the weights u and v
+ASSUME:  The basering is the n-th Weyl algebra in characteristic 0 and  for all
+@*       1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the
+@*       sequence of variables is given by x(1),...,x(n),D(1),...,D(n), 
+@*       where D(i) is the differential operator belonging to x(i).
+PURPOSE: computes the initial ideal with respect to given weights.
+NOTE:    u and v are understood as weight vectors for x(i) and D(i)
+@*       respectively.
+@*       If s<>0, @code{std} is used for Groebner basis computations,
+@*       otherwise, and by default, @code{slimgb} is used.
+@*       If t<>0, a matrix ordering is used for Groebner basis computations,
+@*       otherwise, and by default, a block ordering is used.
+DISPLAY: If printlevel=1, progress debug messages will be printed,
+@*       if printlevel>=2, all the debug messages will be printed.
+EXAMPLE: example initialIdealW; shows examples
+"
+{
+
+  if (dmodappassumeViolation())
+  {
+    ERROR("Basering is inappropriate: characteristic>0 or qring present");
+  }
+
+  if (!isWeyl())
+  {
+    ERROR("Basering is not a Weyl algebra");
+  }
+
+  int ppl = printlevel - voice +2;
+  def save = basering;
+  int n = nvars(save)/2;
+  int N = 2*n+1;
+  list RL = ringlist(save);
+  int i;
+  int whichengine = 0; // default
+  int methodord   = 0; // default
+  if (size(#)>0)
+  {
+    if (typeof(#[1])=="int" || typeof(#[1])=="number")
+    {
+      whichengine = int(#[1]);
+    }
+    if (size(#)>1)
+    {
+      if (typeof(#[2])=="int" || typeof(#[2])=="number")
+      {
+        methodord = int(#[2]);
+      }
+    }
+  }
+  if (char(save) <> 0) { ERROR("characteristic of basering has to be 0"); }
+  if (isWeyl() == 0)   { ERROR("basering is not a Weyl algebra"); }
+  for (i=1; i<=n; i++)
+  {
+    if (bracket(var(i+n),var(i))<>1)
+    {
+      ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i)));
+    }
+  }
+  // 1. create  homogenized Weyl algebra
+  // 1.1 create ordering
+  intvec uv = u,v,0;
+  list Lord = list(list("a",intvec(1:N)));
+  list C0 = list("C",intvec(0));
+  if (methodord == 0) // default: blockordering
+  {
+    Lord[2] = list("dp",intvec(1:(N-1)));
+    Lord[3] = list("lp",intvec(1));
+    Lord[4] = C0;
+  }
+  else // M() ordering
+  {
+    intmat @Ord[N][N];
+    @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1);
+    for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; }
+    dbprint(ppl,"// the ordering matrix:",@Ord);
+    Lord[2] = list("M",intvec(@Ord));
+    Lord[3] = C0;
+  }
+  // 1.2 the new var
+  list Lvar = RL[2]; Lvar[N] = safeVarName("h","cv");
+  // 1.3 create commutative ring
+  list L@@Dh = RL; L@@Dh = L@@Dh[1..4];
+  L@@Dh[2] = Lvar; L@@Dh[3] = Lord;
+  def @Dh = ring(L@@Dh); kill L@@Dh;
+  setring @Dh;
+  dbprint(ppl,"// the ring @Dh:",@Dh);
+  // 1.4 create non-commutative relations
+  matrix @relD[N][N];
+  for (i=1; i<=n; i++) { @relD[i,n+i] = var(N)^2; }
+  dbprint(ppl,"// nc relations:",@relD);
+  def Dh = nc_algebra(1,@relD);
+  setring Dh; kill @Dh;
+  dbprint(ppl,"// computing in ring",DDh);
+  // 2. Compute the initial ideal
+  ideal I = imap(save,I);
+  I = homog(I,h);
+  // 2.1 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 is " +string(I));
+  I = subst(I,var(N),1); // dehomogenization
+  dbprint(ppl, "I after dehomogenization is " +string(I));
+  // 2.2 the initial form
+  I = inForm(I,uv);
+  setring save;
+  I = imap(Dh,I); kill Dh;
+  // 2.3 the final GB
+  dbprint(ppl, "// starting cosmetic Groebner basis computation with engine: "+string(whichengine));
+  I = engine(I, whichengine);
+  dbprint(ppl,"// finished cosmetic Groebner basis computation");
+  return(I);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring @D = 0,(x,Dx),dp;
+  def D = Weyl();
+  setring D;
+  intvec u = -1; intvec v = 2;
+  ideal I = x^2*Dx^2,x*Dx^4;
+  ideal J = initialIdealW(I,u,v); J;
+}
+
+proc initialMalgrange (poly f,list #)
+"USAGE:  initialMalgrange(f,[,s,t,v]); f poly, s,t optional ints, v opt. intvec
+RETURN:  ring, the Weyl algebra induced by basering, extended by two new vars
+PURPOSE: computes the initial Malgrange ideal of a given poly wrt the weight
+@*       vector (-1,0...,0,1,0,...,0) such that the weight of t is -1 and the
+@*       weight of Dt is 1.
+ASSUME:  The basering is commutative and of characteristic 0.
+NOTE:    Activate the output ring with the @code{setring} command.
+@*       The returned ring contains the ideal \"inF\", being the initial ideal
+@*       of the Malgrange ideal of f.
+@*       Varnames of the basering should not include t and Dt.
+@*       If s<>0, @code{std} is used for Groebner basis computations,
+@*       otherwise, and by default, @code{slimgb} is used.
+@*       If t<>0, a matrix ordering is used for Groebner basis computations,
+@*       otherwise, and by default, a block ordering is used.
+@*       If v is a positive weight vector, v is used for homogenization
+@*       computations, otherwise and by default, no weights are used.
+DISPLAY: If printlevel=1, progress debug messages will be printed,
+@*       if printlevel>=2, all the debug messages will be printed.
+EXAMPLE: example initialMalgrange; shows examples
+"
+{
+
+  if (dmodappassumeViolation())
+  {
+    ERROR("Basering is inappropriate: characteristic>0 or qring present");
+  }
+
+  if (!isCommutative())
+  {
+    ERROR("Basering must be commutative");
+  }
+
+  int ppl = printlevel - voice +2;
+  def save = basering;
+  int n = nvars(save);
+  int i;
+  int whichengine = 0; // default
+  int methodord   = 0; // default
+  intvec u0 = 0;
+  if (size(#)>0)
+  {
+    if (typeof(#[1])=="int" || typeof(#[1])=="number")
+    {
+      whichengine = int(#[1]);
+    }
+    if (size(#)>1)
+    {
+      if (typeof(#[2])=="int" || typeof(#[2])=="number")
+      {
+        methodord = int(#[2]);
+      }
+      if (size(#)>2)
+      {
+        if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1)
+        {
+          u0 = #[3];
+        }
+      }
+    }
+  }
+  if (u0 == 0)
+  {
+    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);
+  list C0 = list("C",intvec(0));
+  // 1. get the weighted degree of f
+  list Lord = list(list("wp",u0),C0);
+  list L@r = RL;
+  L@r[3] = Lord;
+  def r = ring(L@r); kill L@r;
+  setring r;
+  poly f = imap(save,f);
+  int d = deg(f);
+  setring save; kill r;
+  // 2. create homogenized extended Weyl algebra
+  int N = 2*n+3;
+  // 2.1 create names for vars
+  string vart  = safeVarName("t","cv");
+  string varDt = safeVarName("D"+vart,"cv");
+  while (varDt <> "D"+vart)
+  {
+    vart  = safeVarName("@"+vart,"cv");
+    varDt = safeVarName("D"+vart,"cv");
+  }
+  list Lvar,Lvarh;
+  Lvar[1] = vart; Lvar[n+2] = varDt;
+  for (i=1; i<=n; i++)
+  {
+    Lvar[i+1]   = string(var(i));
+    Lvar[i+n+2] = safeVarName("D" + string(var(i)),"cv");
+  }
+  Lvarh = Lvar; Lvarh[N] = safeVarName("h","cv");
+  //  2.2 create ordering
+  intvec uv,@a,weighttx,weightD;
+  uv[1] = -1; uv[n+2] = 1; uv[N] = 0;
+  weighttx = d; weightD = 1;
+  for (i=1; i<=n; i++)
+  {
+    weighttx[i+1] = u0[n-i+1];
+    weightD[i+1]  = d+1-u0[n-i+1];
+  }
+  @a = weighttx,weightD,1;
+  Lord[1] = list("a",@a);
+  if (methodord == 0) // default: block ordering
+  {
+    Lord[2] = list("a",uv);
+    Lord[3] = list("dp",intvec(1:(N-1)));
+    Lord[4] = list("lp",intvec(1));
+    Lord[5] = C0;
+  }
+  else // M() ordering
+  {
+    intmat @Ord[N][N];
+    @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1);
+    for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; }
+    dbprint(ppl,"// weights for ordering: "+string(transpose(@a)));
+    dbprint(ppl,"// the ordering matrix:",@Ord);
+    Lord[2] = list("M",intvec(@Ord));
+    Lord[3] = C0;
+  }
+  // 2.3 create commutative ring
+  list L@@Dh = RL;
+  L@@Dh[2] = Lvarh; L@@Dh[3] = Lord;
+  def @Dh = ring(L@@Dh); kill L@@Dh;
+  setring @Dh;
+  dbprint(ppl,"// the ring @Dh:",@Dh);
+  // var(1)=t, var(2..n+1) = x(1..n), var(n+2)=Dt, var(n+3..2*n+2)=D(1..n),var(2*n+3)=h
+  // 2.4 create non-commutative relations
+  matrix @relD[N][N];
+  for (i=1; i<=n+1; i++) { @relD[i,n+1+i] = var(N)^(d+1); }
+  dbprint(ppl,"// nc relations:",@relD);
+  def Dh = nc_algebra(1,@relD);
+  setring Dh; kill @Dh;
+  dbprint(ppl,"// computing in ring",Dh);
+  // 3. compute the initial ideal
+  poly f = imap(save,f);
+  f = homog(f,h);
+  // 3.1 create the Malgrange ideal
+  ideal I = var(1)-f;
+  for (i=1; i<=n; i++)
+  {
+    I = I, var(n+2+i)+diff(f,var(i+1))*var(n+2);
+  }
+  // 3.2 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 is " +string(I));
+  I = subst(I,var(N),1); // dehomogenization
+  dbprint(ppl, "// I after dehomogenization is " +string(I));
+  // 3.3 the initial form
+  I = inForm(I,uv);
+  // 3.4 create Weyl algebra
+  setring save;
+  Lord = list();
+  Lord[1] = list("dp",intvec(1:(N-1)));
+  Lord[2] = C0;
+  list L@@D = RL;
+  L@@D[2] = Lvar; L@@D[3] = Lord;
+  def @D = ring(L@@D); kill L@@D;
+  setring @D; def D = Weyl(); setring D;
+  ideal I = imap(Dh,I);
+  kill @D,Dh;
+  // 3.5 the final GB
+  dbprint(ppl, "// starting cosmetic Groebner basis computation with engine: "+string(whichengine));
+  I = engine(I, whichengine);
+  dbprint(ppl,"// finished cosmetic Groebner basis computation");
+  ideal inF = I; attrib(inF,"isSB",1);
+  export(inF);
+  return(D);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring r = 0,(x,y),dp;
+  poly f = x^2+y^3+x*y^2;
+  def D = initialMalgrange(f);
+  setring D;
+  inF;
+  setring r;
+  intvec v = 3,2;
+  def D2 = initialMalgrange(f,1,0,1,v);
+  setring D2;
+  inF;
+}
+
+static 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)) ) )
+  {
+    return(1);
+  }
+  return(0);
+}
+
+proc bFactor (poly F)
+"USAGE:  bFactor(f);  f poly
+RETURN:  list
+PURPOSE: computes the roots of irreducible factors of an univariate poly
+NOTE:    The output list consists of two or three entries:
+@*       the roots of f as ideal, their multiplicities as intvec, and,
+@*       if present, a third one being the product of all irreducible factors
+@*       of degree greater than one, given as string.
+DISPLAY: If printlevel=1, progress debug messages will be printed,
+@*       if printlevel>=2, all the debug messages will be printed.
+EXAMPLE: example bFactor; shows examples
+"
+{
+  int ppl = printlevel - voice +2;
+  def save = basering;
+  list L = variables(F);
+  int i = size(L);
+  if (i>1) { ERROR("poly has to be univariate")}
+  if (i == 0)
+  {
+    dbprint(ppl,"// poly is constant");
+    L = list(ideal(0),intvec(0),string(F));
+    return(L);
+  }
+  poly v = L[1];
+  L = ringlist(save); L = L[1..4];
+  L[2] = list(string(v));
+  L[3] = list(list("dp",intvec(1)),list("C",intvec(0)));
+  def @S = ring(L);
+  setring @S;
+  poly ir = 1; poly F = imap(save,F);
+  list L = factorize(F);
+  ideal I = L[1];
+  intvec m = L[2];
+  ideal II; intvec mm;
+  for (i=2; i<=ncols(I); i++)
+  {
+    if (deg(I[i]) > 1)
+    {
+      ir = ir * I[i]^m[i];
+    }
+    else
+    {
+      II[size(II)+1] = I[i];
+      mm[size(II)] = m[i];
+    }
+  }
+  II = normalize(II);
+  II = subst(II,var(1),0);
+  II = -II;
+  if (size(II)>0)
+  {
+    dbprint(ppl,"// found roots");
+    dbprint(ppl-1,string(II));
+  }
+  else
+  {
+    dbprint(ppl,"// found no roots");
+  }  
+  L = list(II,mm);
+  if (ir <> 1)
+  {
+    dbprint(ppl,"// found irreducible factors");
+    ir = cleardenom(ir);
+    dbprint(ppl-1,string(ir));
+    L = L + list(string(ir));
+  }
+  else 
+  {
+    dbprint(ppl,"// no irreducible factors found");
+  }  
+  setring save;
+  L = imap(@S,L);
+  return(L);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring r = 0,(x,y),dp;
+  bFactor((x^2-1)^2);
+  bFactor((x^2+1)^2);  
+  bFactor((y^2+1/2)*(y+9)*(y-7));
+}