source: git/Singular/LIB/bfun.lib @ 858cae

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