source: git/Singular/LIB/bfct.lib @ f7a73ed

spielwiese
Last change on this file since f7a73ed was f7a73ed, checked in by Viktor Levandovskyy <levandov@…>, 16 years ago
*levandov: b-function via Noro's algorithm, both char 0 and prime git-svn-id: file:///usr/local/Singular/svn/trunk@10711 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 17.5 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id: bfct.lib,v 1.1 2008-05-09 15:47:02 levandov Exp $";
3category="Noncommutative";
4info="
5LIBRARY: bfct.lib     M. Noro's Algorithm for Bernstein-Sato polynomial
6AUTHORS: Daniel Andres, wolf.daniel.andres@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 Bernstein-Sato polynomial b(s) in K[s],
11@*      defined to be the monic polynomial, satisfying a functional identity
12@*             L * f^{s+1} = b(s) f^s,   for some operator L in D[s].
13@* Here, D stands for an n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1>
14@* One is interested in the following data:
15@*   global Bernstein polynomial in K[s] and
16@*   the list of all roots of b(s), which are known to be rational, with their multiplicities.
17
18MAIN PROCEDURES:
19
20bfct();
21
22AUXILIARY PROCEDURES:
23
24";
25
26
27LIB "qhmoduli.lib"; // for Max
28LIB "nctools.lib";
29
30//--------------- help procedures ---------------------------------------------------------
31
32proc polylist (poly f)
33"USAGE:  polylist(f);  f a poly
34RETURN:  list of exponents and corresponding terms of f
35PURPOSE: convert a polynomial to a list of exponents and corresponding term
36EXAMPLE: example polylist; shows examples
37"
38{
39  list l;
40  int i = 1;
41  if (f==0) // just for the zero polynomial
42  {
43    l[1] = list(leadexp(f), lead(f));
44  }
45  while (f!=0)
46  {
47    l[i] = list(leadexp(f), lead(f));
48    f = f-lead(f);
49    i++;
50  }
51  return(l);
52}
53example
54{
55  "EXAMPLE:"; echo = 2;
56  ring r = 0,(x),Dp;
57  poly F = x;
58  polylist(F);
59  ring r = 0,(x,y),Dp;
60  poly F = x2y+x*y2;
61  polylist(F);
62}
63
64proc ispositive (intvec v)
65"USAGE:  ispositive(v);  v an intvec
66RETURN:  1 if all components of v are bigger than 0, or 0 otherwise
67PURPOSE: check whether all components of an intvec are positive
68EXAMPLE: example ispositive; shows an example
69"
70{
71  int i;
72  for (i=1; i<=size(v); i++)
73  {
74    if (v[i]<=0)
75    {
76      return(0);
77      break;
78    }
79  }
80  return(1);
81}
82example
83{
84  "EXAMPLE:"; echo = 2;
85  intvec v = 1,2,3;
86  ispositive(v);
87  intvec w = 1,-2,3;
88  ispositive(w);
89}
90
91proc isin (list l, i)
92"USAGE:  isin(l,i);  l a list, i an argument of any type
93RETURN:  1 if i is a member of l, or 0 otherwise
94PURPOSE: check whether the second argument is a member of a list
95EXAMPLE: example isin; shows an example
96"
97{
98  int j;
99  for (j=1; j<=size(l); j++)
100  {
101    if (l[j]==i)
102    {
103      return(1);
104      break;
105    }
106  }
107  return(0);
108}
109example
110{
111  "EXAMPLE:"; echo = 2;
112  list l = 1,2,3;
113  isin(l,4);
114  isin(l,2);
115}
116
117proc scalarprod (intvec v, intvec w)
118"USAGE:  scalarprod(v,w);  v,w intvecs
119RETURN:  an int, the scalarprod of v and w
120PURPOSE: compute the scalarprod of two intvecs
121EXAMPLE: example scalarprod; shows examples
122NOTE:    the arguments must have the same size
123"
124{
125  int i; int sp;
126  if (size(v)!=size(w))
127    {
128      ERROR("non-matching dimensions");
129    }
130  else
131  {
132    for (i=1; i<=size(v);i++)
133    {
134      sp = sp + v[i]*w[i];
135    }
136  }
137  return(sp);
138}
139example
140{
141  "EXAMPLE:"; echo = 2;
142  intvec v = 1,2,3;
143  intvec w = 4,5,6;
144  scalarprod(v,w);
145}
146
147proc InForm (list #)
148"USAGE:  InForm(f,w) or InForm(I,w);  f a poly, I an ideal, w an intvec
149RETURN:  the initial of f or I w.r.t. to the weight w
150PURPOSE: compute the initial form of a poly or ideal w.r.t a given weight
151NOTE:    the size of the weight vector must be equal to the number of variables of the basering
152EXAMPLE: example InForm; shows examples
153"
154{
155  if (size(#)<2)
156  {
157    ERROR("InForm needs two arguments of type  poly/ideal,intvec");
158  }
159  if (typeof(#[1])=="poly" || typeof(#[1])=="ideal")
160  {
161    ideal I = #[1];
162  }
163  else
164  {
165    ERROR("first argument must be of type poly or ideal");
166  }
167  if (typeof(#[2])=="intvec")
168  {
169    def w = #[2];
170  }
171  else
172  {
173    ERROR("second argument must be of type intvec");
174  }
175  if (size(w) != nvars(basering))
176  {
177    ERROR("weight vector has wrong dimension");
178  }
179  int j,i,s,m;
180  list l;
181  poly g;
182  ideal J;
183  for (j=1; j<=ncols(I); j++)
184  {
185    l = polylist(I[j]);
186    m = scalarprod(w,l[1][1]);
187    g = 0;
188    for (i=2; i<=size(l); i++)
189    {
190      s = scalarprod(w,l[i][1]);
191      m = Max(list(s,m));
192    }
193    for (i=1; i<=size(l); i++)
194    {
195      if (scalarprod(w,l[i][1]) == m)
196      {
197        g = g+l[i][2];
198      }
199    }
200    J[j] = g;
201   
202  }
203  if (typeof(#[1])=="poly")
204  {
205    return(J[1]);  // if the input was a poly, return a poly
206  }
207  else
208  {
209    return(J);     // otherwise, return an ideal
210  }
211}
212example
213{
214  "EXAMPLE:"; echo = 2;
215  // todo: create example
216}
217
218//-------------- main procedures -------------------------------------------------------
219
220proc ncsolve (ideal I, list #)
221"USAGE:  ncsolve(I[,s]);  I an ideal, s an optional int
222RETURN:  coefficient vector of a linear combination of 0 in the elements of I
223PURPOSE: compute a linear dependency between the elements of an ideal if such one exists
224NOTE:    If s<>0, @code{std} is used for Groebner basis computations,
225@*       otherwise, and by default, @code{slimgb} is used.
226@*       If printlevel=1, progress debug messages will be printed,
227@*       if printlevel>=2, all the debug messages will be printed.
228EXAMPLE: example ncsolve; shows examples
229"
230{
231  int whichengine = 0; //default
232  if (size(#)>0)
233  {
234    if (typeof(#[1])=="int" || typeof(#[1])=="number")
235    {
236      whichengine = int( #[1]);
237    }
238  }
239  int ppl = printlevel - voice +2;
240  int sI = ncols(I);
241  // check if we are done
242  if (I[sI]==0)
243  {
244    vector v = gen(sI);
245  }
246  else
247  {
248    // 1. introduce undefined coeffs
249    def save = basering;
250    int p = char(save);
251    ring @A = p,(@a(1..sI)),lp;
252    ring @aA = (p,@a(1..sI)), (@z),dp;
253    def @B = save + @aA;
254    setring @B;
255    ideal I = imap(save,I);
256    // 2. form the linear system for the undef coeffs
257    int i;   poly W;  ideal QQ;
258    for (i=1; i<=sI; i++)
259    {
260      W = W + @a(i)*I[i];
261    }
262    while (W!=0)
263    {
264      QQ = QQ,leadcoef(W);
265      W = W - lead(W);
266    }
267    // QQ consists of polynomial expressions in @a(i) of type number
268    setring @A;
269    ideal QQ = imap(@B,QQ);
270    // this QQ is a polynomial ideal, so "solve" the system
271    //    QQ = groebner(QQ);
272    dbprint(ppl, "ncsolve: starting Groebner basis computation with engine:", whichengine);
273    QQ = engine(QQ,whichengine);
274    dbprint(ppl, "QQ after groebner:", QQ);
275    //    "QQ after groebner:"; QQ;
276    if (dim(QQ) == -1)
277    {
278      dbprint(ppl+1, "no solutions by ncsolve");
279      //      "no solutions by ncsolve";
280      // output zeroes
281      setring save;
282      kill @A,@aA,@B;
283      return(v);
284    }
285    // in order to get the numeric values
286    matrix AA = matrix(maxideal(1));
287    attrib(QQ,"isSB",1); // to suppress NF warnings
288    AA = NF(AA,QQ);
289    dbprint(ppl, "AA after NF:",AA);
290    //    "AA after NF:"; print(AA);
291    for(i=1; i<=sI; i++)
292    {
293      AA = subst(AA,var(i),1);
294    }
295    dbprint(ppl, "AA after subst:",AA);
296    //    "AA after subst: "; print(AA);
297    vector v = (module(transpose(AA)))[1];
298    setring save;
299    vector v = imap(@A,v);
300    //    "v before output:";  v;
301    kill @A,@aA,@B;
302  }
303  return(v);
304}
305example
306{
307  "EXAMPLE:"; echo = 2;
308  ring r = 0,x,dp;
309  ideal I = x,2x;
310  ncsolve(I);
311  ideal J = x,x2;
312  ncsolve(J);
313}
314
315proc minpol (poly s, ideal II, list #)
316"USAGE:  minpol(f, I [,t,n]);  f a poly, I an ideal, t,n optional ints
317RETURN:  coefficient vector of the minimal polynomial of f in basering modulo I
318PURPOSE: compute the minimal polynomial of an element in a factor Weyl algebra
319NOTE:    If s<>0, @code{std} is used for Groebner basis computations,
320@*       otherwise, and by default, @code{slimgb} is used.
321@*       If n>0, the algorithm only searches for polynomials of degree n
322@*       otherwise and by default, it searches for all degrees.
323@*       If printlevel=1, progress debug messages will be printed,
324@*       if printlevel>=2, all the debug messages will be printed.
325EXAMPLE: example minpol; shows examples
326"
327{
328  // in # can be the degree of the minimal polynomial if it's known
329  // assume I is given in Groebner basis
330  ideal I = II;
331  attrib(I,"isSB",1); // set attribute for suppressing NF messages
332  int ppl = printlevel-voice+2;
333  int whichengine = 0; // default
334  int validdeg    = 0; // default
335  if (size(#)>0)
336  {
337    if (typeof(#[1])=="int" || typeof(#[1])=="number")
338    {
339      whichengine = int( #[1]);
340    }
341    if (size(#)>1)
342    {
343      if (typeof(#[2])=="int" || typeof(#[2])=="number")
344      {
345        if (#[2]>0)
346        {
347          validdeg = int( #[2]);
348        }
349      }
350    }
351  }
352  int j;
353  vector v;
354  poly p;
355  ideal NI = 1;
356  if (validdeg != 0)
357  {
358    for (j=1; j<=validdeg; j++)
359    {
360      NI[j+1] = NF(s^j,I);
361    }
362    v = ncsolve(NI,whichengine);
363  }
364  else
365  {
366    int i = 1;
367    dbprint(ppl+1,"minpolynomial starts...");
368    dbprint(ppl+1,"with ideal I=", I);
369    string st;
370    // while (1)
371    while (i<=63) // todo: remove this when NF is working properly
372    {
373      dbprint(ppl,"i:"+string(i));
374      NI[i+1] = NF(s^i,I);
375      // look for a solution
376      dbprint(ppl,"ncsolve starts with: "+string(matrix(NI)));
377      v = ncsolve(NI,whichengine);
378      matrix MM[1][nrows(v)]=matrix(v);
379      dbprint(ppl,"ncsolve ready  with: "+string(MM));
380      kill MM;
381      //  "ncsolve ready with"; print(v);
382      if (v!=0)
383      {
384        // a solution:
385        //check for the reality of the solution
386        p = 0;
387        for (j=1; j<=i+1; j++)
388        {
389          p = p + v[j]*NI[j];
390        }
391        if (p!=0)
392        {
393          dbprint(ppl,"ncsolve: bad solution!");
394        }
395        else
396        {
397          dbprint(ppl,"ncsolve: got solution!");
398          // "got solution!";
399          break;
400        }
401      }
402      // no solution:
403      i++;
404    }
405  }
406  dbprint(ppl+1,"minpol finished");
407  return(v);
408}
409example
410{
411  "EXAMPLE:"; echo = 2;
412  // todo: create example
413}
414
415
416proc bfct (poly f, list #)
417"USAGE:  bfct(f [,r,s,t,v]);  f a poly, r,s,t optional ints, v an optional intvec
418RETURN:  list of roots of the Bernstein polynomial bs(f) and its multiplicies
419PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, according to the algorithm by Masayuki Noro
420NOTE:    If r<>0, @code{std} is used for Groebner basis computations,
421@*       otherwise, and by default, @code{slimgb} is used.
422@*       If s<>0, a composed ordering is used for Groebner basis computations,
423@*       otherwise, and by default, a matrix ordering is used.
424@*       If t<>0, the minimal polynomial computation is solely performed over charasteristic 0,
425@*       otherwise and by default, a modular method is used.
426@*       If v is a positive weight vector, v is used for homogenization computations,
427@*       otherwise and by default, no weights are used.
428@*       If printlevel=1, progress debug messages will be printed,
429@*       if printlevel>=2, all the debug messages will be printed.
430EXAMPLE: example bfct; shows examples
431"
432{
433  int ppl = printlevel - voice +2;
434  int i,j;
435  // in # we have three switches:
436  // one for the engine used for Groebner basis computations,
437  // one for  M() ordering or its realization
438  // and one for a modular method when computing the minimal polynomial
439  // in # can also be the optional weight vector
440  def save = basering;
441  int n = nvars(save);
442  int whichengine  = 0; // default
443  int methodord    = 0; // default
444  int methodminpol = 0; // default
445  int validweight  = 0; // default
446  if (size(#)>0)
447  {
448    if (typeof(#[1])=="int" || typeof(#[1])=="number")
449    {
450      whichengine = int(#[1]);
451    }
452    if (size(#)>1)
453    {
454      if (typeof(#[2])=="int" || typeof(#[2])=="number")
455      {
456        methodord = int(#[2]);
457      }
458      if (size(#)>2)
459      {
460         if (typeof(#[3])=="int" || typeof(#[3])=="number")
461         {
462           methodminpol = int(#[3]);
463         }
464         if (size(#)>3)
465         {
466           if (typeof(#[4])=="intvec" && size(#[4])==n && ispositive(#[4])==1)
467           {
468             validweight = 1;
469           }
470         }
471      }
472    }
473  }
474  intvec u0;
475  if (validweight==1)
476  {
477    u0 = #[4];
478  }
479  else // no valid weight vector specified, using no weight
480  {
481    for (i=1; i<=n; i++)
482    {
483      u0[i] = 1;
484    }
485  }
486  ring r = 0,(x(1..n)),wp(u0);
487  poly f = fetch(save,f);
488  int d = deg(f); // weighted degree
489  // setting up the weighted homogenized Weyl-Algebra:
490  intvec u = d;
491  for (i=1; i<=n; i++)
492  {
493    u[i+1] = u0[n-i+1];
494  }
495  intvec v = 1;
496  for (i=1; i<=n; i++)
497  {
498    v[i+1] = d+1-u0[n-i+1];
499  }
500  intvec uv = u,v; // length 2*n + 2
501  intvec w = -1; w[n+2] = 1; w[2*n+2] = 0;
502  // try with the matrix ord
503  if (methodord == 0) // default
504  {
505    intmat @Ord[2*n+3][2*n+3];  // the ordering matrix
506    //   for (i=1; i<=2*n+2; i++)
507    //   {
508    //     @Ord[1,i] = uv[i];
509    //     @Ord[2,i] = w[i];
510    //     @Ord[3,i] = 1;
511    //   }
512    //   for (j=1; j<=2*n-1;j++)
513    //   {
514    //     @Ord[3+j,2*n+3 - j] = -1;
515    //   }
516    //   @Ord[1,2*n+3] = 1;
517    // new idea: use 1st row as the a() and the full square matrix
518    // is appended as M() block
519    intvec toa = uv,1;
520    for (i=1; i<=2*n+2; i++)
521    {
522      @Ord[1,i] = w[i];
523      @Ord[2,i] = 1;
524    }
525    for (j=1; j<=2*n+1;j++)
526    {
527      @Ord[2+j,2*n+3 - j] = -1;
528    }
529    dbprint(ppl,"the ordering matrix is",@Ord);
530    //    ring Dh = 0,(t,x(n..1),dt,d(n..1),h),M(@Ord);
531    ring Dh = 0,(t,x(n..1),dt,d(n..1),h),(a(toa),M(@Ord));
532  }
533  else
534  {
535    ring Dh = 0,(t,x(n..1),dt,d(n..1),h),(a(uv,1),a(w,0),dp(2*n+2),lp(1));
536  }
537  matrix @rel[2*n+3][2*n+3]; // non-commutative relations
538  for (i=1; i<=n+1; i++)
539  {
540    @rel[i,n+1+i] = h^(u[i]+v[i]);
541  }
542  def DDh = nc_algebra(1,@rel);
543  setring DDh;
544  poly f = imap(r,f);
545  kill r;
546  f = homog(f,h);
547  ideal I = t-f;  // defining the Malgrange ideal
548  for (i=1; i<=n; i++)
549  {
550    I = I, d(i)+diff(f,x(i))*dt;
551  }
552  dbprint(ppl, "starting Groebner basis computation with engine:", whichengine);
553  I = engine(I, whichengine);
554  dbprint(ppl, "finished Groebner basis computation");
555  I = subst(I,h,1); // dehomogenizantion
556  ring D = 0,(t,x(n..1),dt,d(n..1)),dp;
557  def DD = Weyl(D);
558  setring DD;
559  ideal I = imap(DDh,I);
560  kill Dh, DDh;
561  ideal inI = InForm(I,w); // we are only interested in the initial form w.r.t. w
562  dbprint(ppl, "starting cosmetic Groebner basis computation with engine:", whichengine);
563  inI = engine(inI, whichengine); // just for the NF
564  dbprint(ppl, "finished cosmetic Groebner basis computation");
565  dbprint(ppl, "inI is:", string(matrix(inI)));
566  poly s = t*dt;
567  vector b;
568  // try it modular
569  if (methodminpol == 0) //default
570  {
571    int lb = 2;         // todo: define better bounds
572    int ub = 2147483647;
573    i=1;
574    list usedprimes;
575    list useddegs;
576    while (b == 0)
577    {
578      dbprint(ppl,"number of run in the loop: "+string(i));
579      int q = prime(random(lb,ub));
580      if (isin(usedprimes,q)==0) // if q was not already used
581      {
582        usedprimes = usedprimes,q;
583        dbprint(ppl,"used prime is: "+string(q));
584        // change the characteristic of DD
585        ring Dq = q,(t,x(n..1),dt,d(n..1)),dp;
586        def DDq = Weyl(Dq);
587        setring DDq;
588        map phi = DD,maxideal(1);
589        poly s = phi(s);
590        ideal inI = phi(inI);
591        int sinI1 = size(inI);
592        inI = simplify(inI,2);
593        int sinI2 = size(inI);
594        if (sinI1 == sinI2) // if they differ, q was unlucky
595        {
596          vector bq = minpol(s,inI,whichengine);
597          int degb = nrows(bq) - 1;
598          dbprint(ppl,"found modular solution "+string(bq));
599          dbprint(ppl,"of degree "+string(degb));
600          setring DD;
601          kill Dq,DDq;
602          // if q was lucky, then degb is the degree of the
603          // minimal polynomial in char zero
604          if (isin(useddegs,degb)==0) // if degb was not already used
605          {
606            useddegs = useddegs,degb;
607            b = minpol(s,inI,whichengine,degb);
608          }
609        }
610        i++;
611      }
612    }
613  }
614  else // non modular
615  {
616    b = minpol(s,inI,whichengine);
617  }
618  ring @S = 0,(@s),dp;
619  vector b = imap(DD,b);
620  kill DD;
621  poly p = 0;
622  for (i=1; i<=nrows(b); i++)
623  {
624    p = p + b[i]*(@s)^(i-1);
625  }
626  p = subst(p,@s,-@s-1);
627  // the rest of this proc is nicked from bernsteinBM from dmod.lib
628  list P = factorize(p);//with constants and multiplicities
629  ideal bs; intvec m;   //the Bernstein polynomial is monic, so we are not interested in constants
630  for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
631  {
632    bs[i-1] = P[1][i];
633    m[i-1]  = P[2][i];
634  }
635  bs =  normalize(bs);
636  bs = -subst(bs,@s,0);
637  setring save;
638  ideal bs = imap(@S,bs);
639  kill @S;
640  list BS = bs,m;
641  return(BS);
642}
643example
644{
645  "EXAMPLE:"; echo = 2;
646  ring r = 0,(x),dp;
647  poly F = x;
648  bfct(F);
649  ring r = 0,(x,y),dp;
650  poly F = x2*y+x*y2;
651  bfct(F);
652}
653
654static proc engine(ideal I, int i)
655{
656  /* std - slimgb mix */
657  ideal J;
658  if (i==0)
659  {
660    J = slimgb(I);
661  }
662  else
663  {
664    // without options -> strange! (ringlist?)
665    option(redSB);
666    option(redTail);
667    J = std(I);
668  }
669  return(J);
670}
671
672static proc test_slim_vs_std()
673{
674  LIB "bfct.lib";
675  //  ring r=0,(x,y),dp;
676  //  poly f = x^4+x*y^7+y^8;
677  //  poly f = x^5+x*y^4+y^5;
678  ring r=0,(x1,x2,x3,x4,x5,x6),dp;
679  poly f = (x1*x2)^3 + (x3*x4)^2 + (x5*x6)^2;
680  option(prot);
681  printlevel = 3;
682  list Lslim = bfct(f);
683}
Note: See TracBrowser for help on using the repository browser.