Index: Singular/LIB/bfct.lib
===================================================================
--- Singular/LIB/bfct.lib	(revision 7f3ad4194490497ddcc8b5400f0f38119943ffb9)
+++ Singular/LIB/bfct.lib	(revision 2653d3c2809f0b96d6da699e54066fdb07623b92)
@@ -1,4 +1,4 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: bfct.lib,v 1.11 2009-01-14 16:07:03 Singular Exp $";
+version="$Id: bfct.lib,v 1.12 2009-01-15 19:19:09 levandov Exp $";
 category="Noncommutative";
 info="
@@ -16,5 +16,5 @@
 @*      - the list of its roots, which are known to be rational
 @*      - the multiplicities of the roots.
-@*      References: Saito, Strurmfels, Takayama: Groebner Deformations of Hypergeometric
+@*      References: Saito, Sturmfels, Takayama: Groebner Deformations of Hypergeometric
 @*      Differential Equations (2000), Noro: An Efficient Modular Algorithm for Computing
 @*      the Global b-function, (2002).
@@ -29,11 +29,12 @@
 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
+pIntersectSyz(f,I[,p,s,t]); intersection of the ideal I with the subalgebra K[f] for a 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
 
 AUXILIARY PROCEDURES:
 
-allPositive(v);  checks whether all entries of an intvec are positive
+allPositive(v);  checks whether all entries of an intvec are positive 
 scalarProd(v,w); computes the standard scalar product of two intvecs
 vec2poly(v[,i]); constructs an univariate poly with given coefficients
@@ -79,5 +80,5 @@
 EXAMPLE: example gradedWeyl; shows examples
 NOTE:    u[i] is the weight of x(i), v[i] the weight of D(i).
-@*       u+v has to be a non-negative intvec.
+@*       u+v has to be a non-negative intvec. 
 "
 {
@@ -88,5 +89,5 @@
   {
     ERROR("weight vectors have wrong dimension");
-  }
+  } 
   intvec uv,gr;
   uv = u+v;
@@ -97,9 +98,9 @@
       if (uv[i]==0)
       {
-        gr[i] = 0;
+	gr[i] = 0;
       }
       else
       {
-        gr[i] = 1;
+	gr[i] = 1;
       }
     }
