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