Changeset f52e64 in git for Singular


Ignore:
Timestamp:
Oct 11, 2010, 9:53:40 PM (14 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'fc741b6502fd8a97288eaa3eba6e5220f3c3df87')
Children:
92f670a7754669a5b094209294c8d2af0499650b
Parents:
d5abcfb5a8d5981beafb238eb7ab2390e274e781
Message:
*levandov: combined set of fixes around release

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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/dmodapp.lib

    rd5abcf rf52e64  
    1010
    1111OVERVIEW: Let K be a field of characteristic 0, R = K[x1,...,xN] and
    12 @* D be the Weyl algebra in variables x1,...,xN,d1,...,dN.
    13 @* In this library there are the following procedures for algebraic D-modules:
    14 @*
     12 D be the Weyl algebra in variables x1,...,xN,d1,...,dN.
     13 In this library there are the following procedures for algebraic D-modules:
     14
    1515@* - given a cyclic representation D/I of a holonomic module and a polynomial
    16 @* F in R, it is proved that the localization of D/I with respect to the mult.
    17 @* closed set of all powers of F is a holonomic D-module. Thus we aim to compute
    18 @* its cyclic representaion D/L for an ideal L in D. The procedures for the
    19 @* localization are DLoc, SDLoc and DLoc0.
    20 @*
     16 F in R, it is proved that the localization of D/I with respect to the mult.
     17 closed set of all powers of F is a holonomic D-module. Thus we aim to compute
     18 its cyclic representaion D/L for an ideal L in D. The procedures for the
     19 localization are DLoc, SDLoc and DLoc0.
     20
    2121@* - annihilator in D of a given polynomial F from R as well as
    22 @* of a given rational function G/F from Quot(R). These can be computed via
    23 @* procedures annPoly resp. annRat.
    24 @*
     22 of a given rational function G/F from Quot(R). These can be computed via
     23 procedures annPoly resp. annRat.
     24
    2525@* - Groebner bases with respect to weights (according to (SST), given an
    26 @* arbitrary integer vector containing weights for variables, one computes the
    27 @* homogenization of a given ideal relative to this vector, then one computes a
    28 @* Groebner basis and returns the dehomogenization of the result), initial
    29 @* forms and initial ideals in Weyl algebras with respect to a given weight
    30 @* vector can be computed with GBWeight, inForm, initialMalgrange and
    31 @* initialIdealW.
    32 @*
     26 arbitrary integer vector containing weights for variables, one computes the
     27 homogenization of a given ideal relative to this vector, then one computes a
     28 Groebner basis and returns the dehomogenization of the result), initial
     29 forms and initial ideals in Weyl algebras with respect to a given weight
     30 vector can be computed with GBWeight, inForm, initialMalgrange and
     31 initialIdealW.
     32
    3333@* - restriction and integration of a holonomic module D/I. Suppose I
    34 @* annihilates a function F(x1,...,xn). Our aim is to compute an ideal J
    35 @* directly from I, which annihilates
     34 annihilates a function F(x1,...,xn). Our aim is to compute an ideal J
     35 directly from I, which annihilates
    3636@*   - F(0,...,0,xk,...,xn) in case of restriction or
    3737@*   - the integral of F with respect to x1,...,xm in case of integration.
    38 @* The corresponding procedures are restrictionModule, restrictionIdeal,
    39 @* integralModule and integralIdeal.
    40 @*
     38 The corresponding procedures are restrictionModule, restrictionIdeal,
     39 integralModule and integralIdeal.
     40
    4141@* - characteristic varieties defined by ideals in Weyl algebras can be computed
    42 @* with charVariety and charInfo.
    43 @*
     42 with charVariety and charInfo.
     43
    4444@* - appelF1, appelF2 and appelF4 return ideals in parametric Weyl algebras,
    45 @* which annihilate corresponding Appel hypergeometric functions.
     45 which annihilate corresponding Appel hypergeometric functions.
    4646
    4747
    4848References:
    4949@* (SST) Saito, Sturmfels, Takayama 'Groebner Deformations of Hypergeometric
    50 @*         Differential Equations', Springer, 2000
     50         Differential Equations', Springer, 2000
    5151@* (OTW) Oaku, Takayama, Walther 'A Localization Algorithm for D-modules',
    52 @*         Journal of Symbolic Computation, 2000
     52         Journal of Symbolic Computation, 2000
    5353@* (OT)  Oaku, Takayama 'Algorithms for D-modules',
    54 @*         Journal of Pure and Applied Algebra, 1998
     54         Journal of Pure and Applied Algebra, 1998
    5555
    5656
    … …  
    100100sortIntvec(v); sorts intvec
    101101
    102 
    103 SEE ALSO: bfun_lib, dmod_lib, dmodvar_lib, gmssing_lib
    104 
    105 
    106102KEYWORDS: D-module; annihilator of polynomial; annihilator of rational function;
    107103D-localization; localization of D-module; D-restriction; restriction of
    108104D-module; D-integration; integration of D-module; characteristic variety;
    109105Appel function; Appel hypergeometric function
     106
     107SEE ALSO: bfun_lib, dmod_lib, dmodvar_lib, gmssing_lib
    110108";
    111109
    112 
    113110/*
    114 CHANGELOG
    115 21.09.10 by DA:
    116 - restructured library for better readability
    117 - new / improved procs:
     111  Changelog
     112  21.09.10 by DA:
     113  - restructured library for better readability
     114  - new / improved procs:
    118115  - toolbox: isInt, intRoots, sortIntvec
    119116  - GB wrt weights: GBWeight, initialIdealW rewritten using GBWeight
    120117  - restriction/integration: restrictionX, integralX where X in {Module, Ideal},
    121     fourier, inverseFourier, deRhamCohom, deRhamCohomIdeal
     118  fourier, inverseFourier, deRhamCohom, deRhamCohomIdeal
    122119  - characteristic variety: charVariety, charInfo
    123 - added keywords for features
    124 - reformated help strings and examples in existing procs
    125 - added SUPPORT in header
    126 - added reference (OT)
    127 
    128 04.10.10 by DA:
    129 - incorporated suggestions by Oleksandr Motsak, among other:
     120  - added keywords for features
     121  - reformated help strings and examples in existing procs
     122  - added SUPPORT in header
     123  - added reference (OT)
     124
     125  04.10.10 by DA:
     126  - incorporated suggestions by Oleksandr Motsak, among other:
    130127  - bugfixes for fl2poly, sortIntvec, annPoly, GBWeight
    131128  - enhanced functionality for deleteGenerator, inForm
     129
     130  11.10.10 by DA:
     131  - procedure bFactor now sorts the roots using new static procedure sortNumberIdeal
    132132*/
    133133
    … …  
    233233  s = s[2..size(s)-1];
    234234  return(s)
    235 }
     235    }
    236236
    237237static proc intLike (def i)
    … …  
    252252
    253253proc engine(def I, int i)
    254 "USAGE:  engine(I,i);  I  ideal/module/matrix, i an int
     254  "USAGE:  engine(I,i);  I  ideal/module/matrix, i an int
    255255RETURN:  the same type as I
    256256PURPOSE: compute the Groebner basis of I with the algorithm, chosen via i
    … …  
    287287
    288288proc poly2list (poly f)
    289 "USAGE:  poly2list(f); f a poly
     289  "USAGE:  poly2list(f); f a poly
    290290RETURN:  list of exponents and corresponding terms of f
    291291PURPOSE: converts a poly to a list of pairs consisting of intvecs (1st entry)
    … …  
    326326
    327327proc fl2poly(list L, string s)
    328 "USAGE:  fl2poly(L,s); L a list, s a string
     328  "USAGE:  fl2poly(L,s); L a list, s a string
    329329RETURN:  poly
    330330PURPOSE: reconstruct a monic polynomial in one variable from its factorization
    … …  
    375375
    376376proc insertGenerator (list #)
    377 "USAGE:  insertGenerator(id,p[,k]);
     377  "USAGE:  insertGenerator(id,p[,k]);
    378378@*       id an ideal/module, p a poly/vector, k an optional int
    379379RETURN:  of the same type as id
    … …  
    456456
    457457proc deleteGenerator (def id, int k)
    458 "USAGE:   deleteGenerator(id,k);  id an ideal/module, k an int
     458  "USAGE:   deleteGenerator(id,k);  id an ideal/module, k an int
    459459RETURN:   of the same type as id
    460460PURPOSE:  deletes the k-th generator from the first argument and returns
    … …  
    501501}
    502502
     503static proc sortNumberIdeal (ideal I)
     504// sorts ideal of constant polys (ie numbers), returns according permutation
     505{
     506  int i;
     507  int nI = ncols(I);
     508  intvec dI;
     509  for (i=nI; i>0; i--)
     510  {
     511    dI[i] = int(denominator(leadcoef(I[i])));
     512  }
     513  int lcmI = lcm(dI);
     514  for (i=nI; i>0; i--)
     515  {
     516    dI[i] = int(lcmI*leadcoef(I[i]));
     517  }
     518  intvec perm = sortIntvec(dI)[2];
     519  return(perm);
     520}
     521example
     522{
     523  "EXAMPLE:"; echo = 2;
     524  ring r = 0,s,dp;
     525  ideal I = -9/20,-11/20,-23/20,-19/20,-1,-13/10,-27/20,-13/20,-21/20,-17/20,
     526    -11/10,-9/10,-7/10; // roots of BS poly of reiffen(4,5)
     527  intvec v = sortNumberIdeal(I); v;
     528  I[v];
     529}
     530
    503531proc bFactor (poly F)
    504 "USAGE:  bFactor(f);  f poly
     532  "USAGE:  bFactor(f);  f poly
    505533RETURN:  list of ideal and intvec and possibly a string
    506534PURPOSE: tries to compute the roots of a univariate poly f
    … …  
    554582  II = subst(II,var(1),0);
    555583  II = -II;
     584  intvec perm = sortNumberIdeal(II);
     585  II = II[perm];
     586  mm = mm[perm];
    556587  if (size(II)>0)
    557588  {
    … …  
    591622
    592623proc isInt (number n)
    593 "USAGE:  isInt(n); n a number
     624  "USAGE:  isInt(n); n a number
    594625RETURN:  int, 1 if n is an integer or 0 otherwise
    595626PURPOSE: check whether given object of type number is actually an int
    … …  
    619650
    620651proc intRoots (list l)
    621 "USAGE:  isInt(L); L a list
     652  "USAGE:  isInt(L); L a list
    622653RETURN:  list
    623654PURPOSE: extracts integer roots from a list given in @code{bFactor} format
    … …  
    678709
    679710proc sortIntvec (intvec v)
    680 "USAGE:  sortIntvec(v); v an intvec
     711  "USAGE:  sortIntvec(v); v an intvec
    681712RETURN:  list of two intvecs
    682713PURPOSE: sorts an intvec
    … …  
    760791
    761792proc isFsat(ideal I, poly F)
    762 "USAGE:  isFsat(I, F);  I an ideal, F a poly
     793  "USAGE:  isFsat(I, F);  I an ideal, F a poly
    763794RETURN:  int, 1  if I is F-saturated and 0 otherwise
    764795PURPOSE: checks whether the ideal I is F-saturated
    … …  
    799830
    800831proc annRat(poly g, poly f)
    801 "USAGE:   annRat(g,f);  f, g polynomials
     832  "USAGE:   annRat(g,f);  f, g polynomials
    802833RETURN:   ring (a Weyl algebra) containing an ideal 'LD'
    803834PURPOSE:  compute the annihilator of the rational function g/f in the
    … …  
    916947  // Now, compare with the output of Macaulay2:
    917948  ideal tst = 3*x*Dx + 2*y*Dy + 1, y^3*Dy^2 - x^2*Dy^2 + 6*y^2*Dy + 6*y,
    918 9*y^2*Dx^2*Dy-4*y*Dy^3+27*y*Dx^2+2*Dy^2, 9*y^3*Dx^2-4*y^2*Dy^2+10*y*Dy -10;
    919  option(redSB); option(redTail);
    920  LD = groebner(LD);
    921  tst = groebner(tst);
    922  print(matrix(NF(LD,tst)));  print(matrix(NF(tst,LD)));
    923  // So, these two answers are the same
     949    9*y^2*Dx^2*Dy-4*y*Dy^3+27*y*Dx^2+2*Dy^2, 9*y^3*Dx^2-4*y^2*Dy^2+10*y*Dy -10;
     950  option(redSB); option(redTail);
     951  LD = groebner(LD);
     952  tst = groebner(tst);
     953  print(matrix(NF(LD,tst)));  print(matrix(NF(tst,LD)));
     954  // So, these two answers are the same
    924955}
    925956
    926957proc annPoly(poly f)
    927 "USAGE:   annPoly(f);  f a poly
     958  "USAGE:   annPoly(f);  f a poly
    928959RETURN:   ring (a Weyl algebra) containing an ideal 'LD'
    929960PURPOSE:  compute the complete annihilator ideal of f in the corresponding
    … …  
    9921023
    9931024proc DLoc(ideal I, poly F)
    994 "USAGE:  DLoc(I, f);  I an ideal, f a poly
     1025  "USAGE:  DLoc(I, f);  I an ideal, f a poly
    9951026RETURN:  list of ideal and list
    9961027ASSUME:  the basering is a Weyl algebra
    … …  
    10461077
    10471078proc DLoc0(ideal I, poly F)
    1048 "USAGE:  DLoc0(I, f);  I an ideal, f a poly
     1079  "USAGE:  DLoc0(I, f);  I an ideal, f a poly
    10491080RETURN:  ring (a Weyl algebra) containing an ideal 'LD0' and a list 'BS'
    10501081PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s,
    … …  
    10911122  dbprint(ppl,"// -2-2- attempt the factorization");
    10921123  list PP = factorize(p);          //with constants and multiplicities
    1093   ideal bs; intvec m;             //the Bernstein polynomial is monic, so we are not interested in constants
     1124  ideal bs; intvec m;             //the Bernstein polynomial is monic, so
     1125  // we are not interested in constants
    10941126  for (i=2; i<= size(PP[1]); i++)  //we delete P[1][1] and P[2][1]
    10951127  {
    … …  
    11401172      }
    11411173    }
    1142 //     if ( size(vP)>=2 )
    1143 //     {
    1144 //       vP = vP[2..size(vP)];
    1145 //     }
     1174    //     if ( size(vP)>=2 )
     1175    //     {
     1176    //       vP = vP[2..size(vP)];
     1177    //     }
    11461178    if ( size(vP)==0 )
    11471179    {
    … …  
    11591191      dbprint(ppl-1, sP);
    11601192      //    int sP = minIntRoot(bbs,1);
    1161 //       P =  normalize(P);
    1162 //       bs = -subst(bs,s,0);
     1193      //       P =  normalize(P);
     1194      //       bs = -subst(bs,s,0);
    11631195      if (sP >=0)
    11641196      {
    … …  
    12871319
    12881320proc SDLoc(ideal I, poly F)
    1289 "USAGE:  SDLoc(I, f);  I an ideal, f a poly
     1321  "USAGE:  SDLoc(I, f);  I an ideal, f a poly
    12901322RETURN:  ring (basering extended by a new variable) containing an ideal 'LD'
    12911323PURPOSE: compute a generic presentation of the localization of D/I w.r.t. f^s
    … …  
    15051537
    15061538proc GBWeight (ideal I, intvec u, intvec v, list #)
    1507 "USAGE:  GBWeight(I,u,v [,s,t,w]); 
     1539  "USAGE:  GBWeight(I,u,v [,s,t,w]); 
    15081540@*       I ideal, u,v intvecs, s,t optional ints, w an optional intvec
    15091541RETURN:  ideal, Groebner basis of I w.r.t. the weights u and v
    … …  
    16411673
    16421674proc inForm (def I, intvec w)
    1643 "USAGE:  inForm(I,w);  I ideal or poly, w intvec
     1675  "USAGE:  inForm(I,w);  I ideal or poly, w intvec
    16441676RETURN:  ideal, generated by initial forms of generators of I w.r.t. w, or
    16451677@*       poly, initial form of input poly w.r.t. w
    … …  
    17181750
    17191751proc initialIdealW(ideal I, intvec u, intvec v, list #)
    1720 "USAGE:  initialIdealW(I,u,v [,s,t,w]);
     1752  "USAGE:  initialIdealW(I,u,v [,s,t,w]);
    17211753@*       I ideal, u,v intvecs, s,t optional ints, w an optional intvec
    17221754RETURN:  ideal, GB of initial ideal of the input ideal wrt the weights u and v
    … …  
    17771809
    17781810proc initialMalgrange (poly f,list #)
    1779 "USAGE:  initialMalgrange(f,[,a,b,v]); f poly, a,b optional ints, v opt. intvec
     1811  "USAGE:  initialMalgrange(f,[,a,b,v]); f poly, a,b optional ints, v opt. intvec
    17801812RETURN:  ring, Weyl algebra induced by basering, extended by two new vars t,Dt
    17811813PURPOSE: computes the initial Malgrange ideal of a given polynomial w.r.t. the
    … …  
    21432175
    21442176proc restrictionModule (ideal I, intvec w, list #)
    2145 "USAGE:  restrictionModule(I,w,[,eng,m,G]);
     2177  "USAGE:  restrictionModule(I,w,[,eng,m,G]);
    21462178@*       I ideal, w intvec, eng and m optional ints, G optional ideal
    21472179RETURN:  ring (a Weyl algebra) containing a module 'resMod'
    … …  
    22652297
    22662298proc restrictionIdeal (ideal I, intvec w, list #)
    2267 "USAGE:  restrictionIdeal(I,w,[,eng,m,G]);
     2299  "USAGE:  restrictionIdeal(I,w,[,eng,m,G]);
    22682300@*       I ideal, w intvec, eng and m optional ints, G optional ideal
    22692301RETURN:  ring (a Weyl algebra) containing an ideal 'resIdeal'
    … …  
    23172349
    23182350proc fourier (ideal I, list #)
    2319 "USAGE:  fourier(I[,v]); I an ideal, v an optional intvec
     2351  "USAGE:  fourier(I[,v]); I an ideal, v an optional intvec
    23202352RETURN:  ideal
    23212353PURPOSE: computes the Fourier transform of an ideal in a Weyl algebra
    … …  
    23862418
    23872419proc inverseFourier (ideal I, list #)
    2388 "USAGE:  inverseFourier(I[,v]); I an ideal, v an optional intvec
     2420  "USAGE:  inverseFourier(I[,v]); I an ideal, v an optional intvec
    23892421RETURN:  ideal
    23902422PURPOSE: computes the inverse Fourier transform of an ideal in a Weyl algebra
    … …  
    24552487
    24562488proc integralModule (ideal I, intvec w, list #)
    2457 "USAGE:  integralModule(I,w,[,eng,m,G]);
     2489  "USAGE:  integralModule(I,w,[,eng,m,G]);
    24582490@*       I ideal, w intvec, eng and m optional ints, G optional ideal
    24592491RETURN:  ring (a Weyl algebra) containing a module 'intMod'
    … …  
    25692601
    25702602proc integralIdeal (ideal I, intvec w, list #)
    2571 "USAGE:  integralIdeal(I,w,[,eng,m,G]);
     2603  "USAGE:  integralIdeal(I,w,[,eng,m,G]);
    25722604@*       I ideal, w intvec, eng and m optional ints, G optional ideal
    25732605RETURN:  ring (a Weyl algebra) containing an ideal 'intIdeal'
    … …  
    26162648
    26172649proc deRhamCohomIdeal (ideal I, list #)
    2618 "USAGE:  deRhamCohomIdeal (I[,w,eng,k,G]);
     2650  "USAGE:  deRhamCohomIdeal (I[,w,eng,k,G]);
    26192651@*       I ideal, w optional intvec, eng and k optional ints, G optional ideal
    26202652RETURN:  ideal
    … …  
    27792811
    27802812proc deRhamCohom (poly f, list #)
    2781 "USAGE:  deRhamCohom(f[,eng,m]);  f poly, eng and m optional ints
     2813  "USAGE:  deRhamCohom(f[,eng,m]);  f poly, eng and m optional ints
    27822814RETURN:  ring (a Weyl Algebra) containing an ideal 'DR'
    27832815ASSUME:  Basering is a commutative and over a field of characteristic 0.
    … …  
    28902922
    28912923proc appelF1()
    2892 "USAGE:  appelF1();
     2924  "USAGE:  appelF1();
    28932925RETURN:  ring (a parametric Weyl algebra) containing an ideal 'IAppel1'
    28942926PURPOSE: defines the ideal in a parametric Weyl algebra,
    … …  
    29212953
    29222954proc appelF2()
    2923 "USAGE:  appelF2();
     2955  "USAGE:  appelF2();
    29242956RETURN:  ring (a parametric Weyl algebra) containing an ideal 'IAppel2'
    29252957PURPOSE: defines the ideal in a parametric Weyl algebra,
    … …  
    29512983
    29522984proc appelF4()
    2953 "USAGE:  appelF4();
     2985  "USAGE:  appelF4();
    29542986RETURN:  ring (a parametric Weyl algebra) containing an ideal 'IAppel4'
    29552987PURPOSE: defines the ideal in a parametric Weyl algebra,
    … …  
    29843016
    29853017proc charVariety(ideal I, list #)
    2986 "USAGE:  charVariety(I [,eng]); I an ideal, eng an optional int
     3018  "USAGE:  charVariety(I [,eng]); I an ideal, eng an optional int
    29873019RETURN:  ring (commutative) containing an ideal 'charVar'
    29883020PURPOSE: computes an ideal whose zero set is the characteristic variety of I in
    … …  
    30483080
    30493081proc charInfo(ideal I)
    3050 "USAGE:  charInfo(I);  I an ideal
     3082  "USAGE:  charInfo(I);  I an ideal
    30513083RETURN:  ring (commut.) containing ideals 'charVar','singLoc' and list 'primDec'
    30523084PURPOSE: computes characteristic variety of I (in the sense of D-module theory),
    … …  
    31063138
    31073139/*
    3108 static proc exCusp()
    3109 {
     3140  static proc exCusp()
     3141  {
    31103142  "EXAMPLE:"; echo = 2;
    31113143  ring r = 0,(x,y,Dx,Dy),dp;
    … …  
    31233155  setring R;
    31243156  DLoc(I,F);
    3125 }
    3126 
    3127 static proc exWalther1()
    3128 {
     3157  }
     3158
     3159  static proc exWalther1()
     3160  {
    31293161  // p.18 Rem 3.10
    31303162  ring r = 0,(x,Dx),dp;
    … …  
    31433175  LD0;
    31443176  BS;
    3145 }
    3146 
    3147 static proc exWalther2()
    3148 {
     3177  }
     3178
     3179  static proc exWalther2()
     3180  {
    31493181  // p.19 Rem 3.10 cont'd
    31503182  ring r = 0,(x,Dx),dp;
    … …  
    31663198  setring R;
    31673199  DLoc(I,F);
    3168 }
    3169 
    3170 static proc exWalther3()
    3171 {
     3200  }
     3201
     3202  static proc exWalther3()
     3203  {
    31723204  // can check with annFs too :-)
    31733205  // p.21 Ex 3.15
    … …  
    31973229  setring R;
    31983230  DLoc(I,F);
    3199 }
    3200 
    3201 static proc ex_annRat()
    3202 {
     3231  }
     3232
     3233  static proc ex_annRat()
     3234  {
    32033235  // more complicated example for annRat
    32043236  ring r = 0,(x,y,z),dp;
    … …  
    32073239  def A = annRat(g,f);
    32083240  setring A;
    3209 }
     3241  }
    32103242*/

  • Singular/LIB/dmodvar.lib

    rd5abcf rf52e64  
    55LIBRARY: dmodvar.lib     Algebraic D-modules for varieties
    66
    7 AUTHORS: Daniel Andres,          daniel.andres@math.rwth-aachen.de
    8        Viktor Levandovskyy,    levandov@math.rwth-aachen.de
    9        Jorge Martin-Morales,   jorge@unizar.es
    10 
    11 OVERVIEW:
    12 Theory: Let K be a field of characteristic 0. Given a polynomial ring R = K[x_1,...,x_n] and
    13 @* a set of polynomial f_1,..., f_r in R, define F = f_1 * ... * f_r and F^s:=f_1^s_1*...*f_r^s_r
    14 @* for symbolic discrete (that is shiftable) variables s_1,..., s_r.
    15 @* The module R[1/F]*F^s has a structure of a D<S>-module, where
    16 D<S> := D(R) tensored with S over K, where
    17 @*   - D(R) is an n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>
    18 @*   - S is the universal enveloping algebra of gl_r, generated by s_{ij}, where s_{ii}=s_i.
     7AUTHORS: Daniel Andres,        daniel.andres@math.rwth-aachen.de
     8@*       Viktor Levandovskyy,  levandov@math.rwth-aachen.de
     9@*       Jorge Martin-Morales, jorge@unizar.es
     10
     11OVERVIEW: Let K be a field of characteristic 0. Given a polynomial ring R = K[x_1,...,x_n]
     12and polynomials f_1,...,f_r in R, define F = f_1*...*f_r and F^s = f_1^s_1*...*f_r^s_r
     13for symbolic discrete (that is shiftable) variables s_1,..., s_r.
     14The module R[1/F]*F^s has the structure of a D<S>-module, where D<S> = D(R)
     15tensored with S over K, where
     16@* - D(R) is the n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j + 1>
     17@* - S is the universal enveloping algebra of gl_r, generated by s_i = s_{ii}.
    1918@* One is interested in the following data:
    20 @*   - the left ideal Ann F^s in D<S>, usually denoted by LD in the output
    21 @*   - global Bernstein polynomial in one variable s = s_1 + ...+ s_r, denoted by bs,
    22 @*   - its minimal integer root s0, the list of all roots of bs, which are known
    23 @*     to be negative rational numbers, with their multiplicities, which is denoted by BS
    24 @*   - an r-tuple of operators in D<S>, denoted by PS, such that the functional equality
    25 @*     sum(k=1 to k=r) P_k*f_k*F^s = bs*F^s holds in R[1/F]*F^s.
     19@* - the left ideal Ann F^s in D<S>, usually denoted by LD in the output
     20@* - global Bernstein polynomial in one variable s = s_1+...+s_r, denoted by bs,
     21@* - its minimal integer root s0, the list of all roots of bs, which are known to be
     22    negative rational numbers, with their multiplicities, which is denoted by BS
     23@* - an r-tuple of operators in D<S>, denoted by PS, such that the functional equality
     24     sum(k=1 to k=r) P_k*f_k*F^s = bs*F^s holds in R[1/F]*F^s.
    2625
    2726References:
    28   (BMS06) Budur, Mustata, Saito: Bernstein-Sato polynomials of arbitrary varieties (2006).
    29   (ALM09) Andres, Levandovskyy, Martin-Morales : Principal Intersection and Bernstein-Sato Polynomial of an Affine Variety (2009).
     27   (BMS06) Budur, Mustata, Saito: Bernstein-Sato polynomials of arbitrary varieties (2006).
     28@* (ALM09) Andres, Levandovskyy, Martin-Morales: Principal Intersection and Bernstein-Sato
     29Polynomial of an Affine Variety (2009).
     30
    3031
    3132PROCEDURES:
    32 bfctVarIn(F[,L]);     compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using initial ideal approach
    33 bfctVarAnn(F[,L]);  compute the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using Sannfs approach
    34 SannfsVar(F[,O,e]); compute the annihilator of F^s in the ring D<S>
     33bfctVarIn(F[,L]);       computes the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using initial ideal approach
     34bfctVarAnn(F[,L]);      computes the roots of the Bernstein-Sato polynomial b(s) of the variety V(F) using Sannfs approach
     35SannfsVar(F[,O,e]);     computes the annihilator of F^s in the ring D<S>
     36makeMalgrange(F[,ORD]); creates the Malgrange ideal, associated with F = F[1],..,F[P]
    3537
    3638SEE ALSO: bfun_lib, dmod_lib, dmodapp_lib, gmssing_lib
    … …  
    3941Weyl algebra; parametric annihilator for variety; Budur-Mustata-Saito approach; initial ideal approach
    4042";
    41 //AUXILIARY PROCEDURES:
    42 //makeIF(F[,ORD]);    create the Malgrange ideal, associated with F = F[1],..,F[P]
    43 
    44 
     43
     44/*
    4545// Static procs:
    46 //coDim(I);           compute the codimension of the leading ideal of I
     46// coDim(I);           compute the codimension of the leading ideal of I
    4747// dmodvarAssumeViolation()
    4848// ORDstr2list (ORD, NN)
    4949// smallGenCoDim(I,k)
     50*/
     51
     52/*
     53  CHANGELOG
     54  11.10.10 by DA:
     55  - reformated help strings
     56  - simplified code
     57  - add and use of safeVarName
     58  - renamed makeIF to makeMalgrange
     59*/
    5060
    5161
    … …  
    5969proc testdmodvarlib ()
    6070{
    61   "AUXILIARY PROCEDURES:";
    62   example makeIF;
    63   "MAIN PROCEDURES:";
     71  example makeMalgrange;
    6472  example bfctVarIn;
    6573  example bfctVarAnn;
    6674  example SannfsVar;
    6775}
    68 
    6976//   example coDim;
    7077
    … …  
    7380static proc dmodvarAssumeViolation()
    7481{
    75   // returns Boolean : yes/no [for assume violation]
    76   // char K = 0
    77   // no qring
     82  // char K = 0, no qring
    7883  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
    7984  {
    80     //    "ERROR: no qring is allowed";
    81     return(1);
    82   }
    83   return(0);
    84 }
     85    ERROR("Basering is inappropriate: characteristic>0 or qring present");
     86  }
     87  return();
     88}
     89
     90static proc safeVarName (string s, string cv)
     91// assumes 's' to be a valid variable name
     92// returns valid var name string @@..@s
     93{
     94  string S;
     95  if (cv == "v")  { S = "," + "," + varstr(basering)  + ","; }
     96  if (cv == "c")  { S = "," + "," + charstr(basering) + ","; }
     97  if (cv == "cv") { S = "," + charstr(basering) + "," + varstr(basering) + ","; }
     98  s = "," + s + ",";
     99  while (find(S,s) <> 0)
     100  {
     101    s[1] = "@";
     102    s = "," + s;
     103  }
     104  s = s[2..size(s)-1];
     105  return(s)
     106    }
    85107
    86108// da: in smallGenCoDim(), rewritten using mstd business
    87109static proc coDim (ideal I)
    88 "USAGE:  coDim (I);  I an ideal
     110  "USAGE:  coDim (I);  I an ideal
    89111RETURN:  int
    90112PURPOSE: computes the codimension of the ideal generated by the leading monomials
    91 @*       of the given generators of the ideal. This is also the codimension of
    92 @*       the ideal if it is represented by a standard basis.
     113   of the given generators of the ideal. This is also the codimension of
     114   the ideal if it is represented by a standard basis.
    93115NOTE:    The codimension of an ideal I means the number of variables minus the
    94 @*       Krull dimension of the basering modulo I.
    95 EXAMPLE: example SannfsVar; shows examples
     116   Krull dimension of the basering modulo I.
     117EXAMPLE: example coDim; shows examples
    96118"
    97119{
    … …  
    122144
    123145proc SannfsVar (ideal F, list #)
    124 "USAGE:  SannfsVar(F [,ORD,eng]);  F an ideal, ORD an optional string, eng an optional int
    125 RETURN:  ring
    126 PURPOSE: compute the D<S>-module structure of D<S>*f^s where f = F[1]*..*F[P]
    127 and D<S> is the Weyl algebra D tensored with K<S>=U(gl_P), according to the
    128 generalized algorithm by Briancon and Maisonobe for affine varieties.
    129 NOTE:    activate this ring with the @code{setring} command.
    130 @*       In the ring D<S>, the ideal LD is the needed D<S>-module structure.
    131 @*       The value of ORD must be an elimination ordering in D<Dt,S> for Dt
    132 @*       written in the string form, otherwise the result may have no meaning.
    133 @*       By default ORD = '(a(1..(P)..1),a(1..(P+P^2)..1),dp)'.
    134 @*       If eng<>0, @code{std} is used for Groebner basis computations,
    135 @*       otherwise, and by default @code{slimgb} is used.
    136 @*       If printlevel=1, progress debug messages will be printed,
    137 @*       if printlevel>=2, all the debug messages will be printed.
     146  "USAGE:  SannfsVar(F [,ORD,eng]);  F an ideal, ORD an optional string, eng an optional int
     147RETURN:  ring (Weyl algebra tensored with U(gl_P)), containing an ideal LD
     148PURPOSE: compute the D<S>-module structure of D<S>*f^s where f = F[1]*...*F[P]
     149   and D<S> is the Weyl algebra D tensored with K<S>=U(gl_P), according to the
     150   generalized algorithm by Briancon and Maisonobe for affine varieties
     151ASSUME:  The basering is commutative and over a field of characteristic 0.
     152NOTE:    Activate the output ring D<S> with the @code{setring} command.
     153   In the ring D<S>, the ideal LD is the needed D<S>-module structure.
     154@* The value of ORD must be an elimination ordering in D<Dt,S> for Dt
     155   written in the string form, otherwise the result may have no meaning.
     156   By default ORD = '(a(1..(P)..1),a(1..(P+P^2)..1),dp)'.
     157@* If eng<>0, @code{std} is used for Groebner basis computations,
     158   otherwise, and by default @code{slimgb} is used.
     159DISPLAY: If printlevel=1, progress debug messages will be printed,
     160@* if printlevel>=2, all the debug messages will be printed.
    138161EXAMPLE: example SannfsVar; shows examples
    139162"
    140163{
    141   if (dmodvarAssumeViolation())
    142   {
    143     ERROR("Basering is inappropriate: characteristic>0 or qring present");
    144   }
     164  dmodvarAssumeViolation();
    145165  if (!isCommutative())
    146166  {
    … …  
    173193          eng = int(#[2]);
    174194        }
    175         else
    176         {
    177           eng = 0;
    178         }
    179       }
    180       else
    181       {
    182         // no second arg
    183         eng = 0;
    184195      }
    185196    }
    … …  
    215226  for (i=1; i<=P; i++)
    216227  {
    217     RName[i] = "Dt("+string(i)+")";
     228    RName[i] = safeVarName("Dt("+string(i)+")","cv");
    218229    for (j=1; j<=P; j++)
    219230    {
    220       st = "s("+string(i)+")("+string(j)+")";
    221       RName[P+(i-1)*P+j] = st;
    222     }
    223   }
    224   for(i=1; i<=N; i++)
    225   {
    226     for(j=1; j<=size(RName); j++)
    227     {
    228       if (Name[i] == RName[j])
    229       {
    230         ERROR("Variable names should not include Dt(i), s(i)(j)");
    231       }
     231      RName[P+(i-1)*P+j] = safeVarName("s("+string(i)+")("+string(j)+")","cv");
    232232    }
    233233  }
    … …  
    236236  for(i=1; i<=N; i++)
    237237  {
    238     DName[i] = "D"+Name[i];  //concat
     238    DName[i] = safeVarName("D"+Name[i],"cv");  //concat
    239239  }
    240240  list NName = RName + Name + DName;
    … …  
    258258      {
    259259        //[sij,Dtk] = djk*Dti
    260         @D[k,P+(i-1)*P+j] = (j==k)*Dt(i);
     260        //         @D[k,P+(i-1)*P+j] = (j==k)*Dt(i);
     261        @D[k,P+(i-1)*P+j] = (j==k)*var(i);
    261262        for (l=1; l<=P; l++)
    262263        {
    … …  
    264265          {
    265266            //[sij,skl] = djk*sil - dil*skj
    266             @D[P+(i-1)*P+j,P+(k-1)*P+l] = -(j==k)*s(i)(l) + (i==l)*s(k)(j);
     267            //             @D[P+(i-1)*P+j,P+(k-1)*P+l] = -(j==k)*s(i)(l) + (i==l)*s(k)(j);
     268            @D[P+(i-1)*P+j,P+(k-1)*P+l] = -(j==k)*var(i*P+l) + (i==l)*var(k*P+j);
    267269          }
    268270        }
    … …  
    288290    for (j=1; j<=P; j++)
    289291    {
    290       I[(i-1)*P+j] = Dt(i)*F[j] + s(i)(j);
     292      //       I[(i-1)*P+j] = Dt(i)*F[j] + s(i)(j);
     293      I[(i-1)*P+j] = var(i)*F[j] + var(i*P+j);
    291294    }
    292295  }
    … …  
    297300    for (j=1; j<=P; j++)
    298301    {
    299       q = Dt(j);
     302      //       q = Dt(j);
     303      q = var(j);
    300304      q = q*diff(F[j],var(P+P^2+i));
    301305      p = p + q;
    … …  
    304308  }
    305309  // -------- the ideal I is ready ----------
    306   dbprint(ppl,"// -1-2- starting the elimination of "+string(Dt(1..P))+" in @R");
     310  dbprint(ppl,"// -1-2- starting the elimination of Dt(i) in @R");
    307311  dbprint(ppl-1, I);
    308312  ideal J = engine(I,eng);
    … …  
    363367  def A = SannfsVar(F);
    364368  setring A;
     369  A;
    365370  LD;
    366371}
    367372
    368373proc bfctVarAnn (ideal F, list #)
    369 "USAGE:  bfctVarAnn(F[,gid,eng]); F an ideal, gid,eng optional ints
     374  "USAGE:  bfctVarAnn(F[,gid,eng]); F an ideal, gid,eng optional ints
    370375RETURN:  list of an ideal and an intvec
    371376PURPOSE: computes the roots of the Bernstein-Sato polynomial and their multiplicities
    372 @*       for an affine algebraic variety defined by F = F[1],..,F[r].
    373 ASSUME:  The basering is a commutative polynomial ring in char 0.
    374 BACKGROUND: In this proc, the annihilator of f^s in D[s] is computed and then a
    375 @*       system of linear equations is solved by linear reductions in order to
    376 @*       find the minimal polynomial of S = s(1)(1) + ... + s(P)(P)
    377 NOTE:    In the output list, the ideal contains all the roots and the intvec their multiplicities.
    378 @*       If gid<>0, the ideal is used as given. Otherwise, and by default, a
    379 @*       heuristically better suited generating set is used.
    380 @*       If eng<>0, @code{std} is used for GB computations,
    381 @*       otherwise, and by default, @code{slimgb} is used.
     377   for an affine algebraic variety defined by F = F[1],..,F[r].
     378ASSUME:  The basering is commutative and over a field in char 0.
     379NOTE:    In the output list, the ideal contains all the roots and
     380   the intvec their multiplicities.
     381@* If gid<>0, the ideal is used as given. Otherwise, and by default, a
     382   heuristically better suited generating set is used.
     383@* If eng<>0, @code{std} is used for GB computations,
     384   otherwise, and by default, @code{slimgb} is used.
     385@* Computational remark: The time of computation can be very different depending
     386   on the chosen generators of F, although the result is always the same.
     387@* Further note that in this proc, the annihilator of f^s in D[s] is computed and
     388   then a system of linear equations is solved by linear reductions in order to
     389   find the minimal polynomial of S = s(1)(1) + ... + s(P)(P).
     390   The resulted is shifted by 1-codim(Var(F)) following (BMS06).
    382391DISPLAY: If printlevel=1, progress debug messages will be printed,
    383 @*       if printlevel=2, all the debug messages will be printed.
    384 COMPUTATIONAL REMARK: The time of computation can be very different depending
    385 @*       on the chosen generators of F, although the result is always the same.
     392@* if printlevel=2, all the debug messages will be printed.
    386393EXAMPLE: example bfctVarAnn; shows examples
    387394"
    388395{
    389   if (dmodvarAssumeViolation())
    390   {
    391     ERROR("Basering is inappropriate: characteristic>0 or qring present");
    392   }
     396  dmodvarAssumeViolation();
    393397  if (!isCommutative())
    394398  {
    … …  
    420424  def @R2 = SannfsVar(F,eng);
    421425  printlevel = printlevel-1;
     426  int sF = size(F); // no 0 in F
    422427  setring @R2;
    423428  // we are in D[s] and LD is a std of SannfsVar(F)
    … …  
    433438  poly S;
    434439  int i;
    435   for (i=1; i<=size(F); i++)
    436   {
    437     S = S + s(i)(i);
     440  for (i=1; i<=sF; i++)
     441  {
     442    //     S = S + s(i)(i);
     443    S = S + var((i-1)*sF+i);
    438444  }
    439445  dbprint(ppl,"// -4-1- computing the minimal polynomial of S");
    440   //dbprint(ppl-1,"S = "+string(S));
    441   module M = pIntersect(S,K);
     446  dbprint(ppl-1,"S = "+string(S));
     447  vector M = pIntersect(S,K);
    442448  dbprint(ppl,"// -4-2- the minimal polynomial has been computed");
    443   //dbprint(ppl-1,M);
    444449  ring @R3 = 0,s,dp;
    445   dbprint(ppl,"// -5-1- the ring @R3(s) is ready");
    446   dbprint(ppl-1,@R3);
    447   ideal M = imap(@R2,M);
    448   //kill @R2;
    449   poly p;
    450   for (i=1; i<=size(M); i++)
    451   {
    452     p = p + M[i]*s^(i-1);
    453   }
    454   dbprint(ppl,"// -5-2- factorization of the minimal polynomial");
    455   list P = factorize(p);          //with constants and multiplicities
    456   dbprint(ppl-1,P);               //the Bernstein polynomial is monic,
    457   ideal bs; intvec m;             //so we are not interested in constants
    458   for (i=2; i<=ncols(P[1]); i++)   //and that is why we delete P[1][1] and P[2][1]
    459   {
    460     bs[i-1] = P[1][i];
    461     m[i-1]  = P[2][i];
    462   }
    463   // convert factors to a list of their roots and multiplicities
    464   bs =  normalize(bs);
    465   bs = -subst(bs,s,0);
     450  vector M = imap(@R2,M);
     451  poly p = vec2poly(M);
     452  dbprint(ppl-1,p);
     453  dbprint(ppl,"// -5-1- codimension of the variety");
     454  dbprint(ppl-1,cd);
     455  dbprint(ppl,"// -5-2- shifting BS(s)=minpoly(s-codim+1)");
     456  p = subst(p,var(1),var(1)-cd+1);
     457  dbprint(ppl-1,p);
     458  dbprint(ppl,"// -5-3- factorization of the minimal polynomial");
     459  list BS = bFactor(p);
    466460  setring save;
    467 //   ideal GF = imap(@R2,GF);
    468 //   attrib(GF,"isSB",1);
    469   kill @R2;
    470   dbprint(ppl,"// -5-3- codimension of the variety");
    471 //   int cd = coDim(GF);
    472   dbprint(ppl-1,cd);
    473   ideal bs = imap(@R3,bs);
    474   dbprint(ppl,"// -5-4- shifting BS(s)=minpoly(s-codim+1)");
    475   for (i=1; i<=ncols(bs); i++)
    476   {
    477     bs[i] = bs[i] + cd - 1;
    478   }
    479   kill @R3;
    480   list BS = bs,m;
     461  list BS = imap(@R3,BS);
     462  kill @R2,@R3;
    481463  return(BS);
    482464}
    … …  
    489471}
    490472
    491 static proc makeIF (ideal F, list #)
    492 "USAGE:  makeIF(F [,ORD]);  F an ideal, ORD an optional string
    493 RETURN:  ring
    494 PURPOSE: create the ideal by Malgrange associated with F = F[1],..,F[P].
    495 NOTE:    activate this ring with the @code{setring} command. In this ring,
    496 @*       - the ideal IF is the ideal by Malgrange corresponding to F.
    497 @*       The value of ORD must be an arbitrary ordering in K<_t,_x,_Dt,_Dx>
    498 @*       written in the string form. By default ORD = 'dp'.
    499 @*      If printlevel=1, progress debug messages will be printed,
    500 @*       if printlevel>=2, all the debug messages will be printed.
    501 EXAMPLE: example makeIF; shows examples
     473proc makeMalgrange (ideal F, list #)
     474  "USAGE:  makeMalgrange(F [,ORD]);  F an ideal, ORD an optional string
     475RETURN:  ring (Weyl algebra) containing an ideal IF
     476PURPOSE: create the ideal by Malgrange associated with F = F[1],...,F[P].
     477NOTE:    Activate the output ring with the @code{setring} command. In this ring,
     478  the ideal IF is the ideal by Malgrange corresponding to F.
     479@* The value of ORD must be an arbitrary ordering in K<_t,_x,_Dt,_Dx>
     480   written in the string form. By default ORD = 'dp'.
     481DISPLAY: If printlevel=1, progress debug messages will be printed,
     482@* if printlevel>=2, all the debug messages will be printed.
     483EXAMPLE: example makeMalgrange; shows examples
    502484"
    503485{
    … …  
    528510  for (i=1; i<=P; i++)
    529511  {
    530     TName[i] = "t("+string(i)+")";
    531     DTName[i] = "Dt("+string(i)+")";
    532   }
    533   for (i=1; i<=N; i++)
    534   {
    535     for (j=1; j<=P; j++)
    536     {
    537       if (Name[i] == TName[j])
    538       {
    539         ERROR("Variable names should not include t(i)");
    540       }
    541     }
     512    TName[i] = safeVarName("t("+string(i)+")","cv");
     513    DTName[i] = safeVarName("Dt("+string(i)+")","cv");
    542514  }
    543515  //now, create the names for new vars
    … …  
    545517  for (i=1; i<=N; i++)
    546518  {
    547     DName[i] = "D"+Name[i];  //concat
     519    DName[i] = safeVarName("D"+Name[i],"cv");  //concat
    548520  }
    549521  list NName = TName + Name + DTName + DName;
    … …  
    556528  def @R@ = ring(L);
    557529  setring @R@;
    558   matrix @D[Nnew][Nnew];
    559   for (i=1; i<=N+P; i++)
    560   {
    561     @D[i,i+N+P]=1;
    562   }
    563   def @R = nc_algebra(1,@D);
     530  def @R = Weyl();
    564531  setring @R;
    565532  kill @R@;
    … …  
    571538  for (j=1; j<=P; j++)
    572539  {
    573     I[j] = t(j) - F[j];
     540    //     I[j] = t(j) - F[j];
     541    I[j] = var(j) - F[j];
    574542  }
    575543  poly p,q;
    … …  
    579547    for (j=1; j<=P; j++)
    580548    {
    581       q = Dt(j);
     549      //       q = Dt(j);
     550      q = var(P+N+j);
    582551      q = diff(F[j],var(P+i))*q;
    583552      p = p + q;
    … …  
    595564  ring R = 0,(x,y,z),Dp;
    596565  ideal I = x^2+y^3, z;
    597   def W = makeIF(I);
     566  def W = makeMalgrange(I);
    598567  setring W;
     568  W;
    599569  IF;
    600570}
    601571
    602572proc bfctVarIn (ideal I, list #)
    603 "USAGE:  bfctVarIn(I [,a,b,c]);  I an ideal, a,b,c optional ints
     573  "USAGE:  bfctVarIn(I [,a,b,c]);  I an ideal, a,b,c optional ints
    604574RETURN:  list of ideal and intvec
    605575PURPOSE: computes the roots of the Bernstein-Sato polynomial and their
    606 @*       multiplicities for an affine algebraic variety defined by I.
    607 ASSUME:  The basering is commutative and of characteristic 0.
    608 @*       Varnames of the basering do not include t(1),...,t(r) and
    609 @*       Dt(1),...,Dt(r), where r is the number of entries of the input ideal.
    610 BACKGROUND:  In this proc, the initial ideal of the multivariate Malgrange ideal
    611 @*       defined by I is computed and then a system of linear equations is solved
    612 @*       by linear reductions following the ideas by Noro.
     576   multiplicities for an affine algebraic variety defined by I.
     577ASSUME:  The basering is commutative and over a field of characteristic 0.
     578@* Varnames of the basering do not include t(1),...,t(r) and
     579   Dt(1),...,Dt(r), where r is the number of entries of the input ideal.
    613580NOTE:    In the output list, say L,
    614 @*       - L[1] of type ideal contains all the rational roots of a b-function,
    615 @*       - L[2] of type intvec contains the multiplicities of above roots,
    616 @*       - optional L[3] of type string is the part of b-function without
    617 @*        rational roots.
    618 @*       Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and
    619 @*       L[3] is 1 (for nonzero constant) or 0 (for zero b-function).
    620 @*       If a<>0, the ideal is used as given. Otherwise, and by default, a
    621 @*       heuristically better suited generating set is used to reduce computation
    622 @*       time.
    623 @*       If b<>0, @code{std} is used for GB computations in characteristic 0,
    624 @*       otherwise, and by default, @code{slimgb} is used.
    625 @*       If c<>0, a matrix ordering is used for GB computations, otherwise,
    626 @*       and by default, a block ordering is used.
     581@* - L[1] of type ideal contains all the rational roots of a b-function,
     582@* - L[2] of type intvec contains the multiplicities of above roots,
     583@* - optional L[3] of type string is the part of b-function without rational roots.
     584@* Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and
     585   L[3] is 1 (for nonzero constant) or 0 (for zero b-function).
     586@* If a<>0, the ideal is used as given. Otherwise, and by default, a
     587   heuristically better suited generating set is used to reduce computation time.
     588@* If b<>0, @code{std} is used for GB computations in characteristic 0,
     589   otherwise, and by default, @code{slimgb} is used.
     590@* If c<>0, a matrix ordering is used for GB computations, otherwise,
     591   and by default, a block ordering is used.
     592@* Further note, that in this proc, the initial ideal of the multivariate Malgrange
     593   ideal defined by I is computed and then a system of linear equations is solved
     594   by linear reductions following the ideas by Noro.
     595   The result is shifted by 1-codim(Var(F)) following (BMS06).
    627596DISPLAY: If printlevel=1, progress debug messages will be printed,
    628 @*       if printlevel>=2, all the debug messages will be printed.
     597@* if printlevel>=2, all the debug messages will be printed.
    629598EXAMPLE: example bfctVarIn; shows examples
    630599"
    631600{
    632   if (dmodvarAssumeViolation())
    633   {
    634     ERROR("Basering is inappropriate: characteristic>0 or qring present");
    635   }
     601  dmodvarAssumeViolation();
    636602  if (!isCommutative())
    637603  {
    … …  
    674640  // step 1: setting up the multivariate Malgrange ideal
    675641  int r = size(I);
    676   def D = makeIF(I);
     642  def D = makeMalgrange(I);
    677643  setring D;
    678644  dbprint(ppl-1,"// Computing in " + string(n+r) + "-th Weyl algebra:", D);
    … …  
    692658  {
    693659    ideal B = L[1];
    694     for (i=1; i<=ncols(B); i++)
    695     {
    696       B[i] = -B[i]+c-r-1;
    697     }
    698     L[1] = B;
     660    ideal BB;
     661    int nB = ncols(B);
     662    for (i=nB; i>0; i--)
     663    {
     664      BB[i] = -B[nB-i+1]+c-r-1;
     665    }
     666    L[1] = BB;
    699667  }
    700668  else // should never get here: BS poly has non-rational roots
    … …  
    722690static proc smallGenCoDim (ideal I, int Iasgiven)
    723691{
    724   // call from K[x]
    725   // returns list L
     692  // call from K[x], returns list L
    726693  // L[1]=I or L[1]=smaller generating set of I
    727694  // L[2]=codimension(I)
    728   int ppl = printlevel - voice + 3;
     695  int ppl = printlevel - voice + 2;
    729696  int n = nvars(basering);
    730697  int c;
    … …  
    791758}
    792759
    793 
     760/*
    794761// Some more examples
    795762
    796763static proc TXcups()
    797764{
    798   "EXAMPLE:"; echo = 2;
    799   //TX tangent space of X=V(x^2+y^3)
    800   ring R = 0,(x0,x1,y0,y1),Dp;
    801   ideal F = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
    802   printlevel = 0;
    803   //ORD = "(a(1,1),a(1,1,1,1,1,1),dp)";
    804   //eng = 0;
    805   def A = SannfsVar(F);
    806   setring A;
    807   LD;
     765"EXAMPLE:"; echo = 2;
     766//TX tangent space of X=V(x^2+y^3)
     767ring R = 0,(x0,x1,y0,y1),Dp;
     768ideal F = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
     769printlevel = 0;
     770//ORD = "(a(1,1),a(1,1,1,1,1,1),dp)";
     771//eng = 0;
     772def A = SannfsVar(F);
     773setring A;
     774LD;
    808775}
    809776
    810777static proc ex47()
    811778{
    812   ring r7 = 0,(x0,x1,y0,y1),dp;
    813   ideal I = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
    814   bfctVarIn(I);
    815   // second ex - too big
    816   ideal J = x0^4+y0^5, 4*x0^3*x1+5*y0^4*y1;
    817   bfctVarIn(J);
     779ring r7 = 0,(x0,x1,y0,y1),dp;
     780ideal I = x0^2+y0^3, 2*x0*x1+3*y0^2*y1;
     781bfctVarIn(I);
     782// second ex - too big
     783ideal J = x0^4+y0^5, 4*x0^3*x1+5*y0^4*y1;
     784bfctVarIn(J);
    818785}
    819786
    820787static proc ex48()
    821788{
    822   ring r8 = 0,(x1,x2,x3),dp;
    823   ideal I = x1^3-x2*x3, x2^2-x1*x3, x3^2-x1^2*x2;
    824   bfctVarIn(I);
     789ring r8 = 0,(x1,x2,x3),dp;
     790ideal I = x1^3-x2*x3, x2^2-x1*x3, x3^2-x1^2*x2;
     791bfctVarIn(I);
    825792}
    826793
    827794static proc ex49 ()
    828795{
    829   ring r9 = 0,(z1,z2,z3,z4),dp;
    830   ideal I = z3^2-z2*z4, z2^2*z3-z1*z4, z2^3-z1*z3;
    831   bfctVarIn(I);
     796ring r9 = 0,(z1,z2,z3,z4),dp;
     797ideal I = z3^2-z2*z4, z2^2*z3-z1*z4, z2^3-z1*z3;
     798bfctVarIn(I);
    832799}
    833800
    834801static proc ex410()
    835802{
    836   LIB "toric.lib";
    837   ring r = 0,(z(1..7)),dp;
    838   intmat A[3][7];
    839   A = 6,4,2,0,3,1,0,0,1,2,3,0,1,0,0,0,0,0,1,1,2;
    840   ideal I = toric_ideal(A,"pt");
    841   I = std(I);
     803LIB "toric.lib";
     804ring r = 0,(z(1..7)),dp;
     805intmat A[3][7];
     806A = 6,4,2,0,3,1,0,0,1,2,3,0,1,0,0,0,0,0,1,1,2;
     807ideal I = toric_ideal(A,"pt");
     808I = std(I);
    842809//ideal I = z(6)^2-z(3)*z(7), z(5)*z(6)-z(2)*z(7), z(5)^2-z(1)*z(7),
    843810//  z(4)*z(5)-z(3)*z(6), z(3)*z(5)-z(2)*z(6), z(2)*z(5)-z(1)*z(6),
    844811//  z(3)^2-z(2)*z(4), z(2)*z(3)-z(1)*z(4), z(2)^2-z(1)*z(3);
    845   bfctVarIn(I,1); // no result yet
    846 }
     812bfctVarIn(I,1); // no result yet
     813}
     814*/

  • Singular/LIB/ncfactor.lib

    rd5abcf rf52e64  
    979979    result = list(list(1,h));
    980980  }//only the trivial factorization could be found
    981   return(result);
     981  //now, refine the possible redundant list
     982  return( refineFactList(result) );
    982983}//proc facFirstWeyl
    983984example
    … …  
    15951596    }//Nontrivial factorization
    15961597  }//recursively factorize factors
    1597   return(result);
     1598  //now, refine the possible redundant list
     1599  return( refineFactList(result) );
    15981600}//facFirstShift
    15991601example
    … …  
    18411843
    18421844*/
     1845
     1846/* more things from Martin Lee to fix:
     1847
     1848ring R = 0,(x,s),dp;
     1849def r = nc_algebra(1,s);
     1850setring(r);
     1851poly h = (s2*x+x)*s;
     1852h= h* (x+s);
     1853def l= facFirstShift(h);
     1854l; // contained doubled entries
     1855
     1856
     1857ring R = 0,(x,s),dp;
     1858def r = nc_algebra(1,-1);
     1859setring(r);
     1860poly h = (s2*x+x)*s;
     1861h= h* (x+s);
     1862def l= facFirstWeyl(h);
     1863l; // contained doubled entries: not anymore, fixed!
     1864
     1865*/
     1866
     1867static proc refineFactList(list L)
     1868{
     1869  // assume: list L is an output of factorization proc
     1870  // doing: remove doubled entries
     1871  int s = size(L); int sm;
     1872  int i,j,k,cnt;
     1873  list M, U, A, B;
     1874  A = L;
     1875  k = 0;
     1876  cnt  = 1;
     1877  for (i=1; i<=s; i++)
     1878  {
     1879    if (size(A[i]) != 0)
     1880    {
     1881      M = A[i];
     1882      //      "probing with"; M; i;
     1883      B[cnt] = M; cnt++;
     1884      for (j=i+1; j<=s; j++)
     1885      {
     1886        //      if ( (size(L[j]) == sm) && (isEqualList(M,L[j])) )
     1887        if ( isEqualList(M,A[j]) )
     1888        {
     1889          k++;
     1890        // U consists of intvecs with equal pairs
     1891          U[k] = intvec(i,j);
     1892          A[j] = 0;
     1893        }
     1894      }
     1895    }
     1896  }
     1897  kill A,U,M;
     1898  return(B);
     1899}
     1900example
     1901{
     1902  "EXAMPLE:";echo=2;
     1903  ring R = 0,(x,s),dp;
     1904  def r = nc_algebra(1,1);
     1905  setring(r);
     1906  list l,m;
     1907  l = list(1,s2+1,x,s,x+s);
     1908  m = l,list(1,s,x,s,x),l;
     1909  refineFactList(m);
     1910}
     1911
     1912static proc isEqualList(list L, list M)
     1913{
     1914  // int boolean: 1=yes, 0 =no : test whether two lists are identical
     1915  int s = size(L);
     1916  if (size(M)!=s) { return(0); }
     1917  int j=1;
     1918  while ( (L[j]==M[j]) && (j<s) )
     1919  {
     1920    j++;
     1921  }
     1922  if (L[j]==M[j])
     1923  {
     1924    return(1);
     1925  }
     1926  return(0);
     1927}
     1928example
     1929{
     1930  "EXAMPLE:";echo=2;
     1931  ring R = 0,(x,s),dp;
     1932  def r = nc_algebra(1,1);
     1933  setring(r);
     1934  list l,m;
     1935  l = list(1,s2+1,x,s,x+s);
     1936  m = l;
     1937  isEqualList(m,l);
     1938}
Note: See TracChangeset for help on using the changeset viewer.