Changeset ab3db62 in git


Ignore:
Timestamp:
Jul 7, 2009, 10:46:51 PM (14 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'spielwiese', 'f6c3dc58b0df4bd712574325fe76d0626174ad97')
Children:
4d10648d5056326751bd51ace1aada9cd81a993f
Parents:
642e65dac97b59e7fbbbcd4fa0e29e65a4ae0853
Message:
*levandov: changes by daniel andres, new functions, doc enhancements


git-svn-id: file:///usr/local/Singular/svn/trunk@11958 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular/LIB
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/bfun.lib

    r642e65 rab3db62  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="$Id: bfun.lib,v 1.9 2009-04-15 11:09:27 seelisch Exp $";
     2version="$Id: bfun.lib,v 1.10 2009-07-07 20:46:51 levandov Exp $";
    33category="Noncommutative";
    44info="
     
    2727@*   K-vector space generated by all initial forms w.r.t (-w,w) of elements of I.
    2828@*   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. 
     29@*   the intersection of J with the PID K[s] is called the b-function of I w.r.t.  w.
    3030@*   Unlike Bernstein-Sato polynomial, general b-function with respect to
    3131@*   arbitrary weights need not have rational roots at all. However, b-function
     
    6565LIB "dmodapp.lib";  // for initialIdealW etc
    6666LIB "nctools.lib";  // for isWeyl etc
     67LIB "presolve.lib"; // for valvars
    6768
    6869///////////////////////////////////////////////////////////////////////////////
     
    9091//--------------- auxiliary procedures ----------------------------------------
    9192
     93/*
    9294static proc gradedWeyl (intvec u,intvec v)
    9395"USAGE:  gradedWeyl(u,v); u,v intvecs
     
    154156  setring G; G;
    155157}
     158*/
    156159
    157160static proc safeVarName (string s)
     
    930933"
    931934{
    932   // assume I is given in Groebner basis
     935  // assume I is given as Groebner basis
    933936  if (attrib(I,"isSB") <> 1)
    934937  {
     
    12501253}
    12511254
    1252 static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodpIntersect, int pIntersectchar, int modengine, intvec u0)
     1255static proc addRoot(number q, list L)
     1256{
     1257  // add root to list in bFactor format
     1258  int i,qInL;
     1259  ideal I = L[1];
     1260  intvec v = L[2];
     1261  list LL;
     1262  if (v == 0)
     1263  {
     1264    I = poly(q);
     1265    v = 1;
     1266    LL = I,v;
     1267  }
     1268  else
     1269  {
     1270    for (i=1; i<=ncols(I); i++)
     1271    {
     1272      if (I[i] == q)
     1273      {
     1274        qInL = i;
     1275        break;
     1276      }
     1277    }
     1278    if (qInL)
     1279    {
     1280      v[qInL] = v[qInL] + 1;
     1281    }
     1282    else
     1283    {
     1284      I = q,I;
     1285      v = 1,v;
     1286    }
     1287  }
     1288  LL = I,v;
     1289  if (size(L) == 3) // irreducible factor
     1290  {
     1291    if (L[3] <> "0" && L[3] <> "1")
     1292    {
     1293      LL = LL + list(L[3]);
     1294    }
     1295  }
     1296  return(LL);
     1297}
     1298
     1299static proc bfctengine (poly f, int inorann, int whichengine, int addPD, int stdsum, int methodord, int methodpIntersect, int pIntersectchar, int modengine, intvec u0)
    12531300{
    12541301  int ppl = printlevel - voice +2;
     
    12691316    poly f = imap(save,f);
    12701317  }
    1271   if (size(variables(f)) == 0) // f is constant
    1272   {
    1273     return(list(ideal(0),intvec(0)));
    1274   } 
    12751318  if (inorann == 0) // bfct using initial ideal
    12761319  {
     1320    // list L = ringlist(basering);
     1321    intvec iv = valvars(f)[1]; // heuristacally better ordering for initialMalgrange
     1322    list varL = L[2];
     1323    varL = varL[iv];
     1324    L[2] = varL;
     1325    if (u0 <> 0)
     1326    {
     1327      u0 = u0[iv];
     1328    }
     1329    def newr = ring(L);
     1330    kill varL,iv,L;
     1331    setring newr;
     1332    poly f = imap(save,f);
    12771333    def D = initialMalgrange(f,whichengine,methodord,u0);
    12781334    setring D;
     
    12831339  else // bfct using Ann(f^s)
    12841340  {
    1285     def D = SannfsBFCT(f,whichengine);
     1341    def D = SannfsBFCT(f,addPD,whichengine,stdsum);
    12861342    setring D;
    12871343    ideal J = LD;
    12881344    kill LD;
    1289     poly s = var(2*n+1);
     1345    poly s = var(1);
    12901346  }
    12911347  vector b;
     
    13541410  }
    13551411  if (inorann == 0) { list l = listofroots(b,1); }
    1356   else              { list l = listofroots(b,0); }
     1412  else // bfctAnn
     1413  {
     1414    list l = listofroots(b,0);
     1415    if (addPD)
     1416    {
     1417      l = addRoot(-1,l);
     1418    }
     1419  }
    13571420  setring save;
    13581421  list l = imap(D,l);
     
    14131476    }
    14141477  }
    1415   list b = bfctengine(f,0,whichengine,methodord,0,0,0,u0);
     1478  list b = bfctengine(f,0,whichengine,0,0,methodord,0,0,0,u0);
    14161479  return(b);
    14171480}
     
    15011564    }
    15021565  }
    1503   list b = bfctengine(f,0,whichengine,methodord,1,pIntersectchar,modengine,u0);
     1566  list b = bfctengine(f,0,whichengine,0,0,methodord,1,pIntersectchar,modengine,u0);
    15041567  return(b);
    15051568}
     
    17711834}
    17721835
     1836
     1837
    17731838proc bfctAnn (poly f, list #)
    1774 "USAGE:  bfctAnn(f [,r]);  f a poly, r an optional int
     1839"USAGE:  bfctAnn(f [,a,b,c]);  f a poly, a, b, c optional ints
    17751840RETURN:  list of ideal and intvec
    17761841PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the
    17771842@*       hypersurface defined by f.
    17781843ASSUME:  The basering is commutative and of characteristic 0.
    1779 BACKGROUND: In this proc, ann(f^s) is computed and then a system of linear
     1844BACKGROUND: In this proc, Ann(f^s) is computed and then a system of linear
    17801845@*       equations is solved by linear reductions.
    17811846NOTE:    In the output list, the ideal contains all the roots and the intvec
    1782 @*       their multiplicities.
    1783 @*       If r<>0, @code{std} is used for GB computations,
    1784 @*       otherwise, and by default, @code{slimgb} is used.
     1847@*  their multiplicities.
     1848@*  If a<>0, only f is appended to Ann(f^s), otherwise, and by default,
     1849@*  f and all its partial derivatives are appended.
     1850@*  If b<>0, @code{std} is used for GB computations, otherwise, and by
     1851@*  default, @code{slimgb} is used.
     1852@*  If c<>0, @code{std} is used for Groebner basis computations of ideals
     1853@*  <I+J> when I is already a Groebner basis of <I>.
     1854@*  Otherwise, and by default the engine determined by the switch b is used.
     1855@*  Note that in the case c<>0, the choice for b will be overwritten only
     1856@*  for the types of ideals mentioned above.
     1857@*  This means that if b<>0, specifying c has no effect.
    17851858DISPLAY: If printlevel=1, progress debug messages will be printed,
    17861859@*       if printlevel>=2, all the debug messages will be printed.
     
    17901863  def save = basering;
    17911864  int ppl = printlevel - voice + 2;
     1865  int addPD       = 1; // default
    17921866  int whichengine = 0; // default
     1867  int stdsum      = 0; // default
    17931868  if (size(#)>0)
    17941869  {
    17951870    if (typeof(#[1])=="int" || typeof(#[1])=="number")
    17961871    {
    1797       whichengine = int(#[1]);
    1798     }
    1799   }
    1800   list b = bfctengine(f,1,whichengine,0,0,0,0,0);
     1872      addPD = int(#[1]);
     1873    }
     1874    if (size(#)>1)
     1875    {
     1876      if (typeof(#[2])=="int" || typeof(#[2])=="number")
     1877      {
     1878        whichengine = int(#[2]);
     1879      }
     1880      if (size(#)>2)
     1881      {
     1882        if (typeof(#[3])=="int" || typeof(#[3])=="number")
     1883        {
     1884          stdsum = int(#[3]);
     1885        }
     1886      }
     1887    }
     1888  }
     1889  list b = bfctengine(f,1,whichengine,addPD,stdsum,0,0,0,0,0);
    18011890  return(b);
    18021891}
     
    18071896  poly f = x^2+y^3+x*y^2;
    18081897  bfctAnn(f);
     1898  def R = reiffen(4,5); setring R;
     1899  RC; // the Reiffen curve in 4,5
     1900  bfctAnn(RC,0,1);
    18091901}
    18101902
    18111903/*
    1812 //static proc hardexamples ()
     1904static proc hardexamples ()
    18131905{
    18141906  //  some hard examples
     
    18351927ideal I1 = 28191*y^2+14628*x*Dy, 24865*x^2+24072*x*Dx+17756*Dy^2;
    18361928std(I1);
    1837 
    1838 
    18391929}
    18401930*/
  • Singular/LIB/dmod.lib

    r642e65 rab3db62  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="$Id: dmod.lib,v 1.44 2009-07-07 16:06:02 levandov Exp $";
     2version="$Id: dmod.lib,v 1.45 2009-07-07 20:46:50 levandov Exp $";
    33category="Noncommutative";
    44info="
     
    34993499}
    35003500
     3501static proc safeVarName (string s, list #)
     3502{
     3503  string S;
     3504  int cv = 1;
     3505  if (size(#)>1)
     3506  {
     3507    if (#[1]=="v") { cv = 0; S = varstr(basering); }
     3508    if (#[1]=="c") { cv = 0; S = charstr(basering); }
     3509  }
     3510  if (cv) { S = charstr(basering) + "," + varstr(basering); }
     3511  S = "," + S + ",";
     3512  s = "," + s + ",";
     3513  while (find(S,s) <> 0)
     3514  {
     3515    s[1] = "@";
     3516    s = "," + s;
     3517  }
     3518  s = s[2..size(s)-1];
     3519  return(s)
     3520}
     3521
    35013522proc SannfsBFCT(poly F, list #)
    3502 "USAGE:  SannfsBFCT(f [,eng]);  f a poly, eng an optional int
     3523"USAGE:  SannfsBFCT(f [,a,b,c]);  f a poly, a,b,c optional ints
    35033524RETURN:  ring
    3504 PURPOSE: compute Ann f^s and Groebner basis of Ann f^s+f in D[s]
    3505 NOTE:    activate the output ring with the @code{setring} command.
    3506 @*  This procedure, unlike SannfsBM, returns a ring with the degrevlex
    3507 @*  ordering in all variables.
    3508 @*  In the ring D[s], the ideal LD is the ideal needed (which is returned as a Groebner basis).
    3509 @*  If eng <>0, @code{std} is used for Groebner basis computations,
     3525PURPOSE: compute a Groebner basis either of Ann(f^s)+<f> or of
     3526@*       Ann(f^s)+<f,f_1,...,f_n> in D[s]
     3527NOTE:    Activate the output ring with the @code{setring} command.
     3528@*  This procedure, unlike SannfsBM, returns the ring D[s] with an anti-
     3529@*  elimination ordering for s.
     3530@*  The output ring contains an ideal @code{LD}, being a Groebner basis
     3531@*  either of  Ann(f^s)+<f>, if a=0 (and by default), or of
     3532@*  Ann(f^s)+<f,f_1,...,f_n>, otherwise.
     3533@*  Here, f_i stands for the i-th  partial derivative of f.
     3534@*  If b<>0, @code{std} is used for Groebner basis computations,
    35103535@*  otherwise, and by default @code{slimgb} is used.
     3536@*  If c<>0, @code{std} is used for Groebner basis computations of
     3537@*  ideals <I+J> when I is already a Groebner basis of <I>.
     3538@*  Otherwise, and by default the engine determined by the switch b is
     3539@*  used. Note that in the case c<>0, the choice for b will be
     3540@*  overwritten only for the types of ideals mentioned above.
     3541@*  This means that if b<>0, specifying c has no effect.
    35113542DISPLAY: If printlevel=1, progress debug messages will be printed,
    3512 @*          if printlevel>=2, all the debug messages will be printed.
     3543@*       if printlevel>=2, all the debug messages will be printed.
    35133544EXAMPLE: example SannfsBFCT; shows examples
    35143545"
    35153546{
    3516   // DEBUG INFO: ordering on the output ring = dp, engine on the sum of ideals is used
    3517   // hence it's the situation for slimgb
    3518 
    3519   int eng = 0;
    3520   if ( size(#)>0 )
    3521   {
    3522     if ( typeof(#[1]) == "int" )
    3523     {
    3524       eng = int(#[1]);
    3525     }
    3526   }
    3527   // returns a list with a ring and an ideal LD in it
     3547  int addPD,eng,stdsum;
     3548  if (size(#)>0)
     3549  {
     3550    if (typeof(#[1])=="int" || typeof(#[1])=="number")
     3551    {
     3552      addPD = int(#[1]);
     3553    }
     3554    if (size(#)>1)
     3555    {
     3556      if (typeof(#[2])=="int" || typeof(#[2])=="number")
     3557      {
     3558        eng = int(#[2]);
     3559      }
     3560      if (size(#)>2)
     3561      {
     3562        if (typeof(#[3])=="int" || typeof(#[3])=="number")
     3563        {
     3564          stdsum = int(#[3]);
     3565        }
     3566      }
     3567    }
     3568  }
    35283569  int ppl = printlevel-voice+2;
    3529   //  printf("plevel :%s, voice: %s",printlevel,voice);
    35303570  def save = basering;
    3531   int N = nvars(basering);
    3532   int Nnew = 2*N+2;
     3571  int N = nvars(save);
     3572  intvec optSave = option(get);
    35333573  int i,j;
    3534   list RL = ringlist(basering);
     3574  list RL = ringlist(save);
     3575  // ----- step 1: compute syzigies
     3576  intvec iv;
     3577  list L,Lord;
     3578  iv = 1:N; Lord[1] = list("dp",iv);
     3579  iv = 0;   Lord[2] = list("C",iv);
     3580  L = RL;
     3581  L[3] = Lord;
     3582  def @RM = ring(L);
     3583  kill L,Lord;
     3584  setring @RM;
     3585  option(redSB);
     3586  option(redTail);
     3587  def RM = makeModElimRing(@RM);
     3588  setring RM;
     3589  poly F = imap(save,F);
     3590  ideal J = jacob(F);
     3591  J = F,J;
     3592  dbprint(ppl,"// -1-1- Starting the computation of syz(F,_Dx(F))");
     3593  dbprint(ppl-1, J);
     3594  module M = syz(J);
     3595  dbprint(ppl,"// -1-2- The module syz(F,_Dx(F)) has been computed");
     3596  dbprint(ppl-1, M);
     3597  dbprint(ppl,"// -1-3- Starting GB computation of syz(F,_Dx(F))");
     3598  M = engine(M,eng);
     3599  dbprint(ppl,"// -1-4- GB computation finished");
     3600  dbprint(ppl-1, M);
     3601  // ----- step 2: compute part of Ann(F^s)
     3602  setring save;
     3603  option(set,optSave);
     3604  module M = imap(RM,M);
     3605  kill optSave,RM;
     3606  // ----- create D[s]
     3607  int Nnew = 2*N+1;
    35353608  list L, Lord;
    3536   L[1] = RL[1]; // char
    3537   L[4] = RL[4]; // char, minpoly
    3538   // check whether vars have admissible names
     3609  // ----- keep char, minpoly
     3610  L[1] = RL[1];
     3611  L[4] = RL[4];
     3612  // ----- create names for new vars
    35393613  list Name  = RL[2];
    3540   list RName;
    3541   RName[1] = "@t";
    3542   RName[2] = "@s";
    3543   list PName;
    3544   int BadName;
    3545   if (size(RL[1])>1)
    3546   {
    3547     PName = RL[1][2];
    3548   }
    3549   //  PName;
    3550   for(i=1;i<=2;i++)
    3551   {
    3552     for(j=1; j<=N; j++)
    3553     {
    3554       if (Name[j] == RName[i])
    3555       {
    3556         BadName = 1;
    3557         break;
    3558       }
    3559     }
    3560     for(j=1; j<=size(PName); j++)
    3561     {
    3562       if (PName[j] == RName[i])
    3563       {
    3564         BadName = 1;
    3565         break;
    3566       }
    3567     }
    3568   }
    3569   if (BadName == 1)
    3570   {
    3571     ERROR("Variable/parameter names should not include @t,@s");
    3572   }
    3573   kill PName,BadName;
    3574   // now, create the names for new vars
     3614  string newVar@s = safeVarName("s");
     3615  if (newVar@s[1] == "@")
     3616  {
     3617    print("Name s already assigned to parameter/ringvar.");
     3618    print("Using " + newVar@s + " instead.")
     3619  }
    35753620  list DName;
    3576   for(i=1;i<=N;i++)
    3577   {
    3578     DName[i] = "D"+Name[i]; // concat
    3579   }
    3580   list NName = RName + DName + Name ;
    3581   L[2]   = NName;
    3582   // Name, Dname will be used further
    3583   kill NName;
     3621  for (i=1; i<=N; i++)
     3622  {
     3623    DName[i] = safeVarName("D" + Name[i]);
     3624  }
     3625  L[2] = list(newVar@s) + Name + DName;
     3626  // ----- create ordering
     3627  // --- anti-elimination ordering for s
     3628  iv = 1;       Lord[1] = list("dp",iv);
     3629  iv = 1:(2*N); Lord[2] = list("dp",iv);
     3630  iv = 0;       Lord[3] = list("C",iv);
     3631  L[3] = Lord;
     3632  // ----- create commutative ring
     3633  def @Ds = ring(L);
     3634  kill L,Lord;
     3635  setring @Ds;
     3636  // ----- create nc relations
     3637   matrix Drel[Nnew][Nnew];
     3638  for (i=1; i<=N; i++)
     3639  {
     3640    Drel[i+1,N+1+i] = 1;
     3641  }
     3642  def Ds = nc_algebra(1,Drel);
     3643  setring Ds;
     3644  kill @Ds;
     3645  dbprint(ppl,"// -2-1- The ring D[s] is ready");
     3646  dbprint(ppl-1, Ds);
     3647  matrix M = imap(save,M);
     3648  vector v = var(1)*gen(1);
     3649  for (i=1; i<=N; i++)
     3650  {
     3651    v = v + var(i+1+N)*gen(i+1); //[s,_Dx]
     3652  }
     3653  ideal J = transpose(M)*v;
     3654  kill M,v;
     3655  dbprint(ppl,"// -2-2- Compute part of Ann(F^s)");
     3656  dbprint(ppl-1, J);
     3657  J = engine(J,eng);
     3658  dbprint(ppl,"// -2-3- GB computation finished");
     3659  dbprint(ppl-1, J);
     3660  // ----- step 3: the full annihilator
     3661  // ----- create D<t,s>
     3662  setring save;
     3663  Nnew = 2*N+2;
     3664  list L, Lord;
     3665  // ----- keep char, minpoly
     3666  L[1] = RL[1];
     3667  L[4] = RL[4]; 
     3668  // ----- create vars
     3669  string newVar@t = safeVarName("t");
     3670  L[2] = list(newVar@t,newVar@s) + DName + Name;
     3671  // ----- create ordering for elimination of t
    35843672  // block ord (lp(2),dp);
    3585   Lord[1] = list("lp",intvec(1,1));
    3586   // continue with dp 1,1,1,1...
    3587   Lord[2] = list("dp",intvec(1:Nnew));
    3588   Lord[3] = list("C",intvec(0));
    3589   L[3]      = Lord;
    3590   // we are done with the list
    3591   def @R@ = ring(L);
    3592   setring @R@;
    3593   matrix @D[Nnew][Nnew];
    3594   @D[1,2] = var(1);
     3673  iv = 1,1;    Lord[1] = list("lp",iv);
     3674  iv = 1:Nnew; Lord[2] = list("dp",iv);
     3675  iv = 0;      Lord[3] = list("C",iv);
     3676  L[3] = Lord;
     3677  def @Dts = ring(L);
     3678  kill RL,L,Lord,Name,DName,newVar@s,newVar@t;
     3679  setring @Dts;
     3680  // ----- create nc relations
     3681  matrix Drel[Nnew][Nnew];
     3682  Drel[1,2] = var(1);
    35953683  for(i=1; i<=N; i++)
    35963684  {
    3597     @D[2+i,N+2+i]=-1;
    3598   }
    3599   //  L[5] = matrix(UpOneMatrix(Nnew));
    3600   //  L[6] = @D;
    3601   def @R = nc_algebra(1,@D);
    3602   setring @R;
    3603   kill @R@;
    3604   dbprint(ppl,"// -1-1- the ring @R(@t,@s,_Dx,_x) is ready");
    3605   dbprint(ppl-1, @R);
    3606   // create the ideal I
     3685    Drel[2+i,N+2+i]=-1;
     3686  }
     3687  def Dts = nc_algebra(1,Drel);
     3688  setring Dts;
     3689  kill @Dts;
     3690  dbprint(ppl,"// -3-1- The ring D<t,s> is ready");
     3691  dbprint(ppl-1, Dts);
     3692  // ----- create the ideal I following BM
    36073693  poly  F = imap(save,F);
    3608   ideal I = var(1)*F+var(2);
     3694  ideal I = var(1)*F + var(2); // = t*F + s
    36093695  poly p;
    36103696  for(i=1; i<=N; i++)
    36113697  {
    3612     p = var(1); // t
    3613     p = diff(F,var(N+2+i))*p;
    3614     I = I, var(2+i) + p;
    3615   }
    3616   // -------- the ideal I is ready ----------
    3617   dbprint(ppl,"// -1-2- starting the elimination of @t in @R");
     3698    p = var(1)*diff(F,var(N+2+i)) + var(2+i); // = t*F_i + D_i
     3699    I[i+1] = p;
     3700  }
     3701  // ----- add already computed part to it
     3702  ideal MM = var(2);      // s
     3703  for (i=1; i<=N; i++)
     3704  {
     3705    MM[1+i] = var(2+N+i); // _x
     3706    MM[1+N+i] = var(2+i); // _Dx
     3707  }
     3708  map Ds2Dts = Ds,MM;
     3709  ideal J = Ds2Dts(J);
     3710  attrib(J,"isSB",1);
     3711  kill MM,Ds2Dts;
     3712  // ----- start the elimination
     3713  dbprint(ppl,"// -3-2- Starting the elimination of t in D<t,s>");
    36183714  dbprint(ppl-1, I);
    3619   ideal J = engine(I,eng);
    3620   ideal K = nselect(J,1);
    3621   dbprint(ppl,"// -1-3- @t is eliminated");
    3622   dbprint(ppl-1, K);  // K is without @t
    3623   K = engine(K,eng);  // std does the job too
    3624   // now, we must change the ordering
    3625   // and create a ring without @t
    3626   //  setring S;
    3627   // ----------- the ring @R3 ------------
    3628   // _Dx,_x,s;  +fast ord !
    3629   // keep: N, i,j,RL
    3630   Nnew = 2*N+1;
    3631   kill Lord, RName;
    3632   list Lord;
    3633   list L=imap(save,L);
    3634   list RL=imap(save,RL);
    3635   L[1] = RL[1];
    3636   L[4] = RL[4];  // char, minpoly
    3637   // check whether vars hava admissible names -> done earlier
    3638   // DName is defined earlier
    3639   list NName = DName + Name + list(string(var(2)));
    3640   L[2] = NName;
    3641   kill NName;
    3642   // just dp
    3643   // block ord (dp(N),dp);
    3644   Lord[1] = list("dp",intvec(1:Nnew));
    3645   Lord[2]= list("C",intvec(0));
    3646   L[3] = Lord;
    3647   // we are done with the list. Now add a Plural part
    3648   def @R2@ = ring(L);
    3649   setring @R2@;
    3650   matrix @D[Nnew][Nnew];
     3715  if (stdsum || eng <> 0)
     3716  {
     3717    I = std(J,I);
     3718  }
     3719  else
     3720  {
     3721    I = J,I;
     3722    I = engine(I,eng); 
     3723  }
     3724  kill J;
     3725  I = nselect(I,1);
     3726  dbprint(ppl,"// -3-3- t is eliminated");
     3727  dbprint(ppl-1, I);  // I is without t
     3728  // ----- step 4: add F
     3729  // ----- back to D[s]
     3730  setring Ds;
     3731  ideal MM = 0,var(1);      // 0,s
    36513732  for (i=1; i<=N; i++)
    36523733  {
    3653     @D[i,N+i]=-1;
    3654   }
    3655   def @R2 = nc_algebra(1,@D);
    3656   setring @R2;
    3657   kill @R2@;
    3658   dbprint(ppl,"//  -2-1- the ring @R2(_Dx,_x,@s) is ready");
    3659   dbprint(ppl-1, @R2);
    3660   ideal MM = maxideal(1);
    3661   MM = 0,var(Nnew),MM;
    3662   map R01 = @R, MM;
    3663   ideal K = R01(K);
    3664   // total cleanup
    3665   poly F = imap(save, F);
    3666   ideal LD = K,F;
    3667   dbprint(ppl,"//  -2-2- start GB computations for Ann f^s + f");
     3734    MM[2+i]   = var(1+N+i); // _Dx
     3735    MM[2+N+i] = var(1+i);   // _x
     3736  }
     3737  map Dts2Ds = Dts, MM;
     3738  ideal LD = Dts2Ds(I);
     3739  kill J,Dts,Dts2Ds,MM;
     3740  dbprint(ppl,"// -4-1- Starting cosmetic Groebner computation");
     3741  LD = engine(LD,eng);
     3742  dbprint(ppl,"// -4-2- Finished cosmetic Groebner computation");
    36683743  dbprint(ppl-1, LD);
    3669   LD = engine(LD,eng);
    3670   dbprint(ppl,"//  -2-3- finished GB computations for Ann f^s + f");
     3744  // ----- use reduction trick as Macaulay2 does: compute b(s)/(s+1) by adding all partial derivations also
     3745  ideal J;
     3746  if (addPD)
     3747  {
     3748    setring @RM;
     3749    poly F = imap(save,F);
     3750    ideal J = jacob(F);
     3751    J = F,J;
     3752    dbprint(ppl,"// -4-2-1- Start GB computation <f, f_i>");
     3753    J = engine(J,eng);
     3754    dbprint(ppl,"// -4-2-2- Finished GB computation <f, f_i>");
     3755    dbprint(ppl-1, J);
     3756    setring Ds;
     3757    J = imap(@RM,J);
     3758    attrib(J,"isSB",1);
     3759    dbprint(ppl,"// -4-3- Start GB computations for Ann f^s + <f, f_i>");
     3760  }
     3761  else
     3762  {
     3763    J = imap(save,F);
     3764    dbprint(ppl,"// -4-3- Start GB computations for Ann f^s + <f>");
     3765  }
     3766  kill @RM;
     3767  // ----- the really hard part
     3768  if (stdsum || eng <> 0)
     3769  {
     3770    LD = std(LD,J);
     3771  }
     3772  else
     3773  {
     3774    LD = LD,J;
     3775    LD = engine(LD,eng); 
     3776  }
     3777  if (addPD) { dbprint(ppl,"// -4-4- Finished GB computations for Ann f^s + <f, f_i>"); }
     3778  else       { dbprint(ppl,"// -4-4- Finished GB computations for Ann f^s + <f>"); }
    36713779  dbprint(ppl-1, LD);
    3672   // make leadcoeffs positive
    3673   for (i=1; i<= ncols(LD); i++)
    3674   {
    3675     if (leadcoef(LD[i]) <0 )
    3676     {
    3677       LD[i] = -LD[i];
    3678     }
    3679   }
    36803780  export LD;
    3681   kill @R;
    3682   return(@R2);
     3781  return(Ds);
    36833782}
    36843783example
     
    36873786  ring r = 0,(x,y,z,w),Dp;
    36883787  poly F = x^3+y^3+z^3*w;
    3689   printlevel = 0;
    3690   def A = SannfsBFCT(F); setring A;
    3691   intvec v = 1,2,3,4,5,6,7,8;
    3692   // are there polynomials, depending on @s only?
    3693   nselect(LD,v);
     3788  // compute Ann(F^s)+<F> using slimgb only
     3789  def A = SannfsBFCT(F);
     3790  setring A; A;
     3791  LD;
     3792  // the Bernstein-Sato poly of F:
     3793  vec2poly(pIntersect(s,LD));
    36943794  // a fancier example:
    36953795  def R = reiffen(4,5); setring R;
    3696   v = 1,2,3,4;
    36973796  RC; // the Reiffen curve in 4,5
    3698   def B = SannfsBFCT(RC);
     3797  // compute Ann(RC^s)+<RC,diff(RC,x),diff(RC,y)>
     3798  // using std for GB computations of ideals <I+J>
     3799  // where I is already a GB of <I>
     3800  // and slimgb for other ideals
     3801  def B = SannfsBFCT(RC,1,0,1);
    36993802  setring B;
    3700   // Are there polynomials, depending on @s only?
    3701   nselect(LD,v);
    3702   // It is not the case. Are there leading terms in @s only?
    3703   nselect(lead(LD),v);
     3803  // the Bernstein-Sato poly of RC:
     3804  (s-1)*vec2poly(pIntersect(s,LD));
    37043805}
    37053806
  • Singular/LIB/dmodapp.lib

    r642e65 rab3db62  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="$Id: dmodapp.lib,v 1.30 2009-07-07 16:17:06 levandov Exp $";
     2version="$Id: dmodapp.lib,v 1.31 2009-07-07 20:46:51 levandov Exp $";
    33category="Noncommutative";
    44info="
     
    6363LIB "sing.lib";
    6464LIB "primdec.lib";
    65 LIB "dmod.lib"; // loads e.g. nctools.lib
    66 LIB "bfun.lib"; //formerly LIB "bfct.lib";
     65LIB "dmod.lib";    // loads e.g. nctools.lib
     66LIB "bfun.lib";    //formerly LIB "bfct.lib";
    6767LIB "nctools.lib"; // for isWeyl etc
    6868LIB "gkdim.lib";
     
    14631463
    14641464proc initialIdealW (ideal I, intvec u, intvec v, list #)
    1465 "USAGE:  initialIdealW(I,u,v [,s,t]);  I ideal, u,v intvecs, s,t optional ints
     1465"USAGE:  initialIdealW(I,u,v [,s,t,w]);  I ideal, u,v intvecs, s,t optional ints,
     1466@*       w an optional intvec
    14661467RETURN:  ideal, GB of initial ideal of the input ideal w.r.t. the weights u and v
    14671468ASSUME:  The basering is the n-th Weyl algebra in characteristic 0 and  for all
     
    14701471@*       where D(i) is the differential operator belonging to x(i).
    14711472PURPOSE: computes the initial ideal with respect to given weights.
    1472 NOTE:    u and v are understood as weight vectors for x(i) and D(i)
     1473NOTE:    u and v are understood as weight vectors for x(1..n) and D(1..n)
    14731474@*       respectively.
    14741475@*       If s<>0, @code{std} is used for Groebner basis computations,
     
    14761477@*       If t<>0, a matrix ordering is used for Groebner basis computations,
    14771478@*       otherwise, and by default, a block ordering is used.
     1479@*       If w consist of 2n strictly positive entries, w is used for weighted
     1480@*       homogenization, otherwise, and by default, no weights are used.
    14781481DISPLAY: If printlevel=1, progress debug messages will be printed,
    14791482@*       if printlevel>=2, all the debug messages will be printed.
     
    14811484"
    14821485{
    1483 
    14841486  if (dmodappassumeViolation())
    14851487  {
    14861488    ERROR("Basering is inappropriate: characteristic>0 or qring present");
    14871489  }
    1488 
    14891490  if (!isWeyl())
    14901491  {
    1491     ERROR("Basering is not a Weyl algebra");
    1492   }
    1493 
     1492    ERROR("Basering is not a Weyl algebra.");
     1493  }
    14941494  int ppl = printlevel - voice +2;
    14951495  def save = basering;
     
    14981498  list RL = ringlist(save);
    14991499  int i;
    1500   int whichengine = 0; // default
    1501   int methodord   = 0; // default
     1500  int whichengine = 0;           // default
     1501  int methodord   = 0;           // default
     1502  intvec homogweights = 1:(2*n); // default
    15021503  if (size(#)>0)
    15031504  {
     
    15121513        methodord = int(#[2]);
    15131514      }
    1514     }
    1515   }
    1516   if (char(save) <> 0) { ERROR("characteristic of basering has to be 0"); }
    1517   if (isWeyl() == 0)   { ERROR("basering is not a Weyl algebra"); }
     1515      if (size(#)>2)
     1516      {
     1517        if (typeof(#[3])=="intvec")
     1518        {
     1519          if (size(#[3])==2*n && allPositive(#[3])==1)
     1520          {
     1521            homogweights = #[3];
     1522          }
     1523          else
     1524          {
     1525            print("// Homogenization weight vector must consist of positive entries and be");
     1526            print("// of size " + string(n) + ". Using no weights.");
     1527          }
     1528        }
     1529      }
     1530    }
     1531  }
    15181532  for (i=1; i<=n; i++)
    15191533  {
    15201534    if (bracket(var(i+n),var(i))<>1)
    15211535    {
    1522       ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i)));
    1523     }
    1524   }
    1525  
     1536      ERROR(string(var(i+n)) + " is not a differential operator for " + string(var(i)));
     1537    }
     1538  } 
    15261539  // 1. create  homogenized Weyl algebra
    15271540  // 1.1 create ordering
    15281541  intvec uv = u,v,0;
    1529   list Lord = list(list("a",intvec(1:N)));
     1542  homogweights = homogweights,1;
     1543  list Lord = list(list("a",homogweights));
    15301544  list C0 = list("C",intvec(0));
    15311545  if (methodord == 0) // default: blockordering
    15321546  {
    1533     Lord[2] = list("dp",intvec(1:(N-1)));
    1534     Lord[3] = list("lp",intvec(1));
    1535     Lord[4] = C0;
     1547    Lord[2] = list("a",uv);
     1548    Lord[3] = list("dp",intvec(1:(N-1)));
     1549    Lord[4] = list("lp",intvec(1));
     1550    Lord[5] = C0;
    15361551  }
    15371552  else // M() ordering
     
    15401555    @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1);
    15411556    for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; }
    1542     dbprint(ppl,"// the ordering matrix:",@Ord);
     1557    dbprint(ppl-1,"// the ordering matrix:",@Ord);
    15431558    Lord[2] = list("M",intvec(@Ord));
    15441559    Lord[3] = C0;
     
    15511566  def @Dh = ring(L@@Dh); kill L@@Dh;
    15521567  setring @Dh;
    1553   dbprint(ppl,"// the ring @Dh:",@Dh);
     1568  dbprint(ppl-1,"// the ring @Dh:",@Dh);
    15541569  // 1.4 create non-commutative relations
    15551570  matrix @relD[N][N];
    1556   for (i=1; i<=n; i++) { @relD[i,n+i] = var(N)^2; }
    1557   dbprint(ppl,"// nc relations:",@relD);
     1571  for (i=1; i<=n; i++) { @relD[i,n+i] = var(N)^(homogweights[i]+homogweights[n+i]); }
    15581572  def Dh = nc_algebra(1,@relD);
    1559   setring Dh; kill @Dh;
    1560   dbprint(ppl,"// computing in ring",Dh);
     1573  setring Dh; kill @Dh;   
     1574  dbprint(ppl-1,"// computing in ring",Dh);
    15611575  // 2. Compute the initial ideal
    15621576  ideal I = imap(save,I);
     
    15661580  I = engine(I, whichengine);
    15671581  dbprint(ppl, "// finished Groebner basis computation");
    1568   dbprint(ppl, "// I before dehomogenization is " +string(I));
     1582  dbprint(ppl-1, "// I before dehomogenization is " +string(I));
    15691583  I = subst(I,var(N),1); // dehomogenization
    1570   dbprint(ppl, "I after dehomogenization is " +string(I));
     1584  dbprint(ppl-1, "I after dehomogenization is " +string(I));
    15711585  // 2.2 the initial form
    15721586  I = inForm(I,uv);
     
    15911605
    15921606proc initialMalgrange (poly f,list #)
    1593 "USAGE:  initialMalgrange(f,[,s,t,v]); f poly, s,t optional ints, v opt. intvec
    1594 RETURN:  ring, the Weyl algebra induced by basering, extended by two new vars
     1607"USAGE:  initialMalgrange(f,[,a,b,v]); f poly, a,b optional ints, v opt. intvec
     1608RETURN:  ring, Weyl algebra induced by basering, extended by two new vars t,Dt
    15951609PURPOSE: computes the initial Malgrange ideal of a given polynomial w.r.t. the weight
    15961610@*       vector (-1,0...,0,1,0,...,0) such that the weight of t is -1 and the
     
    16011615@*       of the Malgrange ideal of f.
    16021616@*       Varnames of the basering should not include t and Dt.
    1603 @*       If s<>0, @code{std} is used for Groebner basis computations,
     1617@*       If a<>0, @code{std} is used for Groebner basis computations,
    16041618@*       otherwise, and by default, @code{slimgb} is used.
    1605 @*       If t<>0, a matrix ordering is used for Groebner basis computations,
     1619@*       If b<>0, a matrix ordering is used for Groebner basis computations,
    16061620@*       otherwise, and by default, a block ordering is used.
    1607 @*       If v is a positive weight vector, v is used for homogenization
    1608 @*       computations, otherwise and by default, no weights are used.
     1621@*       If a positive weight vector v is given, the weight
     1622@*       (d,v[1],...,v[n],1,d+1-v[1],...,d+1-v[n]) is used for homogenization
     1623@*       computations, where d denotes the weighted degree of f.
     1624@*       Otherwise and by default, v is set to (1,...,1). See Noro, 2002.
    16091625DISPLAY: If printlevel=1, progress debug messages will be printed,
    16101626@*       if printlevel>=2, all the debug messages will be printed.
     
    16121628"
    16131629{
    1614 
    16151630  if (dmodappassumeViolation())
    16161631  {
    16171632    ERROR("Basering is inappropriate: characteristic>0 or qring present");
    16181633  }
    1619 
    16201634  if (!isCommutative())
    16211635  {
    1622     ERROR("Basering must be commutative");
    1623   }
    1624 
     1636    ERROR("Basering must be commutative.");
     1637  }
    16251638  int ppl = printlevel - voice +2;
    16261639  def save = basering;
     
    16291642  int whichengine = 0; // default
    16301643  int methodord   = 0; // default
    1631   intvec u0 = 0;
     1644  intvec u0 = 1:n;     // default
    16321645  if (size(#)>0)
    16331646  {
     
    16511664    }
    16521665  }
    1653   if (u0 == 0)
    1654   {
    1655     u0 = 1:n;
    1656   }
    1657   if (char(save) <> 0)      { ERROR("characteristic of basering has to be 0"); }
    16581666  // creating the homogenized extended Weyl algebra
    16591667  list RL = ringlist(save);
    16601668  RL = RL[1..4]; // if basering is commutative nc_algebra
    16611669  list C0 = list("C",intvec(0));
    1662   // 1. get the weighted degree of f
     1670  // 1. create homogenization weights
     1671  // 1.1. get the weighted degree of f
    16631672  list Lord = list(list("wp",u0),C0);
    16641673  list L@r = RL;
    16651674  L@r[3] = Lord;
    1666   def r = ring(L@r); kill L@r;
     1675  def r = ring(L@r); kill L@r,Lord;
    16671676  setring r;
    16681677  poly f = imap(save,f);
    16691678  int d = deg(f);
    16701679  setring save; kill r;
    1671   // 2. create homogenized extended Weyl algebra
    1672   int N = 2*n+3;
     1680  // 1.2 the homogenization weights
     1681  intvec homogweights = d;
     1682  homogweights[n+2] = 1;
     1683  for (i=1; i<=n; i++)
     1684  {
     1685    homogweights[i+1]   = u0[i];
     1686    homogweights[n+2+i] = d+1-u0[i];
     1687  }
     1688  // 2. create extended Weyl algebra
     1689  int N = 2*n+2;
    16731690  // 2.1 create names for vars
    16741691  string vart  = safeVarName("t","cv");
     
    16791696    varDt = safeVarName("D"+vart,"cv");
    16801697  }
    1681   list Lvar,Lvarh;
     1698  list Lvar;
    16821699  Lvar[1] = vart; Lvar[n+2] = varDt;
    16831700  for (i=1; i<=n; i++)
     
    16861703    Lvar[i+n+2] = safeVarName("D" + string(var(i)),"cv");
    16871704  }
    1688   Lvarh = Lvar; Lvarh[N] = safeVarName("h","cv");
    16891705  //  2.2 create ordering
    1690   intvec uv,@a,weighttx,weightD;
    1691   uv[1] = -1; uv[n+2] = 1; uv[N] = 0;
    1692   weighttx = d; weightD = 1;
    1693   for (i=1; i<=n; i++)
    1694   {
    1695     weighttx[i+1] = u0[i];
    1696     weightD[i+1]  = d+1-u0[i];
    1697   }
    1698   @a = weighttx,weightD,1;
    1699   Lord[1] = list("a",@a);
    1700   if (methodord == 0) // default: block ordering
    1701   {
    1702     Lord[2] = list("a",uv);
    1703     Lord[3] = list("dp",intvec(1:(N-1)));
    1704     Lord[4] = list("lp",intvec(1));
    1705     Lord[5] = C0;
    1706   }
    1707   else // M() ordering
    1708   {
    1709     intmat @Ord[N][N];
    1710     @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1);
    1711     for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; }
    1712     dbprint(ppl,"// weights for ordering: "+string(transpose(@a)));
    1713     dbprint(ppl,"// the ordering matrix:",@Ord);
    1714     Lord[2] = list("M",intvec(@Ord));
    1715     Lord[3] = C0;
    1716   }
    1717   // 2.3 create commutative ring
    1718   list L@@Dh = RL;
    1719   L@@Dh[2] = Lvarh; L@@Dh[3] = Lord;
    1720   def @Dh = ring(L@@Dh); kill L@@Dh;
    1721   setring @Dh;
    1722   dbprint(ppl,"// the ring @Dh:",@Dh);
    1723   // 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
    1724   // 2.4 create non-commutative relations
    1725   matrix @relD[N][N];
    1726   for (i=1; i<=n+1; i++) { @relD[i,n+1+i] = var(N)^(d+1); }
    1727   dbprint(ppl,"// nc relations:",@relD);
    1728   def Dh = nc_algebra(1,@relD);
    1729   setring Dh; kill @Dh;
    1730   dbprint(ppl,"// computing in ring",Dh);
     1706  list Lord = list(list("dp",intvec(1:N)),C0);
     1707  // 2.3 create the (n+1)-th Weyl algebra
     1708  list L@D = RL; L@D[2] = Lvar; L@D[3] = Lord;
     1709  def @D = ring(L@D); kill L@D;
     1710  setring @D;
     1711  def D = Weyl();
     1712  setring D; kill @D;
     1713  dbprint(ppl,"// the (n+1)-th Weyl algebra :" ,D);
    17311714  // 3. compute the initial ideal
     1715  // 3.1 create the Malgrange ideal
    17321716  poly f = imap(save,f);
    1733   f = homog(f,h);
    1734   // 3.1 create the Malgrange ideal
    17351717  ideal I = var(1)-f;
    17361718  for (i=1; i<=n; i++)
     
    17381720    I = I, var(n+2+i)+diff(f,var(i+1))*var(n+2);
    17391721  }
    1740   // 3.2 the hard part: Groebner basis computation
    1741   dbprint(ppl, "// starting Groebner basis computation with engine: "+string(whichengine));
    1742   I = engine(I, whichengine);
    1743   dbprint(ppl, "// finished Groebner basis computation");
    1744   dbprint(ppl, "// I before dehomogenization is " +string(I));
    1745   I = subst(I,var(N),1); // dehomogenization
    1746   dbprint(ppl, "// I after dehomogenization is " +string(I));
    1747   // 3.3 the initial form
    1748   I = inForm(I,uv);
    1749   // 3.4 create Weyl algebra
    1750   setring save;
    1751   Lord = list();
    1752   Lord[1] = list("dp",intvec(1:(N-1)));
    1753   Lord[2] = C0;
    1754   list L@@D = RL;
    1755   L@@D[2] = Lvar; L@@D[3] = Lord;
    1756   def @D = ring(L@@D); kill L@@D;
    1757   setring @D; def D = Weyl(); setring D;
    1758   ideal I = imap(Dh,I);
    1759   kill @D,Dh;
    1760   // 3.5 the final GB
    1761   dbprint(ppl, "// starting cosmetic Groebner basis computation with engine: "+string(whichengine));
    1762   I = engine(I, whichengine);
    1763   dbprint(ppl,"// finished cosmetic Groebner basis computation");
     1722  // I = engine(I,whichengine); // todo is it efficient to compute a GB first wrt dp first?
     1723  // 3.2 computie the initial ideal
     1724  intvec w = 1,0:n;
     1725  I = initialIdealW(I,-w,w,whichengine,methodord,homogweights);
    17641726  ideal inF = I; attrib(inF,"isSB",1);
    17651727  export(inF);
     1728  setring save;
    17661729  return(D);
    17671730}
Note: See TracChangeset for help on using the changeset viewer.