Changeset eda414 in git for Singular/LIB/bfun.lib


Ignore:
Timestamp:
Mar 10, 2009, 5:12:52 PM (15 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
af9f1bc5804a0d34be405baa2e1604f18c96b4d4
Parents:
094f80e33c853c56c092642444a1dbea668c09f8
Message:
*levandov: more docu enhancements, small bugfixes


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/bfun.lib

    r094f80 reda414  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="$Id: bfun.lib,v 1.4 2009-03-09 18:34:51 levandov Exp $";
     2version="$Id: bfun.lib,v 1.5 2009-03-10 16:12:52 levandov Exp $";
    33category="Noncommutative";
    44info="
    … …  
    88
    99THEORY: Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R,
    10 @*      one is interested in the global b-Function (also known as Bernstein-Sato
    11 @*      polynomial) b(s) in K[s], defined to be the monic polynomial, satisfying
    12 @*      a functional identity    L * F^{s+1} = b(s) F^s,   
    13 @*      for some operator L in D[s]. Here, D stands for an n-th Weyl algebra
     10@*      one is interested in the global b-function (also known as Bernstein-Sato
     11@*      polynomial) b(s) in K[s], defined to be the monic polynomial of minimal
     12@*      degree, satisfying a functional identity L * F^{s+1} = b(s) F^s,   
     13@*      for some operator L in D[s] (* stands for the action of differential operator)
     14@*      By D one denotes the n-th Weyl algebra
    1415@*      K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>.
    1516@*      One is interested in the following data:
    … …  
    1819@*      - the multiplicities of the roots.
    1920@*
    20 @*    There is a general definition of a b-function of a holonomic ideal [SST]
    21 @*    (that is, an ideal I in a Weyl algebra D, such that D/I is holonomic module)
    22 @*    with respect to the given weight vector. B-function b_w(s) is in K[s].
    23 @*    Unlike Bernstein-Sato polynomial, general B-function with respect to
    24 @*    arbitrary weights need not have rational roots at all.
     21@*   There is a general definition of a b-function of a holonomic ideal [SST]
     22@*   (that is, an ideal I in a Weyl algebra D, such that D/I is holonomic module)
     23@*   with respect to the given weight vector w: For a poly p in D, its initial
     24@*   form w.r.t. (-w,w) is defined as the sum of all terms of p which have
     25@*   maximal weighted total degree where the weight of x_i is -w_i and the weight
     26@*   of d_i is w_i. Let J be the initial ideal of I w.r.t. (-w,w), i.e. the
     27@*   K-vector space generated by all initial forms w.r.t (-w,w) of elements of I.
     28@*   Put s = w_1 x_1 d_1 + ... + w_n x_n d_n. Then the monic generator b_w(s) of
     29@*   the intersection of J with the PID K[s] is called the b-function of I w.r.t.  w.
     30@*   Unlike Bernstein-Sato polynomial, general b-function with respect to
     31@*   arbitrary weights need not have rational roots at all. However, b-function
     32@*  of a holonomic ideal is known to be non-zero as well.
    2533@*
    2634@*      References: [SST] Saito, Sturmfels, Takayama: Groebner Deformations of
    … …  
    3240MAIN PROCEDURES:
    3341
    34 bfct(f[,s,t,v]);            computes the Bernstein-Sato polynomial of poly f
    35 bfctSyz(f[,r,s,t,u,v]);     computes the Bernstein-Sato polynomial of poly f
    36 bfctAnn(f[,s]);             computes the Bernstein-Sato polynomial of poly f
    37 bfctOneGB(f[,s,t]);         computes the Bernstein-Sato polynomial of poly f
    38 bfctIdeal(I,w[,s,t]);       computes the global b-function of ideal wrt weight
    39 pIntersect(f,I[,s]);        intersection of ideal with subalgebra K[f] for poly f
    40 pIntersectSyz(f,I[,p,s,t]); intersection of ideal with subalgebra K[f] for poly f
    41 linReduce(f,I[,s]);         reduces a poly by linear reductions only
    42 linReduceIdeal(I,[s,t]);    reduces generators of ideal by linear reductions only
    43 linSyzSolve(I[,s]);         computes a linear dependency of elements of ideal I
     42bfct(f[,s,t,v]);            compute the BS polynomial of f with linear reductions
     43bfctSyz(f[,r,s,t,u,v]);  compute the BS polynomial of f with syzygy-solver
     44bfctAnn(f[,s]);           compute the BS polynomial of f via Ann f^s + f
     45bfctOneGB(f[,s,t]);     compute the BS polynomial of f by just one GB computation
     46bfctIdeal(I,w[,s,t]);     compute the b-function of ideal wrt weight
     47pIntersect(f,I[,s]);      intersection of ideal with subalgebra K[f] for a poly f
     48pIntersectSyz(f,I[,p,s,t]); intersection of ideal with subalgebra K[f] with syz-solver
     49linReduce(f,I[,s]);         reduce a poly by linear reductions wrt ideal
     50linReduceIdeal(I,[s,t]);  interreduce generators of ideal by linear reductions
     51linSyzSolve(I[,s]);         compute a linear dependency of elements of ideal I
    4452
    4553AUXILIARY PROCEDURES:
    … …  
    911919@*       If p>0 is given, this proc computes the generator of the intersection in
    912920@*       char p first and then only searches for a generator of the obtained
    913 @*       degree in the basering. Otherwise, it searched for all degrees by
    914 @*       computing syzygies. Note that the basering must be free of parameters to
    915 @*       use this modular method.
     921@*       degree in the basering. Otherwise, it searches for all degrees by
     922@*       computing syzygies.
    916923@*       If s<>0, @code{std} is used for Groebner basis computations in char 0,
    917924@*       otherwise, and by default, @code{slimgb} is used.
    … …  
    935942  int solveincharp = 0; // default
    936943  def save = basering;
     944  int n = nvars(save);
    937945  if (size(#)>0)
    938946  {
    … …  
    962970  newNF = NF(s,I);
    963971  NI[2] = newNF;
    964   list l = ringlist(save);
    965   if (typeof(l[1]) == "int") { l[1] = solveincharp; }
    966   else // parameters in basering: problems with mappings from Q(a_i)->Z/p(a_i)
    967   {
    968     if (solveincharp <> 0)
    969     {
    970       print("// basering must be free of parameters to use modular method");
    971       print("// solving in characteristic of basering only...");
    972       solveincharp = 0;
    973     }
    974   }
    975   if (solveincharp<>0)
    976   {
    977     if (size(l)>4)
    978     {
    979       matrix l5 = l[5];
    980       matrix l6 = l[6];
    981     }
    982     def @Rp = ring(l);
     972  list RL = ringlist(save);
     973  if (solveincharp)
     974  {
     975    int psolveincharp = prime(solveincharp);
     976    if (solveincharp <> psolveincharp)
     977    {
     978      print("// " + string(solveincharp) + " is invalid characteristic of ground field.");
     979      print("// " + string(psolveincharp) + " is used.");
     980      solveincharp = psolveincharp;
     981      kill psolveincharp;
     982    }
     983    list RLp = RL[1..4];
     984    if (typeof(RL[1]) == "int") { RLp[1] = solveincharp; }
     985    else { RLp[1][1] = solveincharp; } // parameters
     986    def @Rp = ring(RLp);
    983987    setring @Rp;
    984     list l = ringlist(@Rp);
     988    number c;
    985989    setring save;
    986     if (size(l)>4)
    987     {
     990    int shortSave = short; // workaround for maps Q(a_i) -> Z/p(a_i)
     991    short = 0;
     992    string str;
     993    int badprime;
     994    i=1;
     995    while (badprime == 0 && i<=size(s)) // detect bad primes
     996    {
     997      str = string(denominator(leadcoef(s[i])));
     998      str = "c = " + str + ";";
    988999      setring @Rp;
    989       l[5] = fetch(save,l5);
    990       l[6] = fetch(save,l6);
    991     }
    992     setring @Rp;
    993     def Rp = ring(l);
     1000      execute(str);
     1001      if (c == 0) { badprime = 1; }
     1002      setring save;
     1003      i++;
     1004    }
     1005    str = "poly s = " + string(s) + ";";
     1006    if (size(RL) > 4) // basering is NC-algebra
     1007    {
     1008      string RL5 = "@C = " + string(RL[5]) + ";";
     1009      string RL6 = "@D = " + string(RL[6]) + ";";
     1010      setring @Rp;
     1011      matrix @C[n][n]; matrix @D[n][n];
     1012      execute(RL5); execute(RL6);
     1013      def Rp = nc_algebra(@C,@D);
     1014    }
     1015    else { def Rp = @Rp; }
    9941016    setring Rp;
    9951017    kill @Rp;
    9961018    dbprint(ppl-1,"// solving in ring ", Rp);
     1019    execute(str);
    9971020    vector v;
    998     map phi = save,maxideal(1);
    999     poly s = phi(s);
     1021    number c;
    10001022    ideal NI = 1;
    10011023    setring save;
     1024    i=1;
     1025    while (badprime == 0 && i<=size(I)) // detect bad primes
     1026    {
     1027      str = string(leadcoef(cleardenom(I[i])));
     1028      str = "c = " + str + ";";
     1029      setring Rp;
     1030      execute(str);
     1031      if (c == 0) { badprime = 1; }
     1032      setring save;
     1033      i++;
     1034    }
     1035    if (badprime == 1)
     1036    {
     1037      print("// WARNING: bad prime");
     1038      short = shortSave;
     1039      return(vector(0));
     1040    }
    10021041  }
    10031042  i = 1;
    … …  
    10231062    // look for a solution
    10241063    dbprint(ppl-1,"// linSyzSolve starts with: "+string(matrix(NI)));
    1025     if (solveincharp<>0) // modular method
    1026     {
     1064    if (solveincharp) // modular method
     1065    {
     1066      for (j=1; j<=size(newNF); j++)
     1067      {
     1068        str = string(denominator(leadcoef(s[i])));
     1069        str = "c = " + str + ";";
     1070        setring Rp;
     1071        execute(str);
     1072        if (c == 0)
     1073        {
     1074          print("// WARNING: bad prime");
     1075          setring save;
     1076          short = shortSave;
     1077          return(vector(0));
     1078        }
     1079        setring save;
     1080      }
     1081      str = "NI[" + string(i) + "+1] = " + string(newNF) + ";";
    10271082      setring Rp;
    1028       NI[i+1] = phi(newNF);
     1083      execute(str); // NI[i+1] = [newNF]_{solveincharp}
    10291084      v = linSyzSolve(NI,modengine);
    10301085      if (v!=0) // there is a modular solution
    … …  
    10761131  }
    10771132  dbprint(ppl,"// pIntersectSyz finished");
     1133  if (solveincharp) { short = shortSave; } 
    10781134  return(v);
    10791135}
    … …  
    12401296    {
    12411297      list L = ringlist(D);
    1242       if (typeof(L[1])=="int")
    1243       { 
    1244         int lb = 10000;
    1245         int ub = 536870909;
    1246         i = 1;
    1247         list usedprimes;
    1248         matrix L5 = L[5]; matrix L6 = L[6];
    1249         int sJ = size(J); int sJq;
    1250         while (b == 0)
     1298      int lb = 10000; int ub = 536870909; // bounds for random primes
     1299      list usedprimes;
     1300      int sJ = size(J);
     1301      int sLJq;
     1302      ideal LJ;
     1303      for (i=1; i<=sJ; i++)
     1304      {
     1305        LJ[i] = leadcoef(cleardenom(J[i]));
     1306      }
     1307      int short_save = short; // workaround for map Q(a_i) -> Z/q(a_i)
     1308      short = 0;
     1309      string strLJq = "ideal LJq = " + string(LJ) + ";";
     1310      int nD = nvars(D);
     1311      string L5 = "matrix @C[nD][nD]; @C = " + string(L[5]) + ";";
     1312      string L6 = "matrix @D[nD][nD]; @D = " + string(L[6]) + ";";
     1313      L = L[1..4];
     1314      i = 1;
     1315      while (b == 0)
     1316      {
     1317        dbprint(ppl,"// number of run in the loop: "+string(i));
     1318        int q = prime(random(lb,ub));
     1319        if (findFirst(usedprimes,q)==0) // if q has not been used already
    12511320        {
    1252           dbprint(ppl,"// number of run in the loop: "+string(i));
    1253           int q = prime(random(lb,ub));
    1254           if (findFirst(usedprimes,q)==0) // if q was not already used
     1321          usedprimes = usedprimes,q;
     1322          dbprint(ppl,"// using prime: "+string(q));
     1323          if (typeof(L[1]) == "int") { L[1] = q; }
     1324          else { L[1][1] = q; } // parameters
     1325          def @Rq = ring(L); setring @Rq;
     1326          execute(L5); execute(L6);
     1327          def Rq = nc_algebra(@C,@D); // def Rq = nc_algebra(1,@D);
     1328          setring Rq; kill @Rq;
     1329          execute(strLJq);
     1330          sLJq = size(LJq);
     1331          setring D; kill Rq;
     1332          if (sLJq <> sJ ) // detect unlucky prime
    12551333          {
    1256             usedprimes = usedprimes,q;
    1257             dbprint(ppl,"// using prime: "+string(q));
    1258             L[1] = q;
    1259             def @Rq = ring(L); setring @Rq;
    1260             list lq = ringlist(@Rq);
    1261             lq[5] = fetch(D,L5); lq[6] = fetch(D,L6);
    1262             def Rq = ring(lq);
    1263             setring Rq; kill @Rq;
    1264             ideal Jq = fetch(D,J);
    1265             sJq = size(lead(Jq));
    1266             setring D; kill Rq;
    1267             if (sJq <> sJ)
    1268             {
    1269               dbprint(ppl,"// " +string(q) + " is unlucky");
    1270               b = 0;
    1271             }
    1272             else
    1273             {
    1274               b = pIntersectSyz(s,J,q,whichengine,modengine);
    1275             }
     1334            dbprint(ppl,"// " +string(q) + " is unlucky");
     1335            b = 0;
    12761336          }
    1277           i++;
     1337          else
     1338          {
     1339            b = pIntersectSyz(s,J,q,whichengine,modengine);
     1340          }
    12781341        }
    1279       }
    1280       else //parameters: problem when mapping Q(a_i) -> Z/p(a_i)
    1281       {
    1282         print("// basering must be free of parameters to use modular method");
    1283         print("// solving only in char 0...");
    1284         b = pIntersectSyz(s,J,0,whichengine);
    1285       }
    1286     } 
    1287     else // pIntersectSyz::non-modular else // parameters in basering: problems with mappings
     1342        i++;
     1343      }
     1344      short = short_save;
     1345    }
     1346    else // pIntersectSyz::non-modular
    12881347    {
    12891348      b = pIntersectSyz(s,J,0,whichengine);
    … …  
    13841443@*       If t<>0, the computation of the intersection is solely performed over
    13851444@*       charasteristic 0, otherwise and by default, a modular method is used.
    1386 @*       Note that in this case, the basering must be free of parameters.
    13871445@*       If u<>0 and by default, @code{std} is used for GB computations in
    13881446@*       characteristic >0, otherwise, @code{slimgb} is used.
    … …  
    14601518RETURN:  list of ideal and intvec
    14611519PURPOSE: computes the roots of the global b-function of I wrt the weight (-w,w).
    1462 ASSUME:  The basering is the n-th Weyl algebra in characteristic 0 and  for all
     1520ASSUME:  The basering is the n-th Weyl algebra in characteristic 0 and for all
    14631521@*       1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the
    14641522@*       sequence of variables is given by x(1),...,x(n),D(1),...,D(n),
    14651523@*       where D(i) is the differential operator belonging to x(i).
    1466 @*       Further assume that I is holonomic.
     1524@*       Further we assume that I is holonomic.
    14671525BACKGROUND:  In this proc, the initial ideal of I is computed according to the
    14681526@*       algorithm by Masayuki Noro and then a system of linear equations is
    14691527@*       solved by linear reductions.
    1470 NOTE:    In the output list, the ideal contains all the roots and the intvec
    1471 @*       their multiplicities.
     1528NOTE:    In the output list, say L,
     1529@*  - L[1] of type ideal contains all the rational roots of a b-function,
     1530@*  - L[2] of type intvec contains the multiplicities of above roots,
     1531@*  - optional L[3] of type string is the part of b-function without rational roots.
     1532@*  Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and
     1533@*  L[3] is 1 (for nonzero constant) or 0 (for zero b-function).
    14721534@*       If s<>0, @code{std} is used for GB computations in characteristic 0,
    14731535@*       otherwise, and by default, @code{slimgb} is used.
    … …  
    15101572  if (isH<>1)
    15111573  {
    1512     print("// given ideal is not holonomic");
    1513     print("// setting bound for degree of b-function to 10 and proceeding");
     1574    print("WARNING: given ideal is not holonomic");
     1575    print("... setting bound for degree of b-function to 10 and proceeding");
    15141576    isH = 10;
    15151577  }
    … …  
    15401602  setring D;
    15411603  ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6; I = std(I);
    1542   intvec w1 = 1,1;
    1543   intvec w2 = 1,2;
    1544   intvec w3 = 2,3;
     1604  intvec w1 = 0,1;
     1605  intvec w2 = 2,3;
    15451606  bfctIdeal(I,w1);
    1546   bfctIdeal(I,w2,1);
    1547   bfctIdeal(I,w3,0,1);
     1607  bfctIdeal(I,w2,0,1);
     1608  ideal J = I[size(I)]; // J is not holonomic by construction
     1609  bfctIdeal(J,w1); //  b-function of D/J w.r.t. w1 is non-zero
     1610  bfctIdeal(J,w2); //  b-function of D/J w.r.t. w2 is zero
    15481611}
    15491612
    … …  
    17681831  poly xyzreiffen45 = (xy+z)*(y4+yz4+z5);
    17691832  bfct(xyzreiffen45);
     1833
     1834// sparse ideal as suggested by Alex; gives 1 as std
     1835ideal I1 = 28191*y^2+14628*x*Dy, 24865*x^2+24072*x*Dx+17756*Dy^2;
     1836std(I1);
     1837
     1838
    17701839}
    17711840*/
Note: See TracChangeset for help on using the changeset viewer.