# source:git/Singular/LIB/dmod.lib@3cec5e

spielwiese
Last change on this file since 3cec5e was 3cec5e, checked in by Viktor Levandovskyy <levandov@…>, 15 years ago
*levandov: fixes and doc improvements, bfun instead of bfct_lib git-svn-id: file:///usr/local/Singular/svn/trunk@11383 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100755`
File size: 146.5 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="\$Id: dmod.lib,v 1.35 2009-02-12 20:25:22 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 */
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);
486  for (i=1; i<= ncols(J); i++)
487  {
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
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
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));
1253  int i;
1254  J = subst(J,s,-s-1);
1255  for (i=1; i<= ncols(J); i++)
1256  {
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
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);
1449  J = subst(J,s,-s-1);
1450  for (i=1; i<= ncols(J); i++)
1451  {
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
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
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 1nd column of transp. matrix
1890// this does not happen automatically
1891// for this, do special modulo with emphasis on the 1comp (c,<)
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 i.e. kernel of module homomorphism
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];
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
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
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;
3180  for (i=1; i<= ncols(LD); i++)
3181  {
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;
3386  for (i=1; i<= ncols(LD); i++)
3387  {
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 ann f^s and Groebner basis of ann f^s+f in D[s]
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 a Groebner basis) is the needed ideal.
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  // DEBUG INFO: ordering on the output ring = dp, engine on the sum of ideals is used
3423  // hence it's the situation for slimgb
3424
3425  int eng = 0;
3426  if ( size(#)>0 )
3427  {
3428    if ( typeof(#[1]) == "int" )
3429    {
3430      eng = int(#[1]);
3431    }
3432  }
3433  // returns a list with a ring and an ideal LD in it
3434  int ppl = printlevel-voice+2;
3435  //  printf("plevel :%s, voice: %s",printlevel,voice);
3436  def save = basering;
3437  int N = nvars(basering);
3438  int Nnew = 2*N+2;
3439  int i,j;
3440  string s;
3441  list RL = ringlist(basering);
3442  list L, Lord;
3443  list tmp;
3444  intvec iv;
3445  L[1] = RL[1]; // char
3446  L[4] = RL[4]; // char, minpoly
3447  // check whether vars have admissible names
3448  list Name  = RL[2];
3449  list RName;
3450  RName[1] = "@t";
3451  RName[2] = "@s";
3452  for(i=1;i<=N;i++)
3453  {
3454    for(j=1; j<=size(RName);j++)
3455    {
3456      if (Name[i] == RName[j])
3457      {
3458        ERROR("Variable names should not include @t,@s");
3459      }
3460    }
3461  }
3462  // now, create the names for new vars
3463  list DName;
3464  for(i=1;i<=N;i++)
3465  {
3466    DName[i] = "D"+Name[i]; // concat
3467  }
3468  tmp[1] = "t";
3469  tmp[2] = "s";
3470  list NName = tmp + DName + Name ;
3471  L[2]   = NName;
3472  // Name, Dname will be used further
3473  kill NName;
3474  // block ord (lp(2),dp);
3475  tmp[1]  = "lp"; // string
3476  iv      = 1,1;
3477  tmp[2]  = iv; //intvec
3478  Lord[1] = tmp;
3479  // continue with dp 1,1,1,1...
3480  tmp[1]  = "dp"; // string
3481  s       = "iv=";
3482  for(i=1;i<=Nnew;i++)
3483  {
3484    s = s+"1,";
3485  }
3486  s[size(s)]= ";";
3487  execute(s);
3488  kill s;
3489  tmp[2]    = iv;
3490  Lord[2]   = tmp;
3491  tmp[1]    = "C";
3492  iv        = 0;
3493  tmp[2]    = iv;
3494  Lord[3]   = tmp;
3495  tmp       = 0;
3496  L[3]      = Lord;
3497  // we are done with the list
3498  def @R@ = ring(L);
3499  setring @R@;
3500  matrix @D[Nnew][Nnew];
3501  @D[1,2]=t;
3502  for(i=1; i<=N; i++)
3503  {
3504    @D[2+i,N+2+i]=-1;
3505  }
3506  //  L[5] = matrix(UpOneMatrix(Nnew));
3507  //  L[6] = @D;
3508  def @R = nc_algebra(1,@D);
3509  setring @R;
3510  kill @R@;
3511  dbprint(ppl,"// -1-1- the ring @R(t,s,_Dx,_x) is ready");
3512  dbprint(ppl-1, @R);
3513  // create the ideal I
3514  poly  F = imap(save,F);
3515  ideal I = t*F+s;
3516  poly p;
3517  for(i=1; i<=N; i++)
3518  {
3519    p = t; // t
3520    p = diff(F,var(N+2+i))*p;
3521    I = I, var(2+i) + p;
3522  }
3523  // -------- the ideal I is ready ----------
3524  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
3525  dbprint(ppl-1, I);
3526  ideal J = engine(I,eng);
3527  ideal K = nselect(J,1);
3528  dbprint(ppl,"// -1-3- t is eliminated");
3529  dbprint(ppl-1, K);  // K is without t
3530  K = engine(K,eng);  // std does the job too
3531  // now, we must change the ordering
3532  // and create a ring without t
3533  //  setring S;
3534  // ----------- the ring @R3 ------------
3535  // _Dx,_x,s;  +fast ord !
3536  // keep: N, i,j,s, tmp, RL
3537  Nnew = 2*N+1;
3538  kill Lord, tmp, iv, RName;
3539  list Lord, tmp;
3540  intvec iv;
3541  list L=imap(save,L);
3542  list RL=imap(save,RL);
3543  L[1] = RL[1];
3544  L[4] = RL[4];  // char, minpoly
3545  // check whether vars hava admissible names -> done earlier
3546  // now, create the names for new var
3547  tmp[1] = "s";
3548  // DName is defined earlier
3549  list NName = DName + Name + tmp;
3550  L[2] = NName;
3551  tmp = 0;
3552  // just dp
3553  // block ord (dp(N),dp);
3554  string s = "iv=";
3555  for (i=1; i<=Nnew; i++)
3556  {
3557    s = s+"1,";
3558  }
3559  s[size(s)]=";";
3560  execute(s);
3561  tmp[1] = "dp";  // string
3562  tmp[2] = iv;   // intvec
3563  Lord[1] = tmp;
3564  kill s;
3565  kill NName;
3566  tmp[1]      = "C";
3567  Lord[2]     = tmp;  tmp = 0;
3568  L[3]        = Lord;
3569  // we are done with the list. Now add a Plural part
3570  def @R2@ = ring(L);
3571  setring @R2@;
3572  matrix @D[Nnew][Nnew];
3573  for (i=1; i<=N; i++)
3574  {
3575    @D[i,N+i]=-1;
3576  }
3577  def @R2 = nc_algebra(1,@D);
3578  setring @R2;
3579  kill @R2@;
3580  dbprint(ppl,"//  -2-1- the ring @R2(_Dx,_x,s) is ready");
3581  dbprint(ppl-1, @R2);
3582  ideal MM = maxideal(1);
3583  MM = 0,s,MM;
3584  map R01 = @R, MM;
3585  ideal K = R01(K);
3586  // total cleanup
3587  poly F = imap(save, F);
3588  ideal LD = K,F;
3589  dbprint(ppl,"//  -2-2- start GB computations for Ann f^s + f");
3590  dbprint(ppl-1, LD);
3591  LD = engine(LD,eng);
3592  dbprint(ppl,"//  -2-3- finished GB computations for Ann f^s + f");
3593  dbprint(ppl-1, LD);
3595  for (i=1; i<= ncols(LD); i++)
3596  {
3598    {
3599      LD[i] = -LD[i];
3600    }
3601  }
3602  export LD;
3603  kill @R;
3604  return(@R2);
3605}
3606example
3607{
3608  "EXAMPLE:"; echo = 2;
3609  ring r = 0,(x,y,z,w),Dp;
3610  poly F = x^3+y^3+z^3*w;
3611  printlevel = 0;
3612  def A = SannfsBFCT(F); setring A;
3613  intvec v = 1,2,3,4,5,6,7,8;
3614  // are there polynomials, depending on s only?
3615  nselect(LD,v);
3616  // a fancier example:
3617  def R = reiffen(4,5); setring R;
3618  v = 1,2,3,4;
3619  RC; // the Reiffen curve in 4,5
3620  def B = SannfsBFCT(RC);
3621  setring B;
3622  // are there polynomials, depending on s only?
3623  nselect(LD,v);
3624  // it is not the case. Are there leading monomials in s only?
3626}
3627
3628proc SannfsBFCTstd(poly F, list #)
3629"USAGE:  SannfsBFCTstd(f [,eng]);  f a poly, eng an optional int
3630RETURN:  ring
3631PURPOSE: compute ann f^s and Groebner basis of ann f^s+f in D[s]
3632NOTE:    activate this ring with the @code{setring} command.
3633@*       This procedure, unlike SannfsBM, returns a ring with the degrevlex ordering in all variables.
3634@*       In the ring D[s], the ideal LD (which IS a Groebner basis) is the needed ideal.
3635@*       In this procedure @code{std} is used for Groebner basis computations.
3636DISPLAY: If printlevel=1, progress debug messages will be printed,
3637@*       if printlevel>=2, all the debug messages will be printed.
3638EXAMPLE: example SannfsBFCTstd; shows examples
3639"
3640{
3641  // DEBUG INFO: ordering on the output ring = dp,
3642  // use std(K,F); for reusing the std property of K
3643
3644  int eng = 0;
3645  if ( size(#)>0 )
3646  {
3647    if ( typeof(#[1]) == "int" )
3648    {
3649      eng = int(#[1]);
3650    }
3651  }
3652  // returns a list with a ring and an ideal LD in it
3653  int ppl = printlevel-voice+2;
3654  //  printf("plevel :%s, voice: %s",printlevel,voice);
3655  def save = basering;
3656  int N = nvars(basering);
3657  int Nnew = 2*N+2;
3658  int i,j;
3659  string s;
3660  list RL = ringlist(basering);
3661  list L, Lord;
3662  list tmp;
3663  intvec iv;
3664  L[1] = RL[1]; // char
3665  L[4] = RL[4]; // char, minpoly
3666  // check whether vars have admissible names
3667  list Name  = RL[2];
3668  list RName;
3669  RName[1] = "@t";
3670  RName[2] = "@s";
3671  for(i=1;i<=N;i++)
3672  {
3673    for(j=1; j<=size(RName);j++)
3674    {
3675      if (Name[i] == RName[j])
3676      {
3677        ERROR("Variable names should not include @t,@s");
3678      }
3679    }
3680  }
3681  // now, create the names for new vars
3682  list DName;
3683  for(i=1;i<=N;i++)
3684  {
3685    DName[i] = "D"+Name[i]; // concat
3686  }
3687  tmp[1] = "t";
3688  tmp[2] = "s";
3689  list NName = tmp + DName + Name ;
3690  L[2]   = NName;
3691  // Name, Dname will be used further
3692  kill NName;
3693  // block ord (lp(2),dp);
3694  tmp[1]  = "lp"; // string
3695  iv      = 1,1;
3696  tmp[2]  = iv; //intvec
3697  Lord[1] = tmp;
3698  // continue with dp 1,1,1,1...
3699  tmp[1]  = "dp"; // string
3700  s       = "iv=";
3701  for(i=1;i<=Nnew;i++)
3702  {
3703    s = s+"1,";
3704  }
3705  s[size(s)]= ";";
3706  execute(s);
3707  kill s;
3708  tmp[2]    = iv;
3709  Lord[2]   = tmp;
3710  tmp[1]    = "C";
3711  iv        = 0;
3712  tmp[2]    = iv;
3713  Lord[3]   = tmp;
3714  tmp       = 0;
3715  L[3]      = Lord;
3716  // we are done with the list
3717  def @R@ = ring(L);
3718  setring @R@;
3719  matrix @D[Nnew][Nnew];
3720  @D[1,2]=t;
3721  for(i=1; i<=N; i++)
3722  {
3723    @D[2+i,N+2+i]=-1;
3724  }
3725  //  L[5] = matrix(UpOneMatrix(Nnew));
3726  //  L[6] = @D;
3727  def @R = nc_algebra(1,@D);
3728  setring @R;
3729  kill @R@;
3730  dbprint(ppl,"// -1-1- the ring @R(t,s,_Dx,_x) is ready");
3731  dbprint(ppl-1, @R);
3732  // create the ideal I
3733  poly  F = imap(save,F);
3734  ideal I = t*F+s;
3735  poly p;
3736  for(i=1; i<=N; i++)
3737  {
3738    p = t; // t
3739    p = diff(F,var(N+2+i))*p;
3740    I = I, var(2+i) + p;
3741  }
3742  // -------- the ideal I is ready ----------
3743  dbprint(ppl,"// -1-2- starting the elimination of t in @R");
3744  dbprint(ppl-1, I);
3745  ideal J = engine(I,eng);
3746  ideal K = nselect(J,1);
3747  dbprint(ppl,"// -1-3- t is eliminated");
3748  dbprint(ppl-1, K);  // K is without t
3749  K = engine(K,eng);  // std does the job too
3750  // now, we must change the ordering
3751  // and create a ring without t
3752  //  setring S;
3753  // ----------- the ring @R3 ------------
3754  // _Dx,_x,s;  +fast ord !
3755  // keep: N, i,j,s, tmp, RL
3756  Nnew = 2*N+1;
3757  kill Lord, tmp, iv, RName;
3758  list Lord, tmp;
3759  intvec iv;
3760  list L=imap(save,L);
3761  list RL=imap(save,RL);
3762  L[1] = RL[1];
3763  L[4] = RL[4];  // char, minpoly
3764  // check whether vars hava admissible names -> done earlier
3765  // now, create the names for new var
3766  tmp[1] = "s";
3767  // DName is defined earlier
3768  list NName = DName + Name + tmp;
3769  L[2] = NName;
3770  tmp = 0;
3771  // just dp
3772  // block ord (dp(N),dp);
3773  string s = "iv=";
3774  for (i=1; i<=Nnew; i++)
3775  {
3776    s = s+"1,";
3777  }
3778  s[size(s)]=";";
3779  execute(s);
3780  tmp[1] = "dp";  // string
3781  tmp[2] = iv;   // intvec
3782  Lord[1] = tmp;
3783  kill s;
3784  kill NName;
3785  tmp[1]      = "C";
3786  Lord[2]     = tmp;  tmp = 0;
3787  L[3]        = Lord;
3788  // we are done with the list. Now add a Plural part
3789  def @R2@ = ring(L);
3790  setring @R2@;
3791  matrix @D[Nnew][Nnew];
3792  for (i=1; i<=N; i++)
3793  {
3794    @D[i,N+i]=-1;
3795  }
3796  def @R2 = nc_algebra(1,@D);
3797  setring @R2;
3798  kill @R2@;
3799  dbprint(ppl,"//  -2-1- the ring @R2(_Dx,_x,s) is ready");
3800  dbprint(ppl-1, @R2);
3801  ideal MM = maxideal(1);
3802  MM = 0,s,MM;
3803  map R01 = @R, MM;
3804  ideal K = R01(K);
3805  // total cleanup
3806  poly F = imap(save, F);
3807  //  ideal LD = K,F;
3808  dbprint(ppl,"//  -2-2- start GB computations for Ann f^s + f");
3809  //  dbprint(ppl-1, LD);
3810  ideal LD = std(K,F);
3811  //  LD = engine(LD,eng);
3812  dbprint(ppl,"//  -2-3- finished GB computations for Ann f^s + f");
3813  dbprint(ppl-1, LD);
3815  for (i=1; i<= ncols(LD); i++)
3816  {
3818    {
3819      LD[i] = -LD[i];
3820    }
3821  }
3822  export LD;
3823  kill @R;
3824  return(@R2);
3825}
3826example
3827{
3828  "EXAMPLE:"; echo = 2;
3829  ring r = 0,(x,y,z,w),Dp;
3830  poly F = x^3+y^3+z^3*w;
3831  printlevel = 0;
3832  def A = SannfsBFCT(F); setring A;
3833  intvec v = 1,2,3,4,5,6,7,8;
3834  // are there polynomials, depending on s only?
3835  nselect(LD,v);
3836  // a fancier example:
3837  def R = reiffen(4,5); setring R;
3838  v = 1,2,3,4;
3839  RC; // the Reiffen curve in 4,5
3840  def B = SannfsBFCT(RC);
3841  setring B;
3842  // are there polynomials, depending on s only?
3843  nselect(LD,v);
3844  // it is not the case. Are there leading monomials in s only?
3846}
3847
3848// use a finer ordering
3849
3850proc SannfsLOT(poly F, list #)
3851"USAGE:  SannfsLOT(f [,eng]);  f a poly, eng an optional int
3852RETURN:  ring
3853PURPOSE: 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
3854NOTE:    activate this ring with the @code{setring} command.
3855@*       In the ring D[s], the ideal LD (which is NOT a Groebner basis) is the needed D-module structure.
3856@*       If eng <>0, @code{std} is used for Groebner basis computations,
3857@*       otherwise, and by default @code{slimgb} is used.
3858@*       If printlevel=1, progress debug messages will be printed,
3859@*       if printlevel>=2, all the debug messages will be printed.
3860EXAMPLE: example SannfsLOT; shows examples
3861"
3862{
3863  int eng = 0;
3864  if ( size(#)>0 )
3865  {
3866    if ( typeof(#[1]) == "int" )
3867    {
3868      eng = int(#[1]);
3869    }
3870  }
3871  // returns a list with a ring and an ideal LD in it
3872  int ppl = printlevel-voice+2;
3873  //  printf("plevel :%s, voice: %s",printlevel,voice);
3874  def save = basering;
3875  int N = nvars(basering);
3876  //  int Nnew = 2*(N+2);
3877  int Nnew = 2*(N+1)+1; //removed u,v; added s
3878  int i,j;
3879  string s;
3880  list RL = ringlist(basering);
3881  list L, Lord;
3882  list tmp;
3883  intvec iv;
3884  L[1] = RL[1]; // char
3885  L[4] = RL[4]; // char, minpoly
3886  // check whether vars have admissible names
3887  list Name  = RL[2];
3888  list RName;
3889//   RName[1] = "u";
3890//   RName[2] = "v";
3891  RName[1] = "t";
3892  RName[2] = "Dt";
3893  for(i=1;i<=N;i++)
3894  {
3895    for(j=1; j<=size(RName);j++)
3896    {
3897      if (Name[i] == RName[j])
3898      {
3899        ERROR("Variable names should not include t,Dt");
3900      }
3901    }
3902  }
3903  // now, create the names for new vars
3904//   tmp[1]     = "u";
3905//   tmp[2]     = "v";
3906//   list UName = tmp;
3907  list DName;
3908  for(i=1;i<=N;i++)
3909  {
3910    DName[i] = "D"+Name[i]; // concat
3911  }
3912  tmp    =  0;
3913  tmp[1] = "t";
3914  tmp[2] = "Dt";
3915  list SName ; SName[1] = "s";
3916  //  list NName = tmp + Name + DName + SName;
3917  list NName = tmp + SName + Name + DName;
3918  L[2]   = NName;
3919  tmp    = 0;
3920  // Name, Dname will be used further
3921  //  kill UName;
3922  kill NName;
3923  // block ord (a(1,1),dp);
3924  tmp[1]  = "a"; // string
3925  iv      = 1,1;
3926  tmp[2]  = iv; //intvec
3927  Lord[1] = tmp;
3928  // continue with a(0,0,1)
3929  tmp[1]  = "a"; // string
3930  iv      = 0,0,1;
3931  tmp[2]  = iv; //intvec
3932  Lord[2] = tmp;
3933  // continue with dp 1,1,1,1...
3934  tmp[1]  = "dp"; // string
3935  s       = "iv=";
3936  for(i=1;i<=Nnew;i++)
3937  {
3938    s = s+"1,";
3939  }
3940  s[size(s)]= ";";
3941  execute(s);
3942  tmp[2]    = iv;
3943  Lord[3]   = tmp;
3944  tmp[1]    = "C";
3945  iv        = 0;
3946  tmp[2]    = iv;
3947  Lord[4]   = tmp;
3948  tmp       = 0;
3949  L[3]      = Lord;
3950  // we are done with the list
3951  def @R@ = ring(L);
3952  setring @R@;
3953  matrix @D[Nnew][Nnew];
3954  @D[1,2]=1;
3955  for(i=1; i<=N; i++)
3956  {
3957    @D[3+i,N+3+i]=1;
3958  }
3960  @D[1,3] = -var(1);
3961  @D[2,3] = var(2);
3962  //  @D[N+3,2*(N+2)]=1; old t,Dt stuff
3963  //  L[5] = matrix(UpOneMatrix(Nnew));
3964  //  L[6] = @D;
3965  def @R = nc_algebra(1,@D);
3966  setring @R;
3967  kill @R@;
3968  dbprint(ppl,"// -1-1- the ring @R(t,Dt,s,_x,_Dx) is ready");
3969  dbprint(ppl-1, @R);
3970  // create the ideal I
3971  poly  F = imap(save,F);
3972  //  ideal I = u*F-t,u*v-1;
3973  ideal I = F-t;
3974  poly p;
3975  for(i=1; i<=N; i++)
3976  {
3977    //    p = u*Dt; // u*Dt
3978    p = Dt;
3979    p = diff(F,var(3+i))*p;
3980    I = I, var(N+3+i) + p;
3981  }
3982  //  I = I, var(1)*var(2) + var(Nnew) +1; // reduce it with t-f!!!
3983  // t*Dt + s +1 reduced with t-f gives f*Dt + s
3984  I = I, F*var(2) + var(3);
3985  // -------- the ideal I is ready ----------
3986  dbprint(ppl,"// -1-2- starting the elimination of t,Dt in @R");
3987  dbprint(ppl-1, I);
3988  ideal J = engine(I,eng);
3989  ideal K = nselect(J,1..2);
3990  dbprint(ppl,"// -1-3- t,Dt are eliminated");
3991  dbprint(ppl-1, K);  // K is without t, Dt
3992  K = engine(K,eng);  // std does the job too
3993  // now, we must change the ordering
3994  // and create a ring without t, Dt
3995  setring save;
3996  // ----------- the ring @R3 ------------
3997  // _x, _Dx,s;  elim.ord for _x,_Dx.
3998  // keep: N, i,j,s, tmp, RL
3999  Nnew = 2*N+1;
4000  kill Lord, tmp, iv, RName;
4001  list Lord, tmp;
4002  intvec iv;
4003  L[1] = RL[1];
4004  L[4] = RL[4];  // char, minpoly
4005  // check whether vars hava admissible names -> done earlier
4006  // now, create the names for new var
4007  tmp[1] = "s";
4008  // DName is defined earlier
4009  list NName = Name + DName + tmp;
4010  L[2] = NName;
4011  tmp = 0;
4012  // block ord (dp(N),dp);
4013  // string s is already defined
4014  s = "iv=";
4015  for (i=1; i<=Nnew-1; i++)
4016  {
4017    s = s+"1,";
4018  }
4019  s[size(s)]=";";
4020  execute(s);
4021  tmp[1] = "dp";  // string
4022  tmp[2] = iv;   // intvec
4023  Lord[1] = tmp;
4024  // continue with dp 1,1,1,1...
4025  tmp[1] = "dp";  // string
4026  s[size(s)] = ",";
4027  s = s+"1;";
4028  execute(s);
4029  kill s;
4030  kill NName;
4031  tmp[2]      = iv;
4032  Lord[2]     = tmp;
4033  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
4034  Lord[3]     = tmp;  tmp = 0;
4035  L[3]        = Lord;
4036  // we are done with the list. Now add a Plural part
4037  def @R2@ = ring(L);
4038  setring @R2@;
4039  matrix @D[Nnew][Nnew];
4040  for (i=1; i<=N; i++)
4041  {
4042    @D[i,N+i]=1;
4043  }
4044  def @R2 = nc_algebra(1,@D);
4045  setring @R2;
4046  kill @R2@;
4047  dbprint(ppl,"//  -2-1- the ring @R2(_x,_Dx,s) is ready");
4048  dbprint(ppl-1, @R2);
4049  ideal MM = maxideal(1);
4050  //  MM = 0,s,MM;
4051  MM = 0,0,s,MM[1..size(MM)-1];
4052  map R01 = @R, MM;
4053  ideal K = R01(K);
4054  // total cleanup
4055  ideal LD = K;
4057  for (i=1; i<= ncols(LD); i++)
4058  {
4060    {
4061      LD[i] = -LD[i];
4062    }
4063  }
4064  export LD;
4065  kill @R;
4066  return(@R2);
4067}
4068example
4069{
4070  "EXAMPLE:"; echo = 2;
4071  ring r = 0,(x,y,z),Dp;
4072  poly F = x^3+y^3+z^3;
4073  printlevel = 0;
4074  def A  = SannfsLOT(F);
4075  setring A;
4076  LD;
4077}
4078
4079/*
4080proc SannfsLOTold(poly F, list #)
4081"USAGE:  SannfsLOT(f [,eng]);  f a poly, eng an optional int
4082RETURN:  ring
4083PURPOSE: 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
4084NOTE:    activate this ring with the @code{setring} command.
4085@*       In the ring D[s], the ideal LD (which is NOT a Groebner basis) is the needed D-module structure.
4086@*       If eng <>0, @code{std} is used for Groebner basis computations,
4087@*       otherwise, and by default @code{slimgb} is used.
4088@*       If printlevel=1, progress debug messages will be printed,
4089@*       if printlevel>=2, all the debug messages will be printed.
4090EXAMPLE: example SannfsLOT; shows examples
4091"
4092{
4093  int eng = 0;
4094  if ( size(#)>0 )
4095  {
4096    if ( typeof(#[1]) == "int" )
4097    {
4098      eng = int(#[1]);
4099    }
4100  }
4101  // returns a list with a ring and an ideal LD in it
4102  int ppl = printlevel-voice+2;
4103  //  printf("plevel :%s, voice: %s",printlevel,voice);
4104  def save = basering;
4105  int N = nvars(basering);
4106  //  int Nnew = 2*(N+2);
4107  int Nnew = 2*(N+1)+1; //removed u,v; added s
4108  int i,j;
4109  string s;
4110  list RL = ringlist(basering);
4111  list L, Lord;
4112  list tmp;
4113  intvec iv;
4114  L[1] = RL[1]; // char
4115  L[4] = RL[4]; // char, minpoly
4116  // check whether vars have admissible names
4117  list Name  = RL[2];
4118  list RName;
4119//   RName[1] = "u";
4120//   RName[2] = "v";
4121  RName[1] = "t";
4122  RName[2] = "Dt";
4123  for(i=1;i<=N;i++)
4124  {
4125    for(j=1; j<=size(RName);j++)
4126    {
4127      if (Name[i] == RName[j])
4128      {
4129        ERROR("Variable names should not include t,Dt");
4130      }
4131    }
4132  }
4133  // now, create the names for new vars
4134//   tmp[1]     = "u";
4135//   tmp[2]     = "v";
4136//   list UName = tmp;
4137  list DName;
4138  for(i=1;i<=N;i++)
4139  {
4140    DName[i] = "D"+Name[i]; // concat
4141  }
4142  tmp    =  0;
4143  tmp[1] = "t";
4144  tmp[2] = "Dt";
4145  list SName ; SName[1] = "s";
4146  //  list NName = UName +  tmp + Name + DName;
4147  list NName = tmp + Name + DName + SName;
4148  L[2]   = NName;
4149  tmp    = 0;
4150  // Name, Dname will be used further
4151  //  kill UName;
4152  kill NName;
4153  // block ord (a(1,1),dp);
4154  tmp[1]  = "a"; // string
4155  iv      = 1,1;
4156  tmp[2]  = iv; //intvec
4157  Lord[1] = tmp;
4158  // continue with dp 1,1,1,1...
4159  tmp[1]  = "dp"; // string
4160  s       = "iv=";
4161  for(i=1;i<=Nnew;i++)
4162  {
4163    s = s+"1,";
4164  }
4165  s[size(s)]= ";";
4166  execute(s);
4167  tmp[2]    = iv;
4168  Lord[2]   = tmp;
4169  tmp[1]    = "C";
4170  iv        = 0;
4171  tmp[2]    = iv;
4172  Lord[3]   = tmp;
4173  tmp       = 0;
4174  L[3]      = Lord;
4175  // we are done with the list
4176  def @R@ = ring(L);
4177  setring @R@;
4178  matrix @D[Nnew][Nnew];
4179  @D[1,2]=1;
4180  for(i=1; i<=N; i++)
4181  {
4182    @D[2+i,N+2+i]=1;
4183  }
4185  @D[1,Nnew] = -var(1);
4186  @D[2,Nnew] = var(2);
4187  //  @D[N+3,2*(N+2)]=1; old t,Dt stuff
4188  //  L[5] = matrix(UpOneMatrix(Nnew));
4189  //  L[6] = @D;
4190  def @R = nc_algebra(1,@D);
4191  setring @R;
4192  kill @R@;
4193  dbprint(ppl,"// -1-1- the ring @R(t,Dt,_x,_Dx,s) is ready");
4194  dbprint(ppl-1, @R);
4195  // create the ideal I
4196  poly  F = imap(save,F);
4197  //  ideal I = u*F-t,u*v-1;
4198  ideal I = F-t;
4199  poly p;
4200  for(i=1; i<=N; i++)
4201  {
4202    //    p = u*Dt; // u*Dt
4203    p = Dt;
4204    p = diff(F,var(2+i))*p;
4205    I = I, var(N+2+i) + p;
4206  }
4207  //  I = I, var(1)*var(2) + var(Nnew) +1; // reduce it with t-f!!!
4208  // t*Dt + s +1 reduced with t-f gives f*Dt + s
4209  I = I, F*var(2) + var(Nnew);
4210  // -------- the ideal I is ready ----------
4211  dbprint(ppl,"// -1-2- starting the elimination of t,Dt in @R");
4212  dbprint(ppl-1, I);
4213  ideal J = engine(I,eng);
4214  ideal K = nselect(J,1..2);
4215  dbprint(ppl,"// -1-3- t,Dt are eliminated");
4216  dbprint(ppl-1, K);  // K is without t, Dt
4217  K = engine(K,eng);  // std does the job too
4218  // now, we must change the ordering
4219  // and create a ring without t, Dt
4220  setring save;
4221  // ----------- the ring @R3 ------------
4222  // _x, _Dx,s;  elim.ord for _x,_Dx.
4223  // keep: N, i,j,s, tmp, RL
4224  Nnew = 2*N+1;
4225  kill Lord, tmp, iv, RName;
4226  list Lord, tmp;
4227  intvec iv;
4228  L[1] = RL[1];
4229  L[4] = RL[4];  // char, minpoly
4230  // check whether vars hava admissible names -> done earlier
4231  // now, create the names for new var
4232  tmp[1] = "s";
4233  // DName is defined earlier
4234  list NName = Name + DName + tmp;
4235  L[2] = NName;
4236  tmp = 0;
4237  // block ord (dp(N),dp);
4238  // string s is already defined
4239  s = "iv=";
4240  for (i=1; i<=Nnew-1; i++)
4241  {
4242    s = s+"1,";
4243  }
4244  s[size(s)]=";";
4245  execute(s);
4246  tmp[1] = "dp";  // string
4247  tmp[2] = iv;   // intvec
4248  Lord[1] = tmp;
4249  // continue with dp 1,1,1,1...
4250  tmp[1] = "dp";  // string
4251  s[size(s)] = ",";
4252  s = s+"1;";
4253  execute(s);
4254  kill s;
4255  kill NName;
4256  tmp[2]      = iv;
4257  Lord[2]     = tmp;
4258  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
4259  Lord[3]     = tmp;  tmp = 0;
4260  L[3]        = Lord;
4261  // we are done with the list. Now add a Plural part
4262  def @R2@ = ring(L);
4263  setring @R2@;
4264  matrix @D[Nnew][Nnew];
4265  for (i=1; i<=N; i++)
4266  {
4267    @D[i,N+i]=1;
4268  }
4269  def @R2 = nc_algebra(1,@D);
4270  setring @R2;
4271  kill @R2@;
4272  dbprint(ppl,"//  -2-1- the ring @R2(_x,_Dx,s) is ready");
4273  dbprint(ppl-1, @R2);
4274  ideal MM = maxideal(1);
4275  MM = 0,s,MM;
4276  map R01 = @R, MM;
4277  ideal K = R01(K);
4278  // total cleanup
4279  ideal LD = K;
4281  for (i=1; i<= ncols(LD); i++)
4282  {
4284    {
4285      LD[i] = -LD[i];
4286    }
4287  }
4288  export LD;
4289  kill @R;
4290  return(@R2);
4291}
4292example
4293{
4294  "EXAMPLE:"; echo = 2;
4295  ring r = 0,(x,y,z),Dp;
4296  poly F = x^3+y^3+z^3;
4297  printlevel = 0;
4298  def A  = SannfsLOTold(F);
4299  setring A;
4300  LD;
4301}
4302
4303*/
4304
4305proc annfsLOT(poly F, list #)
4306"USAGE:  annfsLOT(F [,eng]);  F a poly, eng an optional int
4307RETURN:  ring
4308PURPOSE: compute the D-module structure of basering[1/f]*f^s, according to the Levandovskyy's modification of the algorithm by Oaku and Takayama
4309NOTE:    activate this ring with the @code{setring} command. In this ring,
4310@*       - the ideal LD (which is a Groebner basis) is the needed D-module structure,
4311@*         which is obtained by substituting the minimal integer root of a Bernstein
4312@*         polynomial into the s-parametric ideal;
4313@*       - the list BS contains the roots with multiplicities of a Bernstein polynomial of f.
4314@*       If eng <>0, @code{std} is used for Groebner basis computations,
4315@*       otherwise and by default @code{slimgb} is used.
4316@*       If printlevel=1, progress debug messages will be printed,
4317@*       if printlevel>=2, all the debug messages will be printed.
4318EXAMPLE: example annfsLOT; shows examples
4319"
4320{
4321  int eng = 0;
4322  if ( size(#)>0 )
4323  {
4324    if ( typeof(#[1]) == "int" )
4325    {
4326      eng = int(#[1]);
4327    }
4328  }
4329  printlevel=printlevel+1;
4330  def save = basering;
4331  def @A = SannfsLOT(F,eng);
4332  setring @A;
4333  poly F = imap(save,F);
4334  def B  = annfs0(LD,F,eng);
4335  return(B);
4336}
4337example
4338{
4339  "EXAMPLE:"; echo = 2;
4340  ring r = 0,(x,y,z),Dp;
4341  poly F = z*x^2+y^3;
4342  printlevel = 0;
4343  def A  = annfsLOT(F);
4344  setring A;
4345  LD;
4346  BS;
4347}
4348
4349proc annfs0(ideal I, poly F, list #)
4350"USAGE:  annfs0(I, F [,eng]);  I an ideal, F a poly, eng an optional int
4351RETURN:  ring
4352PURPOSE: compute the annihilator ideal of f^s in the Weyl Algebra, based on the
4353output of Sannfs-like procedure
4354NOTE:    activate this ring with the @code{setring} command. In this ring,
4355@*       - the ideal LD (which is a Groebner basis) is the annihilator of f^s,
4356@*       - the list BS contains the roots with multiplicities of a Bernstein polynomial of f.
4357@*       If eng <>0, @code{std} is used for Groebner basis computations,
4358@*       otherwise and by default @code{slimgb} is used.
4359@*       If printlevel=1, progress debug messages will be printed,
4360@*       if printlevel>=2, all the debug messages will be printed.
4361EXAMPLE: example annfs0; shows examples
4362"
4363{
4364  int eng = 0;
4365  if ( size(#)>0 )
4366  {
4367    if ( typeof(#[1]) == "int" )
4368    {
4369      eng = int(#[1]);
4370    }
4371  }
4372  def @R2 = basering;
4373  // we're in D_n[s], where the elim ord for s is set
4374  ideal J = NF(I,std(F));
4376  int i;
4377  for (i=1; i<= ncols(J); i++)
4378  {
4380    {
4381      J[i] = -J[i];
4382    }
4383  }
4384  J = J,F;
4385  ideal M = engine(J,eng);
4386  int Nnew = nvars(@R2);
4387  ideal K2 = nselect(M,1..Nnew-1);
4388  int ppl = printlevel-voice+2;
4389  dbprint(ppl,"// -1-1- _x,_Dx are eliminated in basering");
4390  dbprint(ppl-1, K2);
4391  // the ring @R3 and the search for minimal negative int s
4392  ring @R3 = 0,s,dp;
4393  dbprint(ppl,"// -2-1- the ring @R3 i.e. K[s] is ready");
4394  ideal K3 = imap(@R2,K2);
4395  poly p = K3[1];
4396  dbprint(ppl,"// -2-2- factorization");
4397  //  ideal P = factorize(p,1);  //without constants and multiplicities
4398  //  "--------- b-function factorizes into ---------";   P;
4399  // convert factors to the list of their roots with mults
4400  // assume all factors are linear
4401  //  ideal BS = normalize(P);
4402  //  BS = subst(BS,s,0);
4403  //  BS = -BS;
4404  list P = factorize(p);          //with constants and multiplicities
4405  ideal bs; intvec m;             //the Bernstein polynomial is monic, so we are not interested in constants
4406  for (i=2; i<= size(P[1]); i++)  //we delete P[1][1] and P[2][1]
4407  {
4408    bs[i-1] = P[1][i];
4409    m[i-1]  = P[2][i];
4410  }
4411  int sP = minIntRoot(bs,1);
4412  bs =  normalize(bs);
4413  bs = -subst(bs,s,0);
4414  dbprint(ppl,"// -2-3- minimal integer root found");
4415  dbprint(ppl-1, sP);
4416 //TODO: sort BS!
4417  // --------- substitute s found in the ideal ---------
4418  // --------- going back to @R and substitute ---------
4419  setring @R2;
4420  K2 = subst(I,s,sP);
4421  // create the ordinary Weyl algebra and put the result into it,
4422  // thus creating the ring @R5
4423  // keep: N, i,j,s, tmp, RL
4424  Nnew = Nnew - 1; // former 2*N;
4425  // list RL = ringlist(save);  // is defined earlier
4426  //  kill Lord, tmp, iv;
4427  list L = 0;
4428  list Lord, tmp;
4429  intvec iv;
4430  list RL = ringlist(basering);
4431  L[1] = RL[1];
4432  L[4] = RL[4];  //char, minpoly
4433  // check whether vars have admissible names -> done earlier
4434  // list Name = RL[2]M
4435  // DName is defined earlier
4436  list NName; // = RL[2]; // skip the last var 's'
4437  for (i=1; i<=Nnew; i++)
4438  {
4439    NName[i] =  RL[2][i];
4440  }
4441  L[2] = NName;
4442  // dp ordering;
4443  string s = "iv=";
4444  for (i=1; i<=Nnew; i++)
4445  {
4446    s = s+"1,";
4447  }
4448  s[size(s)] = ";";
4449  execute(s);
4450  tmp     = 0;
4451  tmp[1]  = "dp";  // string
4452  tmp[2]  = iv;  // intvec
4453  Lord[1] = tmp;
4454  kill s;
4455  tmp[1]  = "C";
4456  iv = 0;
4457  tmp[2]  = iv;
4458  Lord[2] = tmp;
4459  tmp     = 0;
4460  L[3]    = Lord;
4461  // we are done with the list
4463  def @R4@ = ring(L);
4464  setring @R4@;
4465  int N = Nnew/2;
4466  matrix @D[Nnew][Nnew];
4467  for (i=1; i<=N; i++)
4468  {
4469    @D[i,N+i]=1;
4470  }
4471  def @R4 = nc_algebra(1,@D);
4472  setring @R4;
4473  kill @R4@;
4474  dbprint(ppl,"// -3-1- the ring @R4 is ready");
4475  dbprint(ppl-1, @R4);
4476  ideal K4 = imap(@R2,K2);
4477  option(redSB);
4478  dbprint(ppl,"// -3-2- the final cosmetic std");
4479  K4 = engine(K4,eng);  // std does the job too
4480  // total cleanup
4481  ideal bs = imap(@R3,bs);
4482  kill @R3;
4483  list BS = bs,m;
4484  export BS;
4485  ideal LD = K4;
4486  export LD;
4487  return(@R4);
4488}
4489example
4490{ "EXAMPLE:"; echo = 2;
4491  ring r = 0,(x,y,z),Dp;
4492  poly F = x^3+y^3+z^3;
4493  printlevel = 0;
4494  def A = SannfsBM(F);   setring A;
4495  // alternatively, one can use SannfsOT or SannfsLOT
4496  LD;
4497  poly F = imap(r,F);
4498  def B  = annfs0(LD,F);  setring B;
4499  LD;
4500  BS;
4501}
4502
4503// proc annfsgms(poly F, list #)
4504// "USAGE:  annfsgms(f [,eng]);  f a poly, eng an optional int
4505// ASSUME:  f has an isolated critical point at 0
4506// RETURN:  ring
4507// PURPOSE: compute the D-module structure of basering[1/f]*f^s
4508// NOTE:    activate this ring with the @code{setring} command. In this ring,
4509// @*       - the ideal LD is the needed D-mod structure,
4510// @*       - the ideal BS is the list of roots of a Bernstein polynomial of f.
4511// @*       If eng <>0, @code{std} is used for Groebner basis computations,
4512// @*       otherwise (and by default) @code{slimgb} is used.
4513// @*       If printlevel=1, progress debug messages will be printed,
4514// @*       if printlevel>=2, all the debug messages will be printed.
4515// EXAMPLE: example annfsgms; shows examples
4516// "
4517// {
4518//   LIB "gmssing.lib";
4519//   int eng = 0;
4520//   if ( size(#)>0 )
4521//   {
4522//     if ( typeof(#[1]) == "int" )
4523//     {
4524//       eng = int(#[1]);
4525//     }
4526//   }
4527//   int ppl = printlevel-voice+2;
4528//   // returns a ring with the ideal LD in it
4529//   def save = basering;
4530//   // compute the Bernstein poly from gmssing.lib
4531//   list RL = ringlist(basering);
4532//   // in the descr. of the ordering, replace "p" by "s"
4533//   list NL = convloc(RL);
4534//   // create a ring with the ordering, converted to local
4535//   def @LR = ring(NL);
4536//   setring @LR;
4537//   poly F  = imap(save, F);
4538//   ideal B = bernstein(F)[1];
4539//   // since B may not contain (s+1) [following gmssing.lib]
4541//   B = B,-1;
4542//   B = simplify(B,2+4); // erase zero and repeated entries
4543//   // find the minimal integer value
4544//   int   S = minIntRoot(B,0);
4545//   dbprint(ppl,"// -0- minimal integer root found");
4546//   dbprint(ppl-1,S);
4547//   setring save;
4548//   int N = nvars(basering);
4549//   int Nnew = 2*(N+2);
4550//   int i,j;
4551//   string s;
4552//   //  list RL = ringlist(basering);
4553//   list L, Lord;
4554//   list tmp;
4555//   intvec iv;
4556//   L[1] = RL[1]; // char
4557//   L[4] = RL[4]; // char, minpoly
4558//   // check whether vars have admissible names
4559//   list Name  = RL[2];
4560//   list RName;
4561//   RName[1] = "u";
4562//   RName[2] = "v";
4563//   RName[3] = "t";
4564//   RName[4] = "Dt";
4565//   for(i=1;i<=N;i++)
4566//   {
4567//     for(j=1; j<=size(RName);j++)
4568//     {
4569//       if (Name[i] == RName[j])
4570//       {
4571//         ERROR("Variable names should not include u,v,t,Dt");
4572//       }
4573//     }
4574//   }
4575//   // now, create the names for new vars
4576//   //  tmp[1]     = "u";  tmp[2]     = "v";  tmp[3]     = "t";  tmp[4]     = "Dt";
4577//   list UName = RName;
4578//   list DName;
4579//   for(i=1;i<=N;i++)
4580//   {
4581//     DName[i] = "D"+Name[i]; // concat
4582//   }
4583//   list NName = UName + Name + DName;
4584//   L[2]       = NName;
4585//   tmp        = 0;
4586//   // Name, Dname will be used further
4587//   kill UName;
4588//   kill NName;
4589//   // block ord (a(1,1),dp);
4590//   tmp[1]  = "a"; // string
4591//   iv      = 1,1;
4592//   tmp[2]  = iv; //intvec
4593//   Lord[1] = tmp;
4594//   // continue with dp 1,1,1,1...
4595//   tmp[1]  = "dp"; // string
4596//   s       = "iv=";
4597//   for(i=1; i<=Nnew; i++) // need really all vars!
4598//   {
4599//     s = s+"1,";
4600//   }
4601//   s[size(s)]= ";";
4602//   execute(s);
4603//   tmp[2]    = iv;
4604//   Lord[2]   = tmp;
4605//   tmp[1]    = "C";
4606//   iv        = 0;
4607//   tmp[2]    = iv;
4608//   Lord[3]   = tmp;
4609//   tmp       = 0;
4610//   L[3]      = Lord;
4611//   // we are done with the list
4612//   def @R = ring(L);
4613//   setring @R;
4614//   matrix @D[Nnew][Nnew];
4615//   @D[3,4] = 1; // t,Dt
4616//   for(i=1; i<=N; i++)
4617//   {
4618//     @D[4+i,4+N+i]=1;
4619//   }
4620//   //  L[5] = matrix(UpOneMatrix(Nnew));
4621//   //  L[6] = @D;
4622//   nc_algebra(1,@D);
4623//   dbprint(ppl,"// -1-1- the ring @R is ready");
4624//   dbprint(ppl-1,@R);
4625//   // create the ideal
4626//   poly F  = imap(save,F);
4627//   ideal I = u*F-t,u*v-1;
4628//   poly p;
4629//   for(i=1; i<=N; i++)
4630//   {
4631//     p = u*Dt; // u*Dt
4632//     p = diff(F,var(4+i))*p;
4633//     I = I, var(N+4+i) + p; // Dx, Dy
4634//   }
4635//   // add the relations between t,Dt and s
4636//   //  I = I, t*Dt+1+S;
4637//   // -------- the ideal I is ready ----------
4638//   dbprint(ppl,"// -1-2- starting the elimination of u,v in @R");
4639//   ideal J = engine(I,eng);
4640//   ideal K = nselect(J,1..2);
4641//   dbprint(ppl,"// -1-3- u,v are eliminated in @R");
4642//   dbprint(ppl-1,K); // without u,v: not yet our answer
4643//   //----- create a ring with elim.ord for t,Dt -------
4644//   setring save;
4645//   // ------------ new ring @R2 ------------------
4646//   // without u,v and with the elim.ord for t,Dt
4647//   // keep: N, i,j,s, tmp, RL
4648//   Nnew = 2*N+2;
4649//   //  list RL = ringlist(save); // is defined earlier
4650//   kill Lord,tmp,iv, RName;
4651//   L = 0;
4652//   list Lord, tmp;
4653//   intvec iv;
4654//   L[1] = RL[1]; // char
4655//   L[4] = RL[4]; // char, minpoly
4656//   // check whether vars have admissible names -> done earlier
4657//   //  list Name  = RL[2];
4658//   list RName;
4659//   RName[1] = "t";
4660//   RName[2] = "Dt";
4661//   // DName is defined earlier
4662//   list NName = RName + Name + DName;
4663//   L[2]   = NName;
4664//   tmp    = 0;
4665//   // block ord (a(1,1),dp);
4666//   tmp[1]  = "a"; // string
4667//   iv      = 1,1;
4668//   tmp[2]  = iv; //intvec
4669//   Lord[1] = tmp;
4670//   // continue with dp 1,1,1,1...
4671//   tmp[1]  = "dp"; // string
4672//   s       = "iv=";
4673//   for(i=1;i<=Nnew;i++)
4674//   {
4675//     s = s+"1,";
4676//   }
4677//   s[size(s)]= ";";
4678//   execute(s);
4679//   kill s;
4680//   kill NName;
4681//   tmp[2]    = iv;
4682//   Lord[2]   = tmp;
4683//   tmp[1]    = "C";
4684//   iv        = 0;
4685//   tmp[2]    = iv;
4686//   Lord[3]   = tmp;
4687//   tmp       = 0;
4688//   L[3]      = Lord;
4689//   // we are done with the list
4691//   def @R2 = ring(L);
4692//   setring @R2;
4693//   matrix @D[Nnew][Nnew];
4694//   @D[1,2]=1;
4695//   for(i=1; i<=N; i++)
4696//   {
4697//     @D[2+i,2+N+i]=1;
4698//   }
4699//   nc_algebra(1,@D);
4700//   dbprint(ppl,"// -2-1- the ring @R2 is ready");
4701//   dbprint(ppl-1,@R2);
4702//   ideal MM = maxideal(1);
4703//   MM = 0,0,MM;
4704//   map R01 = @R, MM;
4705//   ideal K2 = R01(K);
4706//   // add the relations between t,Dt and s
4707//   //  K2       = K2, t*Dt+1+S;
4708//   poly G = t*Dt+S+1;
4709//   K2 = NF(K2,std(G)),G;
4710//   dbprint(ppl,"// -2-2- starting elimination for t,Dt in @R2");
4711//   ideal J  = engine(K2,eng);
4712//   ideal K  = nselect(J,1..2);
4713//   dbprint(ppl,"// -2-3- t,Dt are eliminated");
4714//   dbprint(ppl-1,K);
4715//   //  "------- produce a final result --------";
4716//   // ----- create the ordinary Weyl algebra and put
4717//   // ----- the result into it
4718//   // ------ create the ring @R5
4719//   // keep: N, i,j,s, tmp, RL
4720//   setring save;
4721//   Nnew = 2*N;
4722//   //  list RL = ringlist(save); // is defined earlier
4723//   kill Lord, tmp, iv;
4724//   //  kill L;
4725//   L = 0;
4726//   list Lord, tmp;
4727//   intvec iv;
4728//   L[1] = RL[1]; // char
4729//   L[4] = RL[4]; // char, minpoly
4730//   // check whether vars have admissible names -> done earlier
4731//   //  list Name  = RL[2];
4732//   // DName is defined earlier
4733//   list NName = Name + DName;
4734//   L[2]   = NName;
4735//   // dp ordering;
4736//   string   s       = "iv=";
4737//   for(i=1;i<=2*N;i++)
4738//   {
4739//     s = s+"1,";
4740//   }
4741//   s[size(s)]= ";";
4742//   execute(s);
4743//   tmp     = 0;
4744//   tmp[1]  = "dp"; // string
4745//   tmp[2]  = iv; //intvec
4746//   Lord[1] = tmp;
4747//   kill s;
4748//   tmp[1]    = "C";
4749//   iv        = 0;
4750//   tmp[2]    = iv;
4751//   Lord[2]   = tmp;
4752//   tmp       = 0;
4753//   L[3]      = Lord;
4754//   // we are done with the list
4756//   def @R5 = ring(L);
4757//   setring @R5;
4758//   matrix @D[Nnew][Nnew];
4759//   for(i=1; i<=N; i++)
4760//   {
4761//     @D[i,N+i]=1;
4762//   }
4763//   nc_algebra(1,@D);
4764//   dbprint(ppl,"// -3-1- the ring @R5 is ready");
4765//   dbprint(ppl-1,@R5);
4766//   ideal K5 = imap(@R2,K);
4767//   option(redSB);
4768//   dbprint(ppl,"// -3-2- the final cosmetic std");
4769//   K5 = engine(K5,eng); // std does the job too
4770//   // total cleanup
4771//   kill @R;
4772//   kill @R2;
4773//   ideal LD = K5;
4774//   ideal BS = imap(@LR,B);
4775//   kill @LR;
4776//   export BS;
4777//   export LD;
4778//   return(@R5);
4779// }
4780// example
4781// {
4782//   "EXAMPLE:"; echo = 2;
4783//   ring r = 0,(x,y,z),Dp;
4784//   poly F = x^2+y^3+z^5;
4785//   def A = annfsgms(F);
4786//   setring A;
4787//   LD;
4788//   print(matrix(BS));
4789// }
4790
4791
4792proc convloc(list @NL)
4793"USAGE:  convloc(L); L a list
4794RETURN:  list
4795PURPOSE: convert a ringlist L into another ringlist,
4796where all the 'p' orderings are replaced with the 's' orderings.
4797ASSUME: L is a result of a ringlist command
4798EXAMPLE: example convloc; shows examples
4799"
4800{
4801  list NL = @NL;
4802  // given a ringlist, returns a new ringlist,
4803  // where all the p-orderings are replaced with s-ord's
4804  int i,j,k,found;
4805  int nblocks = size(NL[3]);
4806  for(i=1; i<=nblocks; i++)
4807  {
4808    for(j=1; j<=size(NL[3][i]); j++)
4809    {
4810      if (typeof(NL[3][i][j]) == "string" )
4811      {
4812        found = 0;
4813        for (k=1; k<=size(NL[3][i][j]); k++)
4814        {
4815          if (NL[3][i][j][k]=="p")
4816          {
4817            NL[3][i][j][k]="s";
4818            found = 1;
4819            //    printf("replaced at %s,%s,%s",i,j,k);
4820          }
4821        }
4822      }
4823    }
4824  }
4825  return(NL);
4826}
4827example
4828{
4829  "EXAMPLE:"; echo = 2;
4830  ring r = 0,(x,y,z),(Dp(2),dp(1));
4831  list L = ringlist(r);
4832  list N = convloc(L);
4833  def rs = ring(N);
4834  setring rs;
4835  rs;
4836}
4837
4838proc annfspecial(ideal I, poly F, int mir, number n)
4839"USAGE:  annfspecial(I,F,mir,n);  I an ideal, F a poly, int mir, number n
4840RETURN:  ideal
4841PURPOSE: compute the annihilator ideal of F^n in the Weyl Algebra for a rational number n
4842ASSUME:  the basering contains 's' as a variable
4843NOTE:    We assume that the basering is D[s],
4844@*          ideal I is the Ann F^s in D[s] (obtained with e.g. SannfsBM, SannfsLOT, SannfsOT)
4845@*          integer 'mir' is the minimal integer root of the Bernstein polynomial of F
4846@*          and the number n is rational.
4847@*       We compute the real annihilator for any rational value of n (both generic and exceptional).
4848@*       The implementation goes along the lines of Saito-Sturmfels-Takayama, Alg. 5.3.15
4849DISPLAY: If printlevel=1, progress debug messages will be printed,
4850@*       if printlevel>=2, all the debug messages will be printed.
4851EXAMPLE: example annfspecial; shows examples
4852"
4853{
4854  int ppl = printlevel-voice+2;
4855  //  int mir =  minIntRoot(L[1],0);
4856  int ns   = varnum("s");
4857  if (!ns)
4858  {
4859    ERROR("Variable s expected in the ideal AnnFs");
4860  }
4861  int d;
4862  ideal P = subst(I,var(ns),n);
4863  if (denominator(n) == 1)
4864  {
4865    // n is fraction-free
4866    d = int(numerator(n));
4867    if ( (!d) && (n!=0))
4868    {
4869      ERROR("no parametric special values are allowed");
4870    }
4871    d = d - mir;
4872    if (d>0)
4873    {
4874      dbprint(ppl,"// -1-1- starting syzygy computations");
4875      matrix J[1][1] = F^d;
4876      dbprint(ppl-1,"// -1-1-1- of the polynomial ideal");
4877      dbprint(ppl-1,J);
4878      matrix K[1][size(I)] = subst(I,var(ns),mir);
4879      dbprint(ppl-1,"// -1-1-2- modulo ideal:");
4880      dbprint(ppl-1, K);
4881      module M = modulo(J,K);
4882      dbprint(ppl-1,"// -1-1-3- getting the result:");
4883      dbprint(ppl-1, M);
4884      P  = P, ideal(M);
4885      dbprint(ppl,"// -1-2- finished syzygy computations");
4886    }
4887    else
4888    {
4889      dbprint(ppl,"// -1-1- d<=0, no syzygy computations needed");
4890      dbprint(ppl-1,"// -1-2- for d =");
4891      dbprint(ppl-1, d);
4892    }
4893  }
4894  // also the else case: d<=0 or n is not an integer
4895  dbprint(ppl,"// -2-1- starting final Groebner basis");
4896  P = groebner(P);
4897  dbprint(ppl,"// -2-2- finished final Groebner basis");
4898  return(P);
4899}
4900example
4901{
4902  "EXAMPLE:"; echo = 2;
4903  ring r = 0,(x,y),dp;
4904  poly F = x3-y2;
4905  def  B = annfs(F);  setring B;
4906  minIntRoot(BS[1],0);
4907  // So, the minimal integer root is -1
4908  setring r;
4909  def  A = SannfsBM(F);
4910  setring A;
4911  poly F = x3-y2;
4912  annfspecial(LD,F,-1,3/4); // generic root
4913  annfspecial(LD,F,-1,-2);  // integer but still generic root
4914  annfspecial(LD,F,-1,1);   // exceptional root
4915}
4916
4917/*
4918//static proc new_ex_annfspecial()
4919{
4920  //another example for annfspecial: x3+y3+z3
4921  ring r = 0,(x,y,z),dp;
4922  poly F =  x3+y3+z3;
4923  def  B = annfs(F);  setring B;
4924  minIntRoot(BS[1],0);
4925  // So, the minimal integer root is -1
4926  setring r;
4927  def  A = SannfsBM(F);
4928  setring A;
4929  poly F =  x3+y3+z3;
4930  annfspecial(LD,F,-1,3/4); // generic root
4931  annfspecial(LD,F,-1,-2);  // integer but still generic root
4932  annfspecial(LD,F,-1,1);   // exceptional root
4933}
4934*/
4935
4936proc minIntRoot(ideal P, int fact)
4937"USAGE:  minIntRoot(P, fact); P an ideal, fact an int
4938RETURN:  int
4939PURPOSE: minimal integer root of a maximal ideal P
4940NOTE:    if fact==1, P is the result of some 'factorize' call,
4941@*       else P is treated as the result of bernstein::gmssing.lib
4942@*       in both cases without constants and multiplicities
4943EXAMPLE: example minIntRoot; shows examples
4944"
4945{
4946  //    ideal P = factorize(p,1);
4947  // or ideal P = bernstein(F)[1];
4948  intvec vP;
4949  number nP;
4950  int i;
4951  if ( fact )
4952  {
4953    // the result comes from "factorize"
4954    P = normalize(P); // now leadcoef = 1
4955    // TODO: hunt for units and kill then !!!
4957    // nP = leadcoef(P[i]-lead(P[i])); // for 1 var only, extract the coeff
4958  }
4959  else
4960  {
4961    // bernstein deletes -1 usually
4962    P = P, -1;
4963  }
4964  // for both situations:
4965  // now we have an ideal of fractions of type "number"
4966  int sP = size(P);
4967  for(i=1; i<=sP; i++)
4968  {
4970    if ( (nP - int(nP)) == 0 )
4971    {
4972      vP = vP,int(nP);
4973    }
4974  }
4975  if ( size(vP)>=2 )
4976  {
4977    vP = vP[2..size(vP)];
4978  }
4979  sP = -Max(-vP);
4980  if (sP == 0)
4981  {
4982    "Warning: zero root!";
4983  }
4984  return(sP);
4985}
4986example
4987{
4988  "EXAMPLE:"; echo = 2;
4989  ring r   = 0,(x,y),ds;
4990  poly f1  = x*y*(x+y);
4991  ideal I1 = bernstein(f1)[1]; // a local Bernstein poly
4992  minIntRoot(I1,0);
4993  poly  f2  = x2-y3;
4994  ideal I2  = bernstein(f2)[1];
4995  minIntRoot(I2,0);
4996  // now we illustrate the behaviour of factorize
4997  // together with a global ordering
4998  ring r2  = 0,x,dp;
4999  poly f3   = 9*(x+2/3)*(x+1)*(x+4/3); //global b-poly of f1=x*y*(x+y)
5000  ideal I3 = factorize(f3,1);
5001  minIntRoot(I3,1);
5002  // and a more interesting situation
5003  ring  s  = 0,(x,y,z),ds;
5004  poly  f  = x3 + y3 + z3;
5005  ideal I  = bernstein(f)[1];
5006  minIntRoot(I,0);
5007}
5008
5009proc isHolonomic(def M)
5010"USAGE:  isHolonomic(M); M an ideal/module/matrix
5011RETURN:  int, 1 if M is holonomic and 0 otherwise
5012ASSUME: basering is a Weyl algebra
5013PURPOSE: check the modules for the property of holonomy
5014NOTE:    M is holonomic, if 2*dim(M) = dim(R), where R is a
5015ground ring; dim stays for Gelfand-Kirillov dimension
5016EXAMPLE: example isHolonomic; shows examples
5017"
5018{
5019  if ( (typeof(M) != "ideal") && (typeof(M) != "module") && (typeof(M) != "matrix") )
5020  {
5021  //  print(typeof(M));
5022    ERROR("an argument of type ideal/module/matrix expected");
5023  }
5024  if (attrib(M,"isSB")!=1)
5025  {
5026    M = std(M);
5027  }
5028  int dimR = gkdim(std(0));
5029  int dimM = gkdim(M);
5030  return( (dimR==2*dimM) );
5031}
5032example
5033{
5034  "EXAMPLE:"; echo = 2;
5035  ring R = 0,(x,y),dp;
5036  poly F = x*y*(x+y);
5037  def A = annfsBM(F,0);
5038  setring A;
5039  LD;
5040  isHolonomic(LD);
5041  ideal I = std(LD[1]);
5042  I;
5043  isHolonomic(I);
5044}
5045
5046proc reiffen(int p, int q)
5047"USAGE:  reiffen(p, q);  int p, int q
5048RETURN:  ring
5049PURPOSE: set up the polynomial, describing a Reiffen curve
5050NOTE:    activate this ring with the @code{setring} command and find the
5051         curve as a polynomial RC
5052@*   a Reiffen curve is defined as F = x^p + y^q + xy^{q-1}, q >= p+1 >= 5
5053ASSUME: q >= p+1 >= 5. Otherwise an error message is returned
5054EXAMPLE: example reiffen; shows examples
5055"
5056{
5057// a Reiffen curve is defined as
5058// F = x^p + y^q +x*y^{q-1}, q \geq p+1 \geq 5
5059
5060  if ( (p<4) || (q<5) || (q-p<1) )
5061  {
5062    ERROR("Some of conditions p>=4, q>=5 or q>=p+1 is not satisfied!");
5063  }
5064  ring @r = 0,(x,y),dp;
5065  poly RC = y^q +x^p + x*y^(q-1);
5066  export RC;
5067  return(@r);
5068}
5069example
5070{
5071  "EXAMPLE:"; echo = 2;
5072  def r = reiffen(4,5);
5073  setring r;
5074  RC;
5075}
5076
5077proc arrange(int p)
5078"USAGE:  arrange(p);  int p
5079RETURN:  ring
5080PURPOSE: set up the polynomial, describing a hyperplane arrangement
5081NOTE:   must be executed in a ring
5082ASSUME: basering is present
5083EXAMPLE: example arrange; shows examples
5084"
5085{
5086  int UseBasering = 0 ;
5087  if (p<2)
5088  {
5089    ERROR("p>=2 is required!");
5090  }
5091  if ( nameof(basering)!="basering" )
5092  {
5093    if (p > nvars(basering))
5094    {
5095      ERROR("too big p");
5096    }
5097    else
5098    {
5099      def @r = basering;
5100      UseBasering = 1;
5101    }
5102  }
5103  else
5104  {
5105    // no basering found
5106    ERROR("no basering found!");
5107    //    ring @r = 0,(x(1..p)),dp;
5108  }
5109  int i,j,sI;
5110  poly  q = 1;
5111  list ar;
5112  ideal tmp;
5113  tmp = ideal(var(1));
5114  ar[1] = tmp;
5115  for (i = 2; i<=p; i++)
5116  {
5117    // add i-nary sums to the product
5118    ar = indAR(ar,i);
5119  }
5120  for (i = 1; i<=p; i++)
5121  {
5122    tmp = ar[i]; sI = size(tmp);
5123    for (j = 1; j<=sI; j++)
5124    {
5125      q = q*tmp[j];
5126    }
5127  }
5128  if (UseBasering)
5129  {
5130    return(q);
5131  }
5132  //  poly AR = q; export AR;
5133  //  return(@r);
5134}
5135example
5136{
5137  "EXAMPLE:"; echo = 2;
5138  ring X = 0,(x,y,z,t),dp;
5139  poly q = arrange(3);
5140  factorize(q,1);
5141}
5142
5143proc checkRoot(poly F, number a, list #)
5144"USAGE:  checkRoot(f,alpha [,S,eng]);  f a poly, alpha a number, S a string , eng an optional int
5145RETURN:  int
5146PURPOSE: check whether a rational is a root of the global Bernstein polynomial of f (and compute its multiplicity)
5147with the algorithm given in S and with the Groebner basis engine given in eng
5148NOTE:    The annihilator of f^s in D[s] is needed and it is computed according to the algorithm by Briancon and Maisonobe
5149@*       The value of a string S can be
5150@*       'alg1' (default) - for the algorithm 1 of J. Martin-Morales (unpublished)
5151@*       'alg2' - for the algorithm 2 of J. Martin-Morales (unpublished)
5152@*       The output of type int is:
5153@*       'alg1': 1 if -alpha is a root of the global Bernstein polynomial and 0 otherwise
5154@*       'alg2':  the multiplicity of -alpha as a root of the global Bernstein polynomial of f.
5155@*                (If -alpha is not a root, the output is 0)
5156@*       If eng <>0, @code{std} is used for Groebner basis computations,
5157@*       otherwise (and by default) @code{slimgb} is used.
5158DISPLAY: If printlevel=1, progress debug messages will be printed,
5159@*          if printlevel>=2, all the debug messages will be printed.
5160EXAMPLE: example checkRoot; shows examples
5161"
5162{
5163  int eng = 0;
5164  int chs = 0; // choice
5165  string algo = "alg1";
5166  string st;
5167  if ( size(#)>0 )
5168  {
5169   if ( typeof(#[1]) == "string" )
5170   {
5171     st = string(#[1]);
5172     if ( (st == "alg1") || (st == "ALG1") || (st == "Alg1") ||(st == "aLG1"))
5173     {
5174       algo = "alg1";
5175       chs  = 1;
5176     }
5177     if ( (st == "alg2") || (st == "ALG2") || (st == "Alg2") || (st == "aLG2"))
5178     {
5179       algo = "alg2";
5180       chs  = 1;
5181     }
5182     if (chs != 1)
5183     {
5184       // incorrect value of S
5185       print("Incorrect algorithm given, proceed with the default alg1 of J. MartÃ­n-Morales");
5186       algo = "alg1";
5187     }
5188     // second arg
5189     if (size(#)>1)
5190     {
5191       // exists 2nd arg
5192       if ( typeof(#[2]) == "int" )
5193       {
5194         // the case: given alg, given engine
5195         eng = int(#[2]);
5196       }
5197       else
5198       {
5199         eng = 0;
5200       }
5201     }
5202     else
5203     {
5204       // no second arg
5205       eng = 0;
5206     }
5207   }
5208   else
5209   {
5210     if ( typeof(#[1]) == "int" )
5211     {
5212       // the case: default alg, engine
5213       eng = int(#[1]);
5214       // algo = "alg1";  //is already set
5215     }
5216     else
5217     {
5218       // incorr. 1st arg
5219       algo = "alg1";
5220     }
5221   }
5222  }
5223  // size(#)=0, i.e. there is no algorithm and no engine given
5224  //  eng = 0; algo = "alg1";  //are already set
5225  // int ppl = printlevel-voice+2;
5226  printlevel=printlevel+1;
5227  def save = basering;
5228  def @A = SannfsBM(F);
5229  setring @A;
5230  poly F = imap(save,F);
5231  number a = imap(save,a);
5232  if ( algo=="alg1")
5233  {
5234    int output = checkRoot1(LD,F,a,eng);
5235  }
5236  else
5237  {
5238    if ( algo=="alg2")
5239    {
5240      int output = checkRoot2(LD,F,a,eng);
5241    }
5242  }
5243  printlevel=printlevel-1;
5244  return(output);
5245}
5246example
5247{
5248  "EXAMPLE:"; echo = 2;
5249  printlevel=0;
5250  ring r = 0,(x,y),Dp;
5251  poly F = x^4+y^5+x*y^4;
5252  checkRoot(F,11/20);    // -11/20 is a root of bf
5253  poly G = x*y;
5254  checkRoot(G,1,"alg2"); // -1 is a root of bg with multiplicity 2
5255}
5256
5257proc checkRoot1(ideal I, poly F, number a, list #)
5258"USAGE:  checkRoot1(I,f,alpha [,eng]);  I an ideal, f a poly, alpha a number, eng an optional int
5259ASSUME:  I is the annihilator of f^s in D[s], f is a polynomial in K[_x]
5260RETURN:  int, 1 if -alpha is a root of the global Bernstein polynomial of f and 0 otherwise
5261PURPOSE: check whether a rational is a root of the global Bernstein polynomial of f
5262NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
5263@*       otherwise (and by default) @code{slimgb} is used.
5264DISPLAY: If printlevel=1, progress debug messages will be printed,
5265@*          if printlevel>=2, all the debug messages will be printed.
5266EXAMPLE: example checkRoot1; shows examples
5267"
5268{
5269  int eng = 0;
5270  if ( size(#)>0 )
5271  {
5272    if ( typeof(#[1]) == "int" )
5273    {
5274      eng = int(#[1]);
5275    }
5276  }
5277  int ppl = printlevel-voice+2;
5278  dbprint(ppl,"// -0-1- starting the procedure checkRoot1");
5279  def save = basering;
5280  int N = nvars(basering);
5281  int Nnew = N-1;
5282  int n = Nnew / 2;
5283  int i;
5284  string s;
5285  list RL = ringlist(basering);
5286  list L, Lord;
5287  list tmp;
5288  intvec iv;
5289  L[1] = RL[1];  // char
5290  L[4] = RL[4];  // char, minpoly
5291  // check whether basering is D[s]=K(_x,_Dx,s)
5292  list Name = RL[2];
5293//   for (i=1; i<=n; i++)
5294//   {
5295//     if ( bracket(var(i+n),var(i))!=1 )
5296//     {
5297//       ERROR("basering should be D[s]=K(_x,_Dx,s)");
5298//     }
5299//   }
5300  if ( Name[N]!="s" )
5301  {
5302    ERROR("the last variable of basering should be s");
5303  }
5304  // now, create the new vars
5305  list NName;
5306  for (i=1; i<=Nnew; i++)
5307  {
5308    NName[i] = Name[i];
5309  }
5310  L[2] = NName;
5311  kill Name,NName;
5312  // block ord (dp);
5313  tmp[1] = "dp"; // string
5314  s = "iv=";
5315  for (i=1; i<=Nnew; i++)
5316  {
5317    s = s+"1,";
5318  }
5319  s[size(s)]=";";
5320  execute(s);
5321  kill s;
5322  tmp[2]  = iv;
5323  Lord[1] = tmp;
5324  tmp[1]  = "C";
5325  iv      = 0;
5326  tmp[2]  = iv;
5327  Lord[2] = tmp;
5328  tmp     = 0;
5329  L[3]    = Lord;
5330  // we are done with the list
5331  def @R@ = ring(L);
5332  setring @R@;
5333  matrix @D[Nnew][Nnew];
5334  for (i=1; i<=n; i++)
5335  {
5336    @D[i,i+n]=1;
5337  }
5338  def @R = nc_algebra(1,@D);
5339  setring @R;
5340  kill @R@;
5341  dbprint(ppl,"// -1-1- the ring @R(_x,_Dx) is ready");
5342  dbprint(ppl-1, S);
5343  // create the ideal K = ann_D[s](f^s)_{s=-alpha} + < f >
5344  setring save;
5345  ideal K = subst(I,s,-a);
5346  dbprint(ppl,"// -1-2- the variable s has been substituted by "+string(-a));
5347  dbprint(ppl-1, K);
5348  K = NF(K,std(F));
5350  for (i=1; i<=ncols(K); i++)
5351  {
5353    {
5354      K[i] = -K[i];
5355    }
5356  }
5357  K = K,F;
5358  // ------------ the ideal K is ready ------------
5359  setring @R;
5360  ideal K = imap(save,K);
5361  dbprint(ppl,"// -1-3- starting the computation of a Groebner basis of K in @R");
5362  dbprint(ppl-1, K);
5363  ideal G = engine(K,eng);
5364  dbprint(ppl,"// -1-4- the Groebner basis has been computed");
5365  dbprint(ppl-1, G);
5366  return(G[1]!=1);
5367}
5368example
5369{
5370  "EXAMPLE:"; echo = 2;
5371  ring r = 0,(x,y),Dp;
5372  poly F = x^4+y^5+x*y^4;
5373  printlevel = 0;
5374  def A = Sannfs(F);
5375  setring A;
5376  poly F = imap(r,F);
5377  checkRoot1(LD,F,31/20);   // -31/20 is not a root of bs
5378  checkRoot1(LD,F,11/20);   // -11/20 is a root of bs
5379}
5380
5381proc checkRoot2 (ideal I, poly F, number a, list #)
5382"USAGE:  checkRoot2(I,f,alpha [,eng]);  I an ideal, f a poly, alpha a number, eng an optional int
5383ASSUME:  I is the annihilator of f^s in D[s], f is a polynomial in K[_x]
5384RETURN:  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
5385PURPOSE: 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]
5386NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
5387@*       otherwise (and by default) @code{slimgb} is used.
5388DISPLAY: If printlevel=1, progress debug messages will be printed,
5389@*          if printlevel>=2, all the debug messages will be printed.
5390EXAMPLE: example checkRoot2; shows examples
5391"
5392{
5393  int eng = 0;
5394  if ( size(#)>0 )
5395  {
5396    if ( typeof(#[1]) == "int" )
5397    {
5398      eng = int(#[1]);
5399    }
5400  }
5401  int ppl = printlevel-voice+2;
5402  dbprint(ppl,"// -0-1- starting the procedure checkRoot2");
5403  def save = basering;
5404  int N = nvars(basering);
5405  int n = (N-1) / 2;
5406  int i;
5407  string s;
5408  list RL = ringlist(basering);
5409  list L, Lord;
5410  list tmp;
5411  intvec iv;
5412  L[1] = RL[1];  // char
5413  L[4] = RL[4];  // char, minpoly
5414  // check whether basering is D[s]=K(_x,_Dx,s)
5415  list Name = RL[2];
5416  for (i=1; i<=n; i++)
5417  {
5418    if ( bracket(var(i+n),var(i))!=1 )
5419    {
5420      ERROR("basering should be D[s]=K(_x,_Dx,s)");
5421    }
5422  }
5423  if ( Name[N]!="s" )
5424  {
5425    ERROR("the last variable of basering should be s");
5426  }
5427  // now, create the new vars
5428  L[2] = Name;
5429  kill Name;
5430  // block ord (dp);
5431  tmp[1] = "dp"; // string
5432  s = "iv=";
5433  for (i=1; i<=N; i++)
5434  {
5435    s = s+"1,";
5436  }
5437  s[size(s)]=";";
5438  execute(s);
5439  kill s;
5440  tmp[2]  = iv;
5441  Lord[1] = tmp;
5442  tmp[1]  = "C";
5443  iv      = 0;
5444  tmp[2]  = iv;
5445  Lord[2] = tmp;
5446  tmp     = 0;
5447  L[3]    = Lord;
5448  // we are done with the list
5449  def @R@ = ring(L);
5450  setring @R@;
5451  matrix @D[N][N];
5452  for (i=1; i<=n; i++)
5453  {
5454    @D[i,i+n]=1;
5455  }
5456  def @R = nc_algebra(1,@D);
5457  setring @R;
5458  kill @R@;
5459  dbprint(ppl,"// -1-1- the ring @R(_x,_Dx,s) is ready");
5460  dbprint(ppl-1, @R);
5461  // now, continue with the algorithm
5462  ideal I = imap(save,I);
5463  poly F = imap(save,F);
5464  number a = imap(save,a);
5465  ideal II = NF(I,std(F));
5467  for (i=1; i<=ncols(II); i++)
5468  {
5470    {
5471      II[i] = -II[i];
5472    }
5473  }
5474  ideal J,G;
5475  int m;  // the output (multiplicity)
5476  dbprint(ppl,"// -2- starting the bucle");
5477  for (i=0; i<=n; i++)  // the multiplicity has to be <= n
5478  {
5479    // create the ideal Ji = ann_D[s](f^s) + < f, (s+alpha)^{i+1} >
5480    // (s+alpha)^i in Ji <==> -alpha is a root with multiplicity >= i
5481    J = II,F,(s+a)^(i+1);
5482    // ------------ the ideal Ji is ready -----------
5483    dbprint(ppl,"// -2-"+string(i+1)+"-1- starting the computation of a Groebner basis of J"+string(i)+" in @R");
5484    dbprint(ppl-1, J);
5485    G = engine(J,eng);
5486    dbprint(ppl,"// -2-"+string(i+1)+"-2- the Groebner basis has been computed");
5487    dbprint(ppl-1, G);
5488    if ( NF((s+a)^i,G)==0 )
5489    {
5490      dbprint(ppl,"// -2-"+string(i+1)+"-3- the number "+string(-a)+" has not multiplicity "+string(i+1));
5491      m = i;
5492      break;
5493    }
5494    dbprint(ppl,"// -2-"+string(i+1)+"-3- the number "+string(-a)+" has multiplicity at least "+string(i+1));
5495  }
5496  dbprint(ppl,"// -3- the bucle has finished");
5497  return(m);
5498}
5499example
5500{
5501  "EXAMPLE:"; echo = 2;
5502  ring r = 0,(x,y,z),Dp;
5503  poly F = x*y*z;
5504  printlevel = 0;
5505  def A = Sannfs(F);
5506  setring A;
5507  poly F = imap(r,F);
5508  checkRoot2(LD,F,1);    // -1 is a root of bs with multiplicity 3
5509  checkRoot2(LD,F,1/3);  // -1/3 is not a root
5510}
5511
5512proc checkFactor(ideal I, poly F, poly q, list #)
5513"USAGE:  checkFactor(I,f,qs [,eng]);  I an ideal, f a poly, qs a poly, eng an optional int
5514ASSUME:  I is the output of Sannfs, SannfsBM, SannfsLOT or SannfsOT,
5515         f is a polynomial in K[_x], qs is a polynomial in K[s]
5516RETURN:  int, 1 if qs is a factor of the global Bernstein polynomial of f and 0 otherwise
5517PURPOSE: 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]
5518NOTE:    If eng <>0, @code{std} is used for Groebner basis computations,
5519@*       otherwise (and by default) @code{slimgb} is used.
5520DISPLAY: If printlevel=1, progress debug messages will be printed,
5521@*          if printlevel>=2, all the debug messages will be printed.
5522EXAMPLE: example checkFactor; shows examples
5523"
5524{
5525  int eng = 0;
5526  if ( size(#)>0 )
5527  {
5528    if ( typeof(#[1]) == "int" )
5529    {
5530      eng = int(#[1]);
5531    }
5532  }
5533  int ppl = printlevel-voice+2;
5534  def @R2 = basering;
5535  int N = nvars(@R2);
5536  int i;
5537  // we're in D_n[s], where the elim ord for s is set
5538  dbprint(ppl,"// -0-1- starting the procedure checkFactor");
5539  dbprint(ppl,"// -1-1- the ring @R2(_x,_Dx,s) is ready");
5540  dbprint(ppl-1, @R2);
5541  // create the ideal J = ann_D[s](f^s) + < f,q >
5542  ideal J = NF(I,std(F));
5544  for (i=1; i<=ncols(J); i++)
5545  {
5547    {
5548      J[i] = -J[i];
5549    }
5550  }
5551  J = J,F,q;
5552  // ------------ the ideal J is ready -----------
5553  dbprint(ppl,"// -1-2- starting the elimination of _x,_Dx in @R2");
5554  dbprint(ppl-1, J);
5555  ideal G = engine(J,eng);
5556  ideal K = nselect(G,1..N-1);
5557  kill J,G;
5558  dbprint(ppl,"// -1-3- _x,_Dx are eliminated");
5559  dbprint(ppl-1, K);
5560  //q is a factor of bs iff K = < q >
5561  //K = normalize(K);
5562  //q = normalize(q);
5563  //return( (K[1]==q) );
5564  return( NF(K[1],std(q))==0 );
5565}
5566example
5567{
5568  "EXAMPLE:"; echo = 2;
5569  ring r = 0,(x,y),Dp;
5570  poly F = x^4+y^5+x*y^4;
5571  printlevel = 0;
5572  def A = Sannfs(F);
5573  setring A;
5574  poly F = imap(r,F);
5575  checkFactor(LD,F,20*s+31);     // -31/20 is not a root of bs
5576  checkFactor(LD,F,20*s+11);     // -11/20 is a root of bs
5577  checkFactor(LD,F,(20*s+11)^2); // the multiplicity of -11/20 is 1
5578}
5579
5580proc varnum(string s)
5581"USAGE:  varnum(s);  string s
5582RETURN:  int
5583PURPOSE: returns the number of the variable with the name s
5584among the variables of basering or 0 if there is no such variable
5585EXAMPLE: example varnum; shows examples
5586"
5587{
5588  int i;
5589  for (i=1; i<= nvars(basering); i++)
5590  {
5591    if ( string(var(i)) == s )
5592    {
5593      return(i);
5594    }
5595  }
5596  return(0);
5597}
5598example
5599{
5600  "EXAMPLE:"; echo = 2;
5601  ring X = 0,(x,y1,z(0),tTa),dp;
5602  varnum("z(0)");
5603  varnum("tTa");
5604  varnum("xyz");
5605}
5606
5607static proc indAR(list L, int n)
5608"USAGE:  indAR(L,n);  list L, int n
5609RETURN:  list
5610PURPOSE: computes arrangement inductively, using L and var(n) as the
5611next variable
5612ASSUME: L has a structure of an arrangement
5613EXAMPLE: example indAR; shows examples
5614"
5615{
5616  if ( (n<2) || (n>nvars(basering)) )
5617  {
5618    ERROR("incorrect n");
5619  }
5620  int sl = size(L);
5621  list K;
5622  ideal tmp;
5623  poly @t = L[sl][1] + var(n); //1 elt
5624  K[sl+1] = ideal(@t);
5625  tmp = L[1]+var(n);
5626  K[1] = tmp; tmp = 0;
5627  int i,j,sI;
5628  ideal I;
5629  for(i=sl; i>=2; i--)
5630  {
5631    I = L[i-1]; sI = size(I);
5632    for(j=1; j<=sI; j++)
5633    {
5634      I[j] = I[j] + var(n);
5635    }
5636    tmp  = L[i],I;
5637    K[i] = tmp;
5638    I = 0; tmp = 0;
5639  }
5640  kill I; kill tmp;
5641  return(K);
5642}
5643example
5644{
5645  "EXAMPLE:"; echo = 2;
5646  ring r = 0,(x,y,z,t,v),dp;
5647  list L;
5648  L[1] = ideal(x);
5649  list K = indAR(L,2);
5650  K;
5651  list M = indAR(K,3);
5652  M;
5653  M = indAR(M,4);
5654  M;
5655}
5656
5657/// ****** EXAMPLES ************
5658
5659/*
5660
5661//static proc exCheckGenericity()
5662{
5663  LIB "control.lib";
5664  ring r = (0,a,b,c),x,dp;
5665  poly p = (x-a)*(x-b)*(x-c);
5666  def A = annfsBM(p);
5667  setring A;
5668  ideal J = slimgb(LD);
5669  matrix T = lift(LD,J);
5670  T = normalize(T);
5671  genericity(T);
5672  // 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)
5673  // genericity: g = a2-ab-ac+b2-bc+c2 =0
5674  // g = (a -(b+c)/2)^2 + (3/4)*(b-c)^2;
5675  // g ==0 <=> a=b=c
5676  // indeed, Ann = (x-a)^2*(x*Dx+3*t*Dt+(-a)*Dx+3)
5677  // --------------------------------------------
5678  // BUT a direct computation shows
5679  // when a=b=c,
5680  // Ann = x*Dx+3*t*Dt+(-a)*Dx+3
5681}
5682
5683//static proc exOT_17()
5684{
5685  // Oaku-Takayama, p.208
5686  ring R = 0,(x,y),dp;
5687  poly F = x^3-y^2; // x^2+x*y+y^2;
5688  option(prot);
5689  option(mem);
5690  //  option(redSB);
5691  def A = annfsOT(F,0);
5692  setring A;
5693  LD;
5694  gkdim(LD); // a holonomic check
5695  //  poly F = x^3-y^2; // = x^7 - y^5; // x^3-y^4; // x^5 - y^4;
5696}
5697
5698//static proc exOT_16()
5699{
5700  // Oaku-Takayama, p.208
5701  ring R = 0,(x),dp;
5702  poly F = x*(1-x);
5703  option(prot);
5704  option(mem);
5705  //  option(redSB);
5706  def A = annfsOT(F,0);
5707  setring A;
5708  LD;
5709  gkdim(LD); // a holonomic check
5710}
5711
5712//static proc ex_bcheck()
5713{
5714  ring R = 0,(x,y),dp;
5715  poly F = x*y*(x+y);
5716  option(prot);
5717  option(mem);
5718  int eng = 0;
5719  //  option(redSB);
5720  def A = annfsOT(F,eng);
5721  setring A;
5722  LD;
5723}
5724
5725//static proc ex_bcheck2()
5726{
5727  ring R = 0,(x,y),dp;
5728  poly F = x*y*(x+y);
5729  int eng = 0;
5730  def A = annfsBM(F,eng);
5731  setring A;
5732  LD;
5733}
5734
5735//static proc ex_BMI()
5736{
5737  // a hard example
5738  ring r = 0,(x,y),Dp;
5739  poly F1 = (x2-y3)*(x3-y2);
5740  poly F2 = (x2-y3)*(xy4+y5+x4);
5741  ideal F = F1,F2;
5742  def A = annfsBMI(F);
5743  setring A;
5744  LD;
5745  BS;
5746}
5747
5748//static proc ex2_BMI()
5749{
5750  // this example was believed to be intractable in 2005 by Gago-Vargas, Castro and Ucha
5751  ring r = 0,(x,y),Dp;
5752  option(prot);
5753  option(mem);
5754  ideal F = x2+y3,x3+y2;
5755  printlevel = 2;
5756  def A = annfsBMI(F);
5757  setring A;
5758  LD;
5759  BS;
5760}
5761
5762//static proc ex_operatorBM()
5763{
5764  ring r = 0,(x,y,z,w),Dp;
5765  poly F = x^3+y^3+z^2*w;
5766  printlevel = 0;
5767  def A = operatorBM(F);
5768  setring A;
5769  F; // the original polynomial itself
5770  LD; // generic annihilator
5771  LD0; // annihilator
5772  bs; // normalized Bernstein poly
5773  BS; // root and multiplicities of the Bernstein poly
5774  PS; // the operator, s.t. PS*F^{s+1} = bs*F^s mod LD
5775  reduce(PS*F-bs,LD); // check the property of PS
5776}
5777
5778*/
Note: See TracBrowser for help on using the repository browser.