@@ -114,5 +115,5 @@
   for (i=1; i<=n; i++)
   {
-    if (gr[i] == 1)
+    if (gr[i] == 1)  
     {
        l2[n+i] = "xi("+string(i)+")";
@@ -229,7 +230,7 @@
 proc linReduceIdeal(ideal I, list #)
 "USAGE:  linReduceIdeal(I [,s,t]); I an ideal, s,t optional ints
-RETURN:  ideal/list, linear reductum (over field) of f by elements from I
-PURPOSE: reduce a poly only by linear reductions (no monomial multiplications)
-NOTE:    If s<>0, a list consisting of the reduced ideal and the coefficient
+RETURN:  ideal/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 
 @*       vectors of the used reductions given as module is returned.
 @*       Otherwise (and by default), only the reduced ideal is returned.
@@ -255,5 +256,5 @@
       if (typeof(#[2])=="int" || typeof(#[2])=="number")
       {
-        redtail = #[2];
+	redtail = #[2];
       }
     }
@@ -273,15 +274,15 @@
       for (i=1; i<=sI; i++)
       {
-        if (I[i] == 0)
-        {
-          j++;
-          J[j] = 0;
-          ordJ[j] = -1;
-          M[j] = gen(i);
-        }
-        else
-        {
-          M[i+sZeros-j] = gen(lJ[2][i-j]+j);
-        }
+	if (I[i] == 0)
+	{
+	  j++;
+	  J[j] = 0;
+	  ordJ[j] = -1;
+	  M[j] = gen(i);
+	}
+	else
+	{
+	  M[i+sZeros-j] = gen(lJ[2][i-j]+j);
+	}
       }
     }
@@ -290,6 +291,6 @@
       for (i=1; i<=sZeros; i++)
       {
-        J[i] = 0;
-        ordJ[i] = -1;
+	J[i] = 0;
+	ordJ[i] = -1;
       }
     }
@@ -327,22 +328,22 @@
       for (j=sZeros+1; j<i; j++)
       {
-        if (lm == 0) { break; }
-        if (ordlm > maxordJ) { break; }
-        if (ordlm == ordJ[j])
-        {
-          if (lm > maxlmJ) { break; }
-          if (lm == lmJ[j])
-          {
-            dbprint(ppl,"reducing " + string(redpoly));
-            dbprint(ppl," with " + string(J[j]));
-            c = leadcoef(redpoly)/leadcoef(J[j]);
-            redpoly = redpoly - c*J[j];
-            dbprint(ppl," to " + string(redpoly));
-            lm = leadmonom(redpoly);
-            ordlm = ord(lm);
-            if (remembercoeffs <> 0) { M[i] = M[i] - c * M[j]; }
-            reduction = 1;
-          }
-        }
+	if (lm == 0) { break; }
+	if (ordlm > maxordJ) { break; }
+	if (ordlm == ordJ[j])
+	{	    
+	  if (lm > maxlmJ) { break; }
+	  if (lm == lmJ[j])
+	  {
+	    dbprint(ppl,"reducing " + string(redpoly));
+	    dbprint(ppl," with " + string(J[j]));
+	    c = leadcoef(redpoly)/leadcoef(J[j]);
+	    redpoly = redpoly - c*J[j];
+	    dbprint(ppl," to " + string(redpoly));
+	    lm = leadmonom(redpoly);
+	    ordlm = ord(lm);
+	    if (remembercoeffs <> 0) { M[i] = M[i] - c * M[j]; }
+	    reduction = 1;
+	  }
+	}
       }
     }
@@ -354,5 +355,5 @@
     lmJ = insertGenerator(lmJ,lm,j);
     ordJ = insertGenerator(ordJ,poly(ordlm),j);
-    if (remembercoeffs <> 0)
+    if (remembercoeffs <> 0) 
     {
       v = M[i];
@@ -373,19 +374,19 @@
       for (j=i+1; j<=sI; j++)
       {
-        for (k=2; k<=size(J[j]); k++) // run over all terms in J[j]
-        {
-          if (ordJ[i] == ord(J[j][k]))
-          {
-            if (lm == normalize(J[j][k]))
-            {
-              c = leadcoef(J[j][k])/leadcoef(J[i]);
-              dbprint(ppl,"reducing " + string(J[j]));
-              dbprint(ppl," with " + string(J[i]));
-              J[j] = J[j] - c*J[i];
-              dbprint(ppl," to " + string(J[j]));
-              if (remembercoeffs <> 0) { M[j] = M[j] - c * M[i]; }
-            }
-          }
-        }
+	for (k=2; k<=size(J[j]); k++) // run over all terms in J[j]
+	{
+	  if (ordJ[i] == ord(J[j][k]))
+	  {    
+	    if (lm == normalize(J[j][k]))
+	    {
+	      c = leadcoef(J[j][k])/leadcoef(J[i]);
+	      dbprint(ppl,"reducing " + string(J[j]));
+	      dbprint(ppl," with " + string(J[i]));
+	      J[j] = J[j] - c*J[i];
+	      dbprint(ppl," to " + string(J[j]));
+	      if (remembercoeffs <> 0) { M[j] = M[j] - c * M[i]; }
+	    }
+	  }
+	}
       }
     }
@@ -410,5 +411,5 @@
 RETURN:  poly/list, linear reductum (over field) of f by elements from I
 PURPOSE: reduce a poly only by linear reductions (no monomial multiplications)
-NOTE:    If s<>0, a list consisting of the reduced poly and the coefficient
+NOTE:    If s<>0, a list consisting of the reduced poly and the coefficient 
 @*       vector of the used reductions is returned, otherwise (and by default)
 @*       only reduced poly is returned.
@@ -433,12 +434,12 @@
       if (typeof(#[2])=="int" || typeof(#[2])=="number")
       {
-        redtail = #[2];
+	redtail = #[2];
       }
       if (size(#)>2)
       {
-        if (typeof(#[3])=="int" || typeof(#[3])=="number")
+	if (typeof(#[3])=="int" || typeof(#[3])=="number")
         {
-          prepareideal = #[3];
-        }
+	  prepareideal = #[3];
+	}
       }
     }
@@ -469,5 +470,5 @@
       for (i=1; i<=sI; i++)
       {
-        M[i] = gen(i);
+	M[i] = gen(i);
       }
     }
@@ -495,15 +496,15 @@
       if (ordf == ordI[i])
       {
-        if (lm == lmI[i])
-        {
-          c = leadcoef(f)/lcI[i];
-          f = f - c*I[i];
-          lm = leadmonom(f);
-          ordf = ord(lm);
-          if (remembercoeffs <> 0)
-          {
-            v = v - c * M[i];
-          }
-        }
+	if (lm == lmI[i])
+	{
+	  c = leadcoef(f)/lcI[i];
+	  f = f - c*I[i];
+	  lm = leadmonom(f);
+	  ordf = ord(lm);
+	  if (remembercoeffs <> 0)
+	  {
+	    v = v - c * M[i];
+	  }
+	}
       }
     }
@@ -518,17 +519,17 @@
       for (j=1; j<=size(f); j++)
       {
-        if (ord(f[j]) == ordI[i])
-        {
-          if (normalize(f[j]) == lmI[i])
-          {
-            c = leadcoef(f[j])/lcI[i];
-            f = f - c*I[i];
-            dbprint(ppl,"reducing poly to ",f);
-            if (remembercoeffs <> 0)
-            {
-              v = v - c * M[i];
-            }
-          }
-        }
+	if (ord(f[j]) == ordI[i])
+	{
+	  if (normalize(f[j]) == lmI[i])
+	  {
+	    c = leadcoef(f[j])/lcI[i];
+	    f = f - c*I[i];
+	    dbprint(ppl,"reducing poly to ",f);
+	    if (remembercoeffs <> 0)
+	    {
+	      v = v - c * M[i];
+	    }
+	  }
+	}
       }
     }
@@ -556,5 +557,5 @@
   f = x3 + y2 + x2 + y + x;
   I = x3 - y3, y3 - x2, x3 - y2, x2 - y, y2-x;
-  list l = linReduce(f, I, 1);
+  list l = linReduce(f, I, 1); 
   l;
   module M = I;
@@ -601,5 +602,5 @@
       if (p <> 0)
       {
-        whichengine = 1;
+	whichengine = 1;
       }
     }
@@ -668,5 +669,5 @@
 proc pIntersect (poly s, ideal I)
 "USAGE:  pIntersect(f, I);  f a poly, I an ideal
-RETURN:  vector, coefficient vector of the monic polynomial
+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.
@@ -684,5 +685,5 @@
     attrib(I,"isSB",1); // set attribute for suppressing NF messages
   }
-  int ppl = printlevel-voice+2;
+  int ppl = printlevel-voice+1;
   // ---case 1: I = basering---
   if (size(I) == 1)
@@ -712,13 +713,13 @@
       if (degs[j] == 0)
       {
-        if (degI[i][j] <> 0)
-        {
-          break;
-        }
+	if (degI[i][j] <> 0)
+	{
+	  break;
+	}
       }
       if (j == n)
       {
-        k++;
-        possdegbounds[k] = Max(degI[i]);
+	k++;
+	possdegbounds[k] = Max(degI[i]);
       }
     }
@@ -749,9 +750,9 @@
       if (tobracket==0) // todo bug in bracket?
       {
-        toNF = 0;
+	toNF = 0;
       }
       else
       {
-        toNF = bracket(tobracket,secNF);
+	toNF = bracket(tobracket,secNF);
       }
       newNF = NF(toNF+oldNF*secNF,I);  // = NF(s^i,I)
@@ -802,5 +803,5 @@
     v = v + m[j,1]*gen(j);
   }
