# source:git/Singular/LIB/bfun.lib@5ec7e3

jengelh-datetimespielwiese
Last change on this file since 5ec7e3 was 5ec7e3, checked in by Hans Schoenemann <hannes@…>, 12 years ago
• Property mode set to `100755`
File size: 54.0 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="\$Id\$";
3category="Noncommutative";
4info="
5LIBRARY: bfun.lib     Algorithms for b-functions and Bernstein-Sato polynomial
6AUTHORS: Daniel Andres, daniel.andres@math.rwth-aachen.de
7@* Viktor Levandovskyy,      levandov@math.rwth-aachen.de
8
9OVERVIEW:
10Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R,
11one is interested in the global b-function (also known as Bernstein-Sato
12polynomial) b(s) in K[s], defined to be the non-zero monic polynomial of minimal
13degree, satisfying a functional identity L * F^{s+1} = b(s) F^s,
14for some operator L in D[s] (* stands for the action of differential operator)@*
15By D one denotes the n-th Weyl algebra
16K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>.
17One is interested in the following data:@*
18  - Bernstein-Sato polynomial b(s) in K[s],@*
19  - the list of its roots, which are known to be rational@*
20  - the multiplicities of the roots.@*
21
22There is a constructive definition of a b-function of a holonomic ideal I in D
23(that is, an ideal I in a Weyl algebra D, such that D/I is holonomic module)
24with respect to the given weight vector w: For a polynomial p in D, its initial
25form w.r.t. (-w,w) is defined as the sum of all terms of p which have
26maximal weighted total degree where the weight of x_i is -w_i and the weight
27of d_i is w_i. Let J be the initial ideal of I w.r.t. (-w,w), i.e. the
28K-vector space generated by all initial forms w.r.t (-w,w) of elements of I.
29Put s = w_1 x_1 d_1 + ... + w_n x_n d_n. Then the monic generator b_w(s) of
30the intersection of J with the PID K[s] is called the b-function of I w.r.t.  w.
31Unlike Bernstein-Sato polynomial, general b-function with respect to
32arbitrary weights need not have rational roots at all. However, b-function
33of a holonomic ideal is known to be non-zero as well.
34
35REFERENCES: [SST] Saito, Sturmfels, Takayama: Groebner Deformations of
36Hypergeometric Differential Equations (2000),@*
37Noro: An Efficient Modular Algorithm for Computing the Global b-function,
38(2002).
39
40
41PROCEDURES:
42bfct(f[,s,t,v]);            compute the BS polynomial of f with linear reductions
43bfctSyz(f[,r,s,t,u,v]);     compute the BS polynomial of f with syzygy-solver
44bfctAnn(f[,s]);             compute the BS polynomial of f via Ann f^s + f
45bfctOneGB(f[,s,t]);         compute the BS polynomial of f by just one GB computation
46bfctIdeal(I,w[,s,t]);       compute the b-function of ideal w.r.t. weight
47pIntersect(f,I[,s]);        intersection of ideal with subalgebra K[f] for a polynomial f
48pIntersectSyz(f,I[,p,s,t]); intersection of ideal with subalgebra K[f] with syz-solver
49linReduce(f,I[,s]);         reduce a polynomial by linear reductions w.r.t. ideal
50linReduceIdeal(I,[s,t]);    interreduce generators of ideal by linear reductions
51linSyzSolve(I[,s]);         compute a linear dependency of elements of ideal I
52
53allPositive(v);             checks whether all entries of an intvec are positive
54scalarProd(v,w);            computes the standard scalar product of two intvecs
55vec2poly(v[,i]);            constructs an univariate polynomial with given coefficients
56
57
59
60KEYWORDS: D-module; global Bernstein-Sato polynomial; Bernstein-Sato polynomial; b-function;
61graded Weyl algebra; initial ideal; initial form; principal intersection; linear interreduction;
62initial ideal approach
63";
64
65
66
67/*
68CHANGELOG
69
7003.03.11:
71- simplified scalarProd
72- fixed bug in bfct when user used vars t,Dt
73- now bFactor is used by bfct, bfctAnn, i.e. the static procs
75- fixed printlevel/debug message issue in bfct, bfctAnn
76- fixed small issue for zero ideal input in linReduceIdeal
77*/
78
79
80
81LIB "qhmoduli.lib"; // for Max
82LIB "dmod.lib";     // for SannfsBFCT etc
83LIB "dmodapp.lib";  // for initialIdealW etc
84LIB "nctools.lib";  // for isWeyl etc
85LIB "presolve.lib"; // for valvars
86
87//--------------- auxiliary procedures ----------------------------------------
88
89/*
90static proc gradedWeyl (intvec u,intvec v)
92RETURN:  a ring, the associated graded ring of the basering w.r.t. u and v
93PURPOSE: compute the associated graded ring of the basering w.r.t. u and v
94ASSUME: basering is a Weyl algebra
96NOTE:    u[i] is the weight of x(i), v[i] the weight of D(i).
97@*       u+v has to be a non-negative intvec.
98"
99{
100  // todo check assumption
101  int i;
102  def save = basering;
103  int n = nvars(save)/2;
104  if (nrows(u)<>n || nrows(v)<>n)
105  {
106    ERROR("weight vectors have wrong dimension");
107  }
108  intvec uv,gr;
109  uv = u+v;
110  for (i=1; i<=n; i++)
111  {
112    if (uv[i]>=0)
113    {
114      if (uv[i]==0)
115      {
116        gr[i] = 0;
117      }
118      else
119      {
120        gr[i] = 1;
121      }
122    }
123    else
124    {
125      ERROR("the sum of the weight vectors has to be a non-negative intvec");
126    }
127  }
128  list l = ringlist(save);
129  list l2 = l[2];
130  matrix l6 = l[6];
131  for (i=1; i<=n; i++)
132  {
133    if (gr[i] == 1)
134    {
135       l2[n+i] = "xi("+string(i)+")";
136       l6[i,n+i] = 0;
137    }
138  }
139  l[2] = l2;
140  l[6] = l6;
141  def G = ring(l);
142  return(G);
143}
144example
145{
146  "EXAMPLE:"; echo = 2;
147  ring @D = 0,(x,y,z,Dx,Dy,Dz),dp;
148  def D = Weyl();
149  setring D;
150  intvec u = -1,-1,1; intvec v = 2,1,1;
152  setring G; G;
153}
154*/
155
156static proc safeVarName (string s)
157{
158  string S = "," + charstr(basering) + "," + varstr(basering) + ",";
159  s = "," + s + ",";
160  while (find(S,s) <> 0)
161  {
162    s[1] = "@";
163    s = "," + s;
164  }
165  s = s[2..size(s)-1];
166  return(s)
167}
168
169proc allPositive (intvec v)
170"USAGE:  allPositive(v);  v an intvec
171RETURN:  int, 1 if all components of v are positive, or 0 otherwise
172PURPOSE: check whether all components of an intvec are positive
173EXAMPLE: example allPositive; shows an example
174"
175{
176  int i;
177  for (i=1; i<=size(v); i++)
178  {
179    if (v[i]<=0)
180    {
181      return(0);
182      break;
183    }
184  }
185  return(1);
186}
187example
188{
189  "EXAMPLE:"; echo = 2;
190  intvec v = 1,2,3;
191  allPositive(v);
192  intvec w = 1,-2,3;
193  allPositive(w);
194}
195
196static proc findFirst (list l, i)
197"USAGE:  findFirst(l,i);  l a list, i an argument of any type
198RETURN:  int, the position of the first appearance of i in l,
199@*       or 0 if i is not a member of l
200PURPOSE: check whether the second argument is a member of a list
201EXAMPLE: example findFirst; shows an example
202"
203{
204  int j;
205  for (j=1; j<=size(l); j++)
206  {
207    if (l[j]==i)
208    {
209      return(j);
210      break;
211    }
212  }
213  return(0);
214}
215//   findFirst(list(1, 2, list(1)),2);       // seems to be a bit buggy,
216//   findFirst(list(1, 2, list(1)),3);       // but works for the purposes
217//   findFirst(list(1, 2, list(1)),list(1)); // of this library
218//   findFirst(l,list(2));
219example
220{
221  "EXAMPLE:"; echo = 2;
222  ring r = 0,x,dp;
223  list l = 1,2,3;
224  findFirst(l,4);
225  findFirst(l,2);
226}
227
228proc scalarProd (intvec v, intvec w)
229"USAGE:  scalarProd(v,w);  v,w intvecs
230RETURN:  int, the standard scalar product of v and w
231PURPOSE: computes the scalar product of two intvecs
232ASSUME:  the arguments are of the same size
233EXAMPLE: example scalarProd; shows examples
234"
235{
236  if (size(v)!=size(w))
237  {
238    ERROR("non-matching dimensions");
239  }
240  else
241  {
242    intvec u = transpose(v)*w;
243  }
244  return(u[1]);
245}
246example
247{
248  "EXAMPLE:"; echo = 2;
249  intvec v = 1,2,3;
250  intvec w = 4,5,6;
251  scalarProd(v,w);
252}
253
254//-------------- main procedures -------------------------------------------------------
255proc linReduceIdeal(ideal I, list #)
256"USAGE:  linReduceIdeal(I [,s,t,u]); I an ideal, s,t,u optional ints
257RETURN:  ideal or list, linear reductum (over field) of I by its elements
258PURPOSE: reduces a list of polys only by linear reductions (no monomial
259@*        multiplications)
260NOTE:    If s<>0, a list consisting of the reduced ideal and the coefficient
261@*       vectors of the used reductions given as module is returned.
262@*       Otherwise (and by default), only the reduced ideal is returned.
263@*       If t<>0 (and by default) all monomials are reduced (if possible),
264@*       otherwise, only leading monomials are reduced.
265@*       If u<>0 (and by default), the ideal is first sorted in increasing order.
266@*       If u is set to 0 and the given ideal is not sorted in the way described,
267@*       the result might not be as expected.
268DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
269@*       if printlevel>=2, all the debug messages will be printed.
270EXAMPLE: example linReduceIdeal; shows examples
271"
272{
273  // #[1] = remembercoeffs
274  // #[2] = redtail
275  // #[3] = sortideal
276  int ppl = printlevel - voice + 2;
277  int remembercoeffs = 0; // default
278  int redtail        = 1; // default
279  int sortideal      = 1; // default
280  if (size(#)>0)
281  {
282    if (typeof(#[1])=="int" || typeof(#[1])=="number")
283    {
284      remembercoeffs = #[1];
285    }
286    if (size(#)>1)
287    {
288      if (typeof(#[2])=="int" || typeof(#[2])=="number")
289      {
290        redtail = #[2];
291      }
292      if (size(#)>2)
293      {
294        if (typeof(#[3])=="int" || typeof(#[3])=="number")
295        {
296          sortideal = #[3];
297        }
298      }
299    }
300  }
301  int sI = ncols(I);
302  int sZeros = sI - size(I);
303  int i,j;
304  ideal J,lmJ,ordJ;
305  list lJ = sort(I);
306  module M; // for the coefficients
307  // step 1: prepare, e.g. sort I
308  if (sortideal <> 0)
309  {
310    if (sZeros > 0) // I contains zeros
311    {
312      if (remembercoeffs <> 0)
313      {
314        j = 0;
315        for (i=1; i<=sI; i++)
316        {
317          if (I[i] == 0)
318          {
319            j++;
320            J[j] = 0;
321            ordJ[j] = -1;
322            M[j] = gen(i);
323          }
324          else
325          {
326            M[i+sZeros-j] = gen(lJ[2][i-j]+j);
327          }
328        }
329      }
330      else
331      {
332        for (i=1; i<=sZeros; i++)
333        {
334          J[i] = 0;
335          ordJ[i] = -1;
336        }
337      }
338      I = J,lJ[1];
339    }
340    else // I contains no zeros
341    {
342      I = lJ[1];
343      if (remembercoeffs <> 0)
344      {
345        for (i=1; i<=size(lJ[1]); i++) { M[i+sZeros] = gen(lJ[2][i]); }
346      }
347    }
348  }
349  else // assume I is already sorted
350  {
351    if (remembercoeffs <> 0)
352    {
353      for (i=1; i<=ncols(I); i++) { M[i] = gen(i);  }
354    }
355  }
356  dbprint(ppl-1,"// initially sorted ideal:", I);
357  if (remembercoeffs <> 0) { dbprint(ppl-1,"// used permutations:", M); }
358  // step 2: reduce leading monomials by linear reductions
359  poly lm,c,redpoly,maxlmJ;
360  J[sZeros+1] = I[sZeros+1];
361  lm = I[sZeros+1];
363  lmJ[sZeros+1] = maxlmJ;
364  int ordlm,reduction;
365  ordJ[sZeros+1] = ord(lm);
366  vector v;
367  int noRedPast;
368  for (i=sZeros+2; i<=sI; i++)
369  {
370    redpoly = I[i];
372    ordlm = ord(lm);
373    if (remembercoeffs <> 0) { v = M[i]; }
374    reduction = 1;
375    while (reduction == 1) // while there was a reduction
376    {
377      noRedPast = i;
378      reduction = 0;
379      for (j=sZeros+1; j<noRedPast; j++)
380      {
381        if (lm == 0) { break; } // nothing more to reduce
382        if (lm > maxlmJ) { break; } //lm is bigger than maximal monomial to reduce with
383        if (ordlm == ordJ[j])
384        {
385          if (lm == lmJ[j])
386          {
387            dbprint(ppl-1,"// reducing " + string(redpoly));
388            dbprint(ppl-1,"//  with " + string(J[j]));
390            redpoly = redpoly - c*J[j];
391            dbprint(ppl-1,"//  to " + string(redpoly));
393            ordlm = ord(lm);
394            if (remembercoeffs <> 0) { M[i] = M[i] - c * M[j]; }
395            noRedPast = j;
396            reduction = 1;
397          }
398        }
399      }
400    }
401    for (j=sZeros+1; j<i; j++)
402    {
403      if (redpoly < J[j]) { break; }
404    }
405    J = insertGenerator(J,redpoly,j);
406    lmJ = insertGenerator(lmJ,lm,j);
407    ordJ = insertGenerator(ordJ,poly(ordlm),j);
408    if (remembercoeffs <> 0)
409    {
410      v = M[i];
411      M = deleteGenerator(M,i);
412      M = insertGenerator(M,v,j);
413    }
414    maxlmJ = lmJ[i];
415  }
416  // step 3: reduce tails by linear reductions as well
417  if (redtail <> 0)
418  {
419    dbprint(ppl,"// finished reducing leading monomials");
420    dbprint(ppl-1,string(J));
421    if (remembercoeffs <> 0) { dbprint(ppl-1,"// used reductions:" + string(M)); }
422    int k;
423    for (i=sZeros+1; i<=sI; i++)
424    {
425      lm = lmJ[i];
426      for (j=i+1; j<=sI; j++)
427      {
428        for (k=2; k<=size(J[j]); k++) // run over all terms in J[j]
429        {
430          if (ordJ[i] == ord(J[j][k]))
431          {
432            if (lm == normalize(J[j][k]))
433            {
435              dbprint(ppl-1,"// reducing " + string(J[j]));
436              dbprint(ppl-1,"//  with " + string(J[i]));
437              J[j] = J[j] - c*J[i];
438              dbprint(ppl-1,"//  to " + string(J[j]));
439              if (remembercoeffs <> 0) { M[j] = M[j] - c * M[i]; }
440            }
441          }
442        }
443      }
444    }
445  }
446  if (sI == sZeros) // if I=0,0,...,0, we now have one too much by construction due to sort
447  {
448    J = J[1..sZeros];
449  }
450  if (remembercoeffs <> 0) { return(list(J,M)); }
451  else { return(J); }
452}
453example
454{
455  "EXAMPLE:"; echo = 2;
456  ring r = 0,(x,y),dp;
457  ideal I = 3,x+9,y4+5x,2y4+7x+2;
458  linReduceIdeal(I);     // reduces tails
459  linReduceIdeal(I,0,0); // no reductions of tails
460  list l = linReduceIdeal(I,1); // reduces tails and shows reductions used
461  l;
462  module M = I;
463  matrix(l[1]) - M*l[2];
464}
465
466proc linReduce(poly f, ideal I, list #)
467"USAGE:  linReduce(f, I [,s,t,u]); f a poly, I an ideal, s,t,u optional ints
468RETURN:  poly or list, linear reductum (over field) of f by elements from I
469PURPOSE: reduce a polynomial only by linear reductions (no monomial multiplications)
470NOTE:    If s<>0, a list consisting of the reduced polynomial and the coefficient
471@*       vector of the used reductions is returned, otherwise (and by default)
472@*       only reduced polynomial is returned.
473@*       If t<>0 (and by default) all monomials are reduced (if possible),
474@*       otherwise, only leading monomials are reduced.
475@*       If u<>0 (and by default), the ideal is linearly pre-reduced, i.e.
476@*       instead of the given ideal, the output of @code{linReduceIdeal} is used.
477@*       If u is set to 0 and the given ideal does not equal the output of
478@*       @code{linReduceIdeal}, the result might not be as expected.
479DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
480@*       if printlevel>=2, all the debug messages will be printed.
481EXAMPLE: example linReduce; shows examples
482"
483{
484  int ppl = printlevel - voice + 2;
485  int remembercoeffs = 0; // default
486  int redtail        = 1; // default
487  int prepareideal   = 1; // default
488  if (size(#)>0)
489  {
490    if (typeof(#[1])=="int" || typeof(#[1])=="number")
491    {
492      remembercoeffs = #[1];
493    }
494    if (size(#)>1)
495    {
496      if (typeof(#[2])=="int" || typeof(#[2])=="number")
497      {
498        redtail = #[2];
499      }
500      if (size(#)>2)
501      {
502        if (typeof(#[3])=="int" || typeof(#[3])=="number")
503        {
504          prepareideal = #[3];
505        }
506      }
507    }
508  }
509  int i,j,k;
510  int sI = ncols(I);
511  // pre-reduce I:
512  module M;
513  if (prepareideal <> 0)
514  {
515    if (remembercoeffs <> 0)
516    {
517      list sortedI = linReduceIdeal(I,1,redtail);
518      I = sortedI[1];
519      M = sortedI[2];
520      dbprint(ppl-1,"// prepared ideal:" +string(I));
521      dbprint(ppl-1,"//  with operations:" +string(M));
522    }
523    else { I = linReduceIdeal(I,0,redtail); }
524  }
525  else
526  {
527    if (remembercoeffs <> 0)
528    {
529      for (i=1; i<=sI; i++) { M[i] = gen(i); }
530    }
531  }
532  ideal lmI,lcI,ordI;
533  for (i=1; i<=sI; i++)
534  {
537    ordI[i] = ord(lmI[i]);
538  }
539  vector v;
540  poly c;
541  // === reduce leading monomials ===
543  int ordf = ord(lm);
544  for (i=sI; i>=1; i--)  // I is sorted: smallest lm's on top
545  {
546    if (lm == 0) { break; }
547    else
548    {
549      if (ordf == ordI[i])
550      {
551        if (lm == lmI[i])
552        {
554          f = f - c*I[i];
556          ordf = ord(lm);
557          if (remembercoeffs <> 0) { v = v - c * M[i]; }
558        }
559      }
560    }
561  }
562  // === reduce tails as well ===
563  if (redtail <> 0)
564  {
565    dbprint(ppl,"// finished reducing leading monomials");
566    dbprint(ppl-1,string(f));
567    if (remembercoeffs <> 0) { dbprint(ppl-1,"// used reductions:" + string(v)); }
568    for (i=1; i<=sI; i++)
569    {
570      dbprint(ppl,"// testing ideal entry "+string(i));
571      for (j=1; j<=size(f); j++)
572      {
573        if (ord(f[j]) == ordI[i])
574        {
575          if (normalize(f[j]) == lmI[i])
576          {
578            f = f - c*I[i];
579            dbprint(ppl-1,"// reducing with " + string(I[i]));
580            dbprint(ppl-1,"//  to " + string(f));
581            if (remembercoeffs <> 0) { v = v - c * M[i]; }
582          }
583        }
584      }
585    }
586  }
587  if (remembercoeffs <> 0)
588  {
589    list l = f,v;
590    return(l);
591  }
592  else { return(f); }
593}
594example
595{
596  "EXAMPLE:"; echo = 2;
597  ring r = 0,(x,y),dp;
598  ideal I = 1,y,xy;
599  poly f = 5xy+7y+3;
600  poly g = 7x+5y+3;
601  linReduce(g,I);     // reduces tails
602  linReduce(g,I,0,0); // no reductions of tails
603  linReduce(f,I,1);   // reduces tails and shows reductions used
604  f = x3+y2+x2+y+x;
605  I = x3-y3, y3-x2,x3-y2,x2-y,y2-x;
606  list l = linReduce(f,I,1);
607  l;
608  module M = I;
609  f - (l[1]-(M*l[2])[1,1]);
610}
611
612proc linSyzSolve (ideal I, list #)
613"USAGE:  linSyzSolve(I[,s]);  I an ideal, s an optional int
614RETURN:  vector, coefficient vector of linear combination of 0 in elements of I
615PURPOSE: compute a linear dependency between the elements of an ideal
616@*       if such one exists
617NOTE:    If s<>0, @code{std} is used for Groebner basis computations,
618@*       otherwise, @code{slimgb} is used.
619@*       By default, @code{slimgb} is used in char 0 and @code{std} in char >0.
620DISPLAY: If printlevel=1, progress debug messages will be printed,
621@*       if printlevel>=2, all the debug messages will be printed.
622EXAMPLE: example linSyzSolve; shows examples
623"
624{
625  int whichengine     = 0; // default
626  int enginespecified = 0; // default
627  if (size(#)>0)
628  {
629    if (typeof(#[1])=="int" || typeof(#[1])=="number")
630    {
631      whichengine = int( #[1]);
632      enginespecified = 1;
633    }
634  }
635  int ppl = printlevel - voice +2;
636  int sI = ncols(I);
637  // check if we are done
638  if (I[sI]==0)
639  {
640    vector v = gen(sI);
641  }
642  else
643  {
644    // ------- 1. introduce undefined coeffs ------------------
645    def save = basering;
646    list RL = ringlist(save);
647    int nv = nvars(save);
648    int np = npars(save);
649    int p = char(save);
650    string cs = "(" + charstr(save) + ")";
651    if (enginespecified == 0)
652    {
653      if (p <> 0)
654      {
655        whichengine = 1;
656      }
657    }
658    int i;
659    list Lvar;
660    for (i=1; i<=sI; i++)
661    {
662      Lvar[i] = safeVarName("@a(" + string(i) + ")");
663    }
664    list L@A = RL[1..4];
665    L@A[2] = Lvar;
666    L@A[3] = list(list("lp",intvec(1:sI)),list("C",intvec(0)));
667    def @A = ring(L@A);
668    L@A[2] = list(safeVarName("@z"));
669    L@A[3][1] = list("dp",intvec(1));
670    if (size(L@A[1])>1)
671    {
672      L@A[1][2] = L@A[1][2] + Lvar;
673      L@A[1][3][size(L@A[1][3])+1] = list("lp",intvec(1:sI));
674    }
675    else
676    {
677      L@A[1] = list(p,Lvar,list(list("lp",intvec(1:sI))),ideal(0));
678    }
679    def @aA = ring(L@A);
680    def @B = save + @aA;
681    setring @B;
682    ideal I = imap(save,I);
683    // ------- 2. form the linear system for the undef coeffs ---
684     poly W;  ideal QQ;
685    for (i=1; i<=sI; i++)
686    {
687      W = W + par(np+i)*I[i];
688    }
689    while (W!=0)
690    {
692      W = W - lead(W);
693    }
694    // QQ consists of polynomial expressions in @a(i) of type number
695    setring @A;
696    ideal QQ = imap(@B,QQ);
697    // ------- 3. this QQ is a polynomial ideal, so "solve" the system -----
698    dbprint(ppl, "// linSyzSolve: starting Groebner basis computation with engine:", whichengine);
699    QQ = engine(QQ,whichengine);
700    dbprint(ppl, "// QQ after engine:", QQ);
701    if (dim(QQ) == -1)
702    {
703      dbprint(ppl+1, "// no solutions by linSyzSolve");
704      // output zeroes
705      setring save;
706      kill @A,@aA,@B;
707      return(v);
708    }
709    // ------- 4. in order to get the numeric values -------
710    matrix AA = matrix(maxideal(1));
711    module MQQ = std(module(QQ));
712    AA = NF(AA,MQQ); // todo: we still receive NF warnings
713    dbprint(ppl, "// AA after NF:",AA);
714    //    "AA after NF:"; print(AA);
715    for(i=1; i<=sI; i++)
716    {
717      AA = subst(AA,var(i),1);
718    }
719    dbprint(ppl, "// AA after subst:",AA);
720    //    "AA after subst: "; print(AA);
721    vector v = (module(transpose(AA)))[1];
722    setring save;
723    vector v = imap(@A,v);
724    kill @A,@aA,@B;
725  }
726  return(v);
727}
728example
729{
730  "EXAMPLE:"; echo = 2;
731  ring r = 0,x,dp;
732  ideal I = x,2x;
733  linSyzSolve(I);
734  ideal J = x,x2;
735  linSyzSolve(J);
736}
737
738proc pIntersect (poly s, ideal I, list #)
739"USAGE:  pIntersect(f, I [,s]);  f a poly, I an ideal, s an optional int
740RETURN:  vector, coefficient vector of the monic polynomial
741PURPOSE: compute the intersection of ideal I with the subalgebra K[f]
742ASSUME:  I is given as Groebner basis, basering is not a qring.
743NOTE:    If the intersection is zero, this proc might not terminate.
744@*       If s>0 is given, it is searched for the generator of the intersection
745@*       only up to degree s. Otherwise (and by default), no bound is assumed.
746DISPLAY: If printlevel=1, progress debug messages will be printed,
747@*       if printlevel>=2, all the debug messages will be printed.
748EXAMPLE: example pIntersect; shows examples
749"
750{
751  def save = basering;
752  int n = nvars(save);
753  list RL = ringlist(save);
754  int i,j,k;
755  if (RL[4]<>0)
756  {
757    ERROR ("basering must not be a qring");
758  }
759  // assume I is given in Groebner basis
760  if (attrib(I,"isSB") <> 1)
761  {
762    print("// WARNING: The input has no SB attribute!");
763    print("// Treating it as if it were a Groebner basis and proceeding...");
764    attrib(I,"isSB",1); // set attribute for suppressing NF messages
765  }
766  int bound = 0; // default
767  if (size(#)>0)
768  {
769    if (typeof(#[1])=="int" || typeof(#[1])=="number")
770    {
771      bound = #[1];
772    }
773  }
774  int ppl = printlevel-voice+2;
775  // ---case 1: I = basering---
776  if (size(I) == 1)
777  {
778    if (simplify(I,2)[1] == 1)
779    {
780      return(gen(1)); // = 1
781    }
782  }
783  // ---case 2: intersection is zero---
785  intvec possdegbounds;
786  list degI;
787  i = 1;
788  while (i <= ncols(I))
789  {
790    if (i == ncols(I)+1) { break; }
792    for (j=1; j<=n; j++)
793    {
794      if (degs[j] == 0)
795      {
796        if (degI[i][j] <> 0) { break; }
797      }
798      if (j == n)
799      {
800        k++;
801        possdegbounds[k] = Max(degI[i]);
802      }
803    }
804    i++;
805  }
806  int degbound = Min(possdegbounds);
807  if (bound>0 && bound<degbound) // given bound is too small
808  {
809    print("// Try a bound of at least " + string(degbound));
810    return(vector(0));
811  }
812  dbprint(ppl,"// lower bound for the degree of the insection is " +string(degbound));
813  if (degbound == 0) // lm(s) does not appear in lm(I)
814  {
815    print("// Intersection is zero");
816    return(vector(0));
817  }
818  // ---case 3: intersection is non-trivial---
819  ideal redNI = 1;
820  vector v;
821  list l,ll;
822  l[1] = vector(0);
823  poly toNF,tobracket,newNF,rednewNF,oldNF,secNF;
824  i = 1;
825  while (1)
826  {
827    if (bound>0 && i>bound) { return(vector(0)); }
828    dbprint(ppl,"// Testing degree "+string(i));
829    if (i>1)
830    {
831      oldNF = newNF;
832      tobracket = s^(i-1) - oldNF;
833      if (tobracket==0) // todo bug in bracket?
834      {
835        toNF = 0;
836      }
837      else
838      {
839        toNF = bracket(tobracket,secNF);
840      }
841      newNF = NF(toNF+oldNF*secNF,I);  // = NF(s^i,I)
842    }
843    else
844    {
845      newNF = NF(s,I);
846      secNF = newNF;
847    }
848    ll = linReduce(newNF,redNI,1);
849    rednewNF = ll[1];
850    l[i+1] = ll[2];
851    dbprint(ppl-1,"// newNF is: "+string(newNF));
852    dbprint(ppl-1,"// rednewNF is: "+string(rednewNF));
853    if (rednewNF != 0) // no linear dependency
854    {
855      redNI[i+1] = rednewNF;
856      i++;
857    }
858    else // there is a linear dependency, hence we are done
859    {
860      dbprint(ppl,"// degree of the generator of the intersection is: "+string(i));
861      break;
862    }
863  }
864  dbprint(ppl-1,"// used linear reductions:", l);
865  // we obtain the coefficients of the generator by the used reductions:
866  list Lvar;
867  for (j=1; j<=i+1; j++)
868  {
869    Lvar[j] = safeVarName("a("+string(j)+")");
870  }
871  list Lord = list("dp",intvec(1:(i+1))),list("C",intvec(0));
872  list L@R = RL[1..4];
873  L@R[2] = Lvar; L@R[3] = Lord;
874  def @R = ring(L@R); setring @R;
875  list l = imap(save,l);
876  ideal C;
877  for (j=1;j<=i+1;j++)
878  {
879    C[j] = 0;
880    for (k=1;k<=j;k++)
881    {
882      C[j] = C[j]+l[j][k]*var(k);
883    }
884  }
885  for (j=i;j>=1;j--)
886  {
887    C[i+1]= subst(C[i+1],var(j),var(j)+C[j]);
888  }
889  matrix m = coeffs(C[i+1],maxideal(1));
890  vector v = gen(i+1);
891  for (j=1;j<=i+1;j++)
892  {
893    v = v + m[j,1]*gen(j);
894  }
895  setring save;
896  v = imap(@R,v);
897  kill @R;
898  return(v);
899}
900example
901{
902  "EXAMPLE:"; echo = 2;
903  ring r = 0,(x,y),dp;
904  poly f = x^2+y^3+x*y^2;
905  def D = initialMalgrange(f);
906  setring D;
907  inF;
908  pIntersect(t*Dt,inF);
909  pIntersect(t*Dt,inF,1);
910}
911
912proc pIntersectSyz (poly s, ideal I, list #)
913"USAGE:  pIntersectSyz(f, I [,p,s,t]);  f poly, I ideal, p,t optial ints, p prime
914RETURN:  vector, coefficient vector of the monic polynomial
915PURPOSE: compute the intersection of an ideal I with the subalgebra K[f]
916ASSUME:  I is given as Groebner basis.
917NOTE:    If the intersection is zero, this procedure might not terminate.
918@*       If p>0 is given, this proc computes the generator of the intersection in
919@*       char p first and then only searches for a generator of the obtained
920@*       degree in the basering. Otherwise, it searches for all degrees by
921@*       computing syzygies.
922@*       If s<>0, @code{std} is used for Groebner basis computations in char 0,
923@*       otherwise, and by default, @code{slimgb} is used.
924@*       If t<>0 and by default, @code{std} is used for Groebner basis
925@*       computations in char >0, otherwise, @code{slimgb} is used.
926DISPLAY: If printlevel=1, progress debug messages will be printed,
927@*       if printlevel>=2, all the debug messages will be printed.
928EXAMPLE: example pIntersectSyz; shows examples
929"
930{
931  // assume I is given as Groebner basis
932  if (attrib(I,"isSB") <> 1)
933  {
934    print("// WARNING: The input has no SB attribute!");
935    print("//  Treating it as if it were a Groebner basis and proceeding...");
936    attrib(I,"isSB",1); // set attribute for suppressing NF messages
937  }
938  int ppl = printlevel-voice+2;
939  int whichengine  = 0; // default
940  int modengine    = 1; // default
941  int solveincharp = 0; // default
942  def save = basering;
943  int n = nvars(save);
944  if (size(#)>0)
945  {
946    if (typeof(#[1])=="int" || typeof(#[1])=="number")
947    {
948      solveincharp = int(#[1]);
949    }
950    if (size(#)>1)
951    {
952      if (typeof(#[2])=="int" || typeof(#[2])=="number")
953      {
954        whichengine = int(#[2]);
955      }
956      if (size(#)>2)
957      {
958        if (typeof(#[3])=="int" || typeof(#[3])=="number")
959        {
960          modengine = int(#[3]);
961        }
962      }
963    }
964  }
965  int i,j;
966  vector v;
967  poly tobracket,toNF,newNF,p;
968  ideal NI = 1;
969  newNF = NF(s,I);
970  NI[2] = newNF;
971  list RL = ringlist(save);
972  if (solveincharp)
973  {
974    int psolveincharp = prime(solveincharp);
975    if (solveincharp <> psolveincharp)
976    {
977      print("// " + string(solveincharp) + " is invalid characteristic of ground field.");
978      print("// " + string(psolveincharp) + " is used.");
979      solveincharp = psolveincharp;
980      kill psolveincharp;
981    }
982    list RLp = RL[1..4];
983    if (typeof(RL[1]) == "int") { RLp[1] = solveincharp; }
984    else { RLp[1][1] = solveincharp; } // parameters
985    def @Rp = ring(RLp);
986    setring @Rp;
987    number c;
988    setring save;
989    int shortSave = short; // workaround for maps Q(a_i) -> Z/p(a_i)
990    short = 0;
991    string str;
993    i=1;
995    {
997      str = "c = " + str + ";";
998      setring @Rp;
999      execute(str);
1000      if (c == 0) { badprime = 1; }
1001      setring save;
1002      i++;
1003    }
1004    str = "poly s = " + string(s) + ";";
1005    if (size(RL) > 4) // basering is NC-algebra
1006    {
1007      string RL5 = "@C = " + string(RL[5]) + ";";
1008      string RL6 = "@D = " + string(RL[6]) + ";";
1009      setring @Rp;
1010      matrix @C[n][n]; matrix @D[n][n];
1011      execute(RL5); execute(RL6);
1012      def Rp = nc_algebra(@C,@D);
1013    }
1014    else { def Rp = @Rp; }
1015    setring Rp;
1016    kill @Rp;
1017    dbprint(ppl-1,"// solving in ring ", Rp);
1018    execute(str);
1019    vector v;
1020    number c;
1021    ideal NI = 1;
1022    setring save;
1023    i=1;
1025    {
1027      str = "c = " + str + ";";
1028      setring Rp;
1029      execute(str);
1030      if (c == 0) { badprime = 1; }
1031      setring save;
1032      i++;
1033    }
1035    {
1037      short = shortSave;
1038      return(vector(0));
1039    }
1040  }
1041  i = 1;
1042  dbprint(ppl,"// pIntersectSyz starts...");
1043  dbprint(ppl-1,"// with ideal I=", I);
1044  while (1)
1045  {
1046    dbprint(ppl,"// testing degree: "+string(i));
1047    if (i>1)
1048    {
1049      tobracket = s^(i-1)-NI[i];
1050      if (tobracket!=0)
1051      {
1052        toNF = bracket(tobracket,NI[2]) + NI[i]*NI[2];
1053      }
1054      else
1055      {
1056        toNF = NI[i]*NI[2];
1057      }
1058      newNF =  NF(toNF,I);
1059      NI[i+1] = newNF;
1060    }
1061    // look for a solution
1062    dbprint(ppl-1,"// linSyzSolve starts with: "+string(matrix(NI)));
1063    if (solveincharp) // modular method
1064    {
1065      for (j=1; j<=size(newNF); j++)
1066      {
1068        str = "c = " + str + ";";
1069        setring Rp;
1070        execute(str);
1071        if (c == 0)
1072        {
1074          setring save;
1075          short = shortSave;
1076          return(vector(0));
1077        }
1078        setring save;
1079      }
1080      str = "NI[" + string(i) + "+1] = " + string(newNF) + ";";
1081      setring Rp;
1082      execute(str); // NI[i+1] = [newNF]_{solveincharp}
1083      v = linSyzSolve(NI,modengine);
1084      if (v!=0) // there is a modular solution
1085      {
1086        dbprint(ppl,"// got solution in char "+string(solveincharp)+" of degree "+string(i));
1087        setring save;
1088        v = linSyzSolve(NI,whichengine);
1089        if (v==0)
1090        {
1091          break;
1092        }
1093      }
1094      else // no modular solution
1095      {
1096        setring save;
1097        v = 0;
1098      }
1099    }
1100    else // non-modular method
1101    {
1102      v = linSyzSolve(NI,whichengine);
1103    }
1104    matrix MM[1][nrows(v)] = matrix(v);
1105    dbprint(ppl-1,"// linSyzSolve ready  with: "+string(MM));
1106    kill MM;
1107    //  "linSyzSolve ready with"; print(v);
1108    if (v!=0)
1109    {
1110      // a solution:
1111      //check for the reality of the solution
1112      p = 0;
1113      for (j=1; j<=i+1; j++)
1114      {
1115        p = p + v[j]*NI[j];
1116      }
1117      if (p!=0)
1118      {
1120      }
1121      else
1122      {
1123        dbprint(ppl,"// linSyzSolve: got solution!");
1124        // "got solution!";
1125        break;
1126      }
1127    }
1128    // no solution:
1129    i++;
1130  }
1131  dbprint(ppl,"// pIntersectSyz finished");
1132  if (solveincharp) { short = shortSave; }
1133  return(v);
1134}
1135example
1136{
1137  "EXAMPLE:"; echo = 2;
1138  ring r = 0,(x,y),dp;
1139  poly f = x^2+y^3+x*y^2;
1140  def D = initialMalgrange(f);
1141  setring D;
1142  inF;
1143  poly s = t*Dt;
1144  pIntersectSyz(s,inF);
1145  int p = prime(20000);
1146  pIntersectSyz(s,inF,p,0,0);
1147}
1148
1149proc vec2poly (list #)
1150"USAGE:  vec2poly(v [,i]);  v a vector or an intvec, i an optional int
1151RETURN:  poly, an univariate polynomial in i-th variable with coefficients given by v
1152PURPOSE: constructs an univariate polynomial in K[var(i)] with given coefficients,
1153@*       such that the coefficient at var(i)^{j-1} is v[j].
1154NOTE:    The optional argument i must be positive, by default i is 1.
1155EXAMPLE: example vec2poly; shows examples
1156"
1157{
1158  def save = basering;
1159  int i,ringvar;
1160  ringvar = 1; // default
1161  if (size(#) > 0)
1162  {
1163    if (typeof(#[1])=="vector" || typeof(#[1])=="intvec")
1164    {
1165      def v = #[1];
1166    }
1167    else
1168    {
1169      ERROR("wrong input: expected vector/intvec expression");
1170    }
1171    if (size(#) > 1)
1172    {
1173      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1174      {
1175        ringvar = int(#[2]);
1176      }
1177    }
1178  }
1179  if (ringvar > nvars(save))
1180  {
1181    ERROR("var out of range");
1182  }
1183  poly p;
1184  for (i=1; i<=nrows(v); i++)
1185  {
1186    p = p + v[i]*(var(ringvar))^(i-1);
1187  }
1188  return(p);
1189}
1190example
1191{
1192  "EXAMPLE:"; echo = 2;
1193  ring r = 0,(x,y),dp;
1194  vector v = gen(1) + 3*gen(3) + 22/9*gen(4);
1195  intvec iv = 3,2,1;
1196  vec2poly(v,2);
1197  vec2poly(iv);
1198}
1199
1200/*
1201// // listofroots and addRoots aren't needed anymore due to some modifications
1202//
1203// static proc listofroots (list #)
1204// {
1205//   def save = basering;
1206//   int n = nvars(save);
1207//   int i;
1208//   poly p;
1209//   if (typeof(#[1])=="vector")
1210//   {
1211//     vector b = #[1];
1212//     for (i=1; i<=nrows(b); i++)
1213//     {
1214//       p = p + b[i]*(var(1))^(i-1);
1215//     }
1216//   }
1217//   else
1218//   {
1219//     p = #[1];
1220//   }
1221//   int substitution = int(#[2]);
1222//   string s = safeVarName("s");
1223//   list RL = ringlist(save); RL = RL[1..4];
1224//   RL[2] = list(s); RL[3] = list(list("dp",intvec(1)),list("C",0));
1225//   def S = ring(RL); setring S;
1226//   ideal J;
1227//   for (i=1; i<=n; i++)
1228//   {
1229//     J[i] = var(1);
1230//   }
1231//   map @m = save,J;
1232//   poly p = @m(p);
1233//   if (substitution == 1)
1234//   {
1235//     p = subst(p,var(1),-var(1)-1);
1236//   }
1237//   // the rest of this proc is nicked from bernsteinBM from dmod.lib
1238//   list P = factorize(p);//with constants and multiplicities
1239//   ideal bs; intvec m;   //the BS polynomial is monic, so we are not interested in constants
1240//   for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
1241//   {
1242//     bs[i-1] = P[1][i];
1243//     m[i-1]  = P[2][i];
1244//   }
1245//   bs =  normalize(bs);
1246//   bs = -subst(bs,var(1),0);
1247//   setring save;
1248//   ideal bs = imap(S,bs);
1249//   kill S;
1250//   list BS = bs,m;
1251//   return(BS);
1252// }
1253//
1254//
1255// static 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*/
1299
1300static proc bfctengine (poly f, int inorann, int whichengine, int addPD, int stdsum, int methodord, int methodpIntersect, int pIntersectchar, int modengine, intvec u0)
1301{
1302  int printlevelsave = printlevel;
1303  printlevel = printlevel + 1;
1304  int ppl = printlevel - voice + 2;
1305  int i;
1306  def save = basering;
1307  int n = nvars(save);
1308  if (isCommutative() == 0) { ERROR("basering must be commutative"); }
1309  if (char(save) <> 0) { ERROR("characteristic of basering has to be 0"); }
1310  list L = ringlist(save);
1311  int qr;
1312  if (L[4] <> 0) // qring
1313  {
1314    print("// basering is qring:");
1315    print("// discarding the quotient and proceeding...");
1316    L[4] = ideal(0);
1317    qr = 1;
1318    def save2 = ring(L); setring save2;
1319    poly f = imap(save,f);
1320  }
1321  if (inorann == 0) // bfct using initial ideal
1322  {
1323    // list L = ringlist(basering);
1324    intvec iv = valvars(f)[1]; // heuristacally better ordering for initialMalgrange
1325    list varL = L[2];
1326    varL = varL[iv];
1327    L[2] = varL;
1328    if (u0 <> 0)
1329    {
1330      u0 = u0[iv];
1331    }
1332    def newr = ring(L);
1333    kill varL,iv,L;
1334    setring newr;
1335    poly f = imap(save,f);
1336    dbprint(ppl,"// starting computation of the initial ideal of the Malgrange ideal...");
1337    def D = initialMalgrange(f,whichengine,methodord,u0);
1338    dbprint(ppl,"// ...done");
1339    setring D;
1340    ideal J = inF;
1341    kill inF;
1342    poly s = var(1)*var(n+2);
1343  }
1344  else // bfct using Ann(f^s)
1345  {
1346    dbprint(ppl,"// starting computation of the s-parametric annihilator...");
1348    dbprint(ppl,"// ...done");
1349    setring D;
1350    ideal J = LD;
1351    kill LD;
1352    poly s = var(1);
1353  }
1354  vector b;
1355  dbprint(ppl,"// starting to intersect with subalgebra...");
1356  // try it modular
1357  if (methodpIntersect <> 0) // pIntersectSyz
1358  {
1359    if (pIntersectchar == 0) // pIntersectSyz::modular
1360    {
1361      list L = ringlist(D);
1362      int lb = 10000; int ub = 536870909; // bounds for random primes
1363      list usedprimes;
1364      int sJ = size(J);
1365      int sLJq;
1366      ideal LJ;
1367      for (i=1; i<=sJ; i++)
1368      {
1370      }
1371      int short_save = short; // workaround for map Q(a_i) -> Z/q(a_i)
1372      short = 0;
1373      string strLJq = "ideal LJq = " + string(LJ) + ";";
1374      int nD = nvars(D);
1375      string L5 = "matrix @C[nD][nD]; @C = " + string(L[5]) + ";";
1376      string L6 = "matrix @D[nD][nD]; @D = " + string(L[6]) + ";";
1377      L = L[1..4];
1378      i = 1;
1379      while (b == 0)
1380      {
1381        dbprint(ppl,"// number of run in the loop: "+string(i));
1382        int q = prime(random(lb,ub));
1383        if (findFirst(usedprimes,q)==0) // if q has not been used already
1384        {
1385          usedprimes = usedprimes,q;
1386          dbprint(ppl,"// using prime: "+string(q));
1387          if (typeof(L[1]) == "int") { L[1] = q; }
1388          else { L[1][1] = q; } // parameters
1389          def @Rq = ring(L); setring @Rq;
1390          execute(L5); execute(L6);
1391          def Rq = nc_algebra(@C,@D); // def Rq = nc_algebra(1,@D);
1392          setring Rq; kill @Rq;
1393          execute(strLJq);
1394          sLJq = size(LJq);
1395          setring D; kill Rq;
1396          if (sLJq <> sJ ) // detect unlucky prime
1397          {
1398            dbprint(ppl,"// " +string(q) + " is unlucky");
1399            b = 0;
1400          }
1401          else
1402          {
1403            b = pIntersectSyz(s,J,q,whichengine,modengine);
1404          }
1405        }
1406        i++;
1407      }
1408      short = short_save;
1409    }
1410    else // pIntersectSyz::non-modular
1411    {
1412      b = pIntersectSyz(s,J,0,whichengine);
1413    }
1414  }
1415  else // pIntersect: linReduce
1416  {
1417    b = pIntersect(s,J);
1418  }
1419  dbprint(ppl,"// ...done"); // with the intersection
1420  poly pb = vec2poly(b);
1421  if (inorann == 0)
1422  {
1423    pb = subst(pb,var(1),-var(1)-1);
1424  }
1425  else // bfctAnn
1426  {
1428    {
1429      pb = pb*(var(1)+1);
1430    }
1431  }
1432  list l = bFactor(pb);
1433  setring save;
1434  list l = imap(D,l);
1435  printlevel = printlevelsave;
1436  return(l);
1437}
1438
1439proc bfct (poly f, list #)
1440"USAGE:  bfct(f [,s,t,v]);  f a poly, s,t optional ints, v an optional intvec
1441RETURN:  list of ideal and intvec
1442PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
1443@*       for the hypersurface defined by f.
1444ASSUME:  The basering is commutative and of characteristic 0.
1445BACKGROUND: In this proc, the initial Malgrange ideal is computed according to
1446@*       the algorithm by Masayuki Noro and then a system of linear equations is
1447@*       solved by linear reductions.
1448NOTE:    In the output list, the ideal contains all the roots
1449@*       and the intvec their multiplicities.
1450@*       If s<>0, @code{std} is used for GB computations,
1451@*       otherwise, and by default, @code{slimgb} is used.
1452@*       If t<>0, a matrix ordering is used for Groebner basis computations,
1453@*       otherwise, and by default, a block ordering is used.
1454@*       If v is a positive weight vector, v is used for homogenization
1455@*       computations, otherwise and by default, no weights are used.
1456DISPLAY: If printlevel=1, progress debug messages will be printed,
1457@*       if printlevel>=2, all the debug messages will be printed.
1458EXAMPLE: example bfct; shows examples
1459"
1460{
1461  int ppl = printlevel - voice +2;
1462  int i;
1463  int n = nvars(basering);
1464  // in # we have two switches:
1465  // one for the engine used for Groebner basis computations,
1466  // one for  M() ordering or its realization
1467  // in # can also be the optional weight vector
1468  int whichengine  = 0; // default
1469  int methodord    = 0; // default
1470  intvec u0        = 0; // default
1471  if (size(#)>0)
1472  {
1473    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1474    {
1475      whichengine = int(#[1]);
1476    }
1477    if (size(#)>1)
1478    {
1479      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1480      {
1481        methodord = int(#[2]);
1482      }
1483      if (size(#)>2)
1484      {
1485        if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1)
1486        {
1487          u0 = #[3];
1488        }
1489      }
1490    }
1491  }
1492  list b = bfctengine(f,0,whichengine,0,0,methodord,0,0,0,u0);
1493  return(b);
1494}
1495example
1496{
1497  "EXAMPLE:"; echo = 2;
1498  ring r = 0,(x,y),dp;
1499  poly f = x^2+y^3+x*y^2;
1500  bfct(f);
1501  intvec v = 3,2;
1502  bfct(f,1,0,v);
1503}
1504
1505proc bfctSyz (poly f, list #)
1506"USAGE:  bfctSyz(f [,r,s,t,u,v]);  f poly, r,s,t,u optional ints, v opt. intvec
1507RETURN:  list of ideal and intvec
1508PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
1509@*       for the hypersurface defined by f
1510ASSUME:  The basering is commutative and of characteristic 0.
1511BACKGROUND: In this proc, the initial Malgrange ideal is computed according to
1512@*       the algorithm by Masayuki Noro and then a system of linear equations is
1513@*       solved by computing syzygies.
1514NOTE:    In the output list, the ideal contains all the roots and the intvec
1515@*       their multiplicities.
1516@*       If r<>0, @code{std} is used for GB computations in characteristic 0,
1517@*       otherwise, and by default, @code{slimgb} is used.
1518@*       If s<>0, a matrix ordering is used for GB computations, otherwise,
1519@*       and by default, a block ordering is used.
1520@*       If t<>0, the computation of the intersection is solely performed over
1521@*       charasteristic 0, otherwise and by default, a modular method is used.
1522@*       If u<>0 and by default, @code{std} is used for GB computations in
1523@*       characteristic >0, otherwise, @code{slimgb} is used.
1524@*       If v is a positive weight vector, v is used for homogenization
1525@*       computations, otherwise and by default, no weights are used.
1526DISPLAY: If printlevel=1, progress debug messages will be printed,
1527@*       if printlevel>=2, all the debug messages will be printed.
1528EXAMPLE: example bfctSyz; shows examples
1529"
1530{
1531  int ppl = printlevel - voice +2;
1532  int i;
1533  // in # we have four switches:
1534  // one for the engine used for Groebner basis computations in char 0,
1535  // one for  M() ordering or its realization
1536  // one for a modular method when computing the intersection
1537  // and one for the engine used for Groebner basis computations in char >0
1538  // in # can also be the optional weight vector
1539  int n = nvars(basering);
1540  int whichengine = 0; // default
1541  int methodord   = 0; // default
1542  int pIntersectchar = 0; // default
1543  int modengine   = 1; // default
1544  intvec u0       = 0; // default
1545  if (size(#)>0)
1546  {
1547    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1548    {
1549      whichengine = int(#[1]);
1550    }
1551    if (size(#)>1)
1552    {
1553      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1554      {
1555        methodord = int(#[2]);
1556      }
1557      if (size(#)>2)
1558      {
1559        if (typeof(#[3])=="int" || typeof(#[3])=="number")
1560        {
1561          pIntersectchar = int(#[3]);
1562        }
1563        if (size(#)>3)
1564        {
1565          if (typeof(#[4])=="int" || typeof(#[4])=="number")
1566          {
1567            modengine = int(#[4]);
1568          }
1569          if (size(#)>4)
1570          {
1571            if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1)
1572            {
1573              u0 = #[5];
1574            }
1575          }
1576        }
1577      }
1578    }
1579  }
1580  list b = bfctengine(f,0,whichengine,0,0,methodord,1,pIntersectchar,modengine,u0);
1581  return(b);
1582}
1583example
1584{
1585  "EXAMPLE:"; echo = 2;
1586  ring r = 0,(x,y),dp;
1587  poly f = x^2+y^3+x*y^2;
1588  bfctSyz(f);
1589  intvec v = 3,2;
1590  bfctSyz(f,0,1,1,0,v);
1591}
1592
1593proc bfctIdeal (ideal I, intvec w, list #)
1594"USAGE:  bfctIdeal(I,w[,s,t]);  I an ideal, w an intvec, s,t optional ints
1595RETURN:  list of ideal and intvec
1596PURPOSE: computes the roots of the global b-function of I w.r.t. the weight (-w,w).
1597ASSUME:  The basering is the n-th Weyl algebra in characteristic 0 and for all
1598@*       1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the
1599@*       sequence of variables is given by x(1),...,x(n),D(1),...,D(n),
1600@*       where D(i) is the differential operator belonging to x(i).
1601@*       Further we assume that I is holonomic.
1602BACKGROUND:  In this proc, the initial ideal of I is computed according to the
1603@*       algorithm by Masayuki Noro and then a system of linear equations is
1604@*       solved by linear reductions.
1605NOTE:    In the output list, say L,
1606@*  - L[1] of type ideal contains all the rational roots of a b-function,
1607@*  - L[2] of type intvec contains the multiplicities of above roots,
1608@*  - optional L[3] of type string is the part of b-function without rational roots.
1609@*  Note, that a b-function of degree 0 is encoded via L[1][1]=0, L[2]=0 and
1610@*  L[3] is 1 (for nonzero constant) or 0 (for zero b-function).
1611@*       If s<>0, @code{std} is used for GB computations in characteristic 0,
1612@*       otherwise, and by default, @code{slimgb} is used.
1613@*       If t<>0, a matrix ordering is used for GB computations, otherwise,
1614@*       and by default, a block ordering is used.
1615DISPLAY: If printlevel=1, progress debug messages will be printed,
1616@*       if printlevel>=2, all the debug messages will be printed.
1617EXAMPLE: example bfctIdeal; shows examples
1618"
1619{
1620  int ppl = printlevel - voice +2;
1621  int i;
1622  def save = basering;
1623  int n = nvars(save)/2;
1624  int whichengine = 0; // default
1625  int methodord   = 0; // default
1626  if (size(#)>0)
1627  {
1628    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1629    {
1630      whichengine = int(#[1]);
1631    }
1632    if (size(#)>1)
1633    {
1634      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1635      {
1636        methodord = int(#[2]);
1637      }
1638    }
1639  }
1640  if (isWeyl()==0) { ERROR("basering is not a Weyl algebra"); }
1641  for (i=1; i<=n; i++)
1642  {
1643    if (bracket(var(i+n),var(i))<>1)
1644    {
1645      ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i)));
1646    }
1647  }
1648  int isH = isHolonomic(I);
1649  if (isH<>1)
1650  {
1651    print("WARNING: given ideal is not holonomic");
1652    print("... setting bound for degree of b-function to 10 and proceeding");
1653    isH = 10;
1654  }
1655  else { isH = 0; } // no degree bound for pIntersect
1656  ideal J = initialIdealW(I,-w,w,whichengine,methodord);
1657  poly s;
1658  for (i=1; i<=n; i++)
1659  {
1660    s = s + w[i]*var(i)*var(n+i);
1661  }
1662  vector b = pIntersect(s,J,isH);
1663  list RL = ringlist(save); RL = RL[1..4];
1664  RL[2] = list(safeVarName("s"));
1665  RL[3] = list(list("dp",intvec(1)),list("C",intvec(0)));
1666  def @S = ring(RL); setring @S;
1667  vector b = imap(save,b);
1668  poly bs = vec2poly(b);
1669  list l = bFactor(bs);
1670  setring save;
1671  list l = imap(@S,l);
1672  return(l);
1673}
1674example
1675{
1676  "EXAMPLE:"; echo = 2;
1677  ring @D = 0,(x,y,Dx,Dy),dp;
1678  def D = Weyl();
1679  setring D;
1680  ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6; I = std(I);
1681  intvec w1 = 0,1;
1682  intvec w2 = 2,3;
1683  bfctIdeal(I,w1);
1684  bfctIdeal(I,w2,0,1);
1685  ideal J = I[size(I)]; // J is not holonomic by construction
1686  bfctIdeal(J,w1); //  b-function of D/J w.r.t. w1 is non-zero
1687  bfctIdeal(J,w2); //  b-function of D/J w.r.t. w2 is zero
1688}
1689
1690proc bfctOneGB (poly f,list #)
1691"USAGE:  bfctOneGB(f [,s,t]);  f a poly, s,t optional ints
1692RETURN:  list of ideal and intvec
1693PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the
1694@*       hypersurface defined by f, using only one GB computation
1695ASSUME:  The basering is commutative and of characteristic 0.
1696BACKGROUND: In this proc, the initial Malgrange ideal is computed based on the
1697@*       algorithm by Masayuki Noro and combined with an elimination ordering.
1698NOTE:    In the output list, the ideal contains all the roots and the intvec
1699@*       their multiplicities.
1700@*       If s<>0, @code{std} is used for the GB computation, otherwise,
1701@*       and by default, @code{slimgb} is used.
1702@*       If t<>0, a matrix ordering is used for GB computations,
1703@*       otherwise, and by default, a block ordering is used.
1704DISPLAY: If printlevel=1, progress debug messages will be printed,
1705@*       if printlevel>=2, all the debug messages will be printed.
1706EXAMPLE: example bfctOneGB; shows examples
1707"
1708{
1709  int ppl = printlevel - voice +2;
1710  if (!isCommutative()) { ERROR("Basering must be commutative"); }
1711  def save = basering;
1712  int n = nvars(save);
1713  if (char(save) <> 0)
1714  {
1715    ERROR("characteristic of basering has to be 0");
1716  }
1717  list L = ringlist(save);
1718  int qr;
1719  if (L[4] <> 0) // qring?
1720  {
1721    print("// basering is qring:");
1722    print("// discarding the quotient and proceeding...");
1723    L[4] = ideal(0);
1724    qr = 1;
1725    def save2 = ring(L);
1726    setring save2;
1727    poly f = imap(save,f);
1728  }
1729  int N = 2*n+4;
1730  int i;
1731  int whichengine = 0; // default
1732  int methodord   = 0; // default
1733  if (size(#)>0)
1734  {
1735    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1736    {
1737      whichengine = int(#[1]);
1738    }
1739    if (size(#)>1)
1740    {
1741      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1742      {
1743        methodord = int(#[2]);
1744      }
1745    }
1746  }
1747  // creating the homogenized extended Weyl algebra
1748  // create names for vars
1749  list Lvar;
1750  Lvar[1]   = safeVarName("t");
1751  Lvar[2]   = safeVarName("s");
1752  Lvar[n+3] = safeVarName("D"+Lvar[1]);
1753  Lvar[N]   = safeVarName("h");
1754  for (i=1; i<=n; i++)
1755  {
1756    Lvar[i+2]   = string(var(i));
1757    Lvar[i+n+3] = safeVarName("D" + string(var(i)));
1758  }
1759  // create ordering
1760  intvec uv = -1; uv[n+3] = 1; uv[N] = 0;
1761  intvec @a = 1:N; @a[2] = 2;
1762  intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0;
1763  list Lord;
1764  Lord[1] = list("a",@a); Lord[2] = list("a",@a2);
1765  if (methodord == 0) // default: block ordering
1766  {
1767    //ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(N-1),lp(1));
1768    Lord[3] = list("a",uv);
1769    Lord[4] = list("dp",intvec(1:(N-1)));
1770    Lord[5] = list("lp",intvec(1));
1771    Lord[6] = list("C",intvec(0));
1772  }
1773  else // M() ordering
1774  {
1775    intmat @Ord[N][N];
1776    @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1);
1777    for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; }
1778    dbprint(ppl,"// weights for ordering: "+string(transpose(@a)));
1779    dbprint(ppl,"// the ordering matrix:",@Ord);
1780    //ring @Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord));
1781    Lord[3] = list("M",intvec(@Ord));
1782    Lord[4] = list("C",intvec(0));
1783  }
1784  // create commutative ring
1785  list L@Dh = ringlist(basering);
1786  L@Dh = L@Dh[1..4]; // if basering is commutative nc_algebra
1787  L@Dh[2] = Lvar; L@Dh[3] = Lord;
1788  def @Dh = ring(L@Dh); setring @Dh;
1789  dbprint(ppl,"// the ring @Dh:",@Dh);
1790  // create non-commutative relations
1791  matrix @relD[N][N];
1792  @relD[1,2] = var(1)*var(N)^2;     //  s*t = t*s + t*h^2
1793  @relD[2,n+3] = var(n+3)*var(N)^2; // Dt*s = s*Dt+Dt*h^2
1794  @relD[1,n+3] = var(N)^2;
1795  for (i=1; i<=n; i++)
1796  {
1797    @relD[i+2,n+3+i] = var(N)^2;
1798  }
1799  dbprint(ppl,"// nc relations:",@relD);
1800  def Dh = nc_algebra(1,@relD);
1801  setring Dh; kill @Dh;
1802  dbprint(ppl,"// computing in ring",Dh);
1803  poly f = imap(save,f);
1804  f = homog(f,h);
1805  // create the Malgrange ideal
1806  ideal I = var(1) - f, var(1)*var(n+3) - var(2);
1807  for (i=1; i<=n; i++)
1808  {
1809    I[3+i] = var(i+n+3)+diff(f,var(i+2))*var(n+3);
1810  }
1811  dbprint(ppl-1, "// the Malgrange ideal: " +string(I));
1812  // the hard part: Groebner basis computation
1813  dbprint(ppl, "// starting Groebner basis computation with engine: "+string(whichengine));
1814  I = engine(I, whichengine);
1815  dbprint(ppl, "// finished Groebner basis computation");
1816  I = subst(I,h,1); // dehomogenization
1817  dbprint(ppl-1,string(I));
1818  // 3.3 the initial form
1819  I = inForm(I,uv);
1820  dbprint(ppl, "// the initial ideal:", string(matrix(I)));
1821  // read off the solution
1822  intvec tonselect = 1;
1823  for (i=3; i<=2*n+4; i++) { tonselect = tonselect,i; }
1824  I = nselect(I,tonselect);
1825  dbprint(ppl, "// generators containing only s:", string(matrix(I)));
1826  I = engine(I, whichengine); // is now a principal ideal;
1827  if (qr == 1) { setring save2; }
1828  else { setring save; }
1829  ideal J; J[2] = var(1);
1830  map @m = Dh,J;
1831  ideal I = @m(I);
1832  poly p = I[1];
1833  p = subst(p,var(1),-var(1)-1);
1834  list l = bFactor(p);
1835  if (qr == 1)
1836  {
1837    setring save;
1838    list l = imap(save2,l);
1839  }
1840  return(l);
1841}
1842example
1843{
1844  "EXAMPLE:"; echo = 2;
1845  ring r = 0,(x,y),dp;
1846  poly f = x^2+y^3+x*y^2;
1847  bfctOneGB(f);
1848  bfctOneGB(f,1,1);
1849}
1850
1851
1852
1853proc bfctAnn (poly f, list #)
1854"USAGE:  bfctAnn(f [,a,b,c]);  f a poly, a, b, c optional ints
1855RETURN:  list of ideal and intvec
1856PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the
1857@*       hypersurface defined by f.
1858ASSUME:  The basering is commutative and of characteristic 0.
1859BACKGROUND: In this proc, Ann(f^s) is computed and then a system of linear
1860@*       equations is solved by linear reductions.
1861NOTE:    In the output list, the ideal contains all the roots and the intvec
1862@*  their multiplicities.
1863@*  If a<>0, only f is appended to Ann(f^s), otherwise, and by default,
1864@*  f and all its partial derivatives are appended.
1865@*  If b<>0, @code{std} is used for GB computations, otherwise, and by
1866@*  default, @code{slimgb} is used.
1867@*  If c<>0, @code{std} is used for Groebner basis computations of ideals
1868@*  <I+J> when I is already a Groebner basis of <I>.
1869@*  Otherwise, and by default the engine determined by the switch b is used.
1870@*  Note that in the case c<>0, the choice for b will be overwritten only
1871@*  for the types of ideals mentioned above.
1872@*  This means that if b<>0, specifying c has no effect.
1873DISPLAY: If printlevel=1, progress debug messages will be printed,
1874@*       if printlevel>=2, all the debug messages will be printed.
1875EXAMPLE: example bfctAnn; shows examples
1876"
1877{
1878  def save = basering;
1879  int ppl = printlevel - voice + 2;
1880  int addPD       = 1; // default
1881  int whichengine = 0; // default
1882  int stdsum      = 0; // default
1883  if (size(#)>0)
1884  {
1885    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1886    {
1888    }
1889    if (size(#)>1)
1890    {
1891      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1892      {
1893        whichengine = int(#[2]);
1894      }
1895      if (size(#)>2)
1896      {
1897        if (typeof(#[3])=="int" || typeof(#[3])=="number")
1898        {
1899          stdsum = int(#[3]);
1900        }
1901      }
1902    }
1903  }
1905  return(b);
1906}
1907example
1908{
1909  "EXAMPLE:"; echo = 2;
1910  ring r = 0,(x,y),dp;
1911  poly f = x^2+y^3+x*y^2;
1912  bfctAnn(f);
1913  def R = reiffen(4,5); setring R;
1914  RC; // the Reiffen curve in 4,5
1915  bfctAnn(RC,0,1);
1916}
1917
1918/*
1919static proc hardexamples ()
1920{
1921  //  some hard examples
1922  ring r1 = 0,(x,y,z,w),dp;
1923  // ab34
1924  poly ab34 = (z3+w4)*(3z2x+4w3y);
1925  bfct(ab34);
1926  // ha3
1927  poly ha3 = xyzw*(x+y)*(x+z)*(x+w)*(y+z+w);
1928  bfct(ha3);
1929  // ha4
1930  poly ha4 = xyzw*(x+y)*(x+z)*(x+w)*(y+z)*(y+w);
1931  bfct(ha4);
1932  // chal4: reiffen(4,5)*reiffen(5,4)
1933  ring r2 = 0,(x,y),dp;
1934  poly chal4 = (x4+xy4+y5)*(x5+x4y+y4);
1935  bfct(chal4);
1936  // (xy+z)*reiffen(4,5)
1937  ring r3 = 0,(x,y,z),dp;
1938  poly xyzreiffen45 = (xy+z)*(y4+yz4+z5);
1939  bfct(xyzreiffen45);
1940
1941// sparse ideal as suggested by Alex; gives 1 as std
1942ideal I1 = 28191*y^2+14628*x*Dy, 24865*x^2+24072*x*Dx+17756*Dy^2;
1943std(I1);
1944}
1945*/
Note: See TracBrowser for help on using the repository browser.