# source:git/Singular/LIB/dmod.lib@8b4657a

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