source:git/Singular/LIB/bfct.lib@5291ef

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