# source:git/Singular/LIB/dmod.lib@565e86

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