source:git/Singular/LIB/bfct.lib@6cd6983

spielwiese
Last change on this file since 6cd6983 was 6cd6983, checked in by Viktor Levandovskyy <levandov@…>, 15 years ago
*levandov: improvements, bugfixes, docu git-svn-id: file:///usr/local/Singular/svn/trunk@11212 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100755`
File size: 32.4 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="\$Id: bfct.lib,v 1.6 2008-12-01 20:58:20 levandov Exp \$";
3category="Noncommutative";
4info="
5LIBRARY: bfct.lib     M. Noro's Algorithm for 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 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-Sato 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(f[,s,t,v]);            compute the global Bernstein-Sato polynomial of a given poly
21bfctsyz(f[,r,s,t,u,v]);     compute the global Bernstein-Sato polynomial of a given poly
22bfctann(f[,s]);             compute the global Bernstein-Sato polynomial of a given poly
23bfctonestep(f[,s,t]);       compute the global Bernstein-Sato polynomial of a given poly
24bfctideal(I,w[,s,t]);       compute the global b-function of a given ideal w.r.t. a given weight
25pintersect(f,I);            compute the intersection of the ideal generated by a given poly with a given ideal
26pintersectsyz(f,I[,p,s,t]); compute the intersection of the ideal generated by a given poly with a given ideal
27linreduce(f,I[,s]);         reduce a poly by linear reductions only
28ncsolve(I[,s]);             find and compute a linear dependency of the elements of an ideal
29
30AUXILIARY PROCEDURES:
31
32ispositive(v);   check whether all entries of an intvec are positive
33isin(l,i);       check whether an element is a member of a list
34scalarprod(v,w); compute the standard scalar product of two intvecs
35vec2poly(v[,i]); convert a coefficient vector to a poly
36
38";
39
40
41LIB "qhmoduli.lib"; // for Max
42LIB "dmodapp.lib";  // for initialideal etc
43
44
45proc testbfctlib ()
46{
47  // tests all procs for consistency
48  "AUXILIARY PROCEDURES:";
49  example ispositive;
50  example isin;
51  example scalarprod;
52  example vec2poly;
53  "MAIN PROCEDURES:";
54  example bfct;
55  example bfctsyz;
56  example bfctann;
57  example bfctonestep;
58  example bfctideal;
59  example pintersect;
60  example pintersectsyz;
61  example linreduce;
62  example ncsolve;
63}
64
65
66//--------------- auxiliary procedures ---------------------------------------------------------
67
68static proc gradedWeyl (intvec u,intvec v)
70RETURN:  a ring, the associated graded ring of the basering w.r.t. u and v
71PURPOSE: compute the associated graded ring of the basering w.r.t. u and v
73NOTE:    u[i] is the weight of x(i), v[i] the weight of D(i).
74@*       u+v has to be a non-negative intvec.
75"
76{
77  int i;
78  def save = basering;
79  int n = nvars(save)/2;
80  if (nrows(u)<>n || nrows(v)<>n)
81  {
82    ERROR("weight vectors have wrong dimension");
83  }
84  intvec uv,gr;
85  uv = u+v;
86  for (i=1; i<=n; i++)
87  {
88    if (uv[i]>=0)
89    {
90      if (uv[i]==0)
91      {
92        gr[i] = 0;
93      }
94      else
95      {
96        gr[i] = 1;
97      }
98    }
99    else
100    {
101      ERROR("the sum of the weight vectors has to be a non-negative intvec");
102    }
103  }
104  list l = ringlist(save);
105  list l2 = l[2];
106  matrix l6 = l[6];
107  for (i=1; i<=n; i++)
108  {
109    if (gr[i] == 1)
110    {
111       l2[n+i] = "xi("+string(i)+")";
112       l6[i,n+i] = 0;
113    }
114  }
115  l[2] = l2;
116  l[6] = l6;
117  def G = ring(l);
118  return(G);
119}
120example
121{
122  "EXAMPLE:"; echo = 2;
123  LIB "bfct.lib";
124  ring @D = 0,(x,y,z,Dx,Dy,Dz),dp;
125  def D = Weyl();
126  setring D;
127  intvec u = -1,-1,1; intvec v = 2,1,1;
129  setring G; G;
130}
131
132
133proc ispositive (intvec v)
134"USAGE:  ispositive(v);  v an intvec
135RETURN:  1 if all components of v are positive, or 0 otherwise
136PURPOSE: check whether all components of an intvec are positive
137EXAMPLE: example ispositive; shows an example
138"
139{
140  int i;
141  for (i=1; i<=size(v); i++)
142  {
143    if (v[i]<=0)
144    {
145      return(0);
146      break;
147    }
148  }
149  return(1);
150}
151example
152{
153  "EXAMPLE:"; echo = 2;
154  intvec v = 1,2,3;
155  ispositive(v);
156  intvec w = 1,-2,3;
157  ispositive(w);
158}
159
160proc isin (list l, i)
161"USAGE:  isin(l,i);  l a list, i an argument of any type
162RETURN:  an int, the position of the first appearance of i in l, or 0 if i is not a member of l
163PURPOSE: check whether the second argument is a member of a list
164EXAMPLE: example isin; shows an example
165"
166{
167  int j;
168  for (j=1; j<=size(l); j++)
169  {
170    if (l[j]==i)
171    {
172      return(j);
173      break;
174    }
175  }
176  return(0);
177}
178example
179{
180  "EXAMPLE:"; echo = 2;
181  ring r = 0,x,dp;
182  list l = 1,2,3;
183  isin(l,4);
184  isin(l,2);
185}
186
187proc scalarprod (intvec v, intvec w)
188"USAGE:  scalarprod(v,w);  v,w intvecs
189RETURN:  an int, the standard scalar product of v and w
190PURPOSE: compute the scalar product of two intvecs
191NOTE:    the arguments must have the same size
192EXAMPLE: example scalarprod; shows examples
193"
194{
195  int i; int sp;
196  if (size(v)!=size(w))
197  {
198    ERROR("non-matching dimensions");
199  }
200  else
201  {
202    for (i=1; i<=size(v);i++)
203    {
204      sp = sp + v[i]*w[i];
205    }
206  }
207  return(sp);
208}
209example
210{
211  "EXAMPLE:"; echo = 2;
212  intvec v = 1,2,3;
213  intvec w = 4,5,6;
214  scalarprod(v,w);
215}
216
217
218//-------------- main procedures -------------------------------------------------------
219
220
221proc linreduce(poly f, ideal I, list #)
222"USAGE:  linreduce(f, I [,s,t]);  f a poly, I an ideal, s,t optional ints
223RETURN:  a poly obtained by linear reductions of the given poly with the given ideal
224PURPOSE: reduce a poly only by linear reductions
225NOTE:    If s<>0, a list consisting of the reduced poly and the coefficient vector of the used
226@*       reductions is returned.
227@*       If t<>0, only leading monomials are reduced, otherwise, and by default, all monomials
228@*       are reduced, if possible.
229EXAMPLE: example linreduce; shows examples
230"
231{
232  int ppl = printlevel - voice + 2;
233  int remembercoeffs = 0; // default
234  int redlm          = 0; // default
235  if (size(#)>0)
236  {
237    if (typeof(#[1])=="int" || typeof(#[1])=="number")
238    {
239      remembercoeffs = #[1];
240    }
241    if (size(#)>1)
242    {
243      if (typeof(#[2])=="int" || typeof(#[2])=="number")
244      {
245        redlm = #[2];
246      }
247    }
248  }
249  int i,j,k;
250  int sI = ncols(I);
251  ideal lmI,lcI;
252  for (i=1; i<=sI; i++)
253  {
256  }
257  vector v;
258  poly c;
259  int reduction = 1;
260  if (redlm == 0)
261  {
262    ideal monomf;
263    for (k=1; k<=size(f); k++)
264    {
265      monomf[k] = normalize(f[k]);
266    }
267    while (reduction == 1) // while there was a reduction
268    {
269      reduction = 0;
270      for (i=sI; i>=1; i--)
271      {
272        dbprint(ppl,"testing ideal entry:",i);
273        for (j=1; j<=size(f); j++)
274        {
275          if (monomf[j] == lmI[i])
276          {
278            f = f - c*I[i];
279            dbprint(ppl,"reducing poly to ",f);
280            monomf = 0;
281            for (k=1; k<=size(f); k++)
282            {
283              monomf[k] = normalize(f[k]);
284            }
285            reduction = 1;
286            if (remembercoeffs <> 0)
287            {
288              v = v - c * gen(i);
289            }
290            break;
291          }
292        }
293        if (reduction == 1)
294        {
295          break;
296        }
297      }
298    }
299  }
300  else // reduce only leading monomials
301  {
303    while (reduction == 1) // while there was a reduction
304    {
305      reduction = 0;
306      for (i=sI;i>=1;i--)
307      {
308        if (lm <> 0 && lm == lmI[i])
309        {
311          f = f - c*I[i];
313          reduction = 1;
314          if (remembercoeffs <> 0)
315          {
316            v = v - c * gen(i);
317          }
318        }
319      }
320    }
321  }
322  if (remembercoeffs <> 0)
323  {
324    list l = f,v;
325    return(l);
326  }
327  else
328  {
329    return(f);
330  }
331}
332example
333{
334  "EXAMPLE:"; echo = 2;
335  ring r = 0,(x,y),dp;
336  ideal I = 1,y,xy;
337  poly f = 5xy+7y+3;
338  poly g = 7x+5y+3;
339  linreduce(g,I);
340  linreduce(g,I,0,1);
341  linreduce(f,I,1);
342}
343
344proc ncsolve (ideal I, list #)
345"USAGE:  ncsolve(I[,s]);  I an ideal, s an optional int
346RETURN:  coefficient vector of a linear combination of 0 in the elements of I
347PURPOSE: compute a linear dependency between the elements of an ideal if such one exists
348NOTE:    If s<>0, @code{std} is used for Groebner basis computations,
349@*       otherwise, @code{slimgb} is used.
350@*       By default, @code{slimgb} is used in char 0 and @code{std} in char >0.
351@*       If printlevel=1, progress debug messages will be printed,
352@*       if printlevel>=2, all the debug messages will be printed.
353EXAMPLE: example ncsolve; shows examples
354"
355{
356  int whichengine     = 0; // default
357  int enginespecified = 0; // default
358  if (size(#)>0)
359  {
360    if (typeof(#[1])=="int" || typeof(#[1])=="number")
361    {
362      whichengine = int( #[1]);
363      enginespecified = 1;
364    }
365  }
366  int ppl = printlevel - voice +2;
367  int sI = ncols(I);
368  // check if we are done
369  if (I[sI]==0)
370  {
371    vector v = gen(sI);
372  }
373  else
374  {
375    // 1. introduce undefined coeffs
376    def save = basering;
377    int p = char(save);
378    if (enginespecified == 0)
379    {
380      if (p <> 0)
381      {
382        whichengine = 1;
383      }
384    }
385    ring @A = p,(@a(1..sI)),lp;
386    ring @aA = (p,@a(1..sI)), (@z),dp;
387    def @B = save + @aA;
388    setring @B;
389    ideal I = imap(save,I);
390    // 2. form the linear system for the undef coeffs
391    int i;   poly W;  ideal QQ;
392    for (i=1; i<=sI; i++)
393    {
394      W = W + @a(i)*I[i];
395    }
396    while (W!=0)
397    {
399      W = W - lead(W);
400    }
401    // QQ consists of polynomial expressions in @a(i) of type number
402    setring @A;
403    ideal QQ = imap(@B,QQ);
404    // 3. this QQ is a polynomial ideal, so "solve" the system
405    dbprint(ppl, "ncsolve: starting Groebner basis computation with engine:", whichengine);
406    QQ = engine(QQ,whichengine);
407    dbprint(ppl, "QQ after engine:", QQ);
408    if (dim(QQ) == -1)
409    {
410      dbprint(ppl+1, "no solutions by ncsolve");
411      // output zeroes
412      setring save;
413      kill @A,@aA,@B;
414      return(v);
415    }
416    // 4. in order to get the numeric values
417    matrix AA = matrix(maxideal(1));
418    module MQQ = std(module(QQ));
419    AA = NF(AA,MQQ); // todo: we still receive NF warnings
420    dbprint(ppl, "AA after NF:",AA);
421    //    "AA after NF:"; print(AA);
422    for(i=1; i<=sI; i++)
423    {
424      AA = subst(AA,var(i),1);
425    }
426    dbprint(ppl, "AA after subst:",AA);
427    //    "AA after subst: "; print(AA);
428    vector v = (module(transpose(AA)))[1];
429    setring save;
430    vector v = imap(@A,v);
431    kill @A,@aA,@B;
432  }
433  return(v);
434}
435example
436{
437  "EXAMPLE:"; echo = 2;
438  ring r = 0,x,dp;
439  ideal I = x,2x;
440  ncsolve(I);
441  ideal J = x,x2;
442  ncsolve(J);
443}
444
445proc pintersect (poly s, ideal I)
446"USAGE:  pintersect(f, I);  f a poly, I an ideal
447RETURN:  coefficient vector of the monic generator of the intersection of the ideal generated by f with I
448PURPOSE: compute the intersection of an ideal with a principal ideal defined by f
449NOTE:    If the intersection is zero, this proc might not terminate.
450@*       I should be given as standard basis.
451@*       If printlevel=1, progress debug messages will be printed,
452@*       if printlevel>=2, all the debug messages will be printed.
453EXAMPLE: example pintersect; shows examples
454"
455{
456  // assume I is given in Groebner basis
457  attrib(I,"isSB",1); // set attribute for suppressing NF messages
458  int ppl = printlevel-voice+2;
459  // ---case 1: I = basering---
460  if (size(I) == 1)
461  {
462    if (simplify(I,1) == ideal(1))
463    {
464      return(gen(2)); // = s
465    }
466  }
467  def save = basering;
468  int n = nvars(save);
469  int i,j,k;
470  // ---case 2: intersection is zero---
472  intvec possdegbounds;
473  list degI;
474  i = 1;
475  while (i <= ncols(I))
476  {
477    if (i == ncols(I)+1)
478    {
479      break;
480    }
482    for (j=1; j<=n; j++)
483    {
484      if (degs[j] == 0)
485      {
486        if (degI[i][j] <> 0)
487        {
488          break;
489        }
490      }
491      if (j == n)
492      {
493        k++;
494        possdegbounds[k] = Max(degI[i]);
495      }
496    }
497    i++;
498  }
499  int degbound = Min(possdegbounds);
500  dbprint(ppl,"a lower bound for the degree of the insection is:");
501  dbprint(ppl,degbound);
502  if (degbound == 0) // lm(s) does not appear in lm(I)
503  {
504    return(vector(0));
505  }
506  // ---case 3: intersection is non-trivial---
507  ideal redNI = 1;
508  vector v;
509  list l,ll;
510  l[1] = vector(0);
511  poly toNF,tobracket,newNF,rednewNF,oldNF,secNF;
512  i = 1;
513  dbprint(ppl+1,"pintersect starts...");
514  while (1)
515  {
516    dbprint(ppl,"testing degree: "+string(i));
517    if (i>1)
518    {
519      oldNF = newNF;
520      tobracket = s^(i-1) - oldNF;
521      if (tobracket==0) // todo bug in bracket?
522      {
523        toNF = 0;
524      }
525      else
526      {
527        toNF = bracket(tobracket,secNF);
528      }
529      newNF = NF(toNF+oldNF*secNF,I);  // = NF(s^i,I)
530    }
531    else
532    {
533      newNF = NF(s,I);
534      secNF = newNF;
535    }
536    ll = linreduce(newNF,redNI,1);
537    rednewNF = ll[1];
538    l[i+1] = ll[2];
539    dbprint(ppl,"newNF is:", newNF);
540    dbprint(ppl,"rednewNF is:", rednewNF);
541    if (rednewNF != 0) // no linear dependency
542    {
543      redNI[i+1] = rednewNF;
544      i++;
545    }
546    else // there is a linear dependency, hence we are done
547    {
548      dbprint(ppl+1,"the degree of the generator of the intersection is:", i);
549      break;
550    }
551  }
552  dbprint(ppl,"used linear reductions:", l);
553  // we obtain the coefficients of the generator of the intersection by the used reductions:
554  ring @R = 0,(a(1..i+1)),dp;
555  setring @R;
556  list l = imap(save,l);
557  ideal C;
558  for (j=1;j<=i+1;j++)
559  {
560    C[j] = 0;
561    for (k=1;k<=j;k++)
562    {
563      C[j] = C[j]+l[j][k]*a(k);
564    }
565  }
566  for (j=i;j>=1;j--)
567  {
568    C[i+1]= subst(C[i+1],a(j),a(j)+C[j]);
569  }
570  matrix m = coeffs(C[i+1],maxideal(1));
571  vector v = gen(i+1);
572  for (j=1;j<=i+1;j++)
573  {
574    v = v + m[j,1]*gen(j);
575  }
576  setring save;
577  v = imap(@R,v);
578  kill @R;
579  dbprint(ppl+1,"pintersect finished");
580  return(v);
581}
582example
583{
584  "EXAMPLE:"; echo = 2;
585  ring r = 0,(x,y),dp;
586  poly f = x^2+y^3+x*y^2;
587  def D = initialmalgrange(f);
588  setring D;
589  inF;
590  pintersect(t*Dt,inF);
591}
592
593proc pintersectsyz (poly s, ideal II, list #)
594"USAGE:  pintersectsyz(f, I [,p,s,t]);  f a poly, I an ideal, p, t optial ints, p a prime number
595RETURN:  coefficient vector of the monic generator of the intersection of the ideal generated by f with I
596PURPOSE: compute the intersection of an ideal with a principal ideal defined by f
597NOTE:    If the intersection is zero, this proc might not terminate.
598@*       I should be given as standard basis.
599@*       If p>0 is given, this proc computes the generator of the intersection in char p first and
600@*       then only searches for a generator of the obtained degree in the basering.
601@*       Otherwise, it searched for all degrees.
602@*       This is done by computing syzygies.
603@*       If s<>0, @code{std} is used for Groebner basis computations in char 0,
604@*       otherwise, and by default, @code{slimgb} is used.
605@*       If t<>0 and by default, @code{std} is used for Groebner basis computations in char >0,
606@*       otherwise, @code{slimgb} is used.
607@*       If printlevel=1, progress debug messages will be printed,
608@*       if printlevel>=2, all the debug messages will be printed.
609EXAMPLE: example pintersectsyz; shows examples
610"
611{
612  // assume I is given in Groebner basis
613  ideal I = II;
614  attrib(I,"isSB",1); // set attribute for suppressing NF messages
615  int ppl = printlevel-voice+2;
616  int whichengine  = 0; // default
617  int modengine    = 1; // default
618  int solveincharp = 0; // default
619  def save = basering;
620  if (size(#)>0)
621  {
622    if (typeof(#[1])=="int" || typeof(#[1])=="number")
623    {
624      solveincharp = int(#[1]);
625    }
626    if (size(#)>1)
627    {
628      if (typeof(#[2])=="int" || typeof(#[2])=="number")
629      {
630        whichengine = int(#[2]);
631      }
632      if (size(#)>2)
633      {
634        if (typeof(#[3])=="int" || typeof(#[3])=="number")
635        {
636          modengine = int(#[3]);
637        }
638      }
639    }
640  }
641  int i,j;
642  vector v;
643  poly tobracket,toNF,newNF,p;
644  ideal NI = 1;
645  newNF = NF(s,I);
646  NI[2] = newNF;
647  if (solveincharp<>0)
648  {
649    list l = ringlist(save);
650    l[1] = solveincharp;
651    matrix l5 = l[5];
652    matrix l6 = l[6];
653    def @Rp = ring(l);
654    setring @Rp;
655    list l = ringlist(@Rp);
656    l[5] = fetch(save,l5);
657    l[6] = fetch(save,l6);
658    def Rp = ring(l);
659    setring Rp;
660    kill @Rp;
661    dbprint(ppl+1,"solving in ring ", Rp);
662    vector v;
663    map phi = save,maxideal(1);
664    poly s = phi(s);
665    ideal NI = 1;
666    setring save;
667  }
668  i = 1;
669  dbprint(ppl+1,"pintersectsyz starts...");
670  dbprint(ppl+1,"with ideal I=", I);
671  while (1)
672  {
673    dbprint(ppl,"i:"+string(i));
674    if (i>1)
675    {
676      tobracket = s^(i-1)-NI[i];
677      if (tobracket!=0)
678      {
679        toNF = bracket(tobracket,NI[2]) + NI[i]*NI[2];
680      }
681      else
682      {
683        toNF = NI[i]*NI[2];
684      }
685      newNF =  NF(toNF,I);
686      NI[i+1] = newNF;
687    }
688    // look for a solution
689    dbprint(ppl,"ncsolve starts with: "+string(matrix(NI)));
690    if (solveincharp<>0) // modular method
691    {
692      setring Rp;
693      NI[i+1] = phi(newNF);
694      v = ncsolve(NI,modengine);
695      if (v!=0) // there is a modular solution
696      {
697        dbprint(ppl,"got solution in char ",solveincharp," of degree " ,i);
698        setring save;
699        v = ncsolve(NI,whichengine);
700        if (v==0)
701        {
702          break;
703        }
704      }
705      else // no modular solution
706      {
707        setring save;
708        v = 0;
709      }
710    }
711    else // non-modular method
712    {
713      v = ncsolve(NI,whichengine);
714    }
715    matrix MM[1][nrows(v)] = matrix(v);
717    kill MM;
718    //  "ncsolve ready with"; print(v);
719    if (v!=0)
720    {
721      // a solution:
722      //check for the reality of the solution
723      p = 0;
724      for (j=1; j<=i+1; j++)
725      {
726        p = p + v[j]*NI[j];
727      }
728      if (p!=0)
729      {
731      }
732      else
733      {
734        dbprint(ppl,"ncsolve: got solution!");
735        // "got solution!";
736        break;
737      }
738    }
739    // no solution:
740    i++;
741  }
742  dbprint(ppl+1,"pintersectsyz finished");
743  return(v);
744}
745example
746{
747  "EXAMPLE:"; echo = 2;
748  ring r = 0,(x,y),dp;
749  poly f = x^2+y^3+x*y^2;
750  def D = initialmalgrange(f);
751  setring D;
752  inF;
753  poly s = t*Dt;
754  pintersectsyz(s,inF);
755  int p = prime(20000);
756  pintersectsyz(s,inF,p,0,0);
757}
758
759proc vec2poly (list #)
760"USAGE:  vec2poly(v [,i]);  v a vector or an intvec, i an optional int
761RETURN:  a poly with coefficient vector v
762PURPOSE: convert a coefficient vector to a poly
763NOTE:    If i>0 is given, the returned poly is an element of K[var(i)],
764@*       otherwise, and by default, @code{i=1} is used.
765@*       The first entry of v is the coefficient of 1.
766EXAMPLE: example vec2poly; shows examples
767"
768{
769  def save = basering;
770  int i,ringvar;
771  ringvar = 1; // default
772  if (size(#) > 0)
773  {
774    if (typeof(#[1])=="vector" || typeof(#[1])=="intvec")
775    {
776      def v = #[1];
777    }
778    else
779    {
780      ERROR("wrong input: expected vector/intvec expression");
781    }
782    if (size(#) > 1)
783    {
784      if (typeof(#[2])=="int" || typeof(#[2])=="number")
785      {
786        ringvar = int(#[2]);
787      }
788    }
789  }
790  if (ringvar > nvars(save))
791  {
792    ERROR("var out of range");
793  }
794  poly p;
795  for (i=1; i<=nrows(v); i++)
796  {
797    p = p + v[i]*(var(ringvar))^(i-1);
798  }
799  return(p);
800}
801example
802{
803  "EXAMPLE:"; echo = 2;
804  ring r = 0,(x,y),dp;
805  vector v = gen(1) + 3*gen(3) + 22/9*gen(4);
806  intvec iv = 3,2,1;
807  vec2poly(v,2);
808  vec2poly(iv);
809}
810
811static proc listofroots (list #)
812{
813  def save = basering;
814  int n = nvars(save);
815  int i;
816  poly p;
817  if (typeof(#[1])=="vector")
818  {
819    vector b = #[1];
820    for (i=1; i<=nrows(b); i++)
821    {
822      p = p + b[i]*(var(1))^(i-1);
823    }
824  }
825  else
826  {
827    p = #[1];
828  }
829  int substitution = int(#[2]);
830  ring S = 0,s,dp;
831  ideal J;
832  for (i=1; i<=n; i++)
833  {
834    J[i] = s;
835  }
836  map @m = save,J;
837  poly p = @m(p);
838  if (substitution == 1)
839  {
840    p = subst(p,s,-s-1);
841  }
842  // the rest of this proc is nicked from bernsteinBM from dmod.lib
843  list P = factorize(p);//with constants and multiplicities
844  ideal bs; intvec m;   //the Bernstein polynomial is monic, so we are not interested in constants
845  for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
846  {
847    bs[i-1] = P[1][i];
848    m[i-1]  = P[2][i];
849  }
850  bs =  normalize(bs);
851  bs = -subst(bs,s,0);
852  setring save;
853  ideal bs = imap(S,bs);
854  kill S;
855  list BS = bs,m;
856  return(BS);
857}
858
859static proc bfctengine (poly f, int inorann, int whichengine, int methodord, int methodpintersect, int pintersectchar, int modengine, intvec u0)
860{
861  int ppl = printlevel - voice +2;
862  int i;
863  def save = basering;
864  int n = nvars(save);
865  if (inorann == 0) // bfct using initial ideal
866  {
867    def D = initialmalgrange(f,whichengine,methodord,1,u0);
868    setring D;
869    ideal J = inF;
870    kill inF;
871    poly s = t*Dt;
872  }
873  else // bfct using Ann(f^s)
874  {
875    def D = SannfsBFCT(f,whichengine);
876    setring D;
877    ideal J = LD;
878    kill LD;
879  }
880  vector b;
881  // try it modular
882  if (methodpintersect <> 0) // pintersectsyz
883  {
884    if (pintersectchar == 0) // pintersectsyz::modular
885    {
886      int lb = 30000;
887      int ub = 536870909;
888      i = 1;
889      list usedprimes;
890      while (b == 0)
891      {
892        dbprint(ppl,"number of run in the loop: "+string(i));
893        int q = prime(random(lb,ub));
894        if (isin(usedprimes,q)==0) // if q was not already used
895        {
896          usedprimes = usedprimes,q;
897          dbprint(ppl,"used prime is: "+string(q));
898          b = pintersectsyz(s,J,q,whichengine,modengine);
899        }
900        i++;
901      }
902    }
903    else // pintersectsyz::non-modular
904    {
905      b = pintersectsyz(s,J,0,whichengine);
906    }
907  }
908  else // pintersect: linreduce
909  {
910    b = pintersect(s,J);
911  }
912  setring save;
913  vector b = imap(D,b);
914  if (inorann == 0)
915  {
916    list l = listofroots(b,1);
917  }
918  else
919  {
920    list l = listofroots(b,0);
921  }
922  return(l);
923}
924
925proc bfct (poly f, list #)
926"USAGE:  bfct(f [,s,t,v]);  f a poly, s,t optional ints, v an optional intvec
927RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and their multiplicies
928PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, according to the algorithm by Masayuki Noro
929NOTE:    In this proc, the initial Malgrange ideal is computed.
930@*       Further, a system of linear equations is solved by linear reductions.
931@*       If s<>0, @code{std} is used for Groebner basis computations,
932@*       otherwise, and by default, @code{slimgb} is used.
933@*       If t<>0, a matrix ordering is used for Groebner basis computations,
934@*       otherwise, and by default, a block ordering is used.
935@*       If v is a positive weight vector, v is used for homogenization computations,
936@*       otherwise and by default, no weights are used.
937@*       If printlevel=1, progress debug messages will be printed,
938@*       if printlevel>=2, all the debug messages will be printed.
939EXAMPLE: example bfct; shows examples
940"
941{
942  int ppl = printlevel - voice +2;
943  int i;
944  int n = nvars(basering);
945  // in # we have two switches:
946  // one for the engine used for Groebner basis computations,
947  // one for  M() ordering or its realization
948  // in # can also be the optional weight vector
949  int whichengine  = 0; // default
950  int methodord    = 0; // default
951  intvec u0        = 0; // default
952  if (size(#)>0)
953  {
954    if (typeof(#[1])=="int" || typeof(#[1])=="number")
955    {
956      whichengine = int(#[1]);
957    }
958    if (size(#)>1)
959    {
960      if (typeof(#[2])=="int" || typeof(#[2])=="number")
961      {
962        methodord = int(#[2]);
963      }
964      if (size(#)>2)
965      {
966        if (typeof(#[3])=="intvec" && size(#[3])==n && ispositive(#[3])==1)
967        {
968          u0 = #[3];
969        }
970      }
971    }
972  }
973  list b = bfctengine(f,0,whichengine,methodord,0,0,0,u0);
974  return(b);
975}
976example
977{
978  "EXAMPLE:"; echo = 2;
979  ring r = 0,(x,y),dp;
980  poly f = x^2+y^3+x*y^2;
981  bfct(f);
982  intvec v = 3,2;
983  bfct(f,1,0,v);
984}
985
986proc bfctsyz (poly f, list #)
987"USAGE:  bfctsyz(f [,r,s,t,u,v]);  f a poly, r,s,t,u optional ints, v an optional intvec
988RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and its multiplicies
989PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, according to the algorithm by Masayuki Noro
990NOTE:    In this proc, the initial Malgrange ideal is computed.
991@*       Further, a system of linear equations is solved by computing syzygies.
992@*       If r<>0, @code{std} is used for Groebner basis computations in characteristic 0,
993@*       otherwise, and by default, @code{slimgb} is used.
994@*       If s<>0, a matrix ordering is used for Groebner basis computations,
995@*       otherwise, and by default, a block ordering is used.
996@*       If t<>0, the computation of the intersection is solely performed over charasteristic 0,
997@*       otherwise and by default, a modular method is used.
998@*       If u<>0 and by default, @code{std} is used for Groebner basis computations in characteristic >0,
999@*       otherwise, @code{slimgb} is used.
1000@*       If v is a positive weight vector, v is used for homogenization computations,
1001@*       otherwise and by default, no weights are used.
1002@*       If printlevel=1, progress debug messages will be printed,
1003@*       if printlevel>=2, all the debug messages will be printed.
1004EXAMPLE: example bfct; shows examples
1005"
1006{
1007  int ppl = printlevel - voice +2;
1008  int i;
1009  // in # we have four switches:
1010  // one for the engine used for Groebner basis computations in char 0,
1011  // one for  M() ordering or its realization
1012  // one for a modular method when computing the intersection
1013  // and one for the engine used for Groebner basis computations in char >0
1014  // in # can also be the optional weight vector
1015  def save = basering;
1016  int n = nvars(save);
1017  int whichengine = 0; // default
1018  int methodord   = 0; // default
1019  int pintersectchar  = 0; // default
1020  int modengine   = 1; // default
1021  intvec u0       = 0; // default
1022  if (size(#)>0)
1023  {
1024    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1025    {
1026      whichengine = int(#[1]);
1027    }
1028    if (size(#)>1)
1029    {
1030      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1031      {
1032        methodord = int(#[2]);
1033      }
1034      if (size(#)>2)
1035      {
1036        if (typeof(#[3])=="int" || typeof(#[3])=="number")
1037        {
1038          pintersectchar = int(#[3]);
1039        }
1040        if (size(#)>3)
1041        {
1042          if (typeof(#[4])=="int" || typeof(#[4])=="number")
1043          {
1044            modengine = int(#[4]);
1045          }
1046          if (size(#)>4)
1047          {
1048            if (typeof(#[5])=="intvec" && size(#[5])==n && ispositive(#[5])==1)
1049            {
1050              u0 = #[5];
1051            }
1052          }
1053        }
1054      }
1055    }
1056  }
1057  list b = bfctengine(f,0,whichengine,methodord,1,pintersectchar,modengine,u0);
1058  return(b);
1059}
1060example
1061{
1062  "EXAMPLE:"; echo = 2;
1063  ring r = 0,(x,y),dp;
1064  poly f = x^2+y^3+x*y^2;
1065  bfctsyz(f);
1066  intvec v = 3,2;
1067  bfctsyz(f,0,1,1,0,v);
1068}
1069
1070proc bfctideal (ideal I, intvec w, list #)
1071"USAGE:  bfctideal(I,w[,s,t]);  I an ideal, w an intvec, s,t optional ints
1072RETURN:  list of roots and their multiplicies of the global b-function of I w.r.t. the weight vector (-w,w)
1073PURPOSE: compute the global b-function of an ideal according to the algorithm by M. Noro
1074NOTE:    Assume, I is an ideal in the n-th Weyl algebra where the sequence of the
1075@*       variables is x(1),...,x(n),D(1),...,D(n).
1076@*       If s<>0, @code{std} is used for Groebner basis computations in characteristic 0,
1077@*       otherwise, and by default, @code{slimgb} is used.
1078@*       If t<>0, a matrix ordering is used for Groebner basis computations,
1079@*       otherwise, and by default, a block ordering is used.
1080@*       If printlevel=1, progress debug messages will be printed,
1081@*       if printlevel>=2, all the debug messages will be printed.
1082EXAMPLE: example bfctideal; shows examples
1083"
1084{
1085  int ppl = printlevel - voice +2;
1086  int i;
1087  def save = basering;
1088  int n = nvars(save)/2;
1089  int whichengine = 0; // default
1090  int methodord   = 0; // default
1091  if (size(#)>0)
1092  {
1093    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1094    {
1095      whichengine = int(#[1]);
1096    }
1097    if (size(#)>1)
1098    {
1099      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1100      {
1101        methodord = int(#[2]);
1102      }
1103    }
1104  }
1105  ideal J = initialideal(I,-w,w,whichengine,methodord);
1106  poly s;
1107  for (i=1; i<=n; i++)
1108  {
1109    s = s + w[i]*var(i)*var(n+i);
1110  }
1111  vector b = pintersect(s,J);
1112  list l = listofroots(b,0);
1113  return(l);
1114}
1115example
1116{
1117  "EXAMPLE:"; echo = 2;
1118  ring @D = 0,(x,y,Dx,Dy),dp;
1119  def D = Weyl();
1120  setring D;
1121  ideal I = 3*x^2*Dy+2*y*Dx,2*x*Dx+3*y*Dy+6;
1122  intvec w1 = 1,1;
1123  intvec w2 = 1,2;
1124  intvec w3 = 2,3;
1125  bfctideal(I,w1);
1126  bfctideal(I,w2,1);
1127  bfctideal(I,w3,0,1);
1128}
1129
1130
1131proc bfctonestep (poly f,list #)
1132"USAGE:  bfctonestep(f [,s,t]);  f a poly, s,t optional ints
1133RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and its multiplicies
1134PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, using only one Groebner basis computation
1135NOTE:    If s<>0, @code{std} is used for the Groebner basis computation, otherwise,
1136@*       and by default, @code{slimgb} is used.
1137@*       If t<>0, a matrix ordering is used for Groebner basis computations,
1138@*       otherwise, and by default, a block ordering is used.
1139@*       If printlevel=1, progress debug messages will be printed,
1140@*       if printlevel>=2, all the debug messages will be printed.
1141EXAMPLE: example bfctonestep; shows examples
1142"
1143{
1144  int ppl = printlevel - voice +2;
1145  def save = basering;
1146  int n = nvars(save);
1147  int i;
1148  int whichengine = 0; // default
1149  int methodord   = 0; // default
1150  if (size(#)>0)
1151  {
1152    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1153    {
1154      whichengine = int(#[1]);
1155    }
1156    if (size(#)>1)
1157    {
1158      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1159      {
1160        methodord = int(#[2]);
1161      }
1162    }
1163  }
1164  def DDh = initialidealengine("bfctonestep", whichengine, methodord, f);
1165  setring DDh;
1166  dbprint(ppl, "the initial ideal:", string(matrix(inF)));
1167  intvec tonselect = 1;
1168  for (i=3; i<=2*n+4; i++)
1169  {
1170    tonselect = tonselect,i;
1171  }
1172  inF = nselect(inF,tonselect);
1173  dbprint(ppl, "generators containing only s:", string(matrix(inF)));
1174  inF = engine(inF, whichengine); // is now a principal ideal;
1175  setring save;
1176  ideal J; J[2] = var(1);
1177  map @m = DDh,J;
1178  ideal inF = @m(inF);
1179  poly p = inF[1];
1180  list l = listofroots(p,1);
1181  return(l);
1182}
1183example
1184{
1185  "EXAMPLE:"; echo = 2;
1186  ring r = 0,(x,y),dp;
1187  poly f = x^2+y^3+x*y^2;
1188  bfctonestep(f);
1189  bfctonestep(f,1,1);
1190}
1191
1192proc bfctann (poly f, list #)
1193"USAGE:  bfctann(f [,r]);  f a poly, r an optional int
1194RETURN:  list of roots of the Bernstein-Sato polynomial bs(f) and their multiplicies
1195PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f
1196NOTE:    In this proc, ann(f^s) is computed.
1197@*       If r<>0, @code{std} is used for Groebner basis computations,
1198@*       otherwise, and by default, @code{slimgb} is used.
1199@*       If printlevel=1, progress debug messages will be printed,
1200@*       if printlevel>=2, all the debug messages will be printed.
1201EXAMPLE: example bfctann; shows examples
1202"
1203{
1204  def save = basering;
1205  int ppl = printlevel - voice + 2;
1206  int whichengine = 0; // default
1207  if (size(#)>0)
1208  {
1209    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1210    {
1211      whichengine = int(#[1]);
1212    }
1213  }
1214  list b = bfctengine(f,1,whichengine,0,1,0,0,0);
1215  return(b);
1216}
1217example
1218{
1219  "EXAMPLE:"; echo = 2;
1220  ring r = 0,(x,y),dp;
1221  poly f = x^2+y^3+x*y^2;
1222  bfctann(f);
1223}
1224
1225static proc hardexamples ()
1226{
1227  //  some hard examples
1228  ring r1 = 0,(x,y,z,w),dp;
1229  // ab34
1230  poly ab34 = (z3+w4)*(3z2x+4w3y);
1231  bfct(ab34);
1232  // ha3
1233  poly ha3 = xyzw*(x+y)*(x+z)*(x+w)*(y+z+w);
1234  bfct(ha3);
1235  // ha4
1236  poly ha4 = xyzw*(x+y)*(x+z)*(x+w)*(y+z)*(y+w);
1237  bfct(ha4);
1238  // chal4: reiffen(4,5)*reiffen(5,4)
1239  ring r2 = 0,(x,y),dp;
1240  poly chal4 = (x4+xy4+y5)*(x5+x4y+y4);
1241  bfct(chal4);
1242  // (xy+z)*reiffen(4,5)
1243  ring r3 = 0,(x,y,z),dp;
1244  poly xyzreiffen45 = (xy+z)*(y4+yz4+z5);
1245  bfct(xyzreiffen45);
1246}
1247
Note: See TracBrowser for help on using the repository browser.