Changeset 6cd6983 in git

Timestamp:

Dec 1, 2008, 9:58:20 PM (15 years ago)

Author:

Viktor Levandovskyy <levandov@…>

Branches:

(u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')

Children:

a6a5d9ad89ba42bbc129bacba382bb6f463aae65

Parents:

8c66fb0e0dc173c2eb5aff1e988308a8527bc045

Message:

*levandov: improvements, bugfixes, docu


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

File:

• Singular/LIB/bfct.lib (modified) (38 diffs)

Singular/LIB/bfct.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: bfct.lib,v 1.5 2008-10-06 17:04:26 Singular Exp $";
+version="$Id: bfct.lib,v 1.6 2008-12-01 20:58:20 levandov Exp $";
 category="Noncommutative";
 info="
@@ -18 +18 @@
 MAIN PROCEDURES:
 
-bfct(f[,s,t,v]);           compute the global Bernstein-Sato polynomial of a given poly
-bfctsyz(f[,r,s,t,u,v]);    compute the global Bernstein-Sato polynomial of a given poly
-bfctonestep(f[,s,t]);      compute the global Bernstein-Sato polynomial of a given poly
-bfctideal(I,w[,s,t]);      compute the global b-function of a given ideal w.r.t. a given weight
-minpol(f,I);               compute the minimal polynomial of the endormorphism in basering modulo ideal given by a poly
-minpolsyz(f,I[,p,s,t]);    compute the minimal polynomial of the endormorphism in basering modulo ideal given by a poly
-linreduce(f,I[,s]);        reduce a poly by linear reductions of its leading term
-ncsolve(I[,s]);            find and compute a linear dependency of the elements of an ideal
+bfct(f[,s,t,v]);            compute the global Bernstein-Sato polynomial of a given poly
+bfctsyz(f[,r,s,t,u,v]);     compute the global Bernstein-Sato polynomial of a given poly
+bfctann(f[,s]);             compute the global Bernstein-Sato polynomial of a given poly
+bfctonestep(f[,s,t]);       compute the global Bernstein-Sato polynomial of a given poly
+bfctideal(I,w[,s,t]);       compute the global b-function of a given ideal w.r.t. a given weight
+pintersect(f,I);            compute the intersection of the ideal generated by a given poly with a given ideal
+pintersectsyz(f,I[,p,s,t]); compute the intersection of the ideal generated by a given poly with a given ideal
+linreduce(f,I[,s]);         reduce a poly by linear reductions only
+ncsolve(I[,s]);             find and compute a linear dependency of the elements of an ideal
 
 AUXILIARY PROCEDURES:
@@ -32 +33 @@
 isin(l,i);       check whether an element is a member of a list 
 scalarprod(v,w); compute the standard scalar product of two intvecs
+vec2poly(v[,i]); convert a coefficient vector to a poly
 
 SEE ALSO: dmod_lib, dmodapp_lib, gmssing_lib
@@ -38 +40 @@
 
 LIB "qhmoduli.lib"; // for Max
-LIB "dmodapp.lib"; // for initialideal etc
+LIB "dmodapp.lib";  // for initialideal etc
 
 
@@ -48 +50 @@
   example isin;
   example scalarprod;
+  example vec2poly;
   "MAIN PROCEDURES:";
   example bfct;
   example bfctsyz;
+  example bfctann;
   example bfctonestep;
   example bfctideal;
-  example minpol;
-  example minpolsyz;
+  example pintersect;
+  example pintersectsyz;
   example linreduce;
   example ncsolve;
@@ -156 +160 @@
 proc isin (list l, i)
 "USAGE:  isin(l,i);  l a list, i an argument of any type
-RETURN:  1 if i is a member of l, or 0 otherwise
+RETURN:  an 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 isin; shows an example
@@ -166 +170 @@
     if (l[j]==i)
     {
-      return(1);
+      return(j);
       break;
     }
@@ -216 +220 @@
 
 proc linreduce(poly f, ideal I, list #)
-"USAGE:  linreduce(f, I [,s]);  f a poly, I an ideal, s an optional int
-RETURN:  a poly obtained by linear reductions of the leading term of the given poly with an ideal
-PURPOSE: reduce a poly only by linear reductions of its leading term
-NOTE:    If s<>0, a list consisting of the reduced poly and the vector of the used 
+"USAGE:  linreduce(f, I [,s,t]);  f a poly, I an ideal, s,t optional ints
+RETURN:  a poly obtained by linear reductions of the given poly with the given ideal
+PURPOSE: reduce a poly only by linear reductions
+NOTE:    If s<>0, a list consisting of the reduced poly and the coefficient vector of the used 
 @*       reductions is returned.
+@*       If t<>0, only leading monomials are reduced, otherwise, and by default, all monomials
+@*       are reduced, if possible.
 EXAMPLE: example linreduce; shows examples
 "
 {
+  int ppl = printlevel - voice + 2;
   int remembercoeffs = 0; // default
+  int redlm          = 0; // default
   if (size(#)>0)
   {
@@ -231 +239 @@
       remembercoeffs = #[1];
     }
-  }
-  int i;
+    if (size(#)>1)
+    {
+      if (typeof(#[2])=="int" || typeof(#[2])=="number")
+      {
+        redlm = #[2];
+      }
+    }
+  }
+  int i,j,k;
   int sI = ncols(I);
   ideal lmI,lcI;
@@ -241 +256 @@
   }
   vector v;
-  poly lm,c;
-  int reduction;
-  lm = leadmonom(f);
-  reduction = 1;
-  while (reduction == 1) // while there was a reduction
-  {
-    reduction = 0;
-    for (i=sI;i>=1;i--)
-    {
-      if (lm <> 0 && lm == lmI[i])
-      {
-        c = leadcoef(f)/lcI[i];
-        f = f - c*I[i];
-        lm = leadmonom(f);
-        reduction = 1;
-        if (remembercoeffs <> 0)
+  poly c;
+  int reduction = 1;
+  if (redlm == 0)
+  {
+    ideal monomf;
+    for (k=1; k<=size(f); k++)
+    {
+      monomf[k] = normalize(f[k]);
+    }
+    while (reduction == 1) // while there was a reduction
+    {
+      reduction = 0;
+      for (i=sI; i>=1; i--)
+      {
+        dbprint(ppl,"testing ideal entry:",i);
+        for (j=1; j<=size(f); j++)
         {
-          v = v - c * gen(i);
+          if (monomf[j] == lmI[i])
+          {
+            c = leadcoef(f[j])/lcI[i];
+            f = f - c*I[i];
+            dbprint(ppl,"reducing poly to ",f);
+            monomf = 0;
+            for (k=1; k<=size(f); k++)
+            {
+              monomf[k] = normalize(f[k]);
+            }
+            reduction = 1;
+            if (remembercoeffs <> 0)
+            {
+              v = v - c * gen(i);
+            }
+            break;
+          }
+        }
+        if (reduction == 1)
+        {
+          break;
+        }
+      }
+    }
+  }
+  else // reduce only leading monomials
+  {
+    poly lm = leadmonom(f);
+    while (reduction == 1) // while there was a reduction
+    {
+      reduction = 0;
+      for (i=sI;i>=1;i--)
+      {
+        if (lm <> 0 && lm == lmI[i])
+        {
+          c = leadcoef(f)/lcI[i];
+          f = f - c*I[i];
+          lm = leadmonom(f);
+          reduction = 1;
+          if (remembercoeffs <> 0)
+          {
+            v = v - c * gen(i);
+          }
         }
       }
@@ -279 +336 @@
   ideal I = 1,y,xy;
   poly f = 5xy+7y+3;
-  poly g = 5y+7x+3;
-  linreduce(f,I);
+  poly g = 7x+5y+3;
   linreduce(g,I);
+  linreduce(g,I,0,1);
   linreduce(f,I,1);
 }
@@ -386 +443 @@
 }
 
-proc minpol (poly s, ideal I)
-"USAGE:  minpol(f, I);  f a poly, I an ideal
-RETURN:  coefficient vector of the minimal polynomial of the endomorphism of basering modulo I defined by f
-PURPOSE: compute the minimal polynomial
-NOTE:    If f does not define an endomorphism, this proc will not terminate.
+proc pintersect (poly s, ideal I)
+"USAGE:  pintersect(f, I);  f a poly, I an ideal
+RETURN:  coefficient vector of the monic generator of the intersection of the ideal generated by f with I
+PURPOSE: compute the intersection of an ideal with a principal ideal defined by f
+NOTE:    If the intersection is zero, this proc might not terminate.
 @*       I should be given as standard basis.
 @*       If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example minpol; shows examples
+EXAMPLE: example pintersect; shows examples
 "
 {
@@ -400 +457 @@
   attrib(I,"isSB",1); // set attribute for suppressing NF messages
   int ppl = printlevel-voice+2;
+  // ---case 1: I = basering---
+  if (size(I) == 1)
+  {
+    if (simplify(I,1) == ideal(1))
+    {
+      return(gen(2)); // = s
+    }
+  }
   def save = basering;
+  int n = nvars(save);
   int i,j,k;
+  // ---case 2: intersection is zero---
+  intvec degs = leadexp(s);
+  intvec possdegbounds;
+  list degI;
+  i = 1;
+  while (i <= ncols(I))
+  {
+    if (i == ncols(I)+1)
+    {
+      break;
+    }
+    degI[i] = leadexp(I[i]);
+    for (j=1; j<=n; j++)
+    {
+      if (degs[j] == 0)
+      {
+        if (degI[i][j] <> 0)
+        {
+          break;
+        }
+      }
+      if (j == n)
+      {
+        k++;
+        possdegbounds[k] = Max(degI[i]);
+      }
+    }
+    i++;
+  }
+  int degbound = Min(possdegbounds);
+  dbprint(ppl,"a lower bound for the degree of the insection is:");
+  dbprint(ppl,degbound);
+  if (degbound == 0) // lm(s) does not appear in lm(I)
+  {
+    return(vector(0));
+  }
+  // ---case 3: intersection is non-trivial---
+  ideal redNI = 1;
   vector v;
   list l,ll;
   l[1] = vector(0);
-  poly toNF, tobracket, newNF, rednewNF;
-  ideal NI = 1;
+  poly toNF,tobracket,newNF,rednewNF,oldNF,secNF;
   i = 1;
-  ideal redNI = 1;
-  newNF = NF(s,I);
-  dbprint(ppl+1,"minpol starts...");
-  dbprint(ppl+1,"with ideal I=", I);
+  dbprint(ppl+1,"pintersect starts...");
   while (1)
   {
@@ -417 +517 @@
     if (i>1)
     {
-      tobracket = s^(i-1)-NI[i];
-      if (tobracket==0) // bracket doesn't like zeros
+      oldNF = newNF;
+      tobracket = s^(i-1) - oldNF;
+      if (tobracket==0) // todo bug in bracket?
       {
         toNF = 0;
@@ -424 +525 @@
       else
       {
-        toNF = bracket(tobracket,NI[2]);
-      }
-      newNF = NF(toNF+NI[i]*NI[2],I);  // = NF(s^i,I)
+        toNF = bracket(tobracket,secNF);
+      }
+      newNF = NF(toNF+oldNF*secNF,I);  // = NF(s^i,I)
+    }
+    else
+    {
+      newNF = NF(s,I);
+      secNF = newNF;
     }
     ll = linreduce(newNF,redNI,1);
     rednewNF = ll[1];
     l[i+1] = ll[2];
-    dbprint(ppl+1,"newNF is:", newNF);
-    dbprint(ppl+1,"rednewNF is:", rednewNF);
+    dbprint(ppl,"newNF is:", newNF);
+    dbprint(ppl,"rednewNF is:", rednewNF);
     if (rednewNF != 0) // no linear dependency
     {
-      NI[i+1] = newNF;
       redNI[i+1] = rednewNF;
-      dbprint(ppl+1,"NI is:", NI);
-      dbprint(ppl+1,"redNI is:", redNI);
       i++;
     }
     else // there is a linear dependency, hence we are done
     {
-      dbprint(ppl+1,"the degree of the minimal polynomial is:", i);
+      dbprint(ppl+1,"the degree of the generator of the intersection is:", i);
       break;
     }
   }
   dbprint(ppl,"used linear reductions:", l);
-  // we obtain the coefficients of the minimal polynomial by the used reductions:
+  // we obtain the coefficients of the generator of the intersection by the used reductions:
   ring @R = 0,(a(1..i+1)),dp;
   setring @R;
@@ -461 +564 @@
     }
   }
-  for(j=i;j>=1;j--)
+  for (j=i;j>=1;j--)
   {
     C[i+1]= subst(C[i+1],a(j),a(j)+C[j]);
@@ -474 +577 @@
   v = imap(@R,v);
   kill @R;
-  dbprint(ppl+1,"minpol finished");
+  dbprint(ppl+1,"pintersect finished");
   return(v);
 }
@@ -480 +583 @@
 {
   "EXAMPLE:"; echo = 2;
-  printlevel = 0;
   ring r = 0,(x,y),dp;
   poly f = x^2+y^3+x*y^2;
@@ -486 +588 @@
   setring D;
   inF;
-  poly s = t*Dt;
-  minpol(s,inF);
-}
-
-proc minpolsyz (poly s, ideal II, list #)
-"USAGE:  minpolsyz(f, I [,p,s,t]);  f a poly, I an ideal, p, t optial ints, p a prime number
-RETURN:  coefficient vector of the minimal polynomial of the endomorphism of basering modulo I defined by f
-PURPOSE: compute the minimal polynomial
-NOTE:    If f does not define an endomorphism, this proc will not terminate.
+  pintersect(t*Dt,inF);
+}
+
+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
+RETURN:  coefficient vector of the monic generator of the intersection of the ideal generated by f with I
+PURPOSE: compute the intersection of an ideal with a principal ideal defined by f
+NOTE:    If the intersection is zero, this proc might not terminate.
 @*       I should be given as standard basis.
-@*       If p>0 is given, the proc computes the minimal polynomial in char p first and
-@*       then only searches for a minimal polynomial of the obtained degree in the basering.
+@*       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.
 @*       This is done by computing syzygies.
@@ -506 +607 @@
 @*       If printlevel=1, progress debug messages will be printed,
 @*       if printlevel>=2, all the debug messages will be printed.
-EXAMPLE: example minpolsyz; shows examples
+EXAMPLE: example pintersectsyz; shows examples
 "
 {
@@ -566 +667 @@
   }
   i = 1;
-  dbprint(ppl+1,"minpolynomial starts...");
+  dbprint(ppl+1,"pintersectsyz starts...");
   dbprint(ppl+1,"with ideal I=", I);
   while (1)
@@ -639 +740 @@
     i++;
   }
-  dbprint(ppl+1,"minpol finished");
+  dbprint(ppl+1,"pintersectsyz finished");
   return(v);
 }
@@ -645 +746 @@
 {
   "EXAMPLE:"; echo = 2;
-  printlevel = 0;
   ring r = 0,(x,y),dp;
   poly f = x^2+y^3+x*y^2;
@@ -652 +752 @@
   inF;
   poly s = t*Dt;
-  minpolsyz(s,inF);
+  pintersectsyz(s,inF);
   int p = prime(20000);
-  minpolsyz(s,inF,p,0,0);
+  pintersectsyz(s,inF,p,0,0);
+}
+
+proc vec2poly (list #)
+"USAGE:  vec2poly(v [,i]);  v a vector or an intvec, i an optional int
+RETURN:  a poly with coefficient vector v
+PURPOSE: convert a coefficient vector to a poly
+NOTE:    If i>0 is given, the returned poly is an element of K[var(i)], 
+@*       otherwise, and by default, @code{i=1} is used.
+@*       The first entry of v is the coefficient of 1.
+EXAMPLE: example vec2poly; shows examples
+"
+{
+  def save = basering;
+  int i,ringvar;
+  ringvar = 1; // default
+  if (size(#) > 0)
+  {
+    if (typeof(#[1])=="vector" || typeof(#[1])=="intvec")
+    {
+      def v = #[1];
+    }
+    else
+    {
+      ERROR("wrong input: expected vector/intvec expression");
+    }
+    if (size(#) > 1)
+    {
+      if (typeof(#[2])=="int" || typeof(#[2])=="number")
+      {
+        ringvar = int(#[2]);
+      }
+    }
+  }
+  if (ringvar > nvars(save))
+  {
+    ERROR("var out of range");
+  }
+  poly p;
+  for (i=1; i<=nrows(v); i++)
+  {
+    p = p + v[i]*(var(ringvar))^(i-1);
+  }
+  return(p);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring r = 0,(x,y),dp;
+  vector v = gen(1) + 3*gen(3) + 22/9*gen(4);
+  intvec iv = 3,2,1;
+  vec2poly(v,2);
+  vec2poly(iv);
 }
 
@@ -705 +857 @@
 }
 
-static proc bfctengine (poly f, int whichengine, int methodord, int methodminpol, int minpolchar, int modengine, intvec u0)
+static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodpintersect, int pintersectchar, int modengine, intvec u0)
 {
   int ppl = printlevel - voice +2;
@@ -711 +863 @@
   def save = basering;
   int n = nvars(save);
-  def DD = initialmalgrange(f,whichengine,methodord,1,u0);
-  setring DD;
-  ideal inI = inF;
-  kill inF;
-  poly s = t*Dt;
+  if (inorann == 0) // bfct using initial ideal
+  {
+    def D = initialmalgrange(f,whichengine,methodord,1,u0);
+    setring D;
+    ideal J = inF;
+    kill inF;
+    poly s = t*Dt;
+  }
+  else // bfct using Ann(f^s)
+  {
+    def D = SannfsBFCT(f,whichengine);
+    setring D;
+    ideal J = LD;
+    kill LD;
+  }
   vector b;
   // try it modular
-  if (methodminpol <> 0) // minpolsyz
-  {
-    if (minpolchar == 0) // minpolsyz::modular
+  if (methodpintersect <> 0) // pintersectsyz
+  {
+    if (pintersectchar == 0) // pintersectsyz::modular
     {
       int lb = 30000;
@@ -734 +896 @@
           usedprimes = usedprimes,q;
           dbprint(ppl,"used prime is: "+string(q));
-          b = minpolsyz(s,inI,q,whichengine,modengine);
+          b = pintersectsyz(s,J,q,whichengine,modengine);
         }
         i++;
       }
     }
-    else // minpolsyz::non-modular 
-    {
-      b = minpolsyz(s,inI,0,whichengine);
-    }
-  }
-  else // minpol: linreduce
-  {
-    b = minpol(s,inI);
+    else // pintersectsyz::non-modular 
+    {
+      b = pintersectsyz(s,J,0,whichengine);
+    }
+  }
+  else // pintersect: linreduce
+  {
+    b = pintersect(s,J);
   }
   setring save;
-  vector b = imap(DD,b);
-  list l = listofroots(b,1);
+  vector b = imap(D,b);
+  if (inorann == 0)
+  {
+    list l = listofroots(b,1);
+  }
+  else
+  {
+    list l = listofroots(b,0);
+  }
   return(l);
 }
@@ -758 +927 @@
 RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and their multiplicies
 PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, according to the algorithm by Masayuki Noro
-NOTE:    In this proc, a system of linear equations is solved by linear reductions.
+NOTE:    In this proc, the initial Malgrange ideal is computed.
+@*       Further, a system of linear equations is solved by linear reductions.
 @*       If s<>0, @code{std} is used for Groebner basis computations,
 @*       otherwise, and by default, @code{slimgb} is used. 
@@ -801 +971 @@
     }
   }
-  list b = bfctengine(f,whichengine,methodord,0,0,0,u0);
+  list b = bfctengine(f,0,whichengine,methodord,0,0,0,u0);
   return(b);
 }
@@ -818 +988 @@
 RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and its multiplicies
 PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, according to the algorithm by Masayuki Noro
-NOTE:    In this proc, a system of linear equations is solved by computing syzygies.
+NOTE:    In this proc, the initial Malgrange ideal is computed.
+@*       Further, a system of linear equations is solved by computing syzygies.
 @*       If r<>0, @code{std} is used for Groebner basis computations in characteristic 0,
 @*       otherwise, and by default, @code{slimgb} is used. 
 @*       If s<>0, a matrix ordering is used for Groebner basis computations,
 @*       otherwise, and by default, a block ordering is used.
-@*       If t<>0, the minimal polynomial computation is solely performed over charasteristic 0,
+@*       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 Groebner basis computations in characteristic >0,
@@ -839 +1010 @@
   // one for the engine used for Groebner basis computations in char 0,
   // one for  M() ordering or its realization
-  // one for a modular method when computing the minimal polynomial
+  // one for a modular method when computing the intersection
   // and one for the engine used for Groebner basis computations in char >0
   // in # can also be the optional weight vector
@@ -846 +1017 @@
   int whichengine = 0; // default
   int methodord   = 0; // default
-  int minpolchar  = 0; // default
+  int pintersectchar  = 0; // default
   int modengine   = 1; // default
   intvec u0       = 0; // default
@@ -865 +1036 @@
         if (typeof(#[3])=="int" || typeof(#[3])=="number")
         {
-          minpolchar = int(#[3]);
+          pintersectchar = int(#[3]);
         }
         if (size(#)>3)
@@ -884 +1055 @@
     }
   }
-  list b = bfctengine(f,whichengine,methodord,1,minpolchar,modengine,u0);
+  list b = bfctengine(f,0,whichengine,methodord,1,pintersectchar,modengine,u0);
   return(b);
 }
@@ -938 +1109 @@
     s = s + w[i]*var(i)*var(n+i);
   }
-  vector b = minpol(s,J);
+  vector b = pintersect(s,J);
   list l = listofroots(b,0);
   return(l);
@@ -994 +1165 @@
   setring DDh;
   dbprint(ppl, "the initial ideal:", string(matrix(inF)));
-  inF = nselect(inF,3..2*n+4);
-  inF = nselect(inF,1);
+  intvec tonselect = 1;
+  for (i=3; i<=2*n+4; i++)
+  {
+    tonselect = tonselect,i;
+  }
+  inF = nselect(inF,tonselect);
   dbprint(ppl, "generators containing only s:", string(matrix(inF)));
   inF = engine(inF, whichengine); // is now a principal ideal;
@@ -1013 +1188 @@
   bfctonestep(f);
   bfctonestep(f,1,1);
+}
+
+proc bfctann (poly f, list #)
+"USAGE:  bfctann(f [,r]);  f a poly, r an optional int
+RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and their multiplicies
+PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f
+NOTE:    In this proc, ann(f^s) is computed.
+@*       If r<>0, @code{std} is used for Groebner basis computations,
+@*       otherwise, and by default, @code{slimgb} is used. 
+@*       If printlevel=1, progress debug messages will be printed,
+@*       if printlevel>=2, all the debug messages will be printed.
+EXAMPLE: example bfctann; shows examples
+"
+{
+  def save = basering;
+  int ppl = printlevel - voice + 2;
+  int whichengine = 0; // default
+  if (size(#)>0)
+  {
+    if (typeof(#[1])=="int" || typeof(#[1])=="number")
+    {
+      whichengine = int(#[1]);
+    }
+  }
+  list b = bfctengine(f,1,whichengine,0,1,0,0,0);
+  return(b);
+}
+example
+{
+  "EXAMPLE:"; echo = 2;
+  ring r = 0,(x,y),dp;
+  poly f = x^2+y^3+x*y^2;
+  bfctann(f);
 }
 
@@ -1037 +1245 @@
   bfct(xyzreiffen45);
 }
+