# source:git/Singular/LIB/dmod.lib@e42755

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