source: git/Singular/LIB/bfun.lib @ 3cec5e

spielwiese
Last change on this file since 3cec5e was 3cec5e, checked in by Viktor Levandovskyy <levandov@…>, 15 years ago
*levandov: fixes and doc improvements, bfun instead of bfct_lib git-svn-id: file:///usr/local/Singular/svn/trunk@11383 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 43.5 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id: bfun.lib,v 1.1 2009-02-12 20:25:22 levandov Exp $";
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
9THEORY: Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R,
10@*      one is interested in the global b-Function (also known as Bernstein-Sato
11@*      polynomial) b(s) in K[s], defined to be the monic polynomial, satisfying a functional
12@*      identity L * F^{s+1} = b(s) F^s,   for some operator L in D[s]. Here, D stands for an
13@*      n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>.
14@*      One is interested in the following data:
15@*      - Bernstein-Sato polynomial b(s) in K[s],
16@*      - the list of its roots, which are known to be rational
17@*      - the multiplicities of the roots.
18@*      References: Saito, Sturmfels, Takayama: Groebner Deformations of Hypergeometric
19@*      Differential Equations (2000), Noro: An Efficient Modular Algorithm for Computing
20@*      the Global b-function, (2002).
21
22
23MAIN PROCEDURES:
24
25bfct(f[,s,t,v]);            computes the Bernstein-Sato polynomial of poly f
26bfctSyz(f[,r,s,t,u,v]);     computes the Bernstein-Sato polynomial of poly f
27bfctAnn(f[,s]);             computes the Bernstein-Sato polynomial of poly f
28bfctOneGB(f[,s,t]);         computes the Bernstein-Sato polynomial of poly f
29bfctIdeal(I,w[,s,t]);       computes the global b-function of ideal I w.r.t. the weight w
30pIntersect(f,I);            intersection of the ideal I with the subalgebra K[f] for a poly f
31pIntersectSyz(f,I[,p,s,t]); intersection of the ideal I with the subalgebra K[f] for a poly f
32linReduce(f,I[,s]);         reduces a poly by linear reductions only
33linReduceIdeal(I,[s,t]);    reduces generators of ideal by linear reductions only
34linSyzSolve(I[,s]);         computes a linear dependency of the elements of ideal I
35
36AUXILIARY PROCEDURES:
37
38allPositive(v);  checks whether all entries of an intvec are positive
39scalarProd(v,w); computes the standard scalar product of two intvecs
40vec2poly(v[,i]); constructs an univariate poly with given coefficients
41
42SEE ALSO: dmod_lib, dmodapp_lib, gmssing_lib
43";
44
45
46LIB "qhmoduli.lib"; // for Max
47LIB "dmod.lib";     // for SannfsBFCT etc
48LIB "dmodapp.lib";  // for initialIdealW etc
49
50
51///////////////////////////////////////////////////////////////////////////////
52// testing for consistency of the library:
53proc testbfunlib ()
54{
55  // tests all procs for consistency
56  "AUXILIARY PROCEDURES:";
57  example allPositive;
58  example scalarProd;
59  example vec2poly;
60  "MAIN PROCEDURES:";
61  example bfct;
62  example bfctSyz;
63  example bfctAnn;
64  example bfctOneGB;
65  example bfctIdeal;
66  example pIntersect;
67  example pIntersectSyz;
68  example linReduce;
69  example linReduceIdeal;
70  example linSyzSolve;
71}
72
73//--------------- auxiliary procedures ----------------------------------------
74
75static proc gradedWeyl (intvec u,intvec v)
76"USAGE:  gradedWeyl(u,v); u,v intvecs
77RETURN:  a ring, the associated graded ring of the basering w.r.t. u and v
78PURPOSE: compute the associated graded ring of the basering w.r.t. u and v
79ASSUME: basering is a Weyl algebra
80EXAMPLE: example gradedWeyl; shows examples
81NOTE:    u[i] is the weight of x(i), v[i] the weight of D(i).
82@*       u+v has to be a non-negative intvec.
83"
84{
85  int i;
86  def save = basering;
87  int n = nvars(save)/2;
88  if (nrows(u)<>n || nrows(v)<>n)
89  {
90    ERROR("weight vectors have wrong dimension");
91  }
92  intvec uv,gr;
93  uv = u+v;
94  for (i=1; i<=n; i++)
95  {
96    if (uv[i]>=0)
97    {
98      if (uv[i]==0)
99      {
100        gr[i] = 0;
101      }
102      else
103      {
104        gr[i] = 1;
105      }
106    }
107    else
108    {
109      ERROR("the sum of the weight vectors has to be a non-negative intvec");
110    }
111  }
112  list l = ringlist(save);
113  list l2 = l[2];
114  matrix l6 = l[6];
115  for (i=1; i<=n; i++)
116  {
117    if (gr[i] == 1) 
118    {
119       l2[n+i] = "xi("+string(i)+")";
120       l6[i,n+i] = 0;
121    }
122  }
123  l[2] = l2;
124  l[6] = l6;
125  def G = ring(l);
126  return(G);
127}
128example
129{
130  "EXAMPLE:"; echo = 2;
131  ring @D = 0,(x,y,z,Dx,Dy,Dz),dp;
132  def D = Weyl();
133  setring D;
134  intvec u = -1,-1,1; intvec v = 2,1,1;
135  def G = gradedWeyl(u,v);
136  setring G; G;
137}
138
139
140proc allPositive (intvec v)
141"USAGE:  allPositive(v);  v an intvec
142RETURN:  int, 1 if all components of v are positive, or 0 otherwise
143PURPOSE: check whether all components of an intvec are positive
144EXAMPLE: example allPositive; shows an example
145"
146{
147  int i;
148  for (i=1; i<=size(v); i++)
149  {
150    if (v[i]<=0)
151    {
152      return(0);
153      break;
154    }
155  }
156  return(1);
157}
158example
159{
160  "EXAMPLE:"; echo = 2;
161  intvec v = 1,2,3;
162  allPositive(v);
163  intvec w = 1,-2,3;
164  allPositive(w);
165}
166
167static proc findFirst (list l, i)
168"USAGE:  findFirst(l,i);  l a list, i an argument of any type
169RETURN:  int, the position of the first appearance of i in l, or 0 if i is not a member of l
170PURPOSE: check whether the second argument is a member of a list
171EXAMPLE: example findFirst; shows an example
172"
173{
174  int j;
175  for (j=1; j<=size(l); j++)
176  {
177    if (l[j]==i)
178    {
179      return(j);
180      break;
181    }
182  }
183  return(0);
184}
185//   findFirst(list(1, 2, list(1)),2);       // seems to be a bit buggy,
186//   findFirst(list(1, 2, list(1)),3);       // but works for the purposes
187//   findFirst(list(1, 2, list(1)),list(1)); // of this library
188//   findFirst(l,list(2));
189example
190{
191  "EXAMPLE:"; echo = 2;
192  ring r = 0,x,dp;
193  list l = 1,2,3;
194  findFirst(l,4);
195  findFirst(l,2);
196}
197
198proc scalarProd (intvec v, intvec w)
199"USAGE:  scalarProd(v,w);  v,w intvecs
200RETURN:  int, the standard scalar product of v and w
201PURPOSE: computes the scalar product of two intvecs
202ASSUME:  the arguments are of the same size
203EXAMPLE: example scalarProd; shows examples
204"
205{
206  int i; int sp;
207  if (size(v)!=size(w))
208  {
209    ERROR("non-matching dimensions");
210  }
211  else
212  {
213    for (i=1; i<=size(v);i++)
214    {
215      sp = sp + v[i]*w[i];
216    }
217  }
218  return(sp);
219}
220example
221{
222  "EXAMPLE:"; echo = 2;
223  intvec v = 1,2,3;
224  intvec w = 4,5,6;
225  scalarProd(v,w);
226}
227
228//-------------- main procedures -------------------------------------------------------
229proc linReduceIdeal(ideal I, list #)
230"USAGE:  linReduceIdeal(I [,s,t,u]); I an ideal, s,t,u optional ints
231RETURN:  ideal or list, linear reductum (over field) of I by its elements
232PURPOSE: reduce a list of polys only by linear reductions (no monomial multiplications)
233NOTE:    If s<>0, a list consisting of the reduced ideal and the coefficient
234@*       vectors of the used reductions given as module is returned.
235@*       Otherwise (and by default), only the reduced ideal is returned.
236@*       If t<>0 (and by default) all monomials are reduced (if possible),
237@*       otherwise, only leading monomials are reduced.
238@*       If u<>0 (and by default), the ideal is first sorted in increasing order.
239@*       If u is set to 0 and the given ideal is not sorted in the way described,
240@*       the result might not be as expected.
241DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
242@*       if printlevel>=2, all the debug messages will be printed.
243EXAMPLE: example linReduceIdeal; shows examples
244"
245{
246  // #[1] = remembercoeffs
247  // #[2] = redtail
248  // #[3] = sortideal
249  int ppl = printlevel - voice + 2;
250  int remembercoeffs = 0; // default
251  int redtail        = 1; // default
252  int sortideal      = 1; // default
253  if (size(#)>0)
254  {
255    if (typeof(#[1])=="int" || typeof(#[1])=="number")
256    {
257      remembercoeffs = #[1];
258    }
259    if (size(#)>1)
260    {
261      if (typeof(#[2])=="int" || typeof(#[2])=="number")
262      {
263        redtail = #[2];
264      }
265      if (size(#)>2)
266      {
267        if (typeof(#[3])=="int" || typeof(#[3])=="number")
268        {
269          sortideal = #[3];
270        }
271      }
272    }
273  }
274  int sI = ncols(I);
275  int sZeros = sI - size(I);
276  int i,j;
277  ideal J,lmJ,ordJ;
278  list lJ = sort(I);
279  module M; // for the coefficients
280  // step 1: prepare, e.g. sort I
281  if (sortideal <> 0)
282  {
283    if (sZeros > 0) // I contains zeros
284    {
285      if (remembercoeffs <> 0)
286      {
287        j = 0;
288        for (i=1; i<=sI; i++)
289        {
290          if (I[i] == 0)
291          {
292            j++;
293            J[j] = 0;
294            ordJ[j] = -1;
295            M[j] = gen(i);
296          }
297          else
298          {
299            M[i+sZeros-j] = gen(lJ[2][i-j]+j);
300          }
301        }
302      }
303      else
304      {
305        for (i=1; i<=sZeros; i++)
306        {
307          J[i] = 0;
308          ordJ[i] = -1;
309        }
310      }
311      I = J,lJ[1];
312    }
313    else // I contains no zeros
314    {
315      I = lJ[1];
316      if (remembercoeffs <> 0)
317      {
318        for (i=1; i<=size(lJ[1]); i++) { M[i+sZeros] = gen(lJ[2][i]); }
319      }
320    }
321  }
322  else // assume I is already sorted
323  {
324    if (remembercoeffs <> 0)
325    {
326      for (i=1; i<=ncols(I); i++) { M[i] = gen(i);  }
327    }
328  } 
329  dbprint(ppl-1,"initially sorted ideal:", I);
330  if (remembercoeffs <> 0) { dbprint(ppl-1," used permutations:", M); }
331  // step 2: reduce leading monomials by linear reductions
332  poly lm,c,redpoly,maxlmJ;
333  J[sZeros+1] = I[sZeros+1];
334  lm = I[sZeros+1];
335  maxlmJ = leadmonom(lm);
336  lmJ[sZeros+1] = maxlmJ;
337  int ordlm,reduction;
338  ordJ[sZeros+1] = ord(lm);
339  vector v;
340  int noRedPast;
341  for (i=sZeros+2; i<=sI; i++)
342  {
343    redpoly = I[i];
344    lm = leadmonom(redpoly);
345    ordlm = ord(lm);
346    if (remembercoeffs <> 0) { v = M[i]; }
347    reduction = 1;
348    while (reduction == 1) // while there was a reduction
349    {
350      noRedPast = i;
351      reduction = 0;
352      for (j=sZeros+1; j<noRedPast; j++)
353      {
354        if (lm == 0) { break; } // nothing more to reduce
355        if (lm > maxlmJ) { break; } //lm is bigger than maximal monomial to reduce with
356        if (ordlm == ordJ[j])
357        {
358          if (lm == lmJ[j])
359          {
360            dbprint(ppl-1,"reducing " + string(redpoly));
361            dbprint(ppl-1," with " + string(J[j]));
362            c = leadcoef(redpoly)/leadcoef(J[j]);
363            redpoly = redpoly - c*J[j];
364            dbprint(ppl-1," to " + string(redpoly));
365            lm = leadmonom(redpoly);
366            ordlm = ord(lm);
367            if (remembercoeffs <> 0) { M[i] = M[i] - c * M[j]; }
368            noRedPast = j;
369            reduction = 1;
370          }
371        }
372      }
373    }
374    for (j=sZeros+1; j<i; j++)
375    {
376      if (redpoly < J[j]) { break; }
377    }
378    J = insertGenerator(J,redpoly,j);
379    lmJ = insertGenerator(lmJ,lm,j);
380    ordJ = insertGenerator(ordJ,poly(ordlm),j);
381    if (remembercoeffs <> 0)
382    {
383      v = M[i];
384      M = deleteGenerator(M,i);
385      M = insertGenerator(M,v,j);
386    }
387    maxlmJ = lmJ[i];
388  }
389  // step 3: reduce tails by linear reductions as well
390  if (redtail <> 0)
391  {
392    dbprint(ppl,"finished reducing leading monomials");
393    dbprint(ppl-1,string(J));
394    if (remembercoeffs <> 0) { dbprint(ppl-1,"used reductions:" + string(M)); }
395    int k;
396    for (i=sZeros+1; i<=sI; i++)
397    {
398      lm = lmJ[i];
399      for (j=i+1; j<=sI; j++)
400      {
401        for (k=2; k<=size(J[j]); k++) // run over all terms in J[j]
402        {
403          if (ordJ[i] == ord(J[j][k]))
404          {   
405            if (lm == normalize(J[j][k]))
406            {
407              c = leadcoef(J[j][k])/leadcoef(J[i]);
408              dbprint(ppl-1,"reducing " + string(J[j]));
409              dbprint(ppl-1," with " + string(J[i]));
410              J[j] = J[j] - c*J[i];
411              dbprint(ppl-1," to " + string(J[j]));
412              if (remembercoeffs <> 0) { M[j] = M[j] - c * M[i]; }
413            }
414          }
415        }
416      }
417    }
418  }
419  if (remembercoeffs <> 0) { return(list(J,M)); }
420  else { return(J); }
421}
422example
423{
424  "EXAMPLE:"; echo = 2;
425  ring r = 0,(x,y),dp;
426  ideal I = 3,x+9,y4+5x,2y4+7x+2;
427  linReduceIdeal(I);     // reduces tails
428  linReduceIdeal(I,0,0); // no reductions of tails
429  list l = linReduceIdeal(I,1); // reduces tails and shows reductions used
430  l;
431  module M = I;
432  l[1] - ideal(M*l[2]);
433}
434
435proc linReduce(poly f, ideal I, list #)
436"USAGE:  linReduce(f, I [,s,t,u]); f a poly, I an ideal, s,t,u optional ints
437RETURN:  poly or list, linear reductum (over field) of f by elements from I
438PURPOSE: reduce a poly only by linear reductions (no monomial multiplications)
439NOTE:    If s<>0, a list consisting of the reduced poly and the coefficient
440@*       vector of the used reductions is returned, otherwise (and by default)
441@*       only reduced poly is returned.
442@*       If t<>0 (and by default) all monomials are reduced (if possible),
443@*       otherwise, only leading monomials are reduced.
444@*       If u<>0 (and by default), the ideal is linearly pre-reduced, i.e. instead
445@*       of the given ideal, the output of @code{linReduceIdeal} is used. If u is
446@*       set to 0 and the given ideal does not equal the output of
447@*       @code{linReduceIdeal}, the result might not be as expected.
448DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
449@*       if printlevel>=2, all the debug messages will be printed.
450EXAMPLE: example linReduce; shows examples
451"
452{
453  int ppl = printlevel - voice + 2;
454  int remembercoeffs = 0; // default
455  int redtail        = 1; // default
456  int prepareideal   = 1; // default
457  if (size(#)>0)
458  {
459    if (typeof(#[1])=="int" || typeof(#[1])=="number")
460    {
461      remembercoeffs = #[1];
462    }
463    if (size(#)>1)
464    {
465      if (typeof(#[2])=="int" || typeof(#[2])=="number")
466      {
467        redtail = #[2];
468      }
469      if (size(#)>2)
470      {
471        if (typeof(#[3])=="int" || typeof(#[3])=="number")
472        {
473          prepareideal = #[3];
474        }
475      }
476    }
477  }
478  int i,j,k;
479  int sI = ncols(I);
480  // pre-reduce I:
481  module M;
482  if (prepareideal <> 0)
483  {
484    if (remembercoeffs <> 0)
485    {
486      list sortedI = linReduceIdeal(I,1,redtail);
487      I = sortedI[1];
488      M = sortedI[2];
489      dbprint(ppl-1,"prepared ideal:",I);
490      dbprint(ppl-1," with operations:",M);
491    }
492    else { I = linReduceIdeal(I,0,redtail); }
493  }
494  else
495  {
496    if (remembercoeffs <> 0)
497    {
498      for (i=1; i<=sI; i++) { M[i] = gen(i); }
499    }
500  }
501  ideal lmI,lcI,ordI;
502  for (i=1; i<=sI; i++)
503  {
504    lmI[i] = leadmonom(I[i]);
505    lcI[i] = leadcoef(I[i]);
506    ordI[i] = ord(lmI[i]);
507  }
508  vector v;
509  poly c;
510  // === reduce leading monomials ===
511  poly lm = leadmonom(f);
512  int ordf = ord(lm);
513  for (i=sI; i>=1; i--)  // I is sorted: smallest lm's on top
514  {
515    if (lm == 0) { break; }
516    else
517    {
518      if (ordf == ordI[i])
519      {
520        if (lm == lmI[i])
521        {
522          c = leadcoef(f)/lcI[i];
523          f = f - c*I[i];
524          lm = leadmonom(f);
525          ordf = ord(lm);
526          if (remembercoeffs <> 0) { v = v - c * M[i]; }
527        }
528      }
529    }
530  }
531  // === reduce tails as well ===
532  if (redtail <> 0)
533  {
534    dbprint(ppl,"finished reducing leading monomials");
535    dbprint(ppl-1,string(f));
536    if (remembercoeffs <> 0) { dbprint(ppl-1,"used reductions:" + string(v)); }
537    for (i=1; i<=sI; i++)
538    {
539      dbprint(ppl,"testing ideal entry:",i);
540      for (j=1; j<=size(f); j++)
541      {
542        if (ord(f[j]) == ordI[i])
543        {
544          if (normalize(f[j]) == lmI[i])
545          {
546            c = leadcoef(f[j])/lcI[i];
547            f = f - c*I[i];
548            dbprint(ppl-1,"reducing with " + string(I[i]));
549            dbprint(ppl-1," to " + string(f));
550            if (remembercoeffs <> 0) { v = v - c * M[i]; }
551          }
552        }
553      }
554    }
555  }
556  if (remembercoeffs <> 0)
557  {
558    list l = f,v;
559    return(l);
560  }
561  else { return(f); }
562}
563example
564{
565  "EXAMPLE:"; echo = 2;
566  ring r = 0,(x,y),dp;
567  ideal I = 1,y,xy;
568  poly f = 5xy+7y+3;
569  poly g = 7x+5y+3;
570  linReduce(g,I);     // reduces tails
571  linReduce(g,I,0,0); // no reductions of tails
572  linReduce(f,I,1);   // reduces tails and shows reductions used
573  f = x3+y2+x2+y+x;
574  I = x3-y3, y3-x2,x3-y2,x2-y,y2-x;
575  list l = linReduce(f,I,1);
576  l;
577  module M = I;
578  f - (l[1]-(M*l[2])[1,1]);
579}
580
581proc linSyzSolve (ideal I, list #)
582"USAGE:  linSyzSolve(I[,s]);  I an ideal, s an optional int
583RETURN:  vector, coefficient vector of a linear combination of 0 in the elements of I
584PURPOSE: compute a linear dependency between the elements of an ideal
585@*       if such one exists
586NOTE:    If s<>0, @code{std} is used for Groebner basis computations,
587@*       otherwise, @code{slimgb} is used.
588@*       By default, @code{slimgb} is used in char 0 and @code{std} in char >0.
589DISPLAY: If printlevel=1, progress debug messages will be printed,
590@*       if printlevel>=2, all the debug messages will be printed.
591EXAMPLE: example linSyzSolve; shows examples
592"
593{
594  int whichengine     = 0; // default
595  int enginespecified = 0; // default
596  if (size(#)>0)
597  {
598    if (typeof(#[1])=="int" || typeof(#[1])=="number")
599    {
600      whichengine = int( #[1]);
601      enginespecified = 1;
602    }
603  }
604  int ppl = printlevel - voice +2;
605  int sI = ncols(I);
606  // check if we are done
607  if (I[sI]==0)
608  {
609    vector v = gen(sI);
610  }
611  else
612  {
613    // ------- 1. introduce undefined coeffs ------------------
614    def save = basering;
615    int p = char(save);
616    if (enginespecified == 0)
617    {
618      if (p <> 0)
619      {
620        whichengine = 1;
621      }
622    }
623    ring @A = p,(@a(1..sI)),lp;
624    ring @aA = (p,@a(1..sI)), (@z),dp;
625    // TODO: BUG! WHAT IF THE ORIGINAL RING ALREADY HAS SUCH VARIABLES/PARAMETERS!!!?
626    // TODO: YOU CAN OVERCOME THIS BY MEANS SIMILAR TO "chooseSafeVarName" FROM NEW "matrix.lib"
627    def @B = save + @aA;
628    setring @B;
629    ideal I = imap(save,I);
630    // ------- 2. form the linear system for the undef coeffs ---
631    int i;   poly W;  ideal QQ;
632    for (i=1; i<=sI; i++)
633    {
634      W = W + @a(i)*I[i];
635    }
636    while (W!=0)
637    {
638      QQ = QQ,leadcoef(W);
639      W = W - lead(W);
640    }
641    // QQ consists of polynomial expressions in @a(i) of type number
642    setring @A;
643    ideal QQ = imap(@B,QQ);
644    // ------- 3. this QQ is a polynomial ideal, so "solve" the system -----
645    dbprint(ppl, "linSyzSolve: starting Groebner basis computation with engine:", whichengine);
646    QQ = engine(QQ,whichengine);
647    dbprint(ppl, "QQ after engine:", QQ);
648    if (dim(QQ) == -1)
649    {
650      dbprint(ppl+1, "no solutions by linSyzSolve");
651      // output zeroes
652      setring save;
653      kill @A,@aA,@B;
654      return(v);
655    }
656    // ------- 4. in order to get the numeric values -------
657    matrix AA = matrix(maxideal(1));
658    module MQQ = std(module(QQ));
659    AA = NF(AA,MQQ); // todo: we still receive NF warnings
660    dbprint(ppl, "AA after NF:",AA);
661    //    "AA after NF:"; print(AA);
662    for(i=1; i<=sI; i++)
663    {
664      AA = subst(AA,var(i),1);
665    }
666    dbprint(ppl, "AA after subst:",AA);
667    //    "AA after subst: "; print(AA);
668    vector v = (module(transpose(AA)))[1];
669    setring save;
670    vector v = imap(@A,v);
671    kill @A,@aA,@B;
672  }
673  return(v);
674}
675example
676{
677  "EXAMPLE:"; echo = 2;
678  ring r = 0,x,dp;
679  ideal I = x,2x;
680  linSyzSolve(I);
681  ideal J = x,x2;
682  linSyzSolve(J);
683}
684
685proc pIntersect (poly s, ideal I, list #)
686"USAGE:  pIntersect(f, I [,s]);  f a poly, I an ideal, s an optional int
687RETURN:  vector, coefficient vector of the monic polynomial
688PURPOSE: compute the intersection of ideal I with the subalgebra K[f]
689ASSUME:  I is given as Groebner basis.
690NOTE:    If the intersection is zero, this proc might not terminate.
691@*       If s>0 is given, it is searched for the generator of the intersection
692@*       only up to degree s. Otherwise (and by default), no bound is assumed.
693DISPLAY: If printlevel=1, progress debug messages will be printed,
694@*       if printlevel>=2, all the debug messages will be printed.
695EXAMPLE: example pIntersect; shows examples
696"
697{
698  // assume I is given in Groebner basis
699  if (attrib(I,"isSB") <> 1)
700  {
701    print("WARNING: The input has no SB attribute!");
702    print(" Treating it as if it were a Groebner basis and proceeding...");
703    attrib(I,"isSB",1); // set attribute for suppressing NF messages
704  }
705  int bound = 0; // default
706  if (size(#)>0)
707  {
708    if (typeof(#[1])=="int" || typeof(#[1])=="number")
709    {
710      bound = #[1];
711    }
712  }
713  int ppl = printlevel-voice+1;
714  // ---case 1: I = basering---
715  if (size(I) == 1)
716  {
717    if (simplify(I,2)[1] == 1)
718    {
719      return(gen(2)); // = s
720    }
721  }
722  def save = basering;
723  int n = nvars(save);
724  int i,j,k;
725  // ---case 2: intersection is zero---
726  intvec degs = leadexp(s);
727  intvec possdegbounds;
728  list degI;
729  i = 1;
730  while (i <= ncols(I))
731  {
732    if (i == ncols(I)+1) { break; }
733    degI[i] = leadexp(I[i]);
734    for (j=1; j<=n; j++)
735    {
736      if (degs[j] == 0)
737      {
738        if (degI[i][j] <> 0) { break; }
739      }
740      if (j == n)
741      {
742        k++;
743        possdegbounds[k] = Max(degI[i]);
744      }
745    }
746    i++;
747  }
748  int degbound = Min(possdegbounds);
749  if (bound>0 && bound<degbound) // given bound is too small
750  {
751    print("// Try a bound of at least " + string(degbound));
752    return(vector(0));
753  } 
754  dbprint(ppl,"a lower bound for the degree of the insection is:");
755  dbprint(ppl,degbound);
756  if (degbound == 0) // lm(s) does not appear in lm(I)
757  {
758    print("// Intersection is zero");
759    return(vector(0));
760  }
761  // ---case 3: intersection is non-trivial---
762  ideal redNI = 1;
763  vector v;
764  list l,ll;
765  l[1] = vector(0);
766  poly toNF,tobracket,newNF,rednewNF,oldNF,secNF;
767  i = 1;
768  dbprint(ppl+1,"pIntersect starts...");
769  while (1)
770  {
771    if (bound>0 && i>bound) { return(vector(0)); }
772    dbprint(ppl,"testing degree: "+string(i));
773    if (i>1)
774    {
775      oldNF = newNF;
776      tobracket = s^(i-1) - oldNF;
777      if (tobracket==0) // todo bug in bracket?
778      {
779        toNF = 0;
780      }
781      else
782      {
783        toNF = bracket(tobracket,secNF);
784      }
785      newNF = NF(toNF+oldNF*secNF,I);  // = NF(s^i,I)
786    }
787    else
788    {
789      newNF = NF(s,I);
790      secNF = newNF;
791    }
792    ll = linReduce(newNF,redNI,1);
793    rednewNF = ll[1];
794    l[i+1] = ll[2];
795    dbprint(ppl,"newNF is:", newNF);
796    dbprint(ppl,"rednewNF is:", rednewNF);
797    if (rednewNF != 0) // no linear dependency
798    {
799      redNI[i+1] = rednewNF;
800      i++;
801    }
802    else // there is a linear dependency, hence we are done
803    {
804      dbprint(ppl+1,"the degree of the generator of the intersection is:", i);
805      break;
806    }
807  }
808  dbprint(ppl,"used linear reductions:", l);
809  // we obtain the coefficients of the generator of the intersection by the used reductions:
810  ring @R = 0,(a(1..i+1)),dp;
811  setring @R;
812  list l = imap(save,l);
813  ideal C;
814  for (j=1;j<=i+1;j++)
815  {
816    C[j] = 0;
817    for (k=1;k<=j;k++)
818    {
819      C[j] = C[j]+l[j][k]*a(k);
820    }
821  }
822  for (j=i;j>=1;j--)
823  {
824    C[i+1]= subst(C[i+1],a(j),a(j)+C[j]);
825  }
826  matrix m = coeffs(C[i+1],maxideal(1));
827  vector v = gen(i+1);
828  for (j=1;j<=i+1;j++)
829  {
830    v = v + m[j,1]*gen(j);
831  }
832  setring save;
833  v = imap(@R,v);
834  kill @R;
835  dbprint(ppl+1,"pIntersect finished");
836  return(v);
837}
838example
839{
840  "EXAMPLE:"; echo = 2;
841  ring r = 0,(x,y),dp;
842  poly f = x^2+y^3+x*y^2;
843  def D = initialMalgrange(f);
844  setring D;
845  inF;
846  pIntersect(t*Dt,inF);
847  pIntersect(t*Dt,inF,1);
848}
849
850proc pIntersectSyz (poly s, ideal II, list #)
851"USAGE:  pIntersectSyz(f, I [,p,s,t]);  f a poly, I an ideal, p, t optial ints, p a prime number
852RETURN:  vector, coefficient vector of the monic polynomial
853PURPOSE: compute the intersection of an ideal I with the subalgebra K[f]
854ASSUME:  I is given as Groebner basis.
855NOTE:    If the intersection is zero, this procedure might not terminate.
856@*       If p>0 is given, this proc computes the generator of the intersection in char p first
857@*       and then only searches for a generator of the obtained degree in the basering.
858@*       Otherwise, it searched for all degrees by computing syzygies.
859@*       If s<>0, @code{std} is used for Groebner basis computations in char 0,
860@*       otherwise, and by default, @code{slimgb} is used.
861@*       If t<>0 and by default, @code{std} is used for Groebner basis computations in char >0,
862@*       otherwise, @code{slimgb} is used.
863DISPLAY: If printlevel=1, progress debug messages will be printed,
864@*       if printlevel>=2, all the debug messages will be printed.
865EXAMPLE: example pIntersectSyz; shows examples
866"
867{
868  // assume I is given in Groebner basis
869  ideal I = II;
870  if (attrib(I,"isSB") <> 1)
871  {
872    print("WARNING: The input has no SB attribute!");
873    print(" Treating it as if it were a Groebner basis and proceeding...");
874    attrib(I,"isSB",1); // set attribute for suppressing NF messages
875  }
876  int ppl = printlevel-voice+2;
877  int whichengine  = 0; // default
878  int modengine    = 1; // default
879  int solveincharp = 0; // default
880  def save = basering;
881  if (size(#)>0)
882  {
883    if (typeof(#[1])=="int" || typeof(#[1])=="number")
884    {
885      solveincharp = int(#[1]);
886    }
887    if (size(#)>1)
888    {
889      if (typeof(#[2])=="int" || typeof(#[2])=="number")
890      {
891        whichengine = int(#[2]);
892      }
893      if (size(#)>2)
894      {
895        if (typeof(#[3])=="int" || typeof(#[3])=="number")
896        {
897          modengine = int(#[3]);
898        }
899      }
900    }
901  }
902  int i,j;
903  vector v;
904  poly tobracket,toNF,newNF,p;
905  ideal NI = 1;
906  newNF = NF(s,I);
907  NI[2] = newNF;
908  if (solveincharp<>0)
909  {
910    list l = ringlist(save);
911    l[1] = solveincharp;
912    matrix l5 = l[5];
913    matrix l6 = l[6];
914    def @Rp = ring(l);
915    setring @Rp;
916    list l = ringlist(@Rp);
917    l[5] = fetch(save,l5);
918    l[6] = fetch(save,l6);
919    def Rp = ring(l);
920    setring Rp;
921    kill @Rp;
922    dbprint(ppl+1,"solving in ring ", Rp);
923    vector v;
924    map phi = save,maxideal(1);
925    poly s = phi(s);
926    ideal NI = 1;
927    setring save;
928  }
929  i = 1;
930  dbprint(ppl+1,"pIntersectSyz starts...");
931  dbprint(ppl+1,"with ideal I=", I);
932  while (1)
933  {
934    dbprint(ppl,"i:"+string(i));
935    if (i>1)
936    {
937      tobracket = s^(i-1)-NI[i];
938      if (tobracket!=0)
939      {
940        toNF = bracket(tobracket,NI[2]) + NI[i]*NI[2];
941      }
942      else
943      {
944        toNF = NI[i]*NI[2];
945      }
946      newNF =  NF(toNF,I);
947      NI[i+1] = newNF;
948    }
949    // look for a solution
950    dbprint(ppl,"linSyzSolve starts with: "+string(matrix(NI)));
951    if (solveincharp<>0) // modular method
952    {
953      setring Rp;
954      NI[i+1] = phi(newNF);
955      v = linSyzSolve(NI,modengine);
956      if (v!=0) // there is a modular solution
957      {
958        dbprint(ppl,"got solution in char ",solveincharp," of degree " ,i);
959        setring save;
960        v = linSyzSolve(NI,whichengine);
961        if (v==0)
962        {
963          break;
964        }
965      }
966      else // no modular solution
967      {
968        setring save;
969        v = 0;
970      }
971    }
972    else // non-modular method
973    {
974      v = linSyzSolve(NI,whichengine);
975    }
976    matrix MM[1][nrows(v)] = matrix(v);
977    dbprint(ppl,"linSyzSolve ready  with: "+string(MM));
978    kill MM;
979    //  "linSyzSolve ready with"; print(v);
980    if (v!=0)
981    {
982      // a solution:
983      //check for the reality of the solution
984      p = 0;
985      for (j=1; j<=i+1; j++)
986      {
987        p = p + v[j]*NI[j];
988      }
989      if (p!=0)
990      {
991        dbprint(ppl,"linSyzSolve: bad solution!");
992      }
993      else
994      {
995        dbprint(ppl,"linSyzSolve: got solution!");
996        // "got solution!";
997        break;
998      }
999    }
1000    // no solution:
1001    i++;
1002  }
1003  dbprint(ppl+1,"pIntersectSyz finished");
1004  return(v);
1005}
1006example
1007{
1008  "EXAMPLE:"; echo = 2;
1009  ring r = 0,(x,y),dp;
1010  poly f = x^2+y^3+x*y^2;
1011  def D = initialMalgrange(f);
1012  setring D;
1013  inF;
1014  poly s = t*Dt;
1015  pIntersectSyz(s,inF);
1016  int p = prime(20000);
1017  pIntersectSyz(s,inF,p,0,0);
1018}
1019
1020proc vec2poly (list #)
1021"USAGE:  vec2poly(v [,i]);  v a vector or an intvec, i an optional int
1022RETURN:  poly, an univariate poly in i-th variable with coefficients given by v
1023PURPOSE: constructs an univariate poly in K[var(i)] with given coefficients,
1024@*       such that the coefficient at var(i)^{j-1} is v[j].
1025NOTE:    The optional argument i must be positive, by default i is 1.
1026EXAMPLE: example vec2poly; shows examples
1027"
1028{
1029  def save = basering;
1030  int i,ringvar;
1031  ringvar = 1; // default
1032  if (size(#) > 0)
1033  {
1034    if (typeof(#[1])=="vector" || typeof(#[1])=="intvec")
1035    {
1036      def v = #[1];
1037    }
1038    else
1039    {
1040      ERROR("wrong input: expected vector/intvec expression");
1041    }
1042    if (size(#) > 1)
1043    {
1044      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1045      {
1046        ringvar = int(#[2]);
1047      }
1048    }
1049  }
1050  if (ringvar > nvars(save))
1051  {
1052    ERROR("var out of range");
1053  }
1054  poly p;
1055  for (i=1; i<=nrows(v); i++)
1056  {
1057    p = p + v[i]*(var(ringvar))^(i-1);
1058  }
1059  return(p);
1060}
1061example
1062{
1063  "EXAMPLE:"; echo = 2;
1064  ring r = 0,(x,y),dp;
1065  vector v = gen(1) + 3*gen(3) + 22/9*gen(4);
1066  intvec iv = 3,2,1;
1067  vec2poly(v,2);
1068  vec2poly(iv);
1069}
1070
1071static proc listofroots (list #)
1072{
1073  def save = basering;
1074  int n = nvars(save);
1075  int i;
1076  poly p;
1077  if (typeof(#[1])=="vector")
1078  {
1079    vector b = #[1];
1080    for (i=1; i<=nrows(b); i++)
1081    {
1082      p = p + b[i]*(var(1))^(i-1);
1083    }
1084  }
1085  else
1086  {
1087    p = #[1];
1088  }
1089  int substitution = int(#[2]);
1090  ring S = 0,s,dp;
1091  ideal J;
1092  for (i=1; i<=n; i++)
1093  {
1094    J[i] = s;
1095  }
1096  map @m = save,J;
1097  poly p = @m(p);
1098  if (substitution == 1)
1099  {
1100    p = subst(p,s,-s-1);
1101  }
1102  // the rest of this proc is nicked from bernsteinBM from dmod.lib
1103  list P = factorize(p);//with constants and multiplicities
1104  ideal bs; intvec m;   //the Bernstein polynomial is monic, so we are not interested in constants
1105  for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
1106  {
1107    bs[i-1] = P[1][i];
1108    m[i-1]  = P[2][i];
1109  }
1110  bs =  normalize(bs);
1111  bs = -subst(bs,s,0);
1112  setring save;
1113  ideal bs = imap(S,bs);
1114  kill S;
1115  list BS = bs,m;
1116  return(BS);
1117}
1118
1119static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodpIntersect, int pIntersectchar, int modengine, intvec u0)
1120{
1121  int ppl = printlevel - voice +2;
1122  int i;
1123  def save = basering;
1124  int n = nvars(save);
1125  if (size(variables(f)) == 0) // f is constant
1126  {
1127    return(list(ideal(0),intvec(0)));
1128  } 
1129  if (inorann == 0) // bfct using initial ideal
1130  {
1131    def D = initialMalgrange(f,whichengine,methodord,1,u0);
1132    setring D;
1133    ideal J = inF;
1134    kill inF;
1135    poly s = t*Dt;
1136  }
1137  else // bfct using Ann(f^s)
1138  {
1139    def D = SannfsBFCT(f,whichengine);
1140    setring D;
1141    ideal J = LD;
1142    kill LD;
1143  }
1144  vector b;
1145  // try it modular
1146  if (methodpIntersect <> 0) // pIntersectSyz
1147  {
1148    if (pIntersectchar == 0) // pIntersectSyz::modular
1149    {
1150      int lb = 30000;
1151      int ub = 536870909;
1152      i = 1;
1153      list usedprimes;
1154      while (b == 0)
1155      {
1156        dbprint(ppl,"number of run in the loop: "+string(i));
1157        int q = prime(random(lb,ub));
1158        if (findFirst(usedprimes,q)==0) // if q was not already used
1159        {
1160          usedprimes = usedprimes,q;
1161          dbprint(ppl,"used prime is: "+string(q));
1162          b = pIntersectSyz(s,J,q,whichengine,modengine);
1163        }
1164        i++;
1165      }
1166    }
1167    else // pIntersectSyz::non-modular
1168    {
1169      b = pIntersectSyz(s,J,0,whichengine);
1170    }
1171  }
1172  else // pIntersect: linReduce
1173  {
1174    b = pIntersect(s,J);
1175  }
1176  setring save;
1177  vector b = imap(D,b);
1178  if (inorann == 0)
1179  {
1180    list l = listofroots(b,1);
1181  }
1182  else
1183  {
1184    list l = listofroots(b,0);
1185  }
1186  return(l);
1187}
1188
1189proc bfct (poly f, list #)
1190"USAGE:  bfct(f [,s,t,v]);  f a poly, s,t optional ints, v an optional intvec
1191RETURN:  list of ideal and intvec
1192PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
1193@*       for the hypersurface defined by f.
1194ASSUME:  The basering is a commutative polynomial ring in char 0.
1195BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm
1196@*       by Masayuki Noro and then a system of linear equations is solved by linear reductions.
1197NOTE:    In the output list, the ideal contains all the roots
1198@*       and the intvec their multiplicities.
1199@*       If s<>0, @code{std} is used for GB computations,
1200@*       otherwise, and by default, @code{slimgb} is used.
1201@*       If t<>0, a matrix ordering is used for Groebner basis computations,
1202@*       otherwise, and by default, a block ordering is used.
1203@*       If v is a positive weight vector, v is used for homogenization computations,
1204@*       otherwise and by default, no weights are used.
1205DISPLAY: If printlevel=1, progress debug messages will be printed,
1206@*       if printlevel>=2, all the debug messages will be printed.
1207EXAMPLE: example bfct; shows examples
1208"
1209{
1210  int ppl = printlevel - voice +2;
1211  int i;
1212  int n = nvars(basering);
1213  // in # we have two switches:
1214  // one for the engine used for Groebner basis computations,
1215  // one for  M() ordering or its realization
1216  // in # can also be the optional weight vector
1217  int whichengine  = 0; // default
1218  int methodord    = 0; // default
1219  intvec u0        = 0; // default
1220  if (size(#)>0)
1221  {
1222    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1223    {
1224      whichengine = int(#[1]);
1225    }
1226    if (size(#)>1)
1227    {
1228      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1229      {
1230        methodord = int(#[2]);
1231      }
1232      if (size(#)>2)
1233      {
1234        if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1)
1235        {
1236          u0 = #[3];
1237        }
1238      }
1239    }
1240  }
1241  list b = bfctengine(f,0,whichengine,methodord,0,0,0,u0);
1242  return(b);
1243}
1244example
1245{
1246  "EXAMPLE:"; echo = 2;
1247  ring r = 0,(x,y),dp;
1248  poly f = x^2+y^3+x*y^2;
1249  bfct(f);
1250  intvec v = 3,2;
1251  bfct(f,1,0,v);
1252}
1253
1254proc bfctSyz (poly f, list #)
1255"USAGE:  bfctSyz(f [,r,s,t,u,v]);  f a poly, r,s,t,u optional ints, v an optional intvec
1256RETURN:  list of ideal and intvec
1257PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
1258@*       for the hypersurface defined by f
1259ASSUME:  The basering is a commutative polynomial ring in char 0.
1260BACKGROUND: In this proc, the initial Malgrange ideal is computed according to the algorithm
1261@*       by Masayuki Noro and then a system of linear equations is solved by computing syzygies.
1262NOTE:    In the output list, the ideal contains all the roots
1263@*       and the intvec their multiplicities.
1264@*       If r<>0, @code{std} is used for GB computations in characteristic 0,
1265@*       otherwise, and by default, @code{slimgb} is used.
1266@*       If s<>0, a matrix ordering is used for GB computations, otherwise,
1267@*       and by default, a block ordering is used.
1268@*       If t<>0, the computation of the intersection is solely performed over
1269@*       charasteristic 0, otherwise and by default, a modular method is used.
1270@*       If u<>0 and by default, @code{std} is used for GB computations in
1271@*       characteristic >0, otherwise, @code{slimgb} is used.
1272@*       If v is a positive weight vector, v is used for homogenization
1273@*       computations, otherwise and by default, no weights are used.
1274DISPLAY: If printlevel=1, progress debug messages will be printed,
1275@*       if printlevel>=2, all the debug messages will be printed.
1276EXAMPLE: example bfctSyz; shows examples
1277"
1278{
1279  int ppl = printlevel - voice +2;
1280  int i;
1281  // in # we have four switches:
1282  // one for the engine used for Groebner basis computations in char 0,
1283  // one for  M() ordering or its realization
1284  // one for a modular method when computing the intersection
1285  // and one for the engine used for Groebner basis computations in char >0
1286  // in # can also be the optional weight vector
1287  def save = basering;
1288  int n = nvars(save);
1289  int whichengine = 0; // default
1290  int methodord   = 0; // default
1291  int pIntersectchar  = 0; // default
1292  int modengine   = 1; // default
1293  intvec u0       = 0; // default
1294  if (size(#)>0)
1295  {
1296    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1297    {
1298      whichengine = int(#[1]);
1299    }
1300    if (size(#)>1)
1301    {
1302      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1303      {
1304        methodord = int(#[2]);
1305      }
1306      if (size(#)>2)
1307      {
1308        if (typeof(#[3])=="int" || typeof(#[3])=="number")
1309        {
1310          pIntersectchar = int(#[3]);
1311        }
1312        if (size(#)>3)
1313        {
1314          if (typeof(#[4])=="int" || typeof(#[4])=="number")
1315          {
1316            modengine = int(#[4]);
1317          }
1318          if (size(#)>4)
1319          {
1320            if (typeof(#[5])=="intvec" && size(#[5])==n && allPositive(#[5])==1)
1321            {
1322              u0 = #[5];
1323            }
1324          }
1325        }
1326      }
1327    }
1328  }
1329  list b = bfctengine(f,0,whichengine,methodord,1,pIntersectchar,modengine,u0);
1330  return(b);
1331}
1332example
1333{
1334  "EXAMPLE:"; echo = 2;
1335  ring r = 0,(x,y),dp;
1336  poly f = x^2+y^3+x*y^2;
1337  bfctSyz(f);
1338  intvec v = 3,2;
1339  bfctSyz(f,0,1,1,0,v);
1340}
1341
1342proc bfctIdeal (ideal I, intvec w, list #)
1343"USAGE:  bfctIdeal(I,w[,s,t]);  I an ideal, w an intvec, s,t optional ints
1344RETURN:  list of ideal and intvec
1345PURPOSE: computes the roots of the global b-function of I wrt the weight vector (-w,w).
1346ASSUME:  The basering is the n-th Weyl algebra in char 0, where the sequence of
1347@*       the variables is x(1),...,x(n),D(1),...,D(n).
1348BACKGROUND:  In this proc, the initial ideal of I is computed according to the algorithm by
1349@*       Masayuki Noro and then a system of linear equations is solved by linear reductions.
1350NOTE:    In the output list, the ideal contains all the roots
1351@*       and the intvec their multiplicities.
1352@*       If s<>0, @code{std} is used for GB computations in characteristic 0,
1353@*       otherwise, and by default, @code{slimgb} is used.
1354@*       If t<>0, a matrix ordering is used for GB computations, otherwise,
1355@*       and by default, a block ordering is used.
1356DISPLAY: If printlevel=1, progress debug messages will be printed,
1357@*       if printlevel>=2, all the debug messages will be printed.
1358EXAMPLE: example bfctIdeal; shows examples
1359"
1360{
1361  int ppl = printlevel - voice +2;
1362  int i;
1363  def save = basering;
1364  int n = nvars(save)/2;
1365  int whichengine = 0; // default
1366  int methodord   = 0; // default
1367  if (size(#)>0)
1368  {
1369    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1370    {
1371      whichengine = int(#[1]);
1372    }
1373    if (size(#)>1)
1374    {
1375      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1376      {
1377        methodord = int(#[2]);
1378      }
1379    }
1380  }
1381  ideal J = initialIdealW(I,-w,w,whichengine,methodord);
1382  poly s;
1383  for (i=1; i<=n; i++)
1384  {
1385    s = s + w[i]*var(i)*var(n+i);
1386  }
1387  vector b = pIntersect(s,J);
1388  list l = listofroots(b,0);
1389  return(l);
1390}
1391example
1392{
1393  "EXAMPLE:"; echo = 2;
1394  ring @D = 0,(x,y,Dx,Dy),dp;
1395  def D = Weyl();
1396  setring D;
1397  ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6; I = std(I);
1398  intvec w1 = 1,1;
1399  intvec w2 = 1,2;
1400  intvec w3 = 2,3;
1401  bfctIdeal(I,w1);
1402  bfctIdeal(I,w2,1);
1403  bfctIdeal(I,w3,0,1);
1404}
1405
1406proc bfctOneGB (poly f,list #)
1407"USAGE:  bfctOneGB(f [,s,t]);  f a poly, s,t optional ints
1408RETURN:  list of ideal and intvec
1409PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s) for the
1410@*       hypersurface defined by f, using only one GB computation
1411ASSUME:  The basering is a commutative polynomial ring in char 0.
1412BACKGROUND: In this proc, the initial Malgrange ideal is computed based on the
1413@*       algorithm by Masayuki Noro and combined with an elimination ordering.
1414NOTE:    In the output list, the ideal contains all the roots
1415@*       and the intvec their multiplicities.
1416@*       If s<>0, @code{std} is used for the GB computation, otherwise,
1417@*       and by default, @code{slimgb} is used.
1418@*       If t<>0, a matrix ordering is used for GB computations,
1419@*       otherwise, and by default, a block ordering is used.
1420DISPLAY: If printlevel=1, progress debug messages will be printed,
1421@*       if printlevel>=2, all the debug messages will be printed.
1422EXAMPLE: example bfctOneGB; shows examples
1423"
1424{
1425  int ppl = printlevel - voice +2;
1426  def save = basering;
1427  int n = nvars(save);
1428  int noofvars = 2*n+4;
1429  int i;
1430  int whichengine = 0; // default
1431  int methodord   = 0; // default
1432  if (size(#)>0)
1433  {
1434    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1435    {
1436      whichengine = int(#[1]);
1437    }
1438    if (size(#)>1)
1439    {
1440      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1441      {
1442        methodord = int(#[2]);
1443      }
1444    }
1445  }
1446  intvec uv;
1447  uv[n+3] = 1;
1448  ring r = 0,(x(1..n)),dp;
1449  poly f = fetch(save,f);
1450  uv[1] = -1; uv[noofvars] = 0;
1451  //  for the ordering
1452  intvec @a; @a = 1:noofvars; @a[2] = 2;
1453  intvec @a2 = @a; @a2[2] = 0; @a2[2*n+4] = 0;
1454  if (methodord == 0) // default: block ordering
1455  {
1456    ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),a(uv),dp(noofvars-1),lp(1));
1457  }
1458  else // M() ordering
1459  {
1460    intmat @Ord[noofvars][noofvars];
1461    @Ord[1,1..noofvars] = uv;
1462    @Ord[2,1..noofvars] = 1:(noofvars-1);
1463    for (i=1; i<=noofvars-2; i++)
1464    {
1465      @Ord[2+i,noofvars - i] = -1;
1466    }
1467    dbprint(ppl,"weights for ordering:",transpose(@a));
1468    dbprint(ppl,"the ordering matrix:",@Ord);
1469    ring Dh = 0,(t,s,x(n..1),Dt,D(n..1),h),(a(@a),a(@a2),M(@Ord));
1470  }
1471  dbprint(ppl,"the ring Dh:",Dh);
1472  // non-commutative relations
1473  matrix @relD[noofvars][noofvars];
1474  @relD[1,2] = t*h^2;// s*t = t*s+t*h^2
1475  @relD[2,n+3] = Dt*h^2;// Dt*s = s*Dt+h^2
1476  @relD[1,n+3] = h^2;
1477  for (i=1; i<=n; i++)
1478  {
1479    @relD[i+2,n+3+i] = h^2;
1480  }
1481  dbprint(ppl,"nc relations:",@relD);
1482  def DDh = nc_algebra(1,@relD);
1483  setring DDh;
1484  dbprint(ppl,"computing in ring",DDh);
1485  ideal I;
1486  poly f = imap(r,f);
1487  kill r;
1488  f = homog(f,h);
1489  I = t - f, t*Dt - s;  // defining the Malgrange ideal
1490  for (i=1; i<=n; i++)
1491  {
1492    I = I, D(i)+diff(f,x(i))*Dt;
1493  }
1494  dbprint(ppl, "starting Groebner basis computation with engine:", whichengine);
1495  I = engine(I, whichengine);
1496  dbprint(ppl, "finished Groebner basis computation");
1497  dbprint(ppl, "I before dehomogenization is" ,I);
1498  I = subst(I,h,1); // dehomogenization
1499  dbprint(ppl, "I after dehomogenization is" ,I);
1500  I = inForm(I,uv); // we are only interested in the initial form w.r.t. uv
1501  dbprint(ppl, "the initial ideal:", string(matrix(I)));
1502  intvec tonselect = 1;
1503  for (i=3; i<=2*n+4; i++) { tonselect = tonselect,i; }
1504  I = nselect(I,tonselect);
1505  dbprint(ppl, "generators containing only s:", string(matrix(I)));
1506  I = engine(I, whichengine); // is now a principal ideal;
1507  setring save;
1508  ideal J; J[2] = var(1);
1509  map @m = DDh,J;
1510  ideal I = @m(I);
1511  poly p = I[1];
1512  list l = listofroots(p,1);
1513  return(l);
1514}
1515example
1516{
1517  "EXAMPLE:"; echo = 2;
1518  ring r = 0,(x,y),dp;
1519  poly f = x^2+y^3+x*y^2;
1520  bfctOneGB(f);
1521  bfctOneGB(f,1,1);
1522}
1523
1524proc bfctAnn (poly f, list #)
1525"USAGE:  bfctAnn(f [,r]);  f a poly, r an optional int
1526RETURN:  list of ideal and intvec
1527PURPOSE: computes the roots of the Bernstein-Sato polynomial b(s)
1528@*       for the hypersurface defined by f
1529ASSUME:  The basering is a commutative polynomial ring in char 0.
1530BACKGROUND: In this proc, ann(f^s) is computed and then a system of linear
1531@*       equations is solved by linear reductions.
1532NOTE:    In the output list, the ideal contains all the roots
1533@*       and the intvec their multiplicities.
1534@*       If r<>0, @code{std} is used for GB computations,
1535@*       otherwise, and by default, @code{slimgb} is used.
1536DISPLAY: If printlevel=1, progress debug messages will be printed,
1537@*       if printlevel>=2, all the debug messages will be printed.
1538EXAMPLE: example bfctAnn; shows examples
1539"
1540{
1541  def save = basering;
1542  int ppl = printlevel - voice + 2;
1543  int whichengine = 0; // default
1544  if (size(#)>0)
1545  {
1546    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1547    {
1548      whichengine = int(#[1]);
1549    }
1550  }
1551  list b = bfctengine(f,1,whichengine,0,1,0,0,0);
1552  return(b);
1553}
1554example
1555{
1556  "EXAMPLE:"; echo = 2;
1557  ring r = 0,(x,y),dp;
1558  poly f = x^2+y^3+x*y^2;
1559  bfctAnn(f);
1560}
1561
1562/*
1563//static proc hardexamples ()
1564{
1565  //  some hard examples
1566  ring r1 = 0,(x,y,z,w),dp;
1567  // ab34
1568  poly ab34 = (z3+w4)*(3z2x+4w3y);
1569  bfct(ab34);
1570  // ha3
1571  poly ha3 = xyzw*(x+y)*(x+z)*(x+w)*(y+z+w);
1572  bfct(ha3);
1573  // ha4
1574  poly ha4 = xyzw*(x+y)*(x+z)*(x+w)*(y+z)*(y+w);
1575  bfct(ha4);
1576  // chal4: reiffen(4,5)*reiffen(5,4)
1577  ring r2 = 0,(x,y),dp;
1578  poly chal4 = (x4+xy4+y5)*(x5+x4y+y4);
1579  bfct(chal4);
1580  // (xy+z)*reiffen(4,5)
1581  ring r3 = 0,(x,y,z),dp;
1582  poly xyzreiffen45 = (xy+z)*(y4+yz4+z5);
1583  bfct(xyzreiffen45);
1584}
1585*/
Note: See TracBrowser for help on using the repository browser.