# source:git/Singular/LIB/bfun.lib@3360fb

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