source: git/Singular/LIB/dmod.lib @ 85ba0a

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