-  setring save;
+  setring save; 
   v = imap(@R,v);
   kill @R;
@@ -864,8 +865,8 @@
       if (size(#)>2)
       {
-        if (typeof(#[3])=="int" || typeof(#[3])=="number")
-        {
-          modengine = int(#[3]);
-        }
+	if (typeof(#[3])=="int" || typeof(#[3])=="number")
+	{
+	  modengine = int(#[3]);
+	}
       }
     }
@@ -909,9 +910,9 @@
       if (tobracket!=0)
       {
-        toNF = bracket(tobracket,NI[2]) + NI[i]*NI[2];
+	toNF = bracket(tobracket,NI[2]) + NI[i]*NI[2];
       }
       else
       {
-        toNF = NI[i]*NI[2];
+	toNF = NI[i]*NI[2];
       }
       newNF =  NF(toNF,I);
@@ -927,16 +928,16 @@
       if (v!=0) // there is a modular solution
       {
-        dbprint(ppl,"got solution in char ",solveincharp," of degree " ,i);
-        setring save;
-        v = linSyzSolve(NI,whichengine);
-        if (v==0)
-        {
-          break;
-        }
+	dbprint(ppl,"got solution in char ",solveincharp," of degree " ,i);
+	setring save;
+	v = linSyzSolve(NI,whichengine);
+	if (v==0)
+	{
+	  break;
+	}
       }
       else // no modular solution
       {
-        setring save;
-        v = 0;
+	setring save;
+	v = 0;
       }
     }
@@ -965,5 +966,5 @@
       {
         dbprint(ppl,"linSyzSolve: got solution!");
-        // "got solution!";
+	// "got solution!";
         break;
       }
@@ -1015,5 +1016,5 @@
       if (typeof(#[2])=="int" || typeof(#[2])=="number")
       {
-        ringvar = int(#[2]);
+	ringvar = int(#[2]);
       }
     }
@@ -1094,4 +1095,8 @@
   def save = basering;
   int n = nvars(save);
+  if (size(variables(f)) == 0) // f is constant
+  {
+    return(list(ideal(0),intvec(0)));
+  }  
   if (inorann == 0) // bfct using initial ideal
   {
@@ -1121,16 +1126,16 @@
       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
+	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);
-        }
-        i++;
-      }
-    }
-    else // pIntersectSyz::non-modular
+	  usedprimes = usedprimes,q;
+	  dbprint(ppl,"used prime is: "+string(q));
+	  b = pIntersectSyz(s,J,q,whichengine,modengine);
+	}
+	i++;
+      }
+    }
+    else // pIntersectSyz::non-modular 
     {
       b = pIntersectSyz(s,J,0,whichengine);
@@ -1158,15 +1163,15 @@
 RETURN:  list of ideal and intvec
 PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
-@*       for the hypersurface defined by f.
+@*       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
+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
+NOTE:    In the output list, the ideal contains all the roots 
 @*       and the intvec their multiplicities.
 @*       If s<>0, @code{std} is used for GB computations,
-@*       otherwise, and by default, @code{slimgb} is used.
+@*       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,
+@*       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,
@@ -1199,8 +1204,8 @@
       if (size(#)>2)
       {
-        if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1)
+	if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1)
         {
-          u0 = #[3];
-        }
+	  u0 = #[3];
+	}
       }
     }
@@ -1222,24 +1227,24 @@
 "USAGE:  bfctSyz(f [,r,s,t,u,v]);  f a poly, r,s,t,u optional ints, v an optional 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
+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
+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,
+@*       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,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example bfct; shows examples
+EXAMPLE: example bfctSyz; shows examples
 "
 {
@@ -1275,20 +1280,20 @@
         if (typeof(#[3])=="int" || typeof(#[3])=="number")
         {
-          pIntersectchar = int(#[3]);
-        }
-        if (size(#)>3)
+	  pIntersectchar = int(#[3]);
+	}
+	if (size(#)>3)
         {
-          if (typeof(#[4])=="int" || typeof(#[4])=="number")
+	  if (typeof(#[4])=="int" || typeof(#[4])=="number")
           {
-            modengine = int(#[4]);
-          }
-          if (size(#)>4)
-          {
-            if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1)
-            {
-              u0 = #[5];
-            }
-          }
-        }
+	    modengine = int(#[4]);
+	  }
+	  if (size(#)>4)
+	  {
+	    if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1)
+	    {
+	      u0 = #[5];
+	    }
+	  }
+	}
       }
     }
@@ -1315,13 +1320,13 @@
 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
+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.
+@*       otherwise, and by default, @code{slimgb} is used. 
 @*       If t<>0, a matrix ordering is used for GB computations, otherwise,
 @*       and by default, a block ordering is used.
 DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example bfctideal; shows examples
+EXAMPLE: example bfctIdeal; shows examples
 "
 {
@@ -1374,13 +1379,13 @@
 "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
+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
+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.
+@*       and by default, @code{slimgb} is used. 
 @*       If t<>0, a matrix ordering is used for GB computations,
 @*       otherwise, and by default, a block ordering is used.
@@ -1492,16 +1497,16 @@
 "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)
+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, 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
+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.
+@*       otherwise, and by default, @code{slimgb} is used. 
 DISPLAY: If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example bfctann; shows examples
+EXAMPLE: example bfctAnn; shows examples
 "
 {
