source: git/Singular/LIB/dmod.lib @ 75cc14

spielwiese
Last change on this file since 75cc14 was 75cc14, checked in by Viktor Levandovskyy <levandov@…>, 16 years ago
*levandov: doc changes git-svn-id: file:///usr/local/Singular/svn/trunk@10446 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 118.1 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id: dmod.lib,v 1.23 2007-11-27 14:40:11 levandov Exp $";
3category="Noncommutative";
4info="
5LIBRARY: dmod.lib     Algorithms for algebraic D-modules
6AUTHORS: Viktor Levandovskyy,     levandov@risc.uni-linz.ac.at
7@*       Jorge Martin Morales,    jorge@unizar.es
8
9THEORY: Given a polynomial ring R = K[x_1,...,x_n] and a polynomial F in R,
10@*      one is interested in the R[1/F]-module of rank one, generated by F^s
11@*      for a natural number s.
12@* In fact, the module R[1/F]*F^s has a structure of a D(R)[s]-module, where D(R)
13@* is an n-th Weyl algebra K<x_1,...,x_n,d_1,...,d_n | d_j x_j = x_j d_j +1> and
14@* D(R)[s] = D(R) tensored with K[s] over K.
15@* Constructively, one needs to find a left ideal I = I(F^s) in D(R), such
16@* that K[x_1,...,x_n,1/F]*F^s is isomorphic to D(R)/I as a D(R)-module.
17@* We often write just D for D(R) and D[s] for D(R)[s].
18@* One is interested in the following data:
19@* - Ann F^s = I = I(F^s) in D(R)[s], denoted by LD in the output
20@* - global Bernstein polynomial in K[s], denoted by bs, its minimal integer root s0 and
21@*   the list of all roots of bs, which are rational, with their multiplicities is denoted by BS
22@* - Ann F^s0 = I(F^s0) in D(R), denoted by LD0 in the output (LD0 is a holonomic ideal in D(R))
23@* - Ann^(1) F^s in D(R)[s], denoted by LD1 (logarithmic derivations)
24@* - an operator in D(R)[s], denoted by PS, such that PS*F^(s+1) = bs*F^s holds in K[x_1,...,x_n,1/F^s].
25
26@* We provide the following implementations:
27@* OT) the classical Ann F^s algorithm from Oaku and Takayama (J. Pure
28        Applied Math., 1999),
29@* LOT) Levandovskyy's modification of the Oaku-Takayama algorithm (unpublished)
30@* BM) the Ann F^s algorithm by Briancon and Maisonobe (Remarques sur
31        l'ideal de Bernstein associe a des polynomes, preprint, 2002)
32
33GUIDE:
34@* - Ann F^s = I = I(F^s) = LD in D(R)[s] can be computed by SannfsBM, SannfsOT, SannfsLOT
35@* - Ann^(1) F^s in D(R)[s] can be computed by Sannfslog
36@* - global Bernstein polynomial bs resp. BS in K[s] can be computed by bernsteinBM
37@* - Ann F^s0 = I(F^s0) = LD0 in D(R) can be computed by annfs0, annfsBM, annfsOT, annfsLOT
38@* - all the relevant data (LD, LD0, bs, PS) are computed by operatorBM
39
40MAIN PROCEDURES:
41
42annfs(F[,S,eng]);       compute Ann F^s0 in D and Bernstein poly for a poly F
43annfspecial(I, F, m, n);  compute Ann F^n from Ann F^s for a poly F and a number n
44Sannfs(F[,S,eng]);      compute Ann F^s in D[s] for a poly F
45Sannfslog(F[,eng]);     compute Ann^(1) F^s in D[s] for a poly F
46bernsteinBM(F[,eng]);   compute global Bernstein poly for a poly F (algorithm of Briancon-Maisonobe)
47operatorBM(F[,eng]);    compute Ann F^s, Ann F^s0, BS and PS for a poly F (algorithm of Briancon-Maisonobe)
48annfsParamBM(F[,eng]);  compute the generic Ann F^s (algorithm by Briancon and Maisonobe) and exceptional parametric constellations for a poly F with parametric coefficients
49annfsBMI(F[,eng]);      compute Ann F^s and Bernstein ideal for a poly F=f1*..*fP (multivariate algorithm of Briancon-Maisonobe)
50checkRoot(F,a[,S,eng]); check if a given rational is a root of the global Bernstein polynomial of F and compute its multiplicity
51
52SECONDARY PROCEDURES FOR D-MODULES:
53
54annfsBM(F[,eng]);          compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of Briancon-Maisonobe)
55annfsLOT(F[,eng]);         compute Ann F^s0 in D and Bernstein poly for a poly F (Levandovskyy modification of the Oaku-Takayama algorithm)
56annfsOT(F[,eng]);          compute Ann F^s0 in D and Bernstein poly for a poly F (algorithm of Oaku-Takayama)
57SannfsBM(F[,eng]);         compute Ann F^s in D[s] for a poly F (algorithm of Briancon-Maisonobe)
58SannfsLOT(F[,eng]);        compute Ann F^s in D[s] for a poly F (Levandovskyy modification of the Oaku-Takayama algorithm)
59SannfsOT(F[,eng]);         compute Ann F^s in D[s] for a poly F (algorithm of Oaku-Takayama)
60annfs0(I,F [,eng]);        compute Ann F^s0 in D and Bernstein poly from the known Ann F^s in D[s]
61checkRoot1(I,F,a[,eng]);   check whether a rational is a root of the global Bernstein polynomial of F from the known Ann F^s in D[s]
62checkRoot2(I,F,a[,eng]);   check whether a rational is a root of the global Bernstein polynomial of F and compute its multiplicity from the known Ann F^s in D[s]
63checkFactor(I,F,qs[,eng]); check whether a polynomial qs in K[s] is a factor of the global Bernstein polynomial of F from the known Ann F^s in D[s]
64
65AUXILIARY PROCEDURES:
66
67arrange(p);         create a poly, describing a full hyperplane arrangement
68reiffen(p,q);       create a poly, describing a Reiffen curve
69isHolonomic(M);     check whether a module is holonomic
70convloc(L);         replace global orderings with local in the ringlist L
71minIntRoot(P,fact); minimal integer root among the roots of a maximal ideal P
72varnum(s);          the number of the variable with the name s
73
74SEE ALSO: gmssing_lib
75";
76
77LIB "nctools.lib";
78LIB "elim.lib";
79LIB "qhmoduli.lib"; // for Max
80LIB "gkdim.lib";
81LIB "gmssing.lib";
82LIB "control.lib";  // for genericity
83
84proc tstdmodlib()
85{
86  /* tests all procs for consistency */
87  /* adding the new proc, add it here */
88
89  "MAIN PROCEDURES:";
90  example annfs;
91  example annfs0;
92  example Sannfs;
93  example Sannfslog;
94  example bernsteinBM;
95  example operatorBM;
96  example annfspecial;
97  example annfsParamBM;
98  example annfsBMI;
99  example checkRoot;
100  "SECONDARY D-MOD PROCEDURES:";
101  example annfsBM;
102  example annfsLOT;
103  example annfsOT;
104  example SannfsBM;
105  example SannfsLOT;
106  example SannfsOT;
107  example annfs0;
108  example checkRoot1;
109  example checkRoot2;
110  example checkFactor;
111}
112
113static proc engine(ideal I, int i)
114{
115  /* std - slimgb mix */
116  ideal J;
117  if (i==0)
118  {
119    J = slimgb(I);
120  }
121  else
122  {
123    // without options -> strange! (ringlist?)
124    option(redSB);
125    option(redTail);
126    J = std(I);
127  }
128  return(J);
129}
130
131
132// new top-level procedures
133proc annfs(poly F, list #)
134"USAGE:  annfs(f [,S,eng]);  f a poly, S a string, eng an optional int
135RETURN:  ring
136PURPOSE: compute the D-module structure of basering[1/f]*f^s with the algorithm given in S and with the Groebner basis engine given in 'eng'
137NOTE:    activate the output ring with the @code{setring} command.
138@*       The value of a string S can be
139@*       'bm' (default) - for the algorithm of Briancon and Maisonobe,
140@*       'ot'  - for the algorithm of Oaku and Takayama,
141@*       'lot' - for the Levandovskyy's modification of the algorithm of Oaku and Takayama.
142@*       If eng <>0, @code{std} is used for Groebner basis computations,
143@*       otherwise and by default @code{slimgb} is used.
144@*       In the output ring:
145@*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
146@*       - the list  BS is the list of roots and multiplicities of a Bernstein polynomial of f.
147@*       If @code{printlevel}=1, progress debug messages will be printed,
148@*       if @code{printlevel}>=2, all the debug messages will be printed.
149EXAMPLE: example annfs; shows examples
150"
151{
152  int eng = 0;
153  int chs = 0; // choice
154  string algo = "bm";
155  string st;
156  if ( size(#)>0 )
157  {
158   if ( typeof(#[1]) == "string" )
159   {
160     st = string(#[1]);
161     if ( (st == "BM") || (st == "Bm") || (st == "bM") ||(st == "bm"))
162     {
163       algo = "bm";
164       chs  = 1;
165     }
166     if ( (st == "OT") || (st == "Ot") || (st == "oT") || (st == "ot"))
167     {
168       algo = "ot";
169       chs  = 1;
170     }
171     if ( (st == "LOT") || (st == "lOT") || (st == "Lot") || (st == "lot"))
172     {
173       algo = "lot";
174       chs  = 1;
175     }
176     if (chs != 1)
177     {
178       // incorrect value of S
179       print("Incorrect algorithm given, proceed with the default BM");
180       algo = "bm";
181     }
182     // second arg
183     if (size(#)>1)
184     {
185       // exists 2nd arg
186       if ( typeof(#[2]) == "int" )
187       {
188         // the case: given alg, given engine
189         eng = int(#[2]);
190       }
191       else
192       {
193         eng = 0;
194       }
195     }
196     else
197     {
198       // no second arg
199       eng = 0;
200     }
201   }
202   else
203   {
204     if ( typeof(#[1]) == "int" )
205     {
206     // the case: default alg, engine
207     eng = int(#[1]);
208     //      algo = "bm";  //is already set
209     }
210     else
211     {
212       // incorr. 1st arg
213       algo = "bm";
214     }
215   }
216  }
217
218   // size(#)=0, i.e. there is no algorithm and no engine given
219   //  eng = 0; algo = "bm";  //are already set
220   // int ppl = printlevel-voice+2;
221
222  int old_printlevel = printlevel;
223  printlevel=printlevel+1;
224  def save = basering;
225  if ( algo=="ot")
226  {
227    def @A = annfsOT(F,eng);
228  }
229  else
230  {
231    if ( algo=="lot")
232    {
233      def @A = annfsLOT(F,eng);
234    }
235    else
236    {
237      // bm = default
238      def @A = annfsBM(F,eng);
239    }
240  }
241  printlevel = old_printlevel;
242  return(@A);
243}
244example
245{
246  "EXAMPLE:"; echo = 2;
247  ring r = 0,(x,y,z),Dp;
248  poly F = z*x^2+y^3;
249  printlevel = 0;
250  def A  = annfs(F); // here, the default BM algorithm will be used
251  setring A;
252  LD;
253  BS;
254}
255
256
257proc Sannfs(poly F, list #)
258"USAGE:  Sannfs(f [,S,eng]);  f a poly, S a string, eng an optional int
259RETURN:  ring
260PURPOSE: compute the D-module structure of basering[f^s] with the algorithm given in S and with the Groebner basis engine given in eng
261NOTE:    activate the output ring with the @code{setring} command.
262@*       The value of a string S can be
263@*       'bm' (default) - for the algorithm of  Briancon and Maisonobe,
264@*       'lot' - for the Levandovskyy's modification of the algorithm of Oaku and Takayama,
265@*       'ot'  - for the algorithm of Oaku and Takayama.
266@*       If eng <>0, @code{std} is used for Groebner basis computations,
267@*       otherwise, and by default @code{slimgb} is used.
268@*       In the output ring:
269@*       - the ideal LD is the needed D-module structure.
270@*       If @code{printlevel}=1, progress debug messages will be printed,
271@*       if @code{printlevel}>=2, all the debug messages will be printed.
272EXAMPLE: example Sannfs; shows examples
273"
274{
275  int eng = 0;
276  int chs = 0; // choice
277  string algo = "bm";
278  string st;
279  if ( size(#)>0 )
280  {
281   if ( typeof(#[1]) == "string" )
282   {
283     st = string(#[1]);
284     if ( (st == "BM") || (st == "Bm") || (st == "bM") ||(st == "bm"))
285     {
286       algo = "bm";
287       chs  = 1;
288     }
289     if ( (st == "OT") || (st == "Ot") || (st == "oT") || (st == "ot"))
290     {
291       algo = "ot";
292       chs  = 1;
293     }
294     if ( (st == "LOT") || (st == "lOT") || (st == "Lot") || (st == "lot"))
295     {
296       algo = "lot";
297       chs  = 1;
298     }
299     if (chs != 1)
300     {
301       // incorrect value of S
302       print("Incorrect algorithm given, proceed with the default BM");
303       algo = "bm";
304     }
305     // second arg
306     if (size(#)>1)
307     {
308       // exists 2nd arg
309       if ( typeof(#[2]) == "int" )
310       {
311         // the case: given alg, given engine
312         eng = int(#[2]);
313       }
314       else
315       {
316         eng = 0;
317       }
318     }
319     else
320     {
321       // no second arg
322       eng = 0;
323     }
324   }
325   else
326   {
327     if ( typeof(#[1]) == "int" )
328     {
329     // the case: default alg, engine
330     eng = int(#[1]);
331     //      algo = "bm";  //is already set
332     }
333     else
334     {
335       // incorr. 1st arg
336       algo = "bm";
337     }
338   }
339  }
340  // size(#)=0, i.e. there is no algorithm and no engine given
341  //  eng = 0; algo = "bm";  //are already set
342  // int ppl = printlevel-voice+2;
343  printlevel=printlevel+1;
344  def save = basering;
345  if ( algo=="ot")
346  {
347    def @A = SannfsOT(F,eng);
348  }
349  else
350  {
351    if ( algo=="lot")
352    {
353      def @A = SannfsLOT(F,eng);
354    }
355    else
356    {
357      // bm = default
358      def @A = SannfsBM(F,eng);
359    }
360  }
361  printlevel=printlevel-1;
362  return(@A);
363}
364example
365{
366  "EXAMPLE:"; echo = 2;
367  ring r = 0,(x,y,z),Dp;
368  poly F = x^3+y^3+z^3;
369  printlevel = 0;
370  def A  = Sannfs(F); // here, the default BM algorithm will be used
371  setring A;
372  LD;
373}
374
375proc Sannfslog (poly F, list #)
376"USAGE:  Sannfslog(f [,eng]);  f a poly, eng an optional int
377RETURN:  ring
378PURPOSE: compute the D-module structure of basering[1/f]*f^s, where D is the Weyl algebra
379NOTE:    activate this ring with the @code{setring} command.
380@*       In the ring D[s], the ideal LD1 is generated by the elements in Ann F^s in D[s] coming from logarithmic derivations.
381@*       If eng <>0, @code{std} is used for Groebner basis computations,
382@*       otherwise, and by default @code{slimgb} is used.
383@*       If printlevel=1, progress debug messages will be printed,
384@*       if printlevel>=2, all the debug messages will be printed.
385EXAMPLE: example Sannfslog; shows examples
386"
387{
388  int eng = 0;
389  if ( size(#)>0 )
390  {
391    if ( typeof(#[1]) == "int" )
392    {
393      eng = int(#[1]);
394    }
395  }
396  int ppl = printlevel-voice+2;
397  def save = basering;
398  int N = nvars(basering);
399  int Nnew = 2*N+1;
400  int i;
401  string s;
402  list RL = ringlist(basering);
403  list L, Lord;
404  list tmp;
405  intvec iv;
406  L[1] = RL[1]; // char
407  L[4] = RL[4]; // char, minpoly
408  // check whether vars have admissible names
409  list Name = RL[2];
410  for (i=1; i<=N; i++)
411  {
412    if (Name[i] == "s")
413    {
414      ERROR("Variable names should not include s");
415    }
416  }
417  // the ideal I
418  ideal I = -F, jacob(F);
419  dbprint(ppl,"// -1-1- starting the computation of syz(-F,_Dx(F))");
420  dbprint(ppl-1, I);
421  matrix M = syz(I);
422  M = transpose(M);  // it is more usefull working with columns
423  dbprint(ppl,"// -1-2- the module syz(-F,_Dx(F)) has been computed");
424  dbprint(ppl-1, M);
425  // ------------ the ring @R ------------
426  // _x, _Dx, s;  elim.ord for _x,_Dx.
427  // now, create the names for new vars
428  list DName;
429  for (i=1; i<=N; i++)
430  {
431    DName[i] = "D"+Name[i]; // concat
432  }
433  tmp[1] = "s";
434  list NName;
435  for (i=1; i<=N; i++)
436  {
437    NName[2*i-1] = Name[i];
438    NName[2*i] = DName[i];
439    //NName[2*i-1] = DName[i];
440    //NName[2*i] = Name[i];
441  }
442  NName[Nnew] = tmp[1];
443  L[2] = NName;
444  tmp = 0;
445  // block ord (a(1,1),a(0,0,1,1),...,dp);
446  //list("a",intvec(1,1)), list("a",intvec(0,0,1,1)), ...
447  tmp[1] = "a";  // string
448  for (i=1; i<=N; i++)
449  {
450    iv[2*i-1] = 1;
451    iv[2*i]   = 1;
452    tmp[2]    = iv;  iv = 0;  // intvec
453    Lord[i]   = tmp;
454  }
455  //list("dp",intvec(1,1,1,1,1,...))
456  s = "iv=";
457  for (i=1; i<=Nnew; i++)
458  {
459    s = s+"1,";
460  }
461  s[size(s)]=";";
462  execute(s);
463  kill s;
464  tmp[1] = "dp";  // string
465  tmp[2] = iv;    // intvec
466  Lord[N+1] = tmp;
467  //list("C",intvec(0))
468  tmp[1] = "C";  // string
469  iv = 0;
470  tmp[2] = iv;   // intvec
471  Lord[N+2] = tmp;
472  tmp = 0;
473  L[3]    = Lord;
474  // we are done with the list. Now add a Plural part
475  def @R@ = ring(L);
476  setring @R@;
477  matrix @D[Nnew][Nnew];
478  for (i=1; i<=N; i++)
479  {
480    @D[2*i-1,2*i]=1;
481    //@D[2*i-1,2*i]=-1;
482  }
483  def @R = nc_algebra(1,@D);
484  setring @R;
485  kill @R@;
486  dbprint(ppl,"// -2-1- the ring @R(_x,_Dx,s) is ready");
487  dbprint(ppl-1, @R);
488  matrix M = imap(save,M);
489  // now, create the vector [-s,_Dx]
490  vector v = [-s];  // now s is a variable
491  for (i=1; i<=N; i++)
492  {
493    v = v + var(2*i)*gen(i+1);
494    //v = v + var(2*i-1)*gen(i+1);
495  }
496  ideal J = ideal(M*v);
497  // make leadcoeffs positive
498  for (i=1; i<= ncols(J); i++)
499  {
500    if ( leadcoef(J[i])<0 )
501    {
502      J[i] = -J[i];
503    }
504  }
505  ideal LD1 = J;
506  kill J;
507  export LD1;
508  return(@R);
509}
510example
511{
512  "EXAMPLE:"; echo = 2;
513  ring r = 0,(x,y),Dp;
514  poly F = x^4+y^5+x*y^4;
515  printlevel = 0;
516  def A  = Sannfslog(F);
517  setring A;
518  LD1;
519}
520
521
522// alternative code for SannfsBM, renamed from annfsBM to ALTannfsBM
523// is superfluos - will not be included in the official documentation
524proc ALTannfsBM (poly F, list #)
525"USAGE:  ALTannfsBM(f [,eng]);  f a poly, eng an optional int
526RETURN:  ring
527PURPOSE: compute the annihilator ideal of f^s in D[s], where D is the Weyl Algebra, according to the algorithm by Briancon and Maisonobe
528NOTE:    activate this ring with the @code{setring} command. In this ring,
529@*       - the ideal LD is the annihilator of f^s.
530@*       If eng <>0, @code{std} is used for Groebner basis computations,
531@*       otherwise, and by default @code{slimgb} is used.
532@*       If printlevel=1, progress debug messages will be printed,
533@*       if printlevel>=2, all the debug messages will be printed.
534EXAMPLE: example ALTannfsBM; shows examples
535"
536{
537  int eng = 0;
538  if ( size(#)>0 )
539  {
540    if ( typeof(#[1]) == "int" )
541    {
542      eng = int(#[1]);
543    }
544  }
545  // returns a list with a ring and an ideal LD in it
546  int ppl = printlevel-voice+2;
547  //  printf("plevel :%s, voice: %s",printlevel,voice);
548  def save = basering;
549  int N = nvars(basering);
550  int Nnew = 2*N+2;
551  int i,j;
552  string s;
553  list RL = ringlist(basering);
554  list L, Lord;
555  list tmp;
556  intvec iv;
557  L[1] = RL[1];  //char
558  L[4] = RL[4];  //char, minpoly
559  // check whether vars have admissible names
560  list Name  = RL[2];
561  list RName;
562  RName[1] = "t";
563  RName[2] = "s";
564  for (i=1; i<=N; i++)
565  {
566    for(j=1; j<=size(RName); j++)
567    {
568      if (Name[i] == RName[j])
569      {
570        ERROR("Variable names should not include t,s");
571      }
572    }
573  }
574  // now, create the names for new vars
575  list DName;
576  for (i=1; i<=N; i++)
577  {
578    DName[i] = "D"+Name[i];  //concat
579  }
580  tmp[1] = "t";
581  tmp[2] = "s";
582  list NName = tmp + Name + DName;
583  L[2]   = NName;
584  // Name, Dname will be used further
585  kill NName;
586  // block ord (lp(2),dp);
587  tmp[1]  = "lp"; // string
588  iv      = 1,1;
589  tmp[2]  = iv; //intvec
590  Lord[1] = tmp;
591  // continue with dp 1,1,1,1...
592  tmp[1]  = "dp"; // string
593  s       = "iv=";
594  for (i=1; i<=Nnew; i++)
595  {
596    s = s+"1,";
597  }
598  s[size(s)]= ";";
599  execute(s);
600  kill s;
601  tmp[2]    = iv;
602  Lord[2]   = tmp;
603  tmp[1]    = "C";
604  iv        = 0;
605  tmp[2]    = iv;
606  Lord[3]   = tmp;
607  tmp       = 0;
608  L[3]      = Lord;
609  // we are done with the list
610  def @R@ = ring(L);
611  setring @R@;
612  matrix @D[Nnew][Nnew];
613  @D[1,2]=t;
614  for(i=1; i<=N; i++)
615  {
616    @D[2+i,N+2+i]=1;
617  }
618  // L[5] = matrix(UpOneMatrix(Nnew));
619  // L[6] = @D;
620  def @R = nc_algebra(1,@D);
621  setring @R;
622  kill @R@;
623  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
624  dbprint(ppl-1, @R);
625  // create the ideal I
626  poly  F = imap(save,F);
627  ideal I = t*F+s;
628  poly p;
629  for(i=1; i<=N; i++)
630  {
631    p = t;  //t
632    p = diff(F,var(2+i))*p;
633    I = I, var(N+2+i) + p;
634  }
635  // -------- the ideal I is ready ----------
636  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
637  dbprint(ppl-1, I);
638  ideal J = engine(I,eng);
639  ideal K = nselect(J,1);
640  kill I,J;
641  dbprint(ppl,"// -1-3- t is eliminated");
642  dbprint(ppl-1, K);  //K is without t
643  // create Dn[s], where Dn is the ordinary Weyl Algebra, and put the result into it,
644  // thus creating the ring @R2
645  // keep: N, i,j,s, tmp, RL
646  setring save;
647  Nnew = 2*N+1;
648  // list RL = ringlist(save);  //is defined earlier
649  kill Lord, tmp, iv;
650  L = 0;
651  list Lord, tmp;
652  intvec iv;
653  L[1] = RL[1];
654  L[4] = RL[4];  //char, minpoly
655  // check whether vars have admissible names -> done earlier
656  // list Name = RL[2]
657  // DName is defined earlier
658  tmp[1] = "s";
659  list NName = Name + DName + tmp;
660  L[2] = NName;
661  // dp ordering;
662  string s = "iv=";
663  for (i=1; i<=Nnew; i++)
664  {
665    s = s+"1,";
666  }
667  s[size(s)] = ";";
668  execute(s);
669  kill s;
670  tmp     = 0;
671  tmp[1]  = "dp";  //string
672  tmp[2]  = iv;    //intvec
673  Lord[1] = tmp;
674  tmp[1]  = "C";
675  iv      = 0;
676  tmp[2]  = iv;
677  Lord[2] = tmp;
678  tmp     = 0;
679  L[3]    = Lord;
680  // we are done with the list
681  // Add: Plural part
682  def @R2@ = ring(L);
683  setring @R2@;
684  matrix @D[Nnew][Nnew];
685  for (i=1; i<=N; i++)
686  {
687    @D[i,N+i]=1;
688  }
689  def @R2 = nc_algebra(1,@D);
690  setring @R2;
691  kill @R2@;
692  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready");
693  dbprint(ppl-1, @R2);
694  ideal K = imap(@R,K);
695  option(redSB);
696  //dbprint(ppl,"// -2-2- the final cosmetic std");
697  //K = engine(K,eng);  //std does the job too
698  // total cleanup
699  kill @R;
700  ideal LD = K;
701  export LD;
702  return(@R2);
703}
704example
705{
706  "EXAMPLE:"; echo = 2;
707  ring r = 0,(x,y,z,w),Dp;
708  poly F = x^3+y^3+z^2*w;
709  printlevel = 0;
710  def A = ALTannfsBM(F);
711  setring A;
712  LD;
713}
714
715proc bernsteinBM(poly F, list #)
716"USAGE:  bernsteinBM(f [,eng]);  f a poly, eng an optional int
717RETURN:  list of roots of the Bernstein polynomial b and its multiplicies
718PURPOSE: compute the global Bernstein-Sato polynomial for a hypersurface, defined by f, according to the algorithm by Briancon and Maisonobe
719NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
720@*       otherwise, and by default @code{slimgb} is used.
721@*       If printlevel=1, progress debug messages will be printed,
722@*       if printlevel>=2, all the debug messages will be printed.
723EXAMPLE: example bernsteinBM; shows examples
724"
725{
726  int eng = 0;
727  if ( size(#)>0 )
728  {
729    if ( typeof(#[1]) == "int" )
730    {
731      eng = int(#[1]);
732    }
733  }
734  // returns a list with a ring and an ideal LD in it
735  int ppl = printlevel-voice+2;
736  //  printf("plevel :%s, voice: %s",printlevel,voice);
737  def save = basering;
738  int N = nvars(basering);
739  int Nnew = 2*N+2;
740  int i,j;
741  string s;
742  list RL = ringlist(basering);
743  list L, Lord;
744  list tmp;
745  intvec iv;
746  L[1] = RL[1];  //char
747  L[4] = RL[4];  //char, minpoly
748  // check whether vars have admissible names
749  list Name  = RL[2];
750  list RName;
751  RName[1] = "t";
752  RName[2] = "s";
753  for (i=1; i<=N; i++)
754  {
755    for(j=1; j<=size(RName); j++)
756    {
757      if (Name[i] == RName[j])
758      {
759        ERROR("Variable names should not include t,s");
760      }
761    }
762  }
763  // now, create the names for new vars
764  list DName;
765  for (i=1; i<=N; i++)
766  {
767    DName[i] = "D"+Name[i];  //concat
768  }
769  tmp[1] = "t";
770  tmp[2] = "s";
771  list NName = tmp + Name + DName;
772  L[2]   = NName;
773  // Name, Dname will be used further
774  kill NName;
775  // block ord (lp(2),dp);
776  tmp[1]  = "lp"; // string
777  iv      = 1,1;
778  tmp[2]  = iv; //intvec
779  Lord[1] = tmp;
780  // continue with dp 1,1,1,1...
781  tmp[1]  = "dp"; // string
782  s       = "iv=";
783  for (i=1; i<=Nnew; i++)
784  {
785    s = s+"1,";
786  }
787  s[size(s)]= ";";
788  execute(s);
789  kill s;
790  tmp[2]    = iv;
791  Lord[2]   = tmp;
792  tmp[1]    = "C";
793  iv        = 0;
794  tmp[2]    = iv;
795  Lord[3]   = tmp;
796  tmp       = 0;
797  L[3]      = Lord;
798  // we are done with the list
799  def @R@ = ring(L);
800  setring @R@;
801  matrix @D[Nnew][Nnew];
802  @D[1,2]=t;
803  for(i=1; i<=N; i++)
804  {
805    @D[2+i,N+2+i]=1;
806  }
807  // L[5] = matrix(UpOneMatrix(Nnew));
808  // L[6] = @D;
809  def @R = nc_algebra(1,@D);
810  setring @R;
811  kill @R@;
812  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
813  dbprint(ppl-1, @R);
814  // create the ideal I
815  poly  F = imap(save,F);
816  ideal I = t*F+s;
817  poly p;
818  for(i=1; i<=N; i++)
819  {
820    p = t;  //t
821    p = diff(F,var(2+i))*p;
822    I = I, var(N+2+i) + p;
823  }
824  // -------- the ideal I is ready ----------
825  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
826  dbprint(ppl-1, I);
827  ideal J = engine(I,eng);
828  ideal K = nselect(J,1);
829  kill I,J;
830  dbprint(ppl,"// -1-3- t is eliminated");
831  dbprint(ppl-1, K);  //K is without t
832  // ----------- the ring @R2 ------------
833  // _x, _Dx,s;  elim.ord for _x,_Dx.
834  // keep: N, i,j,s, tmp, RL
835  setring save;
836  Nnew = 2*N+1;
837  kill Lord, tmp, iv, RName;
838  list Lord, tmp;
839  intvec iv;
840  L[1] = RL[1];
841  L[4] = RL[4];  //char, minpoly
842  // check whether vars hava admissible names -> done earlier
843  // now, create the names for new var
844  tmp[1] = "s";
845  // DName is defined earlier
846  list NName = Name + DName + tmp;
847  L[2] = NName;
848  tmp = 0;
849  // block ord (dp(N),dp);
850  string s = "iv=";
851  for (i=1; i<=Nnew-1; i++)
852  {
853    s = s+"1,";
854  }
855  s[size(s)]=";";
856  execute(s);
857  tmp[1] = "dp";  //string
858  tmp[2] = iv;    //intvec
859  Lord[1] = tmp;
860  // continue with dp 1,1,1,1...
861  tmp[1] = "dp";  //string
862  s[size(s)] = ",";
863  s = s+"1;";
864  execute(s);
865  kill s;
866  kill NName;
867  tmp[2]  = iv;
868  Lord[2] = tmp;
869  tmp[1]  = "C";
870  iv      = 0;
871  tmp[2]  = iv;
872  Lord[3] = tmp;
873  tmp     = 0;
874  L[3]    = Lord;
875  // we are done with the list. Now add a Plural part
876  def @R2@ = ring(L);
877  setring @R2@;
878  matrix @D[Nnew][Nnew];
879  for (i=1; i<=N; i++)
880  {
881    @D[i,N+i]=1;
882  }
883  def @R2 = nc_algebra(1,@D);
884  setring @R2;
885  kill @R2@;
886  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready");
887  dbprint(ppl-1, @R2);
888  ideal MM = maxideal(1);
889  MM = 0,s,MM;
890  map R01 = @R, MM;
891  ideal K = R01(K);
892  kill @R, R01;
893  poly F = imap(save,F);
894  K = K,F;
895  dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2");
896  dbprint(ppl-1, K);
897  ideal M = engine(K,eng);
898  ideal K2 = nselect(M,1,Nnew-1);
899  kill K,M;
900  dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2");
901  dbprint(ppl-1, K2);
902  // the ring @R3 and the search for minimal negative int s
903  ring @R3 = 0,s,dp;
904  dbprint(ppl,"// -3-1- the ring @R3(s) is ready");
905  ideal K3 = imap(@R2,K2);
906  kill @R2;
907  poly p = K3[1];
908  dbprint(ppl,"// -3-2- factorization");
909  list P = factorize(p);          //with constants and multiplicities
910  ideal bs; intvec m;             //the Bernstein polynomial is monic, so we are not interested in constants
911  for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
912  {
913    bs[i-1] = P[1][i];
914    m[i-1]  = P[2][i];
915  }
916  // "--------- b-function factorizes into ---------";  P;
917  //int sP = minIntRoot(bs,1);
918  //dbprint(ppl,"// -3-3- minimal integer root found");
919  //dbprint(ppl-1, sP);
920  // convert factors to a list of their roots and multiplicities
921  bs =  normalize(bs);
922  bs = -subst(bs,s,0);
923  setring save;
924  ideal bs = imap(@R3,bs);
925  kill @R3;
926  list BS = bs,m;
927  return(BS);
928}
929example
930{
931  "EXAMPLE:"; echo = 2;
932  ring r = 0,(x,y,z,w),Dp;
933  poly F = x^3+y^3+z^2*w;
934  printlevel = 0;
935  bernsteinBM(F);
936}
937
938// some changes
939proc annfsBM (poly F, list #)
940"USAGE:  annfsBM(f [,eng]);  f a poly, eng an optional int
941RETURN:  ring
942PURPOSE: compute the D-module structure of basering[1/f]*f^s, according
943to the algorithm by Briancon and Maisonobe
944NOTE:    activate this ring with the @code{setring} command. In this ring,
945@*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
946@*         which is obtained by substituting the minimal integer root of a Bernstein
947@*         polynomial into the s-parametric ideal;
948@*       - the list BS is the list of roots and multiplicities of a Bernstein polynomial of f.
949@*       If eng <>0, @code{std} is used for Groebner basis computations,
950@*       otherwise, and by default @code{slimgb} is used.
951@*       If printlevel=1, progress debug messages will be printed,
952@*       if printlevel>=2, all the debug messages will be printed.
953EXAMPLE: example annfsBM; shows examples
954"
955{
956  int eng = 0;
957  if ( size(#)>0 )
958  {
959    if ( typeof(#[1]) == "int" )
960    {
961      eng = int(#[1]);
962    }
963  }
964  // returns a list with a ring and an ideal LD in it
965  int ppl = printlevel-voice+2;
966  //  printf("plevel :%s, voice: %s",printlevel,voice);
967  def save = basering;
968  int N = nvars(basering);
969  int Nnew = 2*N+2;
970  int i,j;
971  string s;
972  list RL = ringlist(basering);
973  list L, Lord;
974  list tmp;
975  intvec iv;
976  L[1] = RL[1];  //char
977  L[4] = RL[4];  //char, minpoly
978  // check whether vars have admissible names
979  list Name  = RL[2];
980  list RName;
981  RName[1] = "t";
982  RName[2] = "s";
983  for (i=1; i<=N; i++)
984  {
985    for(j=1; j<=size(RName); j++)
986    {
987      if (Name[i] == RName[j])
988      {
989        ERROR("Variable names should not include t,s");
990      }
991    }
992  }
993  // now, create the names for new vars
994  list DName;
995  for (i=1; i<=N; i++)
996  {
997    DName[i] = "D"+Name[i];  //concat
998  }
999  tmp[1] = "t";
1000  tmp[2] = "s";
1001  list NName = tmp + Name + DName;
1002  L[2]   = NName;
1003  // Name, Dname will be used further
1004  kill NName;
1005  // block ord (lp(2),dp);
1006  tmp[1]  = "lp"; // string
1007  iv      = 1,1;
1008  tmp[2]  = iv; //intvec
1009  Lord[1] = tmp;
1010  // continue with dp 1,1,1,1...
1011  tmp[1]  = "dp"; // string
1012  s       = "iv=";
1013  for (i=1; i<=Nnew; i++)
1014  {
1015    s = s+"1,";
1016  }
1017  s[size(s)]= ";";
1018  execute(s);
1019  kill s;
1020  tmp[2]    = iv;
1021  Lord[2]   = tmp;
1022  tmp[1]    = "C";
1023  iv        = 0;
1024  tmp[2]    = iv;
1025  Lord[3]   = tmp;
1026  tmp       = 0;
1027  L[3]      = Lord;
1028  // we are done with the list
1029  def @R@ = ring(L);
1030  setring @R@;
1031  matrix @D[Nnew][Nnew];
1032  @D[1,2]=t;
1033  for(i=1; i<=N; i++)
1034  {
1035    @D[2+i,N+2+i]=1;
1036  }
1037  // L[5] = matrix(UpOneMatrix(Nnew));
1038  // L[6] = @D;
1039  def @R = nc_algebra(1,@D);
1040  setring @R;
1041  kill @R@;
1042  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
1043  dbprint(ppl-1, @R);
1044  // create the ideal I
1045  poly  F = imap(save,F);
1046  ideal I = t*F+s;
1047  poly p;
1048  for(i=1; i<=N; i++)
1049  {
1050    p = t;  //t
1051    p = diff(F,var(2+i))*p;
1052    I = I, var(N+2+i) + p;
1053  }
1054  // -------- the ideal I is ready ----------
1055  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
1056  dbprint(ppl-1, I);
1057  ideal J = engine(I,eng);
1058  ideal K = nselect(J,1);
1059  kill I,J;
1060  dbprint(ppl,"// -1-3- t is eliminated");
1061  dbprint(ppl-1, K);  //K is without t
1062  setring save;
1063  // ----------- the ring @R2 ------------
1064  // _x, _Dx,s;  elim.ord for _x,_Dx.
1065  // keep: N, i,j,s, tmp, RL
1066  Nnew = 2*N+1;
1067  kill Lord, tmp, iv, RName;
1068  list Lord, tmp;
1069  intvec iv;
1070  L[1] = RL[1];
1071  L[4] = RL[4];  //char, minpoly
1072  // check whether vars hava admissible names -> done earlier
1073  // now, create the names for new var
1074  tmp[1] = "s";
1075  // DName is defined earlier
1076  list NName = Name + DName + tmp;
1077  L[2] = NName;
1078  tmp = 0;
1079  // block ord (dp(N),dp);
1080  string s = "iv=";
1081  for (i=1; i<=Nnew-1; i++)
1082  {
1083    s = s+"1,";
1084  }
1085  s[size(s)]=";";
1086  execute(s);
1087  tmp[1] = "dp";  //string
1088  tmp[2] = iv;    //intvec
1089  Lord[1] = tmp;
1090  // continue with dp 1,1,1,1...
1091  tmp[1] = "dp";  //string
1092  s[size(s)] = ",";
1093  s = s+"1;";
1094  execute(s);
1095  kill s;
1096  kill NName;
1097  tmp[2]  = iv;
1098  Lord[2] = tmp;
1099  tmp[1]  = "C";
1100  iv      = 0;
1101  tmp[2]  = iv;
1102  Lord[3] = tmp;
1103  tmp     = 0;
1104  L[3]    = Lord;
1105  // we are done with the list. Now add a Plural part
1106  def @R2@ = ring(L);
1107  setring @R2@;
1108  matrix @D[Nnew][Nnew];
1109  for (i=1; i<=N; i++)
1110  {
1111    @D[i,N+i]=1;
1112  }
1113  def @R2 = nc_algebra(1,@D);
1114  setring @R2;
1115  kill @R2@;
1116  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready");
1117  dbprint(ppl-1, @R2);
1118  ideal MM = maxideal(1);
1119  MM = 0,s,MM;
1120  map R01 = @R, MM;
1121  ideal K = R01(K);
1122  poly F = imap(save,F);
1123  K = K,F;
1124  dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2");
1125  dbprint(ppl-1, K);
1126  ideal M = engine(K,eng);
1127  ideal K2 = nselect(M,1,Nnew-1);
1128  kill K,M;
1129  dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2");
1130  dbprint(ppl-1, K2);
1131  // the ring @R3 and the search for minimal negative int s
1132  ring @R3 = 0,s,dp;
1133  dbprint(ppl,"// -3-1- the ring @R3(s) is ready");
1134  ideal K3 = imap(@R2,K2);
1135  poly p = K3[1];
1136  dbprint(ppl,"// -3-2- factorization");
1137  list P = factorize(p);          //with constants and multiplicities
1138  ideal bs; intvec m;             //the Bernstein polynomial is monic, so we are not interested in constants
1139  for (i=2; i<= size(P[1]); i++)  //we ignore P[1][1] (constant) and P[2][1] (its mult.)
1140  {
1141    bs[i-1] = P[1][i];
1142    m[i-1]  = P[2][i];
1143  }
1144  // "--------- b-function factorizes into ---------";  P;
1145  int sP = minIntRoot(bs,1);
1146  dbprint(ppl,"// -3-3- minimal integer root found");
1147  dbprint(ppl-1, sP);
1148  // convert factors to a list of their roots
1149  bs = normalize(bs);
1150  bs = -subst(bs,s,0);
1151  list BS = bs,m;
1152  //TODO: sort BS!
1153  // --------- substitute s found in the ideal ---------
1154  // --------- going back to @R and substitute ---------
1155  setring @R;
1156  ideal K2 = subst(K,s,sP);
1157  kill K;
1158  // create the ordinary Weyl algebra and put the result into it,
1159  // thus creating the ring @R4
1160  // keep: N, i,j,s, tmp, RL
1161  setring save;
1162  Nnew = 2*N;
1163  // list RL = ringlist(save);  //is defined earlier
1164  kill Lord, tmp, iv;
1165  L = 0;
1166  list Lord, tmp;
1167  intvec iv;
1168  L[1] = RL[1];
1169  L[4] = RL[4];  //char, minpoly
1170  // check whether vars have admissible names -> done earlier
1171  // list Name = RL[2]M
1172  // DName is defined earlier
1173  list NName = Name + DName;
1174  L[2] = NName;
1175  // dp ordering;
1176  string s = "iv=";
1177  for (i=1; i<=Nnew; i++)
1178  {
1179    s = s+"1,";
1180  }
1181  s[size(s)] = ";";
1182  execute(s);
1183  kill s;
1184  tmp     = 0;
1185  tmp[1]  = "dp";  //string
1186  tmp[2]  = iv;    //intvec
1187  Lord[1] = tmp;
1188  tmp[1]  = "C";
1189  iv      = 0;
1190  tmp[2]  = iv;
1191  Lord[2] = tmp;
1192  tmp     = 0;
1193  L[3]    = Lord;
1194  // we are done with the list
1195  // Add: Plural part
1196  def @R4@ = ring(L);
1197  setring @R4@;
1198  matrix @D[Nnew][Nnew];
1199  for (i=1; i<=N; i++)
1200  {
1201    @D[i,N+i]=1;
1202  }
1203  def @R4 = nc_algebra(1,@D);
1204  setring @R4;
1205  kill @R4@;
1206  dbprint(ppl,"// -4-1- the ring @R4(_x,_Dx) is ready");
1207  dbprint(ppl-1, @R4);
1208  ideal K4 = imap(@R,K2);
1209  option(redSB);
1210  dbprint(ppl,"// -4-2- the final cosmetic std");
1211  K4 = engine(K4,eng);  //std does the job too
1212  // total cleanup
1213  kill @R;
1214  kill @R2;
1215  list BS = imap(@R3,BS);
1216  export BS;
1217  kill @R3;
1218  ideal LD = K4;
1219  export LD;
1220  return(@R4);
1221}
1222example
1223{
1224  "EXAMPLE:"; echo = 2;
1225  ring r = 0,(x,y,z),Dp;
1226  poly F = z*x^2+y^3;
1227  printlevel = 0;
1228  def A = annfsBM(F);
1229  setring A;
1230  LD;
1231  BS;
1232}
1233
1234proc operatorBM(poly F, list #)
1235"USAGE:  operatorBM(f [,eng]);  f a poly, eng an optional int
1236RETURN:  ring
1237PURPOSE: compute the B-operator and other relevant data for Ann F^s, according to the algorithm by Briancon and Maisonobe
1238NOTE:    activate this ring with the @code{setring} command. In this ring D[s]
1239@*       - the polynomial F is the same as the input,
1240@*       - the ideal LD is the annihilator of f^s in Dn[s],
1241@*       - the ideal LD0 is the needed D-mod structure, where LD0 = LD|s=s0,
1242@*       - the polynomial bs is the global Bernstein polynomial of f in the variable s,
1243@*       - the list BS contains all the roots with multiplicities of the global Bernstein polynomial of f,
1244@*       - the polynomial PS is an operator in Dn[s] such that PS*f^(s+1) = bs*f^s.
1245@*       If eng <>0, @code{std} is used for Groebner basis computations,
1246@*       otherwise and by default @code{slimgb} is used.
1247@*       If printlevel=1, progress debug messages will be printed,
1248@*       if printlevel>=2, all the debug messages will be printed.
1249EXAMPLE: example operatorBM; shows examples
1250"
1251{
1252  int eng = 0;
1253  if ( size(#)>0 )
1254  {
1255    if ( typeof(#[1]) == "int" )
1256    {
1257      eng = int(#[1]);
1258    }
1259  }
1260  // returns a list with a ring and an ideal LD in it
1261  int ppl = printlevel-voice+2;
1262  //  printf("plevel :%s, voice: %s",printlevel,voice);
1263  def save = basering;
1264  int N = nvars(basering);
1265  int Nnew = 2*N+2;
1266  int i,j;
1267  string s;
1268  list RL = ringlist(basering);
1269  list L, Lord;
1270  list tmp;
1271  intvec iv;
1272  L[1] = RL[1];  //char
1273  L[4] = RL[4];  //char, minpoly
1274  // check whether vars have admissible names
1275  list Name  = RL[2];
1276  list RName;
1277  RName[1] = "t";
1278  RName[2] = "s";
1279  for (i=1; i<=N; i++)
1280  {
1281    for(j=1; j<=size(RName); j++)
1282    {
1283      if (Name[i] == RName[j])
1284      {
1285        ERROR("Variable names should not include t,s");
1286      }
1287    }
1288  }
1289  // now, create the names for new vars
1290  list DName;
1291  for (i=1; i<=N; i++)
1292  {
1293    DName[i] = "D"+Name[i];  //concat
1294  }
1295  tmp[1] = "t";
1296  tmp[2] = "s";
1297  list NName = tmp + Name + DName;
1298  L[2]   = NName;
1299  // Name, Dname will be used further
1300  kill NName;
1301  // block ord (lp(2),dp);
1302  tmp[1]  = "lp"; // string
1303  iv      = 1,1;
1304  tmp[2]  = iv; //intvec
1305  Lord[1] = tmp;
1306  // continue with dp 1,1,1,1...
1307  tmp[1]  = "dp"; // string
1308  s       = "iv=";
1309  for (i=1; i<=Nnew; i++)
1310  {
1311    s = s+"1,";
1312  }
1313  s[size(s)]= ";";
1314  execute(s);
1315  kill s;
1316  tmp[2]    = iv;
1317  Lord[2]   = tmp;
1318  tmp[1]    = "C";
1319  iv        = 0;
1320  tmp[2]    = iv;
1321  Lord[3]   = tmp;
1322  tmp       = 0;
1323  L[3]      = Lord;
1324  // we are done with the list
1325  def @R@ = ring(L);
1326  setring @R@;
1327  matrix @D[Nnew][Nnew];
1328  @D[1,2]=t;
1329  for(i=1; i<=N; i++)
1330  {
1331    @D[2+i,N+2+i]=1;
1332  }
1333  // L[5] = matrix(UpOneMatrix(Nnew));
1334  // L[6] = @D;
1335  def @R = nc_algebra(1,@D);
1336  setring @R;
1337  kill @R@;
1338  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
1339  dbprint(ppl-1, @R);
1340  // create the ideal I
1341  poly  F = imap(save,F);
1342  ideal I = t*F+s;
1343  poly p;
1344  for(i=1; i<=N; i++)
1345  {
1346    p = t;  //t
1347    p = diff(F,var(2+i))*p;
1348    I = I, var(N+2+i) + p;
1349  }
1350  // -------- the ideal I is ready ----------
1351  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
1352  dbprint(ppl-1, I);
1353  ideal J = engine(I,eng);
1354  ideal K = nselect(J,1);
1355  kill I,J;
1356  dbprint(ppl,"// -1-3- t is eliminated");
1357  dbprint(ppl-1, K);  //K is without t
1358  setring save;
1359  // ----------- the ring @R2 ------------
1360  // _x, _Dx,s;  elim.ord for _x,_Dx.
1361  // keep: N, i,j,s, tmp, RL
1362  Nnew = 2*N+1;
1363  kill Lord, tmp, iv, RName;
1364  list Lord, tmp;
1365  intvec iv;
1366  L[1] = RL[1];
1367  L[4] = RL[4];  //char, minpoly
1368  // check whether vars hava admissible names -> done earlier
1369  // now, create the names for new var
1370  tmp[1] = "s";
1371  // DName is defined earlier
1372  list NName = Name + DName + tmp;
1373  L[2] = NName;
1374  tmp = 0;
1375  // block ord (dp(N),dp);
1376  string s = "iv=";
1377  for (i=1; i<=Nnew-1; i++)
1378  {
1379    s = s+"1,";
1380  }
1381  s[size(s)]=";";
1382  execute(s);
1383  tmp[1] = "dp";  //string
1384  tmp[2] = iv;    //intvec
1385  Lord[1] = tmp;
1386  // continue with dp 1,1,1,1...
1387  tmp[1] = "dp";  //string
1388  s[size(s)] = ",";
1389  s = s+"1;";
1390  execute(s);
1391  kill s;
1392  kill NName;
1393  tmp[2]  = iv;
1394  Lord[2] = tmp;
1395  tmp[1]  = "C";
1396  iv      = 0;
1397  tmp[2]  = iv;
1398  Lord[3] = tmp;
1399  tmp     = 0;
1400  L[3]    = Lord;
1401  // we are done with the list. Now add a Plural part
1402  def @R2@ = ring(L);
1403  setring @R2@;
1404  matrix @D[Nnew][Nnew];
1405  for (i=1; i<=N; i++)
1406  {
1407    @D[i,N+i]=1;
1408  }
1409  def @R2 = nc_algebra(1,@D);
1410  setring @R2;
1411  kill @R2@;
1412  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,s) is ready");
1413  dbprint(ppl-1, @R2);
1414  ideal MM = maxideal(1);
1415  MM = 0,s,MM;
1416  map R01 = @R, MM;
1417  ideal K = R01(K);
1418  poly F = imap(save,F);
1419  K = K,F;
1420  dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2");
1421  dbprint(ppl-1, K);
1422  ideal M = engine(K,eng);
1423  ideal K2 = nselect(M,1,Nnew-1);
1424  kill K,M;
1425  dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2");
1426  dbprint(ppl-1, K2);
1427  // the ring @R3 and the search for minimal negative int s
1428  ring @R3 = 0,s,dp;
1429  dbprint(ppl,"// -3-1- the ring @R3(s) is ready");
1430  ideal K3 = imap(@R2,K2);
1431  kill @R2;
1432  poly p = K3[1];
1433  dbprint(ppl,"// -3-2- factorization");
1434  list P = factorize(p);          //with constants and multiplicities
1435  ideal bs; intvec m;             //the Bernstein polynomial is monic, so we are not interested in constants
1436  for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
1437  {
1438    bs[i-1] = P[1][i];
1439    m[i-1]  = P[2][i];
1440  }
1441  // "--------- b-function factorizes into ---------";  P;
1442  int sP = minIntRoot(bs,1);
1443  dbprint(ppl,"// -3-3- minimal integer root found");
1444  dbprint(ppl-1, sP);
1445  // convert factors to a list of their roots with multiplicities
1446  bs = normalize(bs);
1447  bs = -subst(bs,s,0);
1448  list BS = bs,m;
1449  //TODO: sort BS!
1450  // --------- substitute s found in the ideal ---------
1451  // --------- going back to @R and substitute ---------
1452  setring @R;
1453  ideal K2 = subst(K,s,sP);
1454  // create Dn[s], where Dn is the ordinary Weyl algebra, and put the result into it,
1455  // thus creating the ring @R4
1456  // keep: N, i,j,s, tmp, RL
1457  setring save;
1458  Nnew = 2*N+1;
1459  // list RL = ringlist(save);  //is defined earlier
1460  kill Lord, tmp, iv;
1461  L = 0;
1462  list Lord, tmp;
1463  intvec iv;
1464  L[1] = RL[1];
1465  L[4] = RL[4];  //char, minpoly
1466  // check whether vars have admissible names -> done earlier
1467  // list Name = RL[2]
1468  // DName is defined earlier
1469  tmp[1] = "s";
1470  list NName = Name + DName + tmp;
1471  L[2] = NName;
1472  // dp ordering;
1473  string s = "iv=";
1474  for (i=1; i<=Nnew; i++)
1475  {
1476    s = s+"1,";
1477  }
1478  s[size(s)] = ";";
1479  execute(s);
1480  kill s;
1481  tmp     = 0;
1482  tmp[1]  = "dp";  //string
1483  tmp[2]  = iv;    //intvec
1484  Lord[1] = tmp;
1485  tmp[1]  = "C";
1486  iv      = 0;
1487  tmp[2]  = iv;
1488  Lord[2] = tmp;
1489  tmp     = 0;
1490  L[3]    = Lord;
1491  // we are done with the list
1492  // Add: Plural part
1493  def @R4@ = ring(L);
1494  setring @R4@;
1495  matrix @D[Nnew][Nnew];
1496  for (i=1; i<=N; i++)
1497  {
1498    @D[i,N+i]=1;
1499  }
1500  def @R4 = nc_algebra(1,@D);
1501  setring @R4;
1502  kill @R4@;
1503  dbprint(ppl,"// -4-1- the ring @R4(_x,_Dx,s) is ready");
1504  dbprint(ppl-1, @R4);
1505  ideal LD0 = imap(@R,K2);
1506  ideal LD  = imap(@R,K);
1507  kill @R;
1508  poly bs = imap(@R3,p);
1509  list BS = imap(@R3,BS);
1510  kill @R3;
1511  bs = normalize(bs);
1512  poly F = imap(save,F);
1513  dbprint(ppl,"// -4-2- starting the computation of PS via lift");
1514//better liftstd, I didn't knot it works also for Plural, liftslimgb?
1515// liftstd may give extra coeffs in the resulting ideal
1516  matrix T = lift(F+LD,bs);
1517  poly PS = T[1,1];
1518  dbprint(ppl,"// -4-3- an operator PS found, PS*f^(s+1) = b(s)*f^s");
1519  dbprint(ppl-1,PS);
1520  option(redSB);
1521  dbprint(ppl,"// -4-4- the final cosmetic std");
1522  LD0 = engine(LD0,eng);  //std does the job too
1523  LD  = engine(LD,eng);
1524  export F,LD,LD0,bs,BS,PS;
1525  return(@R4);
1526}
1527example
1528{
1529  "EXAMPLE:"; echo = 2;
1530  //  ring r = 0,(x,y,z,w),Dp;
1531  //  poly F = x^3+y^3+z^2*w;
1532  ring r = 0,(x,y,z),Dp;
1533  poly F = x^3+y^3+z^3;
1534  printlevel = 0;
1535  def A = operatorBM(F);
1536  setring A;
1537  F; // the original polynomial itself
1538  LD; // generic annihilator
1539  LD0; // annihilator
1540  bs; // normalized Bernstein poly
1541  BS; // root and multiplicities of the Bernstein poly
1542  PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD
1543  reduce(PS*F-bs,LD); // check the property of PS
1544}
1545
1546proc annfsParamBM (poly F, list #)
1547"USAGE:  annfsParamBM(f [,eng]);  f a poly, eng an optional int
1548RETURN:  ring
1549PURPOSE: compute the generic Ann F^s and exceptional parametric constellations of a polynomial with parametric coefficients, according to the algorithm by Briancon and Maisonobe
1550NOTE:    activate this ring with the @code{setring} command. In this ring,
1551@*       - the ideal LD is the D-module structure oa Ann F^s
1552@*       - the ideal Param contains the list of the special parameters.
1553@*       If eng <>0, @code{std} is used for Groebner basis computations,
1554@*       otherwise, and by default @code{slimgb} is used.
1555@*       If printlevel=1, progress debug messages will be printed,
1556@*       if printlevel>=2, all the debug messages will be printed.
1557EXAMPLE: example annfsParamBM; shows examples
1558"
1559{
1560  //PURPOSE: compute the list of all possible Bernstein-Sato polynomials for a polynomial with parametric coefficients, according to the algorithm by Briancon and Maisonobe
1561  // @*       - the list BS is the list of roots and multiplicities of a Bernstein polynomial of f.
1562  // ***** not implented yet ****
1563  int eng = 0;
1564  if ( size(#)>0 )
1565  {
1566    if ( typeof(#[1]) == "int" )
1567    {
1568      eng = int(#[1]);
1569    }
1570  }
1571  // returns a list with a ring and an ideal LD in it
1572  int ppl = printlevel-voice+2;
1573  //  printf("plevel :%s, voice: %s",printlevel,voice);
1574  def save = basering;
1575  int N = nvars(basering);
1576  int Nnew = 2*N+2;
1577  int i,j;
1578  string s;
1579  list RL = ringlist(basering);
1580  list L, Lord;
1581  list tmp;
1582  intvec iv;
1583  L[1] = RL[1];  //char
1584  L[4] = RL[4];  //char, minpoly
1585  // check whether vars have admissible names
1586  list Name  = RL[2];
1587  list RName;
1588  RName[1] = "t";
1589  RName[2] = "s";
1590  for (i=1; i<=N; i++)
1591  {
1592    for(j=1; j<=size(RName); j++)
1593    {
1594      if (Name[i] == RName[j])
1595      {
1596        ERROR("Variable names should not include t,s");
1597      }
1598    }
1599  }
1600  // now, create the names for new vars
1601  list DName;
1602  for (i=1; i<=N; i++)
1603  {
1604    DName[i] = "D"+Name[i];  //concat
1605  }
1606  tmp[1] = "t";
1607  tmp[2] = "s";
1608  list NName = tmp + Name + DName;
1609  L[2]   = NName;
1610  // Name, Dname will be used further
1611  kill NName;
1612  // block ord (lp(2),dp);
1613  tmp[1]  = "lp"; // string
1614  iv      = 1,1;
1615  tmp[2]  = iv; //intvec
1616  Lord[1] = tmp;
1617  // continue with dp 1,1,1,1...
1618  tmp[1]  = "dp"; // string
1619  s       = "iv=";
1620  for (i=1; i<=Nnew; i++)
1621  {
1622    s = s+"1,";
1623  }
1624  s[size(s)]= ";";
1625  execute(s);
1626  kill s;
1627  tmp[2]    = iv;
1628  Lord[2]   = tmp;
1629  tmp[1]    = "C";
1630  iv        = 0;
1631  tmp[2]    = iv;
1632  Lord[3]   = tmp;
1633  tmp       = 0;
1634  L[3]      = Lord;
1635  // we are done with the list
1636  def @R@ = ring(L);
1637  setring @R@;
1638  matrix @D[Nnew][Nnew];
1639  @D[1,2]=t;
1640  for(i=1; i<=N; i++)
1641  {
1642    @D[2+i,N+2+i]=1;
1643  }
1644  // L[5] = matrix(UpOneMatrix(Nnew));
1645  // L[6] = @D;
1646  def @R = nc_algebra(1,@D);
1647  setring @R;
1648  kill @R@;
1649  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
1650  dbprint(ppl-1, @R);
1651  // create the ideal I
1652  poly  F = imap(save,F);
1653  ideal I = t*F+s;
1654  poly p;
1655  for(i=1; i<=N; i++)
1656  {
1657    p = t;  //t
1658    p = diff(F,var(2+i))*p;
1659    I = I, var(N+2+i) + p;
1660  }
1661  // -------- the ideal I is ready ----------
1662  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
1663  dbprint(ppl-1, I);
1664  ideal J = engine(I,eng);
1665  ideal K = nselect(J,1);
1666  dbprint(ppl,"// -1-3- t is eliminated");
1667  dbprint(ppl-1, K);  //K is without t
1668  // ----- looking for special parameters -----
1669  dbprint(ppl,"// -2-1- starting the computation of the transformation matrix (via lift)");
1670  J = normalize(J);
1671  matrix T = lift(I,J);  //try also with liftstd
1672  kill I,J;
1673  dbprint(ppl,"// -2-2- the transformation matrix has been computed");
1674  dbprint(ppl-1, T);  //T is the transformation matrix
1675  dbprint(ppl,"// -2-3- genericity does the job");
1676  list lParam = genericity(T);
1677  int ip = size(lParam);
1678  int cip;
1679  string sParam;
1680  if (sParam[1]=="-") { sParam=""; } //genericity returns "-"
1681  // if no parameters exist in a basering
1682  for (cip=1; cip <= ip; cip++)
1683  {
1684    sParam = sParam + "," +lParam[cip];
1685  }
1686  if (size(sParam) >=2)
1687  {
1688    sParam = sParam[2..size(sParam)]; // removes the 1st colon
1689  }
1690  export sParam;
1691  kill T;
1692  dbprint(ppl,"// -2-4- the special parameters has been computed");
1693  dbprint(ppl, sParam);
1694  // create Dn[s], where Dn is the ordinary Weyl Algebra, and put the result into it,
1695  // thus creating the ring @R2
1696  // keep: N, i,j,s, tmp, RL
1697  setring save;
1698  Nnew = 2*N+1;
1699  // list RL = ringlist(save);  //is defined earlier
1700  kill Lord, tmp, iv;
1701  L = 0;
1702  list Lord, tmp;
1703  intvec iv;
1704  L[1] = RL[1];
1705  L[4] = RL[4];  //char, minpoly
1706  // check whether vars have admissible names -> done earlier
1707  // list Name = RL[2]M
1708  // DName is defined earlier
1709  tmp[1] = "s";
1710  list NName = Name + DName + tmp;
1711  L[2] = NName;
1712  // dp ordering;
1713  string s = "iv=";
1714  for (i=1; i<=Nnew; i++)
1715  {
1716    s = s+"1,";
1717  }
1718  s[size(s)] = ";";
1719  execute(s);
1720  kill s;
1721  tmp     = 0;
1722  tmp[1]  = "dp";  //string
1723  tmp[2]  = iv;    //intvec
1724  Lord[1] = tmp;
1725  tmp[1]  = "C";
1726  iv      = 0;
1727  tmp[2]  = iv;
1728  Lord[2] = tmp;
1729  tmp     = 0;
1730  L[3]    = Lord;
1731  // we are done with the list
1732  // Add: Plural part
1733  def @R2@ = ring(L);
1734  setring @R2@;
1735  matrix @D[Nnew][Nnew];
1736  for (i=1; i<=N; i++)
1737  {
1738    @D[i,N+i]=1;
1739  }
1740  def @R2 = nc_algebra(1,@D);
1741  setring @R2;
1742  kill @R2@;
1743  dbprint(ppl,"// -3-1- the ring @R2(_x,_Dx,s) is ready");
1744  dbprint(ppl-1, @R2);
1745  ideal K = imap(@R,K);
1746  kill @R;
1747  option(redSB);
1748  dbprint(ppl,"// -3-2- the final cosmetic std");
1749  K = engine(K,eng);  //std does the job too
1750  ideal LD = K;
1751  export LD;
1752  if (sParam[1] == ",")
1753  {
1754    sParam = sParam[2..size(sParam)];
1755  }
1756  //  || ((sParam[1] == " ") && (sParam[2] == ",")))
1757  execute("ideal Param ="+sParam+";");
1758  export Param;
1759  kill sParam;
1760  return(@R2);
1761}
1762example
1763{
1764  "EXAMPLE:"; echo = 2;
1765  ring r = (0,a,b),(x,y),Dp;
1766  poly F = x^2 - (y-a)*(y-b);
1767  printlevel = 0;
1768  def A = annfsParamBM(F);  setring A;
1769  LD;
1770  Param;
1771  setring r;
1772  poly G = x2-(y-a)^2; // try the exceptional value b=a of parameters
1773  def B = annfsParamBM(G); setring B;
1774  LD;
1775  Param;
1776}
1777
1778// *** the following example is nice, but too complicated for the documentation ***
1779//   ring r = (0,a),(x,y,z),Dp;
1780//   poly F = x^4+y^4+z^2+a*x*y*z;
1781//   printlevel = 2; //0
1782//   def A = annfsParamBM(F);
1783//   setring A;
1784//   LD;
1785//   Param;
1786
1787
1788proc annfsBMI(ideal F, list #)
1789"USAGE:  annfsBMI(F [,eng]);  F an ideal, eng an optional int
1790RETURN:  ring
1791PURPOSE: compute the D-module structure of basering[1/f]*f^s where f = F[1]*..*F[P],
1792according to the algorithm by Briancon and Maisonobe.
1793NOTE:    activate this ring with the @code{setring} command. In this ring,
1794@*       - the ideal LD is the needed D-mod structure,
1795@*       - the list BS is the Bernstein ideal of a polynomial f = F[1]*..*F[P].
1796@*       If eng <>0, @code{std} is used for Groebner basis computations,
1797@*       otherwise, and by default @code{slimgb} is used.
1798@*       If printlevel=1, progress debug messages will be printed,
1799@*       if printlevel>=2, all the debug messages will be printed.
1800EXAMPLE: example annfsBMI; shows examples
1801"
1802{
1803  int eng = 0;
1804  if ( size(#)>0 )
1805  {
1806    if ( typeof(#[1]) == "int" )
1807    {
1808      eng = int(#[1]);
1809    }
1810  }
1811  // returns a list with a ring and an ideal LD in it
1812  int ppl = printlevel-voice+2;
1813  //  printf("plevel :%s, voice: %s",printlevel,voice);
1814  def save = basering;
1815  int N = nvars(basering);
1816  int P = size(F);  //if F has some generators which are zero, int P = ncols(I);
1817  int Nnew = 2*N+2*P;
1818  int i,j;
1819  string s;
1820  list RL = ringlist(basering);
1821  list L, Lord;
1822  list tmp;
1823  intvec iv;
1824  L[1] = RL[1];  //char
1825  L[4] = RL[4];  //char, minpoly
1826  // check whether vars have admissible names
1827  list Name  = RL[2];
1828  list RName;
1829  for (j=1; j<=P; j++)
1830  {
1831    RName[j] = "t("+string(j)+")";
1832    RName[j+P] = "s("+string(j)+")";
1833  }
1834  for(i=1; i<=N; i++)
1835  {
1836    for(j=1; j<=size(RName); j++)
1837    {
1838      if (Name[i] == RName[j])
1839      { ERROR("Variable names should not include t(i),s(i)"); }
1840    }
1841  }
1842  // now, create the names for new vars
1843  list DName;
1844  for(i=1; i<=N; i++)
1845  {
1846    DName[i] = "D"+Name[i];  //concat
1847  }
1848  list NName = RName + Name + DName;
1849  L[2]   = NName;
1850  // Name, Dname will be used further
1851  kill NName;
1852  // block ord (lp(P),dp);
1853  tmp[1] = "lp";  //string
1854  s      = "iv=";
1855  for (i=1; i<=2*P; i++)
1856  {
1857    s = s+"1,";
1858  }
1859  s[size(s)]= ";";
1860  execute(s);
1861  tmp[2] = iv;  //intvec
1862  Lord[1] = tmp;
1863  // continue with dp 1,1,1,1...
1864  tmp[1] = "dp";  //string
1865  s      = "iv=";
1866  for (i=1; i<=Nnew; i++)  //actually i<=2*N
1867  {
1868    s = s+"1,";
1869  }
1870  s[size(s)]= ";";
1871  execute(s);
1872  kill s;
1873  tmp[2]  = iv;
1874  Lord[2] = tmp;
1875  tmp[1]  = "C";
1876  iv      = 0;
1877  tmp[2]  = iv;
1878  Lord[3] = tmp;
1879  tmp     = 0;
1880  L[3]    = Lord;
1881  // we are done with the list
1882  def @R@ = ring(L);
1883  setring @R@;
1884  matrix @D[Nnew][Nnew];
1885  for (i=1; i<=P; i++)
1886  {
1887    @D[i,i+P] = t(i);
1888  }
1889  for(i=1; i<=N; i++)
1890  {
1891    @D[2*P+i,2*P+N+i] = 1;
1892  }
1893  // L[5] = matrix(UpOneMatrix(Nnew));
1894  // L[6] = @D;
1895  def @R = nc_algebra(1,@D);
1896  setring @R;
1897  kill @R@;
1898  dbprint(ppl,"// -1-1- the ring @R(_t,_s,_x,_Dx) is ready");
1899  dbprint(ppl-1, @R);
1900  // create the ideal I
1901  ideal  F = imap(save,F);
1902  ideal I = t(1)*F[1]+s(1);
1903  for (j=2; j<=P; j++)
1904  {
1905    I = I, t(j)*F[j]+s(j);
1906  }
1907  poly p,q;
1908  for (i=1; i<=N; i++)
1909  {
1910    p=0;
1911    for (j=1; j<=P; j++)
1912    {
1913      q = t(j);
1914      q = diff(F[j],var(2*P+i))*q;
1915      p = p + q;
1916    }
1917    I = I, var(2*P+N+i) + p;
1918  }
1919  // -------- the ideal I is ready ----------
1920  dbprint(ppl,"// -1-2- starting the elimination of "+string(t(1..P))+" in @R");
1921  dbprint(ppl-1, I);
1922  ideal J = engine(I,eng);
1923  ideal K = nselect(J,1,P);
1924  kill I,J;
1925  dbprint(ppl,"// -1-3- all t(i) are eliminated");
1926  dbprint(ppl-1, K);  //K is without t(i)
1927  // ----------- the ring @R2 ------------
1928  // _x, _Dx,s;  elim.ord for _x,_Dx.
1929  // keep: N, i,j,s, tmp, RL
1930  setring save;
1931  Nnew = 2*N+P;
1932  kill Lord, tmp, iv, RName;
1933  list Lord, tmp;
1934  intvec iv;
1935  L[1] = RL[1];  //char
1936  L[4] = RL[4];  //char, minpoly
1937  // check whether vars hava admissible names -> done earlier
1938  // now, create the names for new var
1939  for (j=1; j<=P; j++)
1940  {
1941    tmp[j] = "s("+string(j)+")";
1942  }
1943  // DName is defined earlier
1944  list NName = Name + DName + tmp;
1945  L[2] = NName;
1946  tmp = 0;
1947  // block ord (dp(N),dp);
1948  string s = "iv=";
1949  for (i=1; i<=Nnew-P; i++)
1950  {
1951    s = s+"1,";
1952  }
1953  s[size(s)]=";";
1954  execute(s);
1955  tmp[1] = "dp";  //string
1956  tmp[2] = iv;    //intvec
1957  Lord[1] = tmp;
1958  // continue with dp 1,1,1,1...
1959  tmp[1] = "dp";  //string
1960  s[size(s)] = ",";
1961  for (j=1; j<=P; j++)
1962  {
1963    s = s+"1,";
1964  }
1965  s[size(s)]=";";
1966  execute(s);
1967  kill s;
1968  kill NName;
1969  tmp[2]  = iv;
1970  Lord[2] = tmp;
1971  tmp[1]  = "C";
1972  iv      = 0;
1973  tmp[2]  = iv;
1974  Lord[3] = tmp;
1975  tmp     = 0;
1976  L[3]    = Lord;
1977  // we are done with the list. Now add a Plural part
1978  def @R2@ = ring(L);
1979  setring @R2@;
1980  matrix @D[Nnew][Nnew];
1981  for (i=1; i<=N; i++)
1982  {
1983    @D[i,N+i]=1;
1984  }
1985  def @R2 = nc_algebra(1,@D);
1986  setring @R2;
1987  kill @R2@;
1988  dbprint(ppl,"// -2-1- the ring @R2(_x,_Dx,_s) is ready");
1989  dbprint(ppl-1, @R2);
1990//  ideal MM = maxideal(1);
1991//  MM = 0,s,MM;
1992//  map R01 = @R, MM;
1993//  ideal K = R01(K);
1994  ideal F = imap(save,F);  // maybe ideal F = R01(I); ?
1995  ideal K = imap(@R,K);    // maybe ideal K = R01(I); ?
1996  poly f=1;
1997  for (j=1; j<=P; j++)
1998  {
1999    f = f*F[j];
2000  }
2001  K = K,f;       // to compute B (Bernstein-Sato ideal)
2002  //j=2;         // for example
2003  //K = K,F[j];  // to compute Bj (see "On the computation of Bernstein-Sato ideals"; Castro, Ucha)
2004  //K = K,F;     // to compute Bsigma (see "On the computation of Bernstein-Sato ideals"; Castro, Ucha)
2005  dbprint(ppl,"// -2-2- starting the elimination of _x,_Dx in @R2");
2006  dbprint(ppl-1, K);
2007  ideal M = engine(K,eng);
2008  ideal K2 = nselect(M,1,Nnew-P);
2009  kill K,M;
2010  dbprint(ppl,"// -2-3- _x,_Dx are eliminated in @R2");
2011  dbprint(ppl-1, K2);
2012  // the ring @R3 and factorize
2013  ring @R3 = 0,s(1..P),dp;
2014  dbprint(ppl,"// -3-1- the ring @R3(_s) is ready");
2015  ideal K3 = imap(@R2,K2);
2016  if (size(K3)==1)
2017  {
2018    poly p = K3[1];
2019    dbprint(ppl,"// -3-2- factorization");
2020    // Warning: now P is an integer
2021    list Q = factorize(p);         //with constants and multiplicities
2022    ideal bs; intvec m;
2023    for (i=2; i<=size(Q[1]); i++)  //we delete Q[1][1] and Q[2][1]
2024    {
2025      bs[i-1] = Q[1][i];
2026      m[i-1]  = Q[2][i];
2027    }
2028    //  "--------- Q-ideal factorizes into ---------";  list(bs,m);
2029    list BS = bs,m;
2030  }
2031  else
2032  {
2033    // conjecture: the Bernstein ideal is principal
2034    dbprint(ppl,"// -3-2- the Bernstein ideal is not principal");
2035    ideal BS = K3;
2036  }
2037  // create the ring @R4(_x,_Dx,_s) and put the result into it,
2038  // _x, _Dx,s;  ord "dp".
2039  // keep: N, i,j,s, tmp, RL
2040  setring save;
2041  Nnew = 2*N+P;
2042  // list RL = ringlist(save);  //is defined earlier
2043  kill Lord, tmp, iv;
2044  L = 0;
2045  list Lord, tmp;
2046  intvec iv;
2047  L[1] = RL[1];  //char
2048  L[4] = RL[4];  //char, minpoly
2049  // check whether vars hava admissible names -> done earlier
2050  // now, create the names for new var
2051  for (j=1; j<=P; j++)
2052  {
2053    tmp[j] = "s("+string(j)+")";
2054  }
2055  // DName is defined earlier
2056  list NName = Name + DName + tmp;
2057  L[2] = NName;
2058  tmp = 0;
2059  // dp ordering;
2060  string s = "iv=";
2061  for (i=1; i<=Nnew; i++)
2062  {
2063    s = s+"1,";
2064  }
2065  s[size(s)]=";";
2066  execute(s);
2067  kill s;
2068  kill NName;
2069  tmp[1] = "dp";  //string
2070  tmp[2] = iv;    //intvec
2071  Lord[1] = tmp;
2072  tmp[1]  = "C";
2073  iv      = 0;
2074  tmp[2]  = iv;
2075  Lord[2] = tmp;
2076  tmp     = 0;
2077  L[3]    = Lord;
2078  // we are done with the list. Now add a Plural part
2079  def @R4@ = ring(L);
2080  setring @R4@;
2081  matrix @D[Nnew][Nnew];
2082  for (i=1; i<=N; i++)
2083  {
2084    @D[i,N+i]=1;
2085  }
2086  def @R4 = nc_algebra(1,@D);
2087  setring @R4;
2088  kill @R4@;
2089  dbprint(ppl,"// -4-1- the ring @R4i(_x,_Dx,_s) is ready");
2090  dbprint(ppl-1, @R4);
2091  ideal K4 = imap(@R,K);
2092  option(redSB);
2093  dbprint(ppl,"// -4-2- the final cosmetic std");
2094  K4 = engine(K4,eng);  //std does the job too
2095  // total cleanup
2096  kill @R;
2097  kill @R2;
2098  def BS = imap(@R3,BS);
2099  export BS;
2100  kill @R3;
2101  ideal LD = K4;
2102  export LD;
2103  return(@R4);
2104}
2105example
2106{
2107  "EXAMPLE:"; echo = 2;
2108  ring r = 0,(x,y),Dp;
2109  ideal F = x,y,x+y;
2110  printlevel = 0;
2111  def A = annfsBMI(F);
2112  setring A;
2113  LD;
2114  BS;
2115}
2116
2117proc annfsOT(poly F, list #)
2118"USAGE:  annfsOT(f [,eng]);  f a poly, eng an optional int
2119RETURN:  ring
2120PURPOSE: compute the D-module structure of basering[1/f]*f^s, according
2121to the algorithm by Oaku and Takayama
2122NOTE:    activate this ring with the @code{setring} command. In this ring,
2123@*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
2124@*         which is obtained by substituting the minimal integer root of a Bernstein
2125@*         polynomial into the s-parametric ideal;
2126@*       - the list BS contains roots with multiplicities of a Bernstein polynomial of f.
2127@*       If eng <>0, @code{std} is used for Groebner basis computations,
2128@*       otherwise, and by default @code{slimgb} is used.
2129@*       If printlevel=1, progress debug messages will be printed,
2130@*       if printlevel>=2, all the debug messages will be printed.
2131EXAMPLE: example annfsOT; shows examples
2132"
2133{
2134  int eng = 0;
2135  if ( size(#)>0 )
2136  {
2137    if ( typeof(#[1]) == "int" )
2138    {
2139      eng = int(#[1]);
2140    }
2141  }
2142  // returns a list with a ring and an ideal LD in it
2143  int ppl = printlevel-voice+2;
2144  //  printf("plevel :%s, voice: %s",printlevel,voice);
2145  def save = basering;
2146  int N = nvars(basering);
2147  int Nnew = 2*(N+2);
2148  int i,j;
2149  string s;
2150  list RL = ringlist(basering);
2151  list L, Lord;
2152  list tmp;
2153  intvec iv;
2154  L[1] = RL[1]; // char
2155  L[4] = RL[4]; // char, minpoly
2156  // check whether vars have admissible names
2157  list Name  = RL[2];
2158  list RName;
2159  RName[1] = "u";
2160  RName[2] = "v";
2161  RName[3] = "t";
2162  RName[4] = "Dt";
2163  for(i=1;i<=N;i++)
2164  {
2165    for(j=1; j<=size(RName);j++)
2166    {
2167      if (Name[i] == RName[j])
2168      {
2169        ERROR("Variable names should not include u,v,t,Dt");
2170      }
2171    }
2172  }
2173  // now, create the names for new vars
2174  tmp[1]     = "u";
2175  tmp[2]     = "v";
2176  list UName = tmp;
2177  list DName;
2178  for(i=1;i<=N;i++)
2179  {
2180    DName[i] = "D"+Name[i]; // concat
2181  }
2182  tmp    =  0;
2183  tmp[1] = "t";
2184  tmp[2] = "Dt";
2185  list NName = UName +  tmp + Name + DName;
2186  L[2]   = NName;
2187  tmp    = 0;
2188  // Name, Dname will be used further
2189  kill UName;
2190  kill NName;
2191  // block ord (a(1,1),dp);
2192  tmp[1]  = "a"; // string
2193  iv      = 1,1;
2194  tmp[2]  = iv; //intvec
2195  Lord[1] = tmp;
2196  // continue with dp 1,1,1,1...
2197  tmp[1]  = "dp"; // string
2198  s       = "iv=";
2199  for(i=1;i<=Nnew;i++)
2200  {
2201    s = s+"1,";
2202  }
2203  s[size(s)]= ";";
2204  execute(s);
2205  tmp[2]    = iv;
2206  Lord[2]   = tmp;
2207  tmp[1]    = "C";
2208  iv        = 0;
2209  tmp[2]    = iv;
2210  Lord[3]   = tmp;
2211  tmp       = 0;
2212  L[3]      = Lord;
2213  // we are done with the list
2214  def @R@ = ring(L);
2215  setring @R@;
2216  matrix @D[Nnew][Nnew];
2217  @D[3,4]=1;
2218  for(i=1; i<=N; i++)
2219  {
2220    @D[4+i,N+4+i]=1;
2221  }
2222  //  @D[N+3,2*(N+2)]=1; old t,Dt stuff
2223  //  L[5] = matrix(UpOneMatrix(Nnew));
2224  //  L[6] = @D;
2225  def @R = nc_algebra(1,@D);
2226  setring @R;
2227  kill @R@;
2228  dbprint(ppl,"// -1-1- the ring @R(u,v,t,Dt,_x,_Dx) is ready");
2229  dbprint(ppl-1, @R);
2230  // create the ideal I
2231  poly  F = imap(save,F);
2232  ideal I = u*F-t,u*v-1;
2233  poly p;
2234  for(i=1; i<=N; i++)
2235  {
2236    p = u*Dt; // u*Dt
2237    p = diff(F,var(4+i))*p;
2238    I = I, var(N+4+i) + p;
2239  }
2240  // -------- the ideal I is ready ----------
2241  dbprint(ppl,"// -1-2- starting the elimination of u,v in @R");
2242  dbprint(ppl-1, I);
2243  ideal J = engine(I,eng);
2244  ideal K = nselect(J,1,2);
2245  dbprint(ppl,"// -1-3- u,v are eliminated");
2246  dbprint(ppl-1, K);  // K is without u,v
2247  setring save;
2248  // ------------ new ring @R2 ------------------
2249  // without u,v and with the elim.ord for t,Dt
2250  // tensored with the K[s]
2251  // keep: N, i,j,s, tmp, RL
2252  Nnew = 2*N+2+1;
2253  //  list RL = ringlist(save); // is defined earlier
2254  L = 0;  //  kill L;
2255  kill Lord, tmp, iv, RName;
2256  list Lord, tmp;
2257  intvec iv;
2258  L[1] = RL[1]; L[4] = RL[4]; // char, minpoly
2259  // check whether vars have admissible names -> done earlier
2260  //  list Name  = RL[2];
2261  list RName;
2262  RName[1] = "t";
2263  RName[2] = "Dt";
2264  // now, create the names for new var (here, s only)
2265  tmp[1]     = "s";
2266  // DName is defined earlier
2267  list NName = RName + Name + DName + tmp;
2268  L[2]   = NName;
2269  tmp    = 0;
2270  // block ord (a(1,1),dp);
2271  tmp[1]  = "a";  iv = 1,1; tmp[2]  = iv; //intvec
2272  Lord[1] = tmp;
2273  // continue with a(1,1,1,1)...
2274  tmp[1]  = "dp"; s  = "iv=";
2275  for(i=1; i<= Nnew; i++)
2276  {
2277    s = s+"1,";
2278  }
2279  s[size(s)]= ";";  execute(s);
2280  kill NName;
2281  tmp[2]    = iv;
2282  Lord[2]   = tmp;
2283  // extra block for s
2284  // tmp[1] = "dp"; iv = 1;
2285  //  s[size(s)]= ",";  s = s + "1,1,1;";  execute(s);  tmp[2]    = iv;
2286  //  Lord[3]   = tmp;
2287  kill s;
2288  tmp[1]    = "C";   iv  = 0; tmp[2] = iv;
2289  Lord[3]   = tmp;   tmp = 0;
2290  L[3]      = Lord;
2291  // we are done with the list. Now, add a Plural part
2292  def @R2@ = ring(L);
2293  setring @R2@;
2294  matrix @D[Nnew][Nnew];
2295  @D[1,2] = 1;
2296  for(i=1; i<=N; i++)
2297  {
2298    @D[2+i,2+N+i] = 1;
2299  }
2300  def @R2 = nc_algebra(1,@D);
2301  setring @R2;
2302  kill @R2@;
2303  dbprint(ppl,"// -2-1- the ring @R2(t,Dt,_x,_Dx,s) is ready");
2304  dbprint(ppl-1, @R2);
2305  ideal MM = maxideal(1);
2306  MM = 0,0,MM;
2307  map R01 = @R, MM;
2308  ideal K = R01(K);
2309  //  ideal K = imap(@R,K); // names of vars are important!
2310  poly G = t*Dt+s+1; // s is a variable here
2311  K = NF(K,std(G)),G;
2312  // -------- the ideal K_(@R2) is ready ----------
2313  dbprint(ppl,"// -2-2- starting the elimination of t,Dt in @R2");
2314  dbprint(ppl-1, K);
2315  ideal M  = engine(K,eng);
2316  ideal K2 = nselect(M,1,2);
2317  dbprint(ppl,"// -2-3- t,Dt are eliminated");
2318  dbprint(ppl-1, K2);
2319  //  dbprint(ppl-1+1," -2-4- std of K2");
2320  //  option(redSB);  option(redTail);  K2 = std(K2);
2321  //  K2; // without t,Dt, and with s
2322  // -------- the ring @R3 ----------
2323  // _x, _Dx, s; elim.ord for _x,_Dx.
2324  // keep: N, i,j,s, tmp, RL
2325  setring save;
2326  Nnew = 2*N+1;
2327  //  list RL = ringlist(save); // is defined earlier
2328  //  kill L;
2329  kill Lord, tmp, iv, RName;
2330  list Lord, tmp;
2331  intvec iv;
2332  L[1] = RL[1]; L[4] = RL[4]; // char, minpoly
2333  // check whether vars have admissible names -> done earlier
2334  //  list Name  = RL[2];
2335  // now, create the names for new var (here, s only)
2336  tmp[1]     = "s";
2337  // DName is defined earlier
2338  list NName = Name + DName + tmp;
2339  L[2]   = NName;
2340  tmp    = 0;
2341  // block ord (a(1,1...),dp);
2342  string  s = "iv=";
2343  for(i=1; i<=Nnew-1; i++)
2344  {
2345    s = s+"1,";
2346  }
2347  s[size(s)]= ";";
2348  execute(s);
2349  tmp[1]  = "a"; // string
2350  tmp[2]  = iv; //intvec
2351  Lord[1] = tmp;
2352  // continue with dp 1,1,1,1...
2353  tmp[1]  = "dp"; // string
2354  s[size(s)]=","; s= s+"1;";
2355  execute(s);
2356  kill s;
2357  kill NName;
2358  tmp[2]    = iv;
2359  Lord[2]   = tmp;
2360  tmp[1]    = "C";  iv  = 0;  tmp[2] = iv;
2361  Lord[3]   = tmp;  tmp = 0;
2362  L[3]      = Lord;
2363  // we are done with the list. Now add a Plural part
2364  def @R3@ = ring(L);
2365  setring @R3@;
2366  matrix @D[Nnew][Nnew];
2367  for(i=1; i<=N; i++)
2368  {
2369    @D[i,N+i]=1;
2370  }
2371  def @R3 = nc_algebra(1,@D);
2372  setring @R3;
2373  kill @R3@;
2374  dbprint(ppl,"// -3-1- the ring @R3(_x,_Dx,s) is ready");
2375  dbprint(ppl-1, @R3);
2376  ideal MM = maxideal(1);
2377  MM = 0,0,MM;
2378  map R12 = @R2, MM;
2379  ideal K = R12(K2);
2380  poly  F = imap(save,F);
2381  K = K,F;
2382  dbprint(ppl,"// -3-2- starting the elimination of _x,_Dx in @R3");
2383  dbprint(ppl-1, K);
2384  ideal M = engine(K,eng);
2385  ideal K3 = nselect(M,1,Nnew-1);
2386  dbprint(ppl,"// -3-3-  _x,_Dx are eliminated in @R3");
2387  dbprint(ppl-1, K3);
2388  // the ring @R4  and the search for minimal negative int s
2389  ring @R4 = 0,(s),dp;
2390  dbprint(ppl,"// -4-1- the ring @R4 is ready");
2391  ideal K4 = imap(@R3,K3);
2392  poly p = K4[1];
2393  dbprint(ppl,"// -4-2- factorization");
2394////  ideal P = factorize(p,1);  // without constants and multiplicities
2395  list P = factorize(p);         // with constants and multiplicities
2396  ideal bs; intvec m;            // the Bernstein polynomial is monic, so we are not interested in constants
2397  for (i=2; i<=size(P[1]); i++)  // we delete P[1][1] and P[2][1]
2398  {
2399    bs[i-1] = P[1][i];
2400    m[i-1]  = P[2][i];
2401  }
2402  //  "------ b-function factorizes into ----------";  P;
2403////  int sP  = minIntRoot(P, 1);
2404  int sP = minIntRoot(bs,1);
2405  dbprint(ppl,"// -4-3- minimal integer root found");
2406  dbprint(ppl-1, sP);
2407  // convert factors to a list of their roots
2408  // assume all factors are linear
2409////  ideal BS = normalize(P);
2410////  BS = subst(BS,s,0);
2411////  BS = -BS;
2412  bs = normalize(bs);
2413  bs = subst(bs,s,0);
2414  bs = -bs;
2415  list BS = bs,m;
2416  // TODO: sort BS!
2417  // ------ substitute s found in the ideal ------
2418  // ------- going back to @R2 and substitute --------
2419  setring @R2;
2420  ideal K3 = subst(K2,s,sP);
2421  // create the ordinary Weyl algebra and put the result into it,
2422  // thus creating the ring @R5
2423  // keep: N, i,j,s, tmp, RL
2424  setring save;
2425  Nnew = 2*N;
2426  //  list RL = ringlist(save); // is defined earlier
2427  kill Lord, tmp, iv;
2428  L = 0;
2429  list Lord, tmp;
2430  intvec iv;
2431  L[1] = RL[1];   L[4] = RL[4]; // char, minpoly
2432  // check whether vars have admissible names -> done earlier
2433  //  list Name  = RL[2];
2434  // DName is defined earlier
2435  list NName = Name + DName;
2436  L[2]   = NName;
2437  // dp ordering;
2438  string   s       = "iv=";
2439  for(i=1;i<=Nnew;i++)
2440  {
2441    s = s+"1,";
2442  }
2443  s[size(s)]= ";";
2444  execute(s);
2445  tmp     = 0;
2446  tmp[1]  = "dp"; // string
2447  tmp[2]  = iv; //intvec
2448  Lord[1] = tmp;
2449  kill s;
2450  tmp[1]    = "C";
2451  iv        = 0;
2452  tmp[2]    = iv;
2453  Lord[2]   = tmp;
2454  tmp       = 0;
2455  L[3]      = Lord;
2456  // we are done with the list
2457  // Add: Plural part
2458  def @R5@ = ring(L);
2459  setring @R5@;
2460  matrix @D[Nnew][Nnew];
2461  for(i=1; i<=N; i++)
2462  {
2463    @D[i,N+i]=1;
2464  }
2465  def @R5 = nc_algebra(1,@D);
2466  setring @R5;
2467  kill @R5@;
2468  dbprint(ppl,"// -5-1- the ring @R5 is ready");
2469  dbprint(ppl-1, @R5);
2470  ideal K5 = imap(@R2,K3);
2471  option(redSB);
2472  dbprint(ppl,"// -5-2- the final cosmetic std");
2473  K5 = engine(K5,eng); // std does the job too
2474  // total cleanup
2475  kill @R;
2476  kill @R2;
2477  kill @R3;
2478////  ideal BS = imap(@R4,BS);
2479  list BS = imap(@R4,BS);
2480  export BS;
2481  ideal LD = K5;
2482  kill @R4;
2483  export LD;
2484  return(@R5);
2485}
2486example
2487{
2488  "EXAMPLE:"; echo = 2;
2489  ring r = 0,(x,y,z),Dp;
2490  poly F = x^2+y^3+z^5;
2491  printlevel = 0;
2492  def A  = annfsOT(F);
2493  setring A;
2494  LD;
2495  BS;
2496}
2497
2498
2499proc SannfsOT(poly F, list #)
2500"USAGE:  SannfsOT(f [,eng]);  f a poly, eng an optional int
2501RETURN:  ring
2502PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the 1st step of the algorithm by Oaku and Takayama in the ring D[s], where D is the Weyl algebra
2503NOTE:    activate this ring with the @code{setring} command.
2504@*       In the ring D[s], the ideal LD (which is NOT a Groebner basis) is the needed D-module structure.
2505@*       If eng <>0, @code{std} is used for Groebner basis computations,
2506@*       otherwise, and by default @code{slimgb} is used.
2507@*       If printlevel=1, progress debug messages will be printed,
2508@*       if printlevel>=2, all the debug messages will be printed.
2509EXAMPLE: example SannfsOT; shows examples
2510"
2511{
2512  int eng = 0;
2513  if ( size(#)>0 )
2514  {
2515    if ( typeof(#[1]) == "int" )
2516    {
2517      eng = int(#[1]);
2518    }
2519  }
2520  // returns a list with a ring and an ideal LD in it
2521  int ppl = printlevel-voice+2;
2522  //  printf("plevel :%s, voice: %s",printlevel,voice);
2523  def save = basering;
2524  int N = nvars(basering);
2525  int Nnew = 2*(N+2);
2526  int i,j;
2527  string s;
2528  list RL = ringlist(basering);
2529  list L, Lord;
2530  list tmp;
2531  intvec iv;
2532  L[1] = RL[1]; // char
2533  L[4] = RL[4]; // char, minpoly
2534  // check whether vars have admissible names
2535  list Name  = RL[2];
2536  list RName;
2537  RName[1] = "u";
2538  RName[2] = "v";
2539  RName[3] = "t";
2540  RName[4] = "Dt";
2541  for(i=1;i<=N;i++)
2542  {
2543    for(j=1; j<=size(RName);j++)
2544    {
2545      if (Name[i] == RName[j])
2546      {
2547        ERROR("Variable names should not include u,v,t,Dt");
2548      }
2549    }
2550  }
2551  // now, create the names for new vars
2552  tmp[1]     = "u";
2553  tmp[2]     = "v";
2554  list UName = tmp;
2555  list DName;
2556  for(i=1;i<=N;i++)
2557  {
2558    DName[i] = "D"+Name[i]; // concat
2559  }
2560  tmp    =  0;
2561  tmp[1] = "t";
2562  tmp[2] = "Dt";
2563  list NName = UName +  tmp + Name + DName;
2564  L[2]   = NName;
2565  tmp    = 0;
2566  // Name, Dname will be used further
2567  kill UName;
2568  kill NName;
2569  // block ord (a(1,1),dp);
2570  tmp[1]  = "a"; // string
2571  iv      = 1,1;
2572  tmp[2]  = iv; //intvec
2573  Lord[1] = tmp;
2574  // continue with dp 1,1,1,1...
2575  tmp[1]  = "dp"; // string
2576  s       = "iv=";
2577  for(i=1;i<=Nnew;i++)
2578  {
2579    s = s+"1,";
2580  }
2581  s[size(s)]= ";";
2582  execute(s);
2583  tmp[2]    = iv;
2584  Lord[2]   = tmp;
2585  tmp[1]    = "C";
2586  iv        = 0;
2587  tmp[2]    = iv;
2588  Lord[3]   = tmp;
2589  tmp       = 0;
2590  L[3]      = Lord;
2591  // we are done with the list
2592  def @R@ = ring(L);
2593  setring @R@;
2594  matrix @D[Nnew][Nnew];
2595  @D[3,4]=1;
2596  for(i=1; i<=N; i++)
2597  {
2598    @D[4+i,N+4+i]=1;
2599  }
2600  //  @D[N+3,2*(N+2)]=1; old t,Dt stuff
2601  //  L[5] = matrix(UpOneMatrix(Nnew));
2602  //  L[6] = @D;
2603  def @R =  nc_algebra(1,@D);
2604  setring @R;
2605  kill @R@;
2606  dbprint(ppl,"// -1-1- the ring @R(u,v,t,Dt,_x,_Dx) is ready");
2607  dbprint(ppl-1, @R);
2608  // create the ideal I
2609  poly  F = imap(save,F);
2610  ideal I = u*F-t,u*v-1;
2611  poly p;
2612  for(i=1; i<=N; i++)
2613  {
2614    p = u*Dt; // u*Dt
2615    p = diff(F,var(4+i))*p;
2616    I = I, var(N+4+i) + p;
2617  }
2618  // -------- the ideal I is ready ----------
2619  dbprint(ppl,"// -1-2- starting the elimination of u,v in @R");
2620  dbprint(ppl-1, I);
2621  ideal J = engine(I,eng);
2622  ideal K = nselect(J,1,2);
2623  dbprint(ppl,"// -1-3- u,v are eliminated");
2624  dbprint(ppl-1, K);  // K is without u,v
2625
2626
2627  setring save;
2628  // ------------ new ring @R2 ------------------
2629  // without u,v and with the elim.ord for t,Dt
2630  // tensored with the K[s]
2631  // keep: N, i,j,s, tmp, RL
2632  Nnew = 2*N+2+1;
2633  //  list RL = ringlist(save); // is defined earlier
2634  L = 0;  //  kill L;
2635  kill Lord, tmp, iv, RName;
2636  list Lord, tmp;
2637  intvec iv;
2638  L[1] = RL[1]; L[4] = RL[4]; // char, minpoly
2639  // check whether vars have admissible names -> done earlier
2640  //  list Name  = RL[2];
2641  list RName;
2642  RName[1] = "t";
2643  RName[2] = "Dt";
2644  // now, create the names for new var (here, s only)
2645  tmp[1]     = "s";
2646  // DName is defined earlier
2647  list NName = RName + Name + DName + tmp;
2648  L[2]   = NName;
2649  tmp    = 0;
2650  // block ord (a(1,1),dp);
2651  tmp[1]  = "a";  iv = 1,1; tmp[2]  = iv; //intvec
2652  Lord[1] = tmp;
2653  // continue with a(1,1,1,1)...
2654  tmp[1]  = "dp"; s  = "iv=";
2655  for(i=1; i<= Nnew; i++)
2656  {
2657    s = s+"1,";
2658  }
2659  s[size(s)]= ";";  execute(s);
2660  kill NName;
2661  tmp[2]    = iv;
2662  Lord[2]   = tmp;
2663  // extra block for s
2664  // tmp[1] = "dp"; iv = 1;
2665  //  s[size(s)]= ",";  s = s + "1,1,1;";  execute(s);  tmp[2]    = iv;
2666  //  Lord[3]   = tmp;
2667  kill s;
2668  tmp[1]    = "C";   iv  = 0; tmp[2] = iv;
2669  Lord[3]   = tmp;   tmp = 0;
2670  L[3]      = Lord;
2671  // we are done with the list. Now, add a Plural part
2672  def @R2@ = ring(L);
2673  setring @R2@;
2674  matrix @D[Nnew][Nnew];
2675  @D[1,2] = 1;
2676  for(i=1; i<=N; i++)
2677  {
2678    @D[2+i,2+N+i] = 1;
2679  }
2680  def @R2 = nc_algebra(1,@D);
2681  setring @R2;
2682  kill @R2@;
2683  dbprint(ppl,"// -2-1- the ring @R2(t,Dt,_x,_Dx,s) is ready");
2684  dbprint(ppl-1, @R2);
2685  ideal MM = maxideal(1);
2686  MM = 0,0,MM;
2687  map R01 = @R, MM;
2688  ideal K = R01(K);
2689  //  ideal K = imap(@R,K); // names of vars are important!
2690  poly G = t*Dt+s+1; // s is a variable here
2691  K = NF(K,std(G)),G;
2692  // -------- the ideal K_(@R2) is ready ----------
2693  dbprint(ppl,"// -2-2- starting the elimination of t,Dt in @R2");
2694  dbprint(ppl-1, K);
2695  ideal M  = engine(K,eng);
2696  ideal K2 = nselect(M,1,2);
2697  dbprint(ppl,"// -2-3- t,Dt are eliminated");
2698  dbprint(ppl-1, K2);
2699  K2 = engine(K2,eng);
2700  setring save;
2701  // ----------- the ring @R3 ------------
2702  // _x, _Dx,s;  elim.ord for _x,_Dx.
2703  // keep: N, i,j,s, tmp, RL
2704  Nnew = 2*N+1;
2705  kill Lord, tmp, iv, RName;
2706  list Lord, tmp;
2707  intvec iv;
2708  L[1] = RL[1];
2709  L[4] = RL[4];  // char, minpoly
2710  // check whether vars hava admissible names -> done earlier
2711  // now, create the names for new var
2712  tmp[1] = "s";
2713  // DName is defined earlier
2714  list NName = Name + DName + tmp;
2715  L[2] = NName;
2716  tmp = 0;
2717  // block ord (dp(N),dp);
2718  string s = "iv=";
2719  for (i=1; i<=Nnew-1; i++)
2720  {
2721    s = s+"1,";
2722  }
2723  s[size(s)]=";";
2724  execute(s);
2725  tmp[1] = "dp";  // string
2726  tmp[2] = iv;   // intvec
2727  Lord[1] = tmp;
2728  // continue with dp 1,1,1,1...
2729  tmp[1] = "dp";  // string
2730  s[size(s)] = ",";
2731  s = s+"1;";
2732  execute(s);
2733  kill s;
2734  kill NName;
2735  tmp[2]      = iv;
2736  Lord[2]     = tmp;
2737  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
2738  Lord[3]     = tmp;  tmp = 0;
2739  L[3]        = Lord;
2740  // we are done with the list. Now add a Plural part
2741  def @R3@ = ring(L);
2742  setring @R3@;
2743  matrix @D[Nnew][Nnew];
2744  for (i=1; i<=N; i++)
2745  {
2746    @D[i,N+i]=1;
2747  }
2748  def @R3 = nc_algebra(1,@D);
2749  setring @R3;
2750  kill @R3@;
2751  dbprint(ppl,"// -3-1- the ring @R3(_x,_Dx,s) is ready");
2752  dbprint(ppl-1, @R3);
2753  ideal MM = maxideal(1);
2754  MM = 0,s,MM;
2755  map R01 = @R2, MM;
2756  ideal K2 = R01(K2);
2757  // total cleanup
2758  ideal LD = K2;
2759  // make leadcoeffs positive
2760  for (i=1; i<= ncols(LD); i++)
2761  {
2762    if (leadcoef(LD[i]) <0 )
2763    {
2764      LD[i] = -LD[i];
2765    }
2766  }
2767  export LD;
2768  kill @R;
2769  kill @R2;
2770  return(@R3);
2771}
2772example
2773{
2774  "EXAMPLE:"; echo = 2;
2775  ring r = 0,(x,y,z),Dp;
2776  poly F = x^3+y^3+z^3;
2777  printlevel = 0;
2778  def A  = SannfsOT(F);
2779  setring A;
2780  LD;
2781}
2782
2783proc SannfsBM(poly F, list #)
2784"USAGE:  SannfsBM(f [,eng]);  f a poly, eng an optional int
2785RETURN:  ring
2786PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the 1st step of the algorithm by Briancon and Maisonobe in the ring D[s], where D is the Weyl algebra
2787NOTE:    activate this ring with the @code{setring} command.
2788@*       In the ring D[s], the ideal LD (which is NOT a Groebner basis) is the needed D-module structure,
2789@*       If eng <>0, @code{std} is used for Groebner basis computations,
2790@*       otherwise, and by default @code{slimgb} is used.
2791@*       If printlevel=1, progress debug messages will be printed,
2792@*       if printlevel>=2, all the debug messages will be printed.
2793EXAMPLE: example SannfsBM; shows examples
2794"
2795{
2796  int eng = 0;
2797  if ( size(#)>0 )
2798  {
2799    if ( typeof(#[1]) == "int" )
2800    {
2801      eng = int(#[1]);
2802    }
2803  }
2804  // returns a list with a ring and an ideal LD in it
2805  int ppl = printlevel-voice+2;
2806  //  printf("plevel :%s, voice: %s",printlevel,voice);
2807  def save = basering;
2808  int N = nvars(basering);
2809  int Nnew = 2*N+2;
2810  int i,j;
2811  string s;
2812  list RL = ringlist(basering);
2813  list L, Lord;
2814  list tmp;
2815  intvec iv;
2816  L[1] = RL[1]; // char
2817  L[4] = RL[4]; // char, minpoly
2818  // check whether vars have admissible names
2819  list Name  = RL[2];
2820  list RName;
2821  RName[1] = "t";
2822  RName[2] = "s";
2823  for(i=1;i<=N;i++)
2824  {
2825    for(j=1; j<=size(RName);j++)
2826    {
2827      if (Name[i] == RName[j])
2828      {
2829        ERROR("Variable names should not include t,s");
2830      }
2831    }
2832  }
2833  // now, create the names for new vars
2834  list DName;
2835  for(i=1;i<=N;i++)
2836  {
2837    DName[i] = "D"+Name[i]; // concat
2838  }
2839  tmp[1] = "t";
2840  tmp[2] = "s";
2841  list NName = tmp + Name + DName;
2842  L[2]   = NName;
2843  // Name, Dname will be used further
2844  kill NName;
2845  // block ord (lp(2),dp);
2846  tmp[1]  = "lp"; // string
2847  iv      = 1,1;
2848  tmp[2]  = iv; //intvec
2849  Lord[1] = tmp;
2850  // continue with dp 1,1,1,1...
2851  tmp[1]  = "dp"; // string
2852  s       = "iv=";
2853  for(i=1;i<=Nnew;i++)
2854  {
2855    s = s+"1,";
2856  }
2857  s[size(s)]= ";";
2858  execute(s);
2859  kill s;
2860  tmp[2]    = iv;
2861  Lord[2]   = tmp;
2862  tmp[1]    = "C";
2863  iv        = 0;
2864  tmp[2]    = iv;
2865  Lord[3]   = tmp;
2866  tmp       = 0;
2867  L[3]      = Lord;
2868  // we are done with the list
2869  def @R@ = ring(L);
2870  setring @R@;
2871  matrix @D[Nnew][Nnew];
2872  @D[1,2]=t;
2873  for(i=1; i<=N; i++)
2874  {
2875    @D[2+i,N+2+i]=1;
2876  }
2877  //  L[5] = matrix(UpOneMatrix(Nnew));
2878  //  L[6] = @D;
2879  def @R = nc_algebra(1,@D);
2880  setring @R;
2881  kill @R@;
2882  dbprint(ppl,"// -1-1- the ring @R(t,s,_x,_Dx) is ready");
2883  dbprint(ppl-1, @R);
2884  // create the ideal I
2885  poly  F = imap(save,F);
2886  ideal I = t*F+s;
2887  poly p;
2888  for(i=1; i<=N; i++)
2889  {
2890    p = t; // t
2891    p = diff(F,var(2+i))*p;
2892    I = I, var(N+2+i) + p;
2893  }
2894  // -------- the ideal I is ready ----------
2895  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
2896  dbprint(ppl-1, I);
2897  ideal J = engine(I,eng);
2898  ideal K = nselect(J,1);
2899  dbprint(ppl,"// -1-3- t is eliminated");
2900  dbprint(ppl-1, K);  // K is without t
2901  K = engine(K,eng);  // std does the job too
2902  // now, we must change the ordering
2903  // and create a ring without t, Dt
2904  //  setring S;
2905  // ----------- the ring @R3 ------------
2906  // _x, _Dx,s;  elim.ord for _x,_Dx.
2907  // keep: N, i,j,s, tmp, RL
2908  Nnew = 2*N+1;
2909  kill Lord, tmp, iv, RName;
2910  list Lord, tmp;
2911  intvec iv;
2912  list L=imap(save,L);
2913  list RL=imap(save,RL);
2914  L[1] = RL[1];
2915  L[4] = RL[4];  // char, minpoly
2916  // check whether vars hava admissible names -> done earlier
2917  // now, create the names for new var
2918  tmp[1] = "s";
2919  // DName is defined earlier
2920  list NName = Name + DName + tmp;
2921  L[2] = NName;
2922  tmp = 0;
2923  // block ord (dp(N),dp);
2924  string s = "iv=";
2925  for (i=1; i<=Nnew-1; i++)
2926  {
2927    s = s+"1,";
2928  }
2929  s[size(s)]=";";
2930  execute(s);
2931  tmp[1] = "dp";  // string
2932  tmp[2] = iv;   // intvec
2933  Lord[1] = tmp;
2934  // continue with dp 1,1,1,1...
2935  tmp[1] = "dp";  // string
2936  s[size(s)] = ",";
2937  s = s+"1;";
2938  execute(s);
2939  kill s;
2940  kill NName;
2941  tmp[2]      = iv;
2942  Lord[2]     = tmp;
2943  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
2944  Lord[3]     = tmp;  tmp = 0;
2945  L[3]        = Lord;
2946  // we are done with the list. Now add a Plural part
2947  def @R2@ = ring(L);
2948  setring @R2@;
2949  matrix @D[Nnew][Nnew];
2950  for (i=1; i<=N; i++)
2951  {
2952    @D[i,N+i]=1;
2953  }
2954  def @R2 = nc_algebra(1,@D);
2955  setring @R2;
2956  kill @R2@;
2957  dbprint(ppl,"//  -2-1- the ring @R2(_x,_Dx,s) is ready");
2958  dbprint(ppl-1, @R2);
2959  ideal MM = maxideal(1);
2960  MM = 0,s,MM;
2961  map R01 = @R, MM;
2962  ideal K = R01(K);
2963  // total cleanup
2964  ideal LD = K;
2965  // make leadcoeffs positive
2966  for (i=1; i<= ncols(LD); i++)
2967  {
2968    if (leadcoef(LD[i]) <0 )
2969    {
2970      LD[i] = -LD[i];
2971    }
2972  }
2973  export LD;
2974  kill @R;
2975  return(@R2);
2976}
2977example
2978{
2979  "EXAMPLE:"; echo = 2;
2980  ring r = 0,(x,y,z),Dp;
2981  poly F = x^3+y^3+z^3;
2982  printlevel = 0;
2983  def A = SannfsBM(F);
2984  setring A;
2985  LD;
2986}
2987
2988proc SannfsLOT(poly F, list #)
2989"USAGE:  SannfsLOT(f [,eng]);  f a poly, eng an optional int
2990RETURN:  ring
2991PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the Levandovskyy's modification of the algorithm by Oaku and Takayama in the ring D[s], where D is the Weyl algebra
2992NOTE:    activate this ring with the @code{setring} command.
2993@*       In the ring D[s], the ideal LD (which is NOT a Groebner basis) is the needed D-module structure.
2994@*       If eng <>0, @code{std} is used for Groebner basis computations,
2995@*       otherwise, and by default @code{slimgb} is used.
2996@*       If printlevel=1, progress debug messages will be printed,
2997@*       if printlevel>=2, all the debug messages will be printed.
2998EXAMPLE: example SannfsLOT; shows examples
2999"
3000{
3001  int eng = 0;
3002  if ( size(#)>0 )
3003  {
3004    if ( typeof(#[1]) == "int" )
3005    {
3006      eng = int(#[1]);
3007    }
3008  }
3009  // returns a list with a ring and an ideal LD in it
3010  int ppl = printlevel-voice+2;
3011  //  printf("plevel :%s, voice: %s",printlevel,voice);
3012  def save = basering;
3013  int N = nvars(basering);
3014  //  int Nnew = 2*(N+2);
3015  int Nnew = 2*(N+1)+1; //removed u,v; added s
3016  int i,j;
3017  string s;
3018  list RL = ringlist(basering);
3019  list L, Lord;
3020  list tmp;
3021  intvec iv;
3022  L[1] = RL[1]; // char
3023  L[4] = RL[4]; // char, minpoly
3024  // check whether vars have admissible names
3025  list Name  = RL[2];
3026  list RName;
3027//   RName[1] = "u";
3028//   RName[2] = "v";
3029  RName[1] = "t";
3030  RName[2] = "Dt";
3031  for(i=1;i<=N;i++)
3032  {
3033    for(j=1; j<=size(RName);j++)
3034    {
3035      if (Name[i] == RName[j])
3036      {
3037        ERROR("Variable names should not include t,Dt");
3038      }
3039    }
3040  }
3041  // now, create the names for new vars
3042//   tmp[1]     = "u";
3043//   tmp[2]     = "v";
3044//   list UName = tmp;
3045  list DName;
3046  for(i=1;i<=N;i++)
3047  {
3048    DName[i] = "D"+Name[i]; // concat
3049  }
3050  tmp    =  0;
3051  tmp[1] = "t";
3052  tmp[2] = "Dt";
3053  list SName ; SName[1] = "s";
3054  //  list NName = UName +  tmp + Name + DName;
3055  list NName = tmp + Name + DName + SName;
3056  L[2]   = NName;
3057  tmp    = 0;
3058  // Name, Dname will be used further
3059  //  kill UName;
3060  kill NName;
3061  // block ord (a(1,1),dp);
3062  tmp[1]  = "a"; // string
3063  iv      = 1,1;
3064  tmp[2]  = iv; //intvec
3065  Lord[1] = tmp;
3066  // continue with dp 1,1,1,1...
3067  tmp[1]  = "dp"; // string
3068  s       = "iv=";
3069  for(i=1;i<=Nnew;i++)
3070  {
3071    s = s+"1,";
3072  }
3073  s[size(s)]= ";";
3074  execute(s);
3075  tmp[2]    = iv;
3076  Lord[2]   = tmp;
3077  tmp[1]    = "C";
3078  iv        = 0;
3079  tmp[2]    = iv;
3080  Lord[3]   = tmp;
3081  tmp       = 0;
3082  L[3]      = Lord;
3083  // we are done with the list
3084  def @R@ = ring(L);
3085  setring @R@;
3086  matrix @D[Nnew][Nnew];
3087  @D[1,2]=1;
3088  for(i=1; i<=N; i++)
3089  {
3090    @D[2+i,N+2+i]=1;
3091  }
3092  // ADD [s,t]=-t, [s,Dt]=Dt
3093  @D[1,Nnew] = -var(1);
3094  @D[2,Nnew] = var(2);
3095  //  @D[N+3,2*(N+2)]=1; old t,Dt stuff
3096  //  L[5] = matrix(UpOneMatrix(Nnew));
3097  //  L[6] = @D;
3098  def @R = nc_algebra(1,@D);
3099  setring @R;
3100  kill @R@;
3101  dbprint(ppl,"// -1-1- the ring @R(t,Dt,_x,_Dx,s) is ready");
3102  dbprint(ppl-1, @R);
3103  // create the ideal I
3104  poly  F = imap(save,F);
3105  //  ideal I = u*F-t,u*v-1;
3106  ideal I = F-t;
3107  poly p;
3108  for(i=1; i<=N; i++)
3109  {
3110    //    p = u*Dt; // u*Dt
3111    p = Dt;
3112    p = diff(F,var(2+i))*p;
3113    I = I, var(N+2+i) + p;
3114  }
3115  //  I = I, var(1)*var(2) + var(Nnew) +1; // reduce it with t-f!!!
3116  // t*Dt + s +1 reduced with t-f gives f*Dt + s
3117  I = I, F*var(2) + var(Nnew);
3118  // -------- the ideal I is ready ----------
3119  dbprint(ppl,"// -1-2- starting the elimination of t,Dt in @R");
3120  dbprint(ppl-1, I);
3121  ideal J = engine(I,eng);
3122  ideal K = nselect(J,1,2);
3123  dbprint(ppl,"// -1-3- t,Dt are eliminated");
3124  dbprint(ppl-1, K);  // K is without t, Dt
3125  K = engine(K,eng);  // std does the job too
3126  // now, we must change the ordering
3127  // and create a ring without t, Dt
3128  setring save;
3129  // ----------- the ring @R3 ------------
3130  // _x, _Dx,s;  elim.ord for _x,_Dx.
3131  // keep: N, i,j,s, tmp, RL
3132  Nnew = 2*N+1;
3133  kill Lord, tmp, iv, RName;
3134  list Lord, tmp;
3135  intvec iv;
3136  L[1] = RL[1];
3137  L[4] = RL[4];  // char, minpoly
3138  // check whether vars hava admissible names -> done earlier
3139  // now, create the names for new var
3140  tmp[1] = "s";
3141  // DName is defined earlier
3142  list NName = Name + DName + tmp;
3143  L[2] = NName;
3144  tmp = 0;
3145  // block ord (dp(N),dp);
3146  // string s is already defined
3147  s = "iv=";
3148  for (i=1; i<=Nnew-1; i++)
3149  {
3150    s = s+"1,";
3151  }
3152  s[size(s)]=";";
3153  execute(s);
3154  tmp[1] = "dp";  // string
3155  tmp[2] = iv;   // intvec
3156  Lord[1] = tmp;
3157  // continue with dp 1,1,1,1...
3158  tmp[1] = "dp";  // string
3159  s[size(s)] = ",";
3160  s = s+"1;";
3161  execute(s);
3162  kill s;
3163  kill NName;
3164  tmp[2]      = iv;
3165  Lord[2]     = tmp;
3166  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
3167  Lord[3]     = tmp;  tmp = 0;
3168  L[3]        = Lord;
3169  // we are done with the list. Now add a Plural part
3170  def @R2@ = ring(L);
3171  setring @R2@;
3172  matrix @D[Nnew][Nnew];
3173  for (i=1; i<=N; i++)
3174  {
3175    @D[i,N+i]=1;
3176  }
3177  def @R2 = nc_algebra(1,@D);
3178  setring @R2;
3179  kill @R2@;
3180  dbprint(ppl,"//  -2-1- the ring @R2(_x,_Dx,s) is ready");
3181  dbprint(ppl-1, @R2);
3182  ideal MM = maxideal(1);
3183  MM = 0,s,MM;
3184  map R01 = @R, MM;
3185  ideal K = R01(K);
3186  // total cleanup
3187  ideal LD = K;
3188  // make leadcoeffs positive
3189  for (i=1; i<= ncols(LD); i++)
3190  {
3191    if (leadcoef(LD[i]) <0 )
3192    {
3193      LD[i] = -LD[i];
3194    }
3195  }
3196  export LD;
3197  kill @R;
3198  return(@R2);
3199}
3200example
3201{
3202  "EXAMPLE:"; echo = 2;
3203  ring r = 0,(x,y,z),Dp;
3204  poly F = x^3+y^3+z^3;
3205  printlevel = 0;
3206  def A  = SannfsLOT(F);
3207  setring A;
3208  LD;
3209}
3210
3211
3212proc annfsLOT(poly F, list #)
3213"USAGE:  annfsLOT(F [,eng]);  F a poly, eng an optional int
3214RETURN:  ring
3215PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the Levandovskyy's modification of the algorithm by Oaku and Takayama
3216NOTE:    activate this ring with the @code{setring} command. In this ring,
3217@*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
3218@*         which is obtained by substituting the minimal integer root of a Bernstein
3219@*         polynomial into the s-parametric ideal;
3220@*       - the list BS contains the roots with multiplicities of a Bernstein polynomial of f.
3221@*       If eng <>0, @code{std} is used for Groebner basis computations,
3222@*       otherwise and by default @code{slimgb} is used.
3223@*       If printlevel=1, progress debug messages will be printed,
3224@*       if printlevel>=2, all the debug messages will be printed.
3225EXAMPLE: example annfsLOT; shows examples
3226"
3227{
3228  int eng = 0;
3229  if ( size(#)>0 )
3230  {
3231    if ( typeof(#[1]) == "int" )
3232    {
3233      eng = int(#[1]);
3234    }
3235  }
3236  printlevel=printlevel+1;
3237  def save = basering;
3238  def @A = SannfsLOT(F,eng);
3239  setring @A;
3240  poly F = imap(save,F);
3241  def B  = annfs0(LD,F,eng);
3242  return(B);
3243}
3244example
3245{
3246  "EXAMPLE:"; echo = 2;
3247  ring r = 0,(x,y,z),Dp;
3248  poly F = z*x^2+y^3;
3249  printlevel = 0;
3250  def A  = annfsLOT(F);
3251  setring A;
3252  LD;
3253  BS;
3254}
3255
3256proc annfs0(ideal I, poly F, list #)
3257"USAGE:  annfs0(I, F [,eng]);  I an ideal, F a poly, eng an optional int
3258RETURN:  ring
3259PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based on the
3260output of procedures SannfsBM, SannfsOT or SannfsLOT
3261NOTE:    activate this ring with the @code{setring} command. In this ring,
3262@*       - the ideal LD (which is a Groebner basis) is the annihilator of f^s,
3263@*       - the list BS contains the roots with multiplicities of a Bernstein polynomial of f.
3264@*       If eng <>0, @code{std} is used for Groebner basis computations,
3265@*       otherwise and by default @code{slimgb} is used.
3266@*       If printlevel=1, progress debug messages will be printed,
3267@*       if printlevel>=2, all the debug messages will be printed.
3268EXAMPLE: example annfs0; shows examples
3269"
3270{
3271  int eng = 0;
3272  if ( size(#)>0 )
3273  {
3274    if ( typeof(#[1]) == "int" )
3275    {
3276      eng = int(#[1]);
3277    }
3278  }
3279  def @R2 = basering;
3280  // we're in D_n[s], where the elim ord for s is set
3281  ideal J = NF(I,std(F));
3282  // make leadcoeffs positive
3283  int i;
3284  for (i=1; i<= ncols(J); i++)
3285  {
3286    if (leadcoef(J[i]) <0 )
3287    {
3288      J[i] = -J[i];
3289    }
3290  }
3291  J = J,F;
3292  ideal M = engine(J,eng);
3293  int Nnew = nvars(@R2);
3294  ideal K2 = nselect(M,1,Nnew-1);
3295  int ppl = printlevel-voice+2;
3296  dbprint(ppl,"// -1-1- _x,_Dx are eliminated in basering");
3297  dbprint(ppl-1, K2);
3298  // the ring @R3 and the search for minimal negative int s
3299  ring @R3 = 0,s,dp;
3300  dbprint(ppl,"// -2-1- the ring @R3 i.e. K[s] is ready");
3301  ideal K3 = imap(@R2,K2);
3302  poly p = K3[1];
3303  dbprint(ppl,"// -2-2- factorization");
3304  //  ideal P = factorize(p,1);  //without constants and multiplicities
3305  //  "--------- b-function factorizes into ---------";   P;
3306  // convert factors to the list of their roots with mults
3307  // assume all factors are linear
3308  //  ideal BS = normalize(P);
3309  //  BS = subst(BS,s,0);
3310  //  BS = -BS;
3311  list P = factorize(p);          //with constants and multiplicities
3312  ideal bs; intvec m;             //the Bernstein polynomial is monic, so we are not interested in constants
3313  for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
3314  {
3315    bs[i-1] = P[1][i];
3316    m[i-1]  = P[2][i];
3317  }
3318  int sP = minIntRoot(bs,1);
3319  bs =  normalize(bs);
3320  bs = -subst(bs,s,0);
3321  dbprint(ppl,"// -2-3- minimal integer root found");
3322  dbprint(ppl-1, sP);
3323 //TODO: sort BS!
3324  // --------- substitute s found in the ideal ---------
3325  // --------- going back to @R and substitute ---------
3326  setring @R2;
3327  K2 = subst(I,s,sP);
3328  // create the ordinary Weyl algebra and put the result into it,
3329  // thus creating the ring @R5
3330  // keep: N, i,j,s, tmp, RL
3331  Nnew = Nnew - 1; // former 2*N;
3332  // list RL = ringlist(save);  // is defined earlier
3333  //  kill Lord, tmp, iv;
3334  list L = 0;
3335  list Lord, tmp;
3336  intvec iv;
3337  list RL = ringlist(basering);
3338  L[1] = RL[1];
3339  L[4] = RL[4];  //char, minpoly
3340  // check whether vars have admissible names -> done earlier
3341  // list Name = RL[2]M
3342  // DName is defined earlier
3343  list NName; // = RL[2]; // skip the last var 's'
3344  for (i=1; i<=Nnew; i++)
3345  {
3346    NName[i] =  RL[2][i];
3347  }
3348  L[2] = NName;
3349  // dp ordering;
3350  string s = "iv=";
3351  for (i=1; i<=Nnew; i++)
3352  {
3353    s = s+"1,";
3354  }
3355  s[size(s)] = ";";
3356  execute(s);
3357  tmp     = 0;
3358  tmp[1]  = "dp";  // string
3359  tmp[2]  = iv;  // intvec
3360  Lord[1] = tmp;
3361  kill s;
3362  tmp[1]  = "C";
3363  iv = 0;
3364  tmp[2]  = iv;
3365  Lord[2] = tmp;
3366  tmp     = 0;
3367  L[3]    = Lord;
3368  // we are done with the list
3369  // Add: Plural part
3370  def @R4@ = ring(L);
3371  setring @R4@;
3372  int N = Nnew/2;
3373  matrix @D[Nnew][Nnew];
3374  for (i=1; i<=N; i++)
3375  {
3376    @D[i,N+i]=1;
3377  }
3378  def @R4 = nc_algebra(1,@D);
3379  setring @R4;
3380  kill @R4@;
3381  dbprint(ppl,"// -3-1- the ring @R4 is ready");
3382  dbprint(ppl-1, @R4);
3383  ideal K4 = imap(@R2,K2);
3384  option(redSB);
3385  dbprint(ppl,"// -3-2- the final cosmetic std");
3386  K4 = engine(K4,eng);  // std does the job too
3387  // total cleanup
3388  ideal bs = imap(@R3,bs);
3389  kill @R3;
3390  list BS = bs,m;
3391  export BS;
3392  ideal LD = K4;
3393  export LD;
3394  return(@R4);
3395}
3396example
3397{ "EXAMPLE:"; echo = 2;
3398  ring r = 0,(x,y,z),Dp;
3399  poly F = x^3+y^3+z^3;
3400  printlevel = 0;
3401  def A = SannfsBM(F);
3402  // alternatively, one can use SannfsOT or SannfsLOT
3403  setring A;
3404  LD;
3405  poly F = imap(r,F);
3406  def B  = annfs0(LD,F);
3407  setring B;
3408  LD;
3409  BS;
3410}
3411
3412// proc annfsgms(poly F, list #)
3413// "USAGE:  annfsgms(f [,eng]);  f a poly, eng an optional int
3414// ASSUME:  f has an isolated critical point at 0
3415// RETURN:  ring
3416// PURPOSE: compute the D-module structure of basering[1/f]*f^s
3417// NOTE:    activate this ring with the @code{setring} command. In this ring,
3418// @*       - the ideal LD is the needed D-mod structure,
3419// @*       - the ideal BS is the list of roots of a Bernstein polynomial of f.
3420// @*       If eng <>0, @code{std} is used for Groebner basis computations,
3421// @*       otherwise (and by default) @code{slimgb} is used.
3422// @*       If printlevel=1, progress debug messages will be printed,
3423// @*       if printlevel>=2, all the debug messages will be printed.
3424// EXAMPLE: example annfsgms; shows examples
3425// "
3426// {
3427//   LIB "gmssing.lib";
3428//   int eng = 0;
3429//   if ( size(#)>0 )
3430//   {
3431//     if ( typeof(#[1]) == "int" )
3432//     {
3433//       eng = int(#[1]);
3434//     }
3435//   }
3436//   int ppl = printlevel-voice+2;
3437//   // returns a ring with the ideal LD in it
3438//   def save = basering;
3439//   // compute the Bernstein poly from gmssing.lib
3440//   list RL = ringlist(basering);
3441//   // in the descr. of the ordering, replace "p" by "s"
3442//   list NL = convloc(RL);
3443//   // create a ring with the ordering, converted to local
3444//   def @LR = ring(NL);
3445//   setring @LR;
3446//   poly F  = imap(save, F);
3447//   ideal B = bernstein(F)[1];
3448//   // since B may not contain (s+1) [following gmssing.lib]
3449//   // add it!
3450//   B = B,-1;
3451//   B = simplify(B,2+4); // erase zero and repeated entries
3452//   // find the minimal integer value
3453//   int   S = minIntRoot(B,0);
3454//   dbprint(ppl,"// -0- minimal integer root found");
3455//   dbprint(ppl-1,S);
3456//   setring save;
3457//   int N = nvars(basering);
3458//   int Nnew = 2*(N+2);
3459//   int i,j;
3460//   string s;
3461//   //  list RL = ringlist(basering);
3462//   list L, Lord;
3463//   list tmp;
3464//   intvec iv;
3465//   L[1] = RL[1]; // char
3466//   L[4] = RL[4]; // char, minpoly
3467//   // check whether vars have admissible names
3468//   list Name  = RL[2];
3469//   list RName;
3470//   RName[1] = "u";
3471//   RName[2] = "v";
3472//   RName[3] = "t";
3473//   RName[4] = "Dt";
3474//   for(i=1;i<=N;i++)
3475//   {
3476//     for(j=1; j<=size(RName);j++)
3477//     {
3478//       if (Name[i] == RName[j])
3479//       {
3480//         ERROR("Variable names should not include u,v,t,Dt");
3481//       }
3482//     }
3483//   }
3484//   // now, create the names for new vars
3485//   //  tmp[1]     = "u";  tmp[2]     = "v";  tmp[3]     = "t";  tmp[4]     = "Dt";
3486//   list UName = RName;
3487//   list DName;
3488//   for(i=1;i<=N;i++)
3489//   {
3490//     DName[i] = "D"+Name[i]; // concat
3491//   }
3492//   list NName = UName + Name + DName;
3493//   L[2]       = NName;
3494//   tmp        = 0;
3495//   // Name, Dname will be used further
3496//   kill UName;
3497//   kill NName;
3498//   // block ord (a(1,1),dp);
3499//   tmp[1]  = "a"; // string
3500//   iv      = 1,1;
3501//   tmp[2]  = iv; //intvec
3502//   Lord[1] = tmp;
3503//   // continue with dp 1,1,1,1...
3504//   tmp[1]  = "dp"; // string
3505//   s       = "iv=";
3506//   for(i=1; i<=Nnew; i++) // need really all vars!
3507//   {
3508//     s = s+"1,";
3509//   }
3510//   s[size(s)]= ";";
3511//   execute(s);
3512//   tmp[2]    = iv;
3513//   Lord[2]   = tmp;
3514//   tmp[1]    = "C";
3515//   iv        = 0;
3516//   tmp[2]    = iv;
3517//   Lord[3]   = tmp;
3518//   tmp       = 0;
3519//   L[3]      = Lord;
3520//   // we are done with the list
3521//   def @R = ring(L);
3522//   setring @R;
3523//   matrix @D[Nnew][Nnew];
3524//   @D[3,4] = 1; // t,Dt
3525//   for(i=1; i<=N; i++)
3526//   {
3527//     @D[4+i,4+N+i]=1;
3528//   }
3529//   //  L[5] = matrix(UpOneMatrix(Nnew));
3530//   //  L[6] = @D;
3531//   nc_algebra(1,@D);
3532//   dbprint(ppl,"// -1-1- the ring @R is ready");
3533//   dbprint(ppl-1,@R);
3534//   // create the ideal
3535//   poly F  = imap(save,F);
3536//   ideal I = u*F-t,u*v-1;
3537//   poly p;
3538//   for(i=1; i<=N; i++)
3539//   {
3540//     p = u*Dt; // u*Dt
3541//     p = diff(F,var(4+i))*p;
3542//     I = I, var(N+4+i) + p; // Dx, Dy
3543//   }
3544//   // add the relations between t,Dt and s
3545//   //  I = I, t*Dt+1+S;
3546//   // -------- the ideal I is ready ----------
3547//   dbprint(ppl,"// -1-2- starting the elimination of u,v in @R");
3548//   ideal J = engine(I,eng);
3549//   ideal K = nselect(J,1,2);
3550//   dbprint(ppl,"// -1-3- u,v are eliminated in @R");
3551//   dbprint(ppl-1,K); // without u,v: not yet our answer
3552//   //----- create a ring with elim.ord for t,Dt -------
3553//   setring save;
3554//   // ------------ new ring @R2 ------------------
3555//   // without u,v and with the elim.ord for t,Dt
3556//   // keep: N, i,j,s, tmp, RL
3557//   Nnew = 2*N+2;
3558//   //  list RL = ringlist(save); // is defined earlier
3559//   kill Lord,tmp,iv, RName;
3560//   L = 0;
3561//   list Lord, tmp;
3562//   intvec iv;
3563//   L[1] = RL[1]; // char
3564//   L[4] = RL[4]; // char, minpoly
3565//   // check whether vars have admissible names -> done earlier
3566//   //  list Name  = RL[2];
3567//   list RName;
3568//   RName[1] = "t";
3569//   RName[2] = "Dt";
3570//   // DName is defined earlier
3571//   list NName = RName + Name + DName;
3572//   L[2]   = NName;
3573//   tmp    = 0;
3574//   // block ord (a(1,1),dp);
3575//   tmp[1]  = "a"; // string
3576//   iv      = 1,1;
3577//   tmp[2]  = iv; //intvec
3578//   Lord[1] = tmp;
3579//   // continue with dp 1,1,1,1...
3580//   tmp[1]  = "dp"; // string
3581//   s       = "iv=";
3582//   for(i=1;i<=Nnew;i++)
3583//   {
3584//     s = s+"1,";
3585//   }
3586//   s[size(s)]= ";";
3587//   execute(s);
3588//   kill s;
3589//   kill NName;
3590//   tmp[2]    = iv;
3591//   Lord[2]   = tmp;
3592//   tmp[1]    = "C";
3593//   iv        = 0;
3594//   tmp[2]    = iv;
3595//   Lord[3]   = tmp;
3596//   tmp       = 0;
3597//   L[3]      = Lord;
3598//   // we are done with the list
3599//   // Add: Plural part
3600//   def @R2 = ring(L);
3601//   setring @R2;
3602//   matrix @D[Nnew][Nnew];
3603//   @D[1,2]=1;
3604//   for(i=1; i<=N; i++)
3605//   {
3606//     @D[2+i,2+N+i]=1;
3607//   }
3608//   nc_algebra(1,@D);
3609//   dbprint(ppl,"// -2-1- the ring @R2 is ready");
3610//   dbprint(ppl-1,@R2);
3611//   ideal MM = maxideal(1);
3612//   MM = 0,0,MM;
3613//   map R01 = @R, MM;
3614//   ideal K2 = R01(K);
3615//   // add the relations between t,Dt and s
3616//   //  K2       = K2, t*Dt+1+S;
3617//   poly G = t*Dt+S+1;
3618//   K2 = NF(K2,std(G)),G;
3619//   dbprint(ppl,"// -2-2- starting elimination for t,Dt in @R2");
3620//   ideal J  = engine(K2,eng);
3621//   ideal K  = nselect(J,1,2);
3622//   dbprint(ppl,"// -2-3- t,Dt are eliminated");
3623//   dbprint(ppl-1,K);
3624//   //  "------- produce a final result --------";
3625//   // ----- create the ordinary Weyl algebra and put
3626//   // ----- the result into it
3627//   // ------ create the ring @R5
3628//   // keep: N, i,j,s, tmp, RL
3629//   setring save;
3630//   Nnew = 2*N;
3631//   //  list RL = ringlist(save); // is defined earlier
3632//   kill Lord, tmp, iv;
3633//   //  kill L;
3634//   L = 0;
3635//   list Lord, tmp;
3636//   intvec iv;
3637//   L[1] = RL[1]; // char
3638//   L[4] = RL[4]; // char, minpoly
3639//   // check whether vars have admissible names -> done earlier
3640//   //  list Name  = RL[2];
3641//   // DName is defined earlier
3642//   list NName = Name + DName;
3643//   L[2]   = NName;
3644//   // dp ordering;
3645//   string   s       = "iv=";
3646//   for(i=1;i<=2*N;i++)
3647//   {
3648//     s = s+"1,";
3649//   }
3650//   s[size(s)]= ";";
3651//   execute(s);
3652//   tmp     = 0;
3653//   tmp[1]  = "dp"; // string
3654//   tmp[2]  = iv; //intvec
3655//   Lord[1] = tmp;
3656//   kill s;
3657//   tmp[1]    = "C";
3658//   iv        = 0;
3659//   tmp[2]    = iv;
3660//   Lord[2]   = tmp;
3661//   tmp       = 0;
3662//   L[3]      = Lord;
3663//   // we are done with the list
3664//   // Add: Plural part
3665//   def @R5 = ring(L);
3666//   setring @R5;
3667//   matrix @D[Nnew][Nnew];
3668//   for(i=1; i<=N; i++)
3669//   {
3670//     @D[i,N+i]=1;
3671//   }
3672//   nc_algebra(1,@D);
3673//   dbprint(ppl,"// -3-1- the ring @R5 is ready");
3674//   dbprint(ppl-1,@R5);
3675//   ideal K5 = imap(@R2,K);
3676//   option(redSB);
3677//   dbprint(ppl,"// -3-2- the final cosmetic std");
3678//   K5 = engine(K5,eng); // std does the job too
3679//   // total cleanup
3680//   kill @R;
3681//   kill @R2;
3682//   ideal LD = K5;
3683//   ideal BS = imap(@LR,B);
3684//   kill @LR;
3685//   export BS;
3686//   export LD;
3687//   return(@R5);
3688// }
3689// example
3690// {
3691//   "EXAMPLE:"; echo = 2;
3692//   ring r = 0,(x,y,z),Dp;
3693//   poly F = x^2+y^3+z^5;
3694//   def A = annfsgms(F);
3695//   setring A;
3696//   LD;
3697//   print(matrix(BS));
3698// }
3699
3700
3701proc convloc(list @NL)
3702"USAGE:  convloc(L); L a list
3703RETURN:  list
3704PURPOSE: convert a ringlist L into another ringlist,
3705where all the 'p' orderings are replaced with the 's' orderings.
3706ASSUME: L is a result of a ringlist command
3707EXAMPLE: example convloc; shows examples
3708"
3709{
3710  list NL = @NL;
3711  // given a ringlist, returns a new ringlist,
3712  // where all the p-orderings are replaced with s-ord's
3713  int i,j,k,found;
3714  int nblocks = size(NL[3]);
3715  for(i=1; i<=nblocks; i++)
3716  {
3717    for(j=1; j<=size(NL[3][i]); j++)
3718    {
3719      if (typeof(NL[3][i][j]) == "string" )
3720      {
3721        found = 0;
3722        for (k=1; k<=size(NL[3][i][j]); k++)
3723        {
3724          if (NL[3][i][j][k]=="p")
3725          {
3726            NL[3][i][j][k]="s";
3727            found = 1;
3728            //    printf("replaced at %s,%s,%s",i,j,k);
3729          }
3730        }
3731      }
3732    }
3733  }
3734  return(NL);
3735}
3736example
3737{
3738  "EXAMPLE:"; echo = 2;
3739  ring r = 0,(x,y,z),(Dp(2),dp(1));
3740  list L = ringlist(r);
3741  list N = convloc(L);
3742  def rs = ring(N);
3743  setring rs;
3744  rs;
3745}
3746
3747proc annfspecial(ideal I, poly F, int mir, number n)
3748"USAGE:  annfspecial(I,F,mir,n);  I an ideal, F a poly, int mir, number n
3749RETURN:  ideal
3750PURPOSE: compute the annihilator ideal of F^n in the Weyl Algebra for a rational number n
3751ASSUME:  the basering contains 's' as a variable
3752NOTE:    We assume that the basering is D[s],
3753@*          ideal I is the Ann F^s in D[s] (obtained with e.g. SannfsBM, SannfsLOT, SannfsOT)
3754@*          integer 'mir' is the minimal integer root of the Bernstein polynomial of F
3755@*          and the number n is rational.
3756@*       We compute the real annihilator for any rational value of n (both generic and exceptional).
3757@*       The implementation goes along the lines of Saito-Sturmfels-Takayama, Alg. 5.3.15
3758@*       If printlevel=1, progress debug messages will be printed,
3759@*       if printlevel>=2, all the debug messages will be printed.
3760EXAMPLE: example annfspecial; shows examples
3761"
3762{
3763  int ppl = printlevel-voice+2;
3764  //  int mir =  minIntRoot(L[1],0);
3765  int ns   = varnum("s");
3766  if (!ns)
3767  {
3768    ERROR("Variable s expected in the ideal AnnFs");
3769  }
3770  int d;
3771  ideal P = subst(I,var(ns),n);
3772  if (denominator(n) == 1)
3773  {
3774    // n is fraction-free
3775    d = int(numerator(n));
3776    if ( (!d) && (n!=0))
3777    {
3778      ERROR("no parametric special values are allowed");
3779    }
3780    d = d - mir;
3781    if (d>0)
3782    {
3783      dbprint(ppl,"// -1-1- starting syzygy computations");
3784      matrix J[1][1] = F^d;
3785      dbprint(ppl-1,"// -1-1-1- of the polynomial ideal");
3786      dbprint(ppl-1,J);
3787      matrix K[1][size(I)] = subst(I,var(ns),mir);
3788      dbprint(ppl-1,"// -1-1-2- modulo ideal:");
3789      dbprint(ppl-1, K);
3790      module M = modulo(J,K);
3791      dbprint(ppl-1,"// -1-1-3- getting the result:");
3792      dbprint(ppl-1, M);
3793      P  = P, ideal(M);
3794      dbprint(ppl,"// -1-2- finished syzygy computations");
3795    }
3796    else
3797    {
3798      dbprint(ppl,"// -1-1- d<=0, no syzygy computations needed");
3799      dbprint(ppl-1,"// -1-2- for d =");
3800      dbprint(ppl-1, d);
3801    }
3802  }
3803  // also the else case: d<=0 or n is not an integer
3804  dbprint(ppl,"// -2-1- starting final Groebner basis");
3805  P = groebner(P);
3806  dbprint(ppl,"// -2-2- finished final Groebner basis");
3807  return(P);
3808}
3809example
3810{
3811  "EXAMPLE:"; echo = 2;
3812  ring r = 0,(x,y),dp;
3813  poly F = x3-y2;
3814  def  B = annfs(F);  setring B;
3815  minIntRoot(BS[1],0);
3816  // So, the minimal integer root is -1
3817  setring r;
3818  def  A = SannfsBM(F);
3819  setring A;
3820  poly F = x3-y2;
3821  annfspecial(LD,F,-1,3/4); // generic root
3822  annfspecial(LD,F,-1,-2);  // integer but still generic root
3823  annfspecial(LD,F,-1,1);   // exceptional root
3824}
3825
3826static proc new_ex_annfspecial()
3827{
3828  //another example for annfspecial: x3+y3+z3
3829  ring r = 0,(x,y,z),dp;
3830  poly F =  x3+y3+z3;
3831  def  B = annfs(F);  setring B;
3832  minIntRoot(BS[1],0);
3833  // So, the minimal integer root is -1
3834  setring r;
3835  def  A = SannfsBM(F);
3836  setring A;
3837  poly F =  x3+y3+z3;
3838  annfspecial(LD,F,-1,3/4); // generic root
3839  annfspecial(LD,F,-1,-2);  // integer but still generic root
3840  annfspecial(LD,F,-1,1);   // exceptional root
3841}
3842
3843proc minIntRoot(ideal P, int fact)
3844"USAGE:  minIntRoot(P, fact); P an ideal, fact an int
3845RETURN:  int
3846PURPOSE: minimal integer root of a maximal ideal P
3847NOTE:    if fact==1, P is the result of some 'factorize' call,
3848@*       else P is treated as the result of bernstein::gmssing.lib
3849@*       in both cases without constants and multiplicities
3850EXAMPLE: example minIntRoot; shows examples
3851"
3852{
3853  //    ideal P = factorize(p,1);
3854  // or ideal P = bernstein(F)[1];
3855  intvec vP;
3856  number nP;
3857  int i;
3858  if ( fact )
3859  {
3860    // the result comes from "factorize"
3861    P = normalize(P); // now leadcoef = 1
3862    // TODO: hunt for units and kill then !!!
3863    P = lead(P)-P;
3864    // nP = leadcoef(P[i]-lead(P[i])); // for 1 var only, extract the coeff
3865  }
3866  else
3867  {
3868    // bernstein deletes -1 usually
3869    P = P, -1;
3870  }
3871  // for both situations:
3872  // now we have an ideal of fractions of type "number"
3873  int sP = size(P);
3874  for(i=1; i<=sP; i++)
3875  {
3876    nP = leadcoef(P[i]);
3877    if ( (nP - int(nP)) == 0 )
3878    {
3879      vP = vP,int(nP);
3880    }
3881  }
3882  if ( size(vP)>=2 )
3883  {
3884    vP = vP[2..size(vP)];
3885  }
3886  sP = -Max(-vP);
3887  if (sP == 0)
3888  {
3889    "Warning: zero root!";
3890  }
3891  return(sP);
3892}
3893example
3894{
3895  "EXAMPLE:"; echo = 2;
3896  ring r   = 0,(x,y),ds;
3897  poly f1  = x*y*(x+y);
3898  ideal I1 = bernstein(f1)[1]; // a local Bernstein poly
3899  minIntRoot(I1,0);
3900  poly  f2  = x2-y3;
3901  ideal I2  = bernstein(f2)[1];
3902  minIntRoot(I2,0);
3903  // now we illustrate the behaviour of factorize
3904  // together with a global ordering
3905  ring r2  = 0,x,dp;
3906  poly f3   = 9*(x+2/3)*(x+1)*(x+4/3); //global b-poly of f1=x*y*(x+y)
3907  ideal I3 = factorize(f3,1);
3908  minIntRoot(I3,1);
3909  // and a more interesting situation
3910  ring  s  = 0,(x,y,z),ds;
3911  poly  f  = x3 + y3 + z3;
3912  ideal I  = bernstein(f)[1];
3913  minIntRoot(I,0);
3914}
3915
3916proc isHolonomic(def M)
3917"USAGE:  isHolonomic(M); M an ideal/module/matrix
3918RETURN:  int, 1 if M is holonomic and 0 otherwise
3919PURPOSE: check the modules for the property of holonomy
3920NOTE:    M is holonomic, if 2*dim(M) = dim(R), where R is a
3921ground ring; dim stays for Gelfand-Kirillov dimension
3922EXAMPLE: example isHolonomic; shows examples
3923"
3924{
3925  if ( (typeof(M) != "ideal") && (typeof(M) != "module") && (typeof(M) != "matrix") )
3926  {
3927  //  print(typeof(M));
3928    ERROR("an argument of type ideal/module/matrix expected");
3929  }
3930  if (attrib(M,"isSB")!=1)
3931  {
3932    M = std(M);
3933  }
3934  int dimR = gkdim(std(0));
3935  int dimM = gkdim(M);
3936  return( (dimR==2*dimM) );
3937}
3938example
3939{
3940  "EXAMPLE:"; echo = 2;
3941  ring R = 0,(x,y),dp;
3942  poly F = x*y*(x+y);
3943  def A = annfsBM(F,0);
3944  setring A;
3945  LD;
3946  isHolonomic(LD);
3947  ideal I = std(LD[1]);
3948  I;
3949  isHolonomic(I);
3950}
3951
3952proc reiffen(int p, int q)
3953"USAGE:  reiffen(p, q);  int p, int q
3954RETURN:  ring
3955PURPOSE: set up the polynomial, describing a Reiffen curve
3956NOTE:    activate this ring with the @code{setring} command and find the
3957         curve as a polynomial RC
3958@*   a Reiffen curve is defined as F = x^p + y^q + xy^{q-1}, q >= p+1 >= 5
3959ASSUME: q >= p+1 >= 5. Otherwise an error message is returned
3960EXAMPLE: example reiffen; shows examples
3961"
3962{
3963// a Reiffen curve is defined as
3964// F = x^p + y^q +x*y^{q-1}, q \geq p+1 \geq 5
3965
3966  if ( (p<4) || (q<5) || (q-p<1) )
3967  {
3968    ERROR("Some of conditions p>=4, q>=5 or q>=p+1 is not satisfied!");
3969  }
3970  ring @r = 0,(x,y),dp;
3971  poly RC = y^q +x^p + x*y^(q-1);
3972  export RC;
3973  return(@r);
3974}
3975example
3976{
3977  "EXAMPLE:"; echo = 2;
3978  def r = reiffen(4,5);
3979  setring r;
3980  RC;
3981}
3982
3983proc arrange(int p)
3984"USAGE:  arrange(p);  int p
3985RETURN:  ring
3986PURPOSE: set up the polynomial, describing a hyperplane arrangement
3987NOTE:   must be executed in a ring
3988ASSUME: basering is present
3989EXAMPLE: example arrange; shows examples
3990"
3991{
3992  int UseBasering = 0 ;
3993  if (p<2)
3994  {
3995    ERROR("p>=2 is required!");
3996  }
3997  if ( nameof(basering)!="basering" )
3998  {
3999    if (p > nvars(basering))
4000    {
4001      ERROR("too big p");
4002    }
4003    else
4004    {
4005      def @r = basering;
4006      UseBasering = 1;
4007    }
4008  }
4009  else
4010  {
4011    // no basering found
4012    ERROR("no basering found!");
4013    //    ring @r = 0,(x(1..p)),dp;
4014  }
4015  int i,j,sI;
4016  poly  q = 1;
4017  list ar;
4018  ideal tmp;
4019  tmp = ideal(var(1));
4020  ar[1] = tmp;
4021  for (i = 2; i<=p; i++)
4022  {
4023    // add i-nary sums to the product
4024    ar = indAR(ar,i);
4025  }
4026  for (i = 1; i<=p; i++)
4027  {
4028    tmp = ar[i]; sI = size(tmp);
4029    for (j = 1; j<=sI; j++)
4030    {
4031      q = q*tmp[j];
4032    }
4033  }
4034  if (UseBasering)
4035  {
4036    return(q);
4037  }
4038  //  poly AR = q; export AR;
4039  //  return(@r);
4040}
4041example
4042{
4043  "EXAMPLE:"; echo = 2;
4044  ring X = 0,(x,y,z,t),dp;
4045  poly q = arrange(3);
4046  factorize(q,1);
4047}
4048
4049proc checkRoot(poly F, number a, list #)
4050"USAGE:  checkRoot(f,alpha [,S,eng]);  f a poly, alpha a number, S a string , eng an optional int
4051RETURN:  int
4052PURPOSE: check whether a rational is a root of the global Bernstein polynomial of f (and compute its multiplicity)
4053with the algorithm given in S and with the Groebner basis engine given in eng
4054NOTE:    The annihilator of f^s in D[s] is needed and it is computed according to the algorithm by Briancon and Maisonobe
4055@*       The value of a string S can be
4056@*       'alg1' (default) - for the algorithm 1 of J. Martin-Morales (unpublished)
4057@*       'alg2' - for the algorithm 2 of J. Martin-Morales (unpublished)
4058@*       The output int is:
4059@*       - if the algorithm 1 is chosen: 1 if -alpha is a root of the global Bernstein polynomial and 0 otherwise
4060@*       - if the algorithm 2 is chosen: the multiplicity of -alpha as a root of the global Bernstein polynomial of f.
4061@*                                       (If -alpha is not a root, the output is 0)
4062@*       If eng <>0, @code{std} is used for Groebner basis computations,
4063@*       otherwise (and by default) @code{slimgb} is used.
4064@*       If printlevel=1, progress debug messages will be printed,
4065@*       if printlevel>=2, all the debug messages will be printed.
4066EXAMPLE: example checkRoot; shows examples
4067"
4068{
4069  int eng = 0;
4070  int chs = 0; // choice
4071  string algo = "alg1";
4072  string st;
4073  if ( size(#)>0 )
4074  {
4075   if ( typeof(#[1]) == "string" )
4076   {
4077     st = string(#[1]);
4078     if ( (st == "alg1") || (st == "ALG1") || (st == "Alg1") ||(st == "aLG1"))
4079     {
4080       algo = "alg1";
4081       chs  = 1;
4082     }
4083     if ( (st == "alg2") || (st == "ALG2") || (st == "Alg2") || (st == "aLG2"))
4084     {
4085       algo = "alg2";
4086       chs  = 1;
4087     }
4088     if (chs != 1)
4089     {
4090       // incorrect value of S
4091       print("Incorrect algorithm given, proceed with the default alg1 of J. Martín-Morales");
4092       algo = "alg1";
4093     }
4094     // second arg
4095     if (size(#)>1)
4096     {
4097       // exists 2nd arg
4098       if ( typeof(#[2]) == "int" )
4099       {
4100         // the case: given alg, given engine
4101         eng = int(#[2]);
4102       }
4103       else
4104       {
4105         eng = 0;
4106       }
4107     }
4108     else
4109     {
4110       // no second arg
4111       eng = 0;
4112     }
4113   }
4114   else
4115   {
4116     if ( typeof(#[1]) == "int" )
4117     {
4118       // the case: default alg, engine
4119       eng = int(#[1]);
4120       // algo = "alg1";  //is already set
4121     }
4122     else
4123     {
4124       // incorr. 1st arg
4125       algo = "alg1";
4126     }
4127   }
4128  }
4129  // size(#)=0, i.e. there is no algorithm and no engine given
4130  //  eng = 0; algo = "alg1";  //are already set
4131  // int ppl = printlevel-voice+2;
4132  printlevel=printlevel+1;
4133  def save = basering;
4134  def @A = SannfsBM(F);
4135  setring @A;
4136  poly F = imap(save,F);
4137  number a = imap(save,a);
4138  if ( algo=="alg1")
4139  {
4140    int output = checkRoot1(LD,F,a,eng);
4141  }
4142  else
4143  {
4144    if ( algo=="alg2")
4145    {
4146      int output = checkRoot2(LD,F,a,eng);
4147    }
4148  }
4149  printlevel=printlevel-1;
4150  return(output);
4151}
4152example
4153{
4154  "EXAMPLE:"; echo = 2;
4155  printlevel=0;
4156  ring r = 0,(x,y),Dp;
4157  poly F = x^4+y^5+x*y^4;
4158  checkRoot(F,11/20);    // -11/20 is a root of bf
4159  poly G = x*y;
4160  checkRoot(G,1,"alg2"); // -1 is a root of bg with multiplicity 2
4161}
4162
4163proc checkRoot1(ideal I, poly F, number a, list #)
4164"USAGE:  checkRoot1(I,f,alpha [,eng]);  I an ideal, f a poly, alpha a number, eng an optional int
4165ASSUME:  I is the annihilator of f^s in D[s], f is a polynomial in K[_x]
4166RETURN:  int, 1 if -alpha is a root of the global Bernstein polynomial of f and 0 otherwise
4167PURPOSE: check whether a rational is a root of the global Bernstein polynomial of f
4168NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
4169@*       otherwise (and by default) @code{slimgb} is used.
4170@*       If printlevel=1, progress debug messages will be printed,
4171@*       if printlevel>=2, all the debug messages will be printed.
4172EXAMPLE: example checkRoot1; shows examples
4173"
4174{
4175  int eng = 0;
4176  if ( size(#)>0 )
4177  {
4178    if ( typeof(#[1]) == "int" )
4179    {
4180      eng = int(#[1]);
4181    }
4182  }
4183  int ppl = printlevel-voice+2;
4184  dbprint(ppl,"// -0-1- starting the procedure checkRoot1");
4185  def save = basering;
4186  int N = nvars(basering);
4187  int Nnew = N-1;
4188  int n = Nnew / 2;
4189  int i;
4190  string s;
4191  list RL = ringlist(basering);
4192  list L, Lord;
4193  list tmp;
4194  intvec iv;
4195  L[1] = RL[1];  // char
4196  L[4] = RL[4];  // char, minpoly
4197  // check whether basering is D[s]=K(_x,_Dx,s)
4198  list Name = RL[2];
4199  for (i=1; i<=n; i++)
4200  {
4201    if ( bracket(var(i+n),var(i))!=1 )
4202    {
4203      ERROR("basering should be D[s]=K(_x,_Dx,s)");
4204    }
4205  }
4206  if ( Name[N]!="s" )
4207  {
4208    ERROR("the last variable of basering should be s");
4209  }
4210  // now, create the new vars
4211  list NName;
4212  for (i=1; i<=Nnew; i++)
4213  {
4214    NName[i] = Name[i];
4215  }
4216  L[2] = NName;
4217  kill Name,NName;
4218  // block ord (dp);
4219  tmp[1] = "dp"; // string
4220  s = "iv=";
4221  for (i=1; i<=Nnew; i++)
4222  {
4223    s = s+"1,";
4224  }
4225  s[size(s)]=";";
4226  execute(s);
4227  kill s;
4228  tmp[2]  = iv;
4229  Lord[1] = tmp;
4230  tmp[1]  = "C";
4231  iv      = 0;
4232  tmp[2]  = iv;
4233  Lord[2] = tmp;
4234  tmp     = 0;
4235  L[3]    = Lord;
4236  // we are done with the list
4237  def @R@ = ring(L);
4238  setring @R@;
4239  matrix @D[Nnew][Nnew];
4240  for (i=1; i<=n; i++)
4241  {
4242    @D[i,i+n]=1;
4243  }
4244  def @R = nc_algebra(1,@D);
4245  setring @R;
4246  kill @R@;
4247  dbprint(ppl,"// -1-1- the ring @R(_x,_Dx) is ready");
4248  dbprint(ppl-1, S);
4249  // create the ideal K = ann_D[s](f^s)_{s=-alpha} + < f >
4250  setring save;
4251  ideal K = subst(I,s,-a);
4252  dbprint(ppl,"// -1-2- the variable s has been substituted by "+string(-a));
4253  dbprint(ppl-1, K);
4254  K = NF(K,std(F));
4255  // make leadcoeffs positive
4256  for (i=1; i<=ncols(K); i++)
4257  {
4258    if ( leadcoef(K[i])<0 )
4259    {
4260      K[i] = -K[i];
4261    }
4262  }
4263  K = K,F;
4264  // ------------ the ideal K is ready ------------
4265  setring @R;
4266  ideal K = imap(save,K);
4267  dbprint(ppl,"// -1-3- starting the computation of a Groebner basis of K in @R");
4268  dbprint(ppl-1, K);
4269  ideal G = engine(K,eng);
4270  dbprint(ppl,"// -1-4- the Groebner basis has been computed");
4271  dbprint(ppl-1, G);
4272  return(G[1]!=1);
4273}
4274example
4275{
4276  "EXAMPLE:"; echo = 2;
4277  ring r = 0,(x,y),Dp;
4278  poly F = x^4+y^5+x*y^4;
4279  printlevel = 0;
4280  def A = Sannfs(F);
4281  setring A;
4282  poly F = imap(r,F);
4283  checkRoot1(LD,F,31/20);   // -31/20 is not a root of bs
4284  checkRoot1(LD,F,11/20);   // -11/20 is a root of bs
4285}
4286
4287proc checkRoot2 (ideal I, poly F, number a, list #)
4288"USAGE:  checkRoot2(I,f,alpha [,eng]);  I an ideal, f a poly, alpha a number, eng an optional int
4289ASSUME:  I is the annihilator of f^s in D[s], f is a polynomial in K[_x]
4290RETURN:  int, the multiplicity of -alpha as a root of the global Bernstein polynomial of f. If -alpha is not a root, the output is 0
4291PURPOSE: check whether a rational is a root of the global Bernstein polynomial of f and compute its multiplicity from the known Ann F^s in D[s]
4292NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
4293@*       otherwise (and by default) @code{slimgb} is used.
4294@*       If printlevel=1, progress debug messages will be printed,
4295@*       if printlevel>=2, all the debug messages will be printed.
4296EXAMPLE: example checkRoot2; shows examples
4297"
4298{
4299  int eng = 0;
4300  if ( size(#)>0 )
4301  {
4302    if ( typeof(#[1]) == "int" )
4303    {
4304      eng = int(#[1]);
4305    }
4306  }
4307  int ppl = printlevel-voice+2;
4308  dbprint(ppl,"// -0-1- starting the procedure checkRoot2");
4309  def save = basering;
4310  int N = nvars(basering);
4311  int n = (N-1) / 2;
4312  int i;
4313  string s;
4314  list RL = ringlist(basering);
4315  list L, Lord;
4316  list tmp;
4317  intvec iv;
4318  L[1] = RL[1];  // char
4319  L[4] = RL[4];  // char, minpoly
4320  // check whether basering is D[s]=K(_x,_Dx,s)
4321  list Name = RL[2];
4322  for (i=1; i<=n; i++)
4323  {
4324    if ( bracket(var(i+n),var(i))!=1 )
4325    {
4326      ERROR("basering should be D[s]=K(_x,_Dx,s)");
4327    }
4328  }
4329  if ( Name[N]!="s" )
4330  {
4331    ERROR("the last variable of basering should be s");
4332  }
4333  // now, create the new vars
4334  L[2] = Name;
4335  kill Name;
4336  // block ord (dp);
4337  tmp[1] = "dp"; // string
4338  s = "iv=";
4339  for (i=1; i<=N; i++)
4340  {
4341    s = s+"1,";
4342  }
4343  s[size(s)]=";";
4344  execute(s);
4345  kill s;
4346  tmp[2]  = iv;
4347  Lord[1] = tmp;
4348  tmp[1]  = "C";
4349  iv      = 0;
4350  tmp[2]  = iv;
4351  Lord[2] = tmp;
4352  tmp     = 0;
4353  L[3]    = Lord;
4354  // we are done with the list
4355  def @R@ = ring(L);
4356  setring @R@;
4357  matrix @D[N][N];
4358  for (i=1; i<=n; i++)
4359  {
4360    @D[i,i+n]=1;
4361  }
4362  def @R = nc_algebra(1,@D);
4363  setring @R;
4364  kill @R@;
4365  dbprint(ppl,"// -1-1- the ring @R(_x,_Dx,s) is ready");
4366  dbprint(ppl-1, @R);
4367  // now, continue with the algorithm
4368  ideal I = imap(save,I);
4369  poly F = imap(save,F);
4370  number a = imap(save,a);
4371  ideal II = NF(I,std(F));
4372  // make leadcoeffs positive
4373  for (i=1; i<=ncols(II); i++)
4374  {
4375    if ( leadcoef(II[i])<0 )
4376    {
4377      II[i] = -II[i];
4378    }
4379  }
4380  ideal J,G;
4381  int m;  // the output (multiplicity)
4382  dbprint(ppl,"// -2- starting the bucle");
4383  for (i=0; i<=n; i++)  // the multiplicity has to be <= n
4384  {
4385    // create the ideal Ji = ann_D[s](f^s) + < f, (s+alpha)^{i+1} >
4386    // (s+alpha)^i in Ji <==> -alpha is a root with multiplicity >= i
4387    J = II,F,(s+a)^(i+1);
4388    // ------------ the ideal Ji is ready -----------
4389    dbprint(ppl,"// -2-"+string(i+1)+"-1- starting the computation of a Groebner basis of J"+string(i)+" in @R");
4390    dbprint(ppl-1, J);
4391    G = engine(J,eng);
4392    dbprint(ppl,"// -2-"+string(i+1)+"-2- the Groebner basis has been computed");
4393    dbprint(ppl-1, G);
4394    if ( NF((s+a)^i,G)==0 )
4395    {
4396      dbprint(ppl,"// -2-"+string(i+1)+"-3- the number "+string(-a)+" has not multiplicity "+string(i+1));
4397      m = i;
4398      break;
4399    }
4400    dbprint(ppl,"// -2-"+string(i+1)+"-3- the number "+string(-a)+" has multiplicity at least "+string(i+1));
4401  }
4402  dbprint(ppl,"// -3- the bucle has finished");
4403  return(m);
4404}
4405example
4406{
4407  "EXAMPLE:"; echo = 2;
4408  ring r = 0,(x,y,z),Dp;
4409  poly F = x*y*z;
4410  printlevel = 0;
4411  def A = Sannfs(F);
4412  setring A;
4413  poly F = imap(r,F);
4414  checkRoot2(LD,F,1);    // -1 is a root of bs with multiplicity 3
4415  checkRoot2(LD,F,1/3);  // -1/3 is not a root
4416}
4417
4418proc checkFactor(ideal I, poly F, poly q, list #)
4419"USAGE:  checkFactor(I,f,qs [,eng]);  I an ideal, f a poly, qs a poly, eng an optional int
4420ASSUME:  I is the output of Sannfs, SannfsBM, SannfsLOT or SannfsOT,
4421         f is a polynomial in K[_x], qs is a polynomial in K[s]
4422RETURN:  int, 1 if qs is a factor of the global Bernstein polynomial of f and 0 otherwise
4423PURPOSE: check whether a rational is a root of the global Bernstein polynomial of f and compute its multiplicity from the known Ann F^s in D[s]
4424NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
4425@*       otherwise (and by default) @code{slimgb} is used.
4426@*       If printlevel=1, progress debug messages will be printed,
4427@*       if printlevel>=2, all the debug messages will be printed.
4428EXAMPLE: example checkFactor; shows examples
4429"
4430{
4431  int eng = 0;
4432  if ( size(#)>0 )
4433  {
4434    if ( typeof(#[1]) == "int" )
4435    {
4436      eng = int(#[1]);
4437    }
4438  }
4439  int ppl = printlevel-voice+2;
4440  def @R2 = basering;
4441  int N = nvars(@R2);
4442  int i;
4443  // we're in D_n[s], where the elim ord for s is set
4444  dbprint(ppl,"// -0-1- starting the procedure checkFactor");
4445  dbprint(ppl,"// -1-1- the ring @R2(_x,_Dx,s) is ready");
4446  dbprint(ppl-1, @R2);
4447  // create the ideal J = ann_D[s](f^s) + < f,q >
4448  ideal J = NF(I,std(F));
4449  // make leadcoeffs positive
4450  for (i=1; i<=ncols(J); i++)
4451  {
4452    if ( leadcoef(J[i])<0 )
4453    {
4454      J[i] = -J[i];
4455    }
4456  }
4457  J = J,F,q;
4458  // ------------ the ideal J is ready -----------
4459  dbprint(ppl,"// -1-2- starting the elimination of _x,_Dx in @R2");
4460  dbprint(ppl-1, J);
4461  ideal G = engine(J,eng);
4462  ideal K = nselect(G,1,N-1);
4463  kill J,G;
4464  dbprint(ppl,"// -1-3- _x,_Dx are eliminated");
4465  dbprint(ppl-1, K);
4466  //q is a factor of bs iff K = < q >
4467  //K = normalize(K);
4468  //q = normalize(q);
4469  //return( (K[1]==q) );
4470  return( NF(K[1],std(q))==0 );
4471}
4472example
4473{
4474  "EXAMPLE:"; echo = 2;
4475  ring r = 0,(x,y),Dp;
4476  poly F = x^4+y^5+x*y^4;
4477  printlevel = 0;
4478  def A = Sannfs(F);
4479  setring A;
4480  poly F = imap(r,F);
4481  checkFactor(LD,F,20*s+31);     // -31/20 is not a root of bs
4482  checkFactor(LD,F,20*s+11);     // -11/20 is a root of bs
4483  checkFactor(LD,F,(20*s+11)^2); // the multiplicity of -11/20 is 1
4484}
4485
4486proc varnum(string s)
4487"USAGE:  varnum(s);  string s
4488RETURN:  int
4489PURPOSE: returns the number of the variable with the name s
4490among the variables of basering or 0 if there is no such variable
4491EXAMPLE: example varnum; shows examples
4492"
4493{
4494  int i;
4495  for (i=1; i<= nvars(basering); i++)
4496  {
4497    if ( string(var(i)) == s )
4498    {
4499      return(i);
4500    }
4501  }
4502  return(0);
4503}
4504example
4505{
4506  "EXAMPLE:"; echo = 2;
4507  ring X = 0,(x,y1,z(0),tTa),dp;
4508  varnum("z(0)");
4509  varnum("tTa");
4510  varnum("xyz");
4511}
4512
4513static proc indAR(list L, int n)
4514"USAGE:  indAR(L,n);  list L, int n
4515RETURN:  list
4516PURPOSE: computes arrangement inductively, using L and var(n) as the
4517next variable
4518ASSUME: L has a structure of an arrangement
4519EXAMPLE: example indAR; shows examples
4520"
4521{
4522  if ( (n<2) || (n>nvars(basering)) )
4523  {
4524    ERROR("incorrect n");
4525  }
4526  int sl = size(L);
4527  list K;
4528  ideal tmp;
4529  poly @t = L[sl][1] + var(n); //1 elt
4530  K[sl+1] = ideal(@t);
4531  tmp = L[1]+var(n);
4532  K[1] = tmp; tmp = 0;
4533  int i,j,sI;
4534  ideal I;
4535  for(i=sl; i>=2; i--)
4536  {
4537    I = L[i-1]; sI = size(I);
4538    for(j=1; j<=sI; j++)
4539    {
4540      I[j] = I[j] + var(n);
4541    }
4542    tmp  = L[i],I;
4543    K[i] = tmp;
4544    I = 0; tmp = 0;
4545  }
4546  kill I; kill tmp;
4547  return(K);
4548}
4549example
4550{
4551  "EXAMPLE:"; echo = 2;
4552  ring r = 0,(x,y,z,t,v),dp;
4553  list L;
4554  L[1] = ideal(x);
4555  list K = indAR(L,2);
4556  K;
4557  list M = indAR(K,3);
4558  M;
4559  M = indAR(M,4);
4560  M;
4561}
4562
4563static proc exCheckGenericity()
4564{
4565  LIB "control.lib";
4566  ring r = (0,a,b,c),x,dp;
4567  poly p = (x-a)*(x-b)*(x-c);
4568  def A = annfsBM(p);
4569  setring A;
4570  ideal J = slimgb(LD);
4571  matrix T = lift(LD,J);
4572  T = normalize(T);
4573  genericity(T);
4574  // Ann =x^3*Dx+3*x^2*t*Dt+(-a-b-c)*x^2*Dx+(-2*a-2*b-2*c)*x*t*Dt+3*x^2+(a*b+a*c+b*c)*x*Dx+(a*b+a*c+b*c)*t*Dt+(-2*a-2*b-2*c)*x+(-a*b*c)*Dx+(a*b+a*c+b*c)
4575  // genericity: g = a2-ab-ac+b2-bc+c2 =0
4576  // g = (a -(b+c)/2)^2 + (3/4)*(b-c)^2;
4577  // g ==0 <=> a=b=c
4578  // indeed, Ann = (x-a)^2*(x*Dx+3*t*Dt+(-a)*Dx+3)
4579  // --------------------------------------------
4580  // BUT a direct computation shows
4581  // when a=b=c,
4582  // Ann = x*Dx+3*t*Dt+(-a)*Dx+3
4583}
4584
4585static proc exOT_17()
4586{
4587  // Oaku-Takayama, p.208
4588  ring R = 0,(x,y),dp;
4589  poly F = x^3-y^2; // x^2+x*y+y^2;
4590  option(prot);
4591  option(mem);
4592  //  option(redSB);
4593  def A = annfsOT(F,0);
4594  setring A;
4595  LD;
4596  gkdim(LD); // a holonomic check
4597  //  poly F = x^3-y^2; // = x^7 - y^5; // x^3-y^4; // x^5 - y^4;
4598}
4599
4600static proc exOT_16()
4601{
4602  // Oaku-Takayama, p.208
4603  ring R = 0,(x),dp;
4604  poly F = x*(1-x);
4605  option(prot);
4606  option(mem);
4607  //  option(redSB);
4608  def A = annfsOT(F,0);
4609  setring A;
4610  LD;
4611  gkdim(LD); // a holonomic check
4612}
4613
4614static proc ex_bcheck()
4615{
4616  ring R = 0,(x,y),dp;
4617  poly F = x*y*(x+y);
4618  option(prot);
4619  option(mem);
4620  int eng = 0;
4621  //  option(redSB);
4622  def A = annfsOT(F,eng);
4623  setring A;
4624  LD;
4625}
4626
4627static proc ex_bcheck2()
4628{
4629  ring R = 0,(x,y),dp;
4630  poly F = x*y*(x+y);
4631  int eng = 0;
4632  def A = annfsBM(F,eng);
4633  setring A;
4634  LD;
4635}
4636
4637static proc ex_BMI()
4638{
4639  // a hard example
4640  ring r = 0,(x,y),Dp;
4641  poly F1 = (x2-y3)*(x3-y2);
4642  poly F2 = (x2-y3)*(xy4+y5+x4);
4643  ideal F = F1,F2;
4644  def A = annfsBMI(F);
4645  setring A;
4646  LD;
4647  BS;
4648}
4649
4650static proc ex2_BMI()
4651{
4652  // this example was believed to be intractable in 2005 by Gago-Vargas, Castro and Ucha
4653  ring r = 0,(x,y),Dp;
4654  option(prot);
4655  option(mem);
4656  ideal F = x2+y3,x3+y2;
4657  printlevel = 2;
4658  def A = annfsBMI(F);
4659  setring A;
4660  LD;
4661  BS;
4662}
4663
4664static proc ex_operatorBM()
4665{
4666  ring r = 0,(x,y,z,w),Dp;
4667  poly F = x^3+y^3+z^2*w;
4668  printlevel = 0;
4669  def A = operatorBM(F);
4670  setring A;
4671  F; // the original polynomial itself
4672  LD; // generic annihilator
4673  LD0; // annihilator
4674  bs; // normalized Bernstein poly
4675  BS; // root and multiplicities of the Bernstein poly
4676  PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD
4677  reduce(PS*F-bs,LD); // check the property of PS
4678}
Note: See TracBrowser for help on using the repository browser.