source: git/Singular/LIB/dmodapp.lib @ ab3db62

spielwiese
Last change on this file since ab3db62 was ab3db62, checked in by Viktor Levandovskyy <levandov@…>, 14 years ago
*levandov: changes by daniel andres, new functions, doc enhancements git-svn-id: file:///usr/local/Singular/svn/trunk@11958 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 49.3 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id: dmodapp.lib,v 1.31 2009-07-07 20:46:51 levandov Exp $";
3category="Noncommutative";
4info="
5LIBRARY: dmodapp.lib     Applications of algebraic D-modules
6AUTHORS: Viktor Levandovskyy,      levandov@math.rwth-aachen.de
7@*             Daniel Andres, daniel.andres@math.rwth-aachen.de
8
9GUIDE: Let K be a field of characteristic 0, R = K[x1,..xN] and
10@* D be the Weyl algebra in variables x1,..xN,d1,..dN.
11@* In this library there are the following procedures for algebraic D-modules:
12@* - localization of a holonomic module D/I with respect to a mult. closed set
13@* of all powers of a given polynomial F from R. Our aim is to compute an
14@* ideal L in D, such that D/L is a presentation of a localized module. Such L
15@* always exists, since such localizations are known to be holonomic and thus
16@* cyclic modules. The procedures for the localization are DLoc,SDLoc and DLoc0.
17@*
18@* - annihilator in D of a given polynomial F from R as well as
19@* of a given rational function G/F from Quot(R). These can be computed via
20@* procedures annPoly resp. annRat.
21@*
22@* - initial form and initial ideals in Weyl algebras with respect to a given
23@* weight vector can be computed with  inForm, initialMalgrange, initialIdealW.
24@*
25@* - appelF1, appelF2 and appelF4 return ideals in parametric Weyl algebras,
26@* which annihilate corresponding Appel hypergeometric functions.
27
28REFERENCES:
29@* (SST) Saito, Sturmfels, Takayama 'Groebner Deformations of Hypergeometric
30@*         Differential Equations', Springer, 2000
31@* (ONW) Oaku, Takayama, Walther 'A Localization Algorithm for D-modules', 2000
32
33MAIN PROCEDURES:
34
35annPoly(f);   annihilator of a polynomial f in the corr. Weyl algebra
36annRat(f,g);  annihilator of a rational function f/g in the corr. Weyl algebra
37DLoc(I,F);     presentation of the localization of D/I w.r.t. f^s
38SDLoc(I, F);  a generic presentation of the localization of D/I w.r.t. f^s
39DLoc0(I, F);  presentation of the localization of D/I w.r.t. f^s, based on SDLoc
40
41initialMalgrange(f[,s,t,v]); Groebner basis of the initial Malgrange ideal for f
42initialIdealW(I,u,v[,s,t]); initial ideal of a given ideal w.r.t. given weights
43inForm(f,w);     initial form of a poly/ideal w.r.t. a given weight
44isFsat(I, F);       check whether the ideal I is F-saturated
45
46AUXILIARY PROCEDURES:
47
48bFactor(F);  computes the roots of irreducible factors of an univariate poly
49appelF1();      create an ideal annihilating Appel F1 function
50appelF2();      create an ideal annihilating Appel F2 function
51appelF4();      create an ideal annihilating Appel F4 function
52engine(I,i);     computes a Groebner basis with the algorithm given by i
53poly2list(f);   decompose a polynomial into a list of terms and exponents
54fl2poly(L,s);  reconstruct a monic univariate polynomial from its factorization
55insertGenerator(id,p[,k]); insert an element into an ideal/module
56deleteGenerator(id,k); delete the k-th element from an ideal/module
57
58
59SEE ALSO: dmod_lib, gmssing_lib, bfun_lib
60";
61
62LIB "poly.lib";
63LIB "sing.lib";
64LIB "primdec.lib";
65LIB "dmod.lib";    // loads e.g. nctools.lib
66LIB "bfun.lib";    //formerly LIB "bfct.lib";
67LIB "nctools.lib"; // for isWeyl etc
68LIB "gkdim.lib";
69
70// todo: complete and include into above list
71// charVariety(I);       compute the characteristic variety of the ideal I
72// charInfo(); ???
73
74
75///////////////////////////////////////////////////////////////////////////////
76// testing for consistency of the library:
77proc testdmodapp()
78{
79  example initialIdealW;
80  example initialMalgrange;
81  example DLoc;
82  example DLoc0;
83  example SDLoc;
84  example inForm;
85  example isFsat;
86  example annRat;
87  example annPoly;
88  example appelF1;
89  example appelF2;
90  example appelF4;
91  example poly2list;
92  example fl2poly;
93  example insertGenerator;
94  example deleteGenerator;
95  example bFactor;
96}
97
98proc inForm (ideal I, intvec w)
99"USAGE:  inForm(I,w);  I ideal, w intvec
100RETURN:  the initial form of I w.r.t. the weight vector w
101PURPOSE: computes the initial form of an ideal w.r.t. a given weight vector
102NOTE:  the size of the weight vector must be equal to the number of variables
103@*     of the basering.
104EXAMPLE: example inForm; shows examples
105"
106{
107  if (size(w) != nvars(basering))
108  {
109    ERROR("weight vector has wrong dimension");
110  }
111  if (I == 0)
112  {
113    return(I);
114  }
115  int j,i,s,m;
116  list l;
117  poly g;
118  ideal J;
119  for (j=1; j<=ncols(I); j++)
120  {
121    l = poly2list(I[j]);
122    m = scalarProd(w,l[1][1]);
123    g = l[1][2];
124    for (i=2; i<=size(l); i++)
125    {
126      s = scalarProd(w,l[i][1]);
127      if (s == m)
128      {
129        g = g + l[i][2];
130      }
131      else
132      {
133        if (s > m)
134        {
135          m = s;
136          g = l[i][2];
137        }
138      }
139    }
140    J[j] = g;
141  }
142  return(J);
143}
144example
145{
146  "EXAMPLE:"; echo = 2;
147  ring @D = 0,(x,y,Dx,Dy),dp;
148  def D = Weyl();
149  setring D;
150  poly F = 3*x^2*Dy+2*y*Dx;
151  poly G = 2*x*Dx+3*y*Dy+6;
152  ideal I = F,G;
153  intvec w1 = -1,-1,1,1;
154  intvec w2 = -1,-2,1,2;
155  intvec w3 = -2,-3,2,3;
156  inForm(I,w1);
157  inForm(I,w2);
158  inForm(I,w3);
159}
160
161/*
162
163proc charVariety(ideal I)
164"USAGE:  charVariety(I);  I an ideal
165RETURN:  ring
166PURPOSE: compute the characteristic variety of a D-module D/I
167STATUS: experimental, todo
168ASSUME: the ground ring is the Weyl algebra with x's before d's
169NOTE:    activate the output ring with the @code{setring} command.
170@*       In the output (in a commutative ring):
171@*       - the ideal CV is the characteristic variety char(I)
172@*       If @code{printlevel}=1, progress debug messages will be printed,
173@*       if @code{printlevel}>=2, all the debug messages will be printed.
174EXAMPLE: example charVariety; shows examples
175"
176{
177  // 1. introduce the weights 0, 1
178  def save = basering;
179  list LL = ringlist(save);
180  list L;
181  int i;
182  for(i=1;i<=4;i++)
183  {
184    L[i] = LL[i];
185  }
186  list OLD = L[3];
187  list NEW; list tmp;
188  tmp[1] = "a";  // string
189  intvec iv;
190  int N = nvars(basering); N = N div 2;
191  for(i=N+1; i<=2*N; i++)
192  {
193    iv[i] = 1;
194  }
195  tmp[2] = iv;
196  NEW[1] = tmp;
197  for (i=2; i<=size(OLD);i++)
198  {
199    NEW[i] = OLD[i-1];
200  }
201  L[3] = NEW;
202  list ncr =ncRelations(save);
203  matrix @C = ncr[1];
204  matrix @D = ncr[2];
205  def @U = ring(L);
206  // 2. create the commutative ring
207  setring save;
208  list CL;
209  for(i=1;i<=4;i++)
210  {
211    CL[i] = L[i];
212  }
213  CL[3] = OLD;
214  def @CU = ring(CL);
215  // comm ring is ready
216  setring @U;
217  // make @U noncommutative
218  matrix @C = imap(save,@C);
219  matrix @D = imap(save,@D);
220  def @@U = nc_algebra(@C,@D);
221  setring @@U; kill @U;
222  // 2. compute Groebner basis
223  ideal I = imap(save,I);
224  //  I = groebner(I);
225  I = slimgb(I); // a bug?
226  setring @CU;
227  ideal CV = imap(@@U,I);
228  //  CV = groebner(CV); // cosmetics
229  CV = slimgb(CV);
230  export CV;
231  return(@CU);
232}
233example
234{
235  "EXAMPLE:"; echo = 2;
236  ring r = 0,(x,y),Dp;
237  poly F = x3-y2;
238  printlevel = 0;
239  def A  = annfs(F);
240  setring A; // Weyl algebra
241  LD; // the ideal
242  def CA = charVariety(LD);
243  setring CA;
244  CV;
245  dim(CV);
246}
247
248/*
249
250// TODO
251
252/*
253proc charInfo(ideal I)
254"USAGE:  charInfo(I);  I an ideal
255RETURN:  ring
256STATUS: experimental, todo
257PURPOSE: compute the characteristic information for I
258ASSUME: the ground ring is the Weyl algebra with x's before d's
259NOTE:    activate the output ring with the @code{setring} command.
260@*       In the output (in a commutative ring):
261@*       - the ideal CV is the characteristic variety char(I)
262@*       - the ideal SL is the singular locus of char(I)
263@*       - the list PD is the primary decomposition of char(I)
264@*       If @code{printlevel}=1, progress debug messages will be printed,
265@*       if @code{printlevel}>=2, all the debug messages will be printed.
266EXAMPLE: example annfs; shows examples
267"
268{
269  def save = basering;
270  def @A = charVariety(I);
271  setring @A;
272  // run slocus
273  // run primdec
274}
275*/
276
277
278proc appelF1()
279"USAGE:  appelF1();
280RETURN:  ring  (and exports an ideal into it)
281PURPOSE: define the ideal in a parametric Weyl algebra,
282@* which annihilates Appel F1 hypergeometric function
283NOTE: the ideal called  IAppel1 is exported to the output ring
284EXAMPLE: example appelF1; shows examples
285"
286{
287  // Appel F1, d = b', SST p.48
288  ring @r = (0,a,b,c,d),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp);
289  matrix @D[4][4];
290  @D[1,3]=1; @D[2,4]=1;
291  def @S = nc_algebra(1,@D);
292  setring @S;
293  ideal IAppel1 =
294    (x*Dx)*(x*Dx+y*Dy+c-1) - x*(x*Dx+y*Dy+a)*(x*Dx+b),
295    (y*Dy)*(x*Dx+y*Dy+c-1) - y*(x*Dx+y*Dy+a)*(y*Dy+d),
296    (x-y)*Dx*Dy - d*Dx + b*Dy;
297  export IAppel1;
298  kill @r;
299  return(@S);
300}
301example
302{
303  "EXAMPLE:"; echo = 2;
304  def A = appelF1();
305  setring A;
306  IAppel1;
307}
308
309proc appelF2()
310"USAGE:  appelF2();
311RETURN:  ring (and exports an ideal into it)
312PURPOSE: define the ideal in a parametric Weyl algebra,
313@* which annihilates Appel F2 hypergeometric function
314NOTE: the ideal called  IAppel2 is exported to the output ring
315EXAMPLE: example appelF2; shows examples
316"
317{
318  // Appel F2, c = b', SST p.85
319  ring @r = (0,a,b,c),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp);
320  matrix @D[4][4];
321  @D[1,3]=1; @D[2,4]=1;
322  def @S = nc_algebra(1,@D);
323  setring @S;
324  ideal IAppel2 =
325    (x*Dx)^2 - x*(x*Dx+y*Dy+a)*(x*Dx+b),
326    (y*Dy)^2 - y*(x*Dx+y*Dy+a)*(y*Dy+c);
327  export IAppel2;
328  kill @r;
329  return(@S);
330}
331example
332{
333  "EXAMPLE:"; echo = 2;
334  def A = appelF2();
335  setring A;
336  IAppel2;
337}
338
339proc appelF4()
340"USAGE:  appelF4();
341RETURN:  ring  (and exports an ideal into it)
342PURPOSE: define the ideal in a parametric Weyl algebra,
343@* which annihilates Appel F4 hypergeometric function
344NOTE: the ideal called  IAppel4 is exported to the output ring
345EXAMPLE: example appelF4; shows examples
346"
347{
348  // Appel F4, d = c', SST, p. 39
349  ring @r = (0,a,b,c,d),(x,y,Dx,Dy),(a(0,0,1,1),a(0,0,1,0),dp);
350  matrix @D[4][4];
351  @D[1,3]=1; @D[2,4]=1;
352  def @S = nc_algebra(1,@D);
353  setring @S;
354  ideal IAppel4 =
355    Dx*(x*Dx+c-1) - (x*Dx+y*Dy+a)*(x*Dx+y*Dy+b),
356    Dy*(y*Dy+d-1) - (x*Dx+y*Dy+a)*(x*Dx+y*Dy+b);
357  export IAppel4;
358  kill @r;
359  return(@S);
360}
361example
362{
363  "EXAMPLE:"; echo = 2;
364  def A = appelF4();
365  setring A;
366  IAppel4;
367}
368
369proc poly2list (poly f)
370"USAGE:  poly2list(f); f a poly
371RETURN:  list of exponents and corresponding terms of f
372PURPOSE: convert a polynomial to a list of exponents and corresponding terms
373EXAMPLE: example poly2list; shows examples
374"
375{
376  list l;
377  int i = 1;
378  if (f == 0) // just for the zero polynomial
379  {
380    l[1] = list(leadexp(f), lead(f));
381  }
382  else { l[size(f)] = list(); } // memory pre-allocation
383  while (f != 0)
384  {
385    l[i] = list(leadexp(f), lead(f));
386    f = f - lead(f);
387    i++;
388  }
389  return(l);
390}
391example
392{
393  "EXAMPLE:"; echo = 2;
394  ring r = 0,x,dp;
395  poly F = x;
396  poly2list(F);
397  ring r2 = 0,(x,y),dp;
398  poly F = x2y+x*y2;
399  poly2list(F);
400}
401
402proc isFsat(ideal I, poly F)
403"USAGE:  isFsat(I, F);  I an ideal, F a poly
404RETURN: int
405PURPOSE: check whether the ideal I is F-saturated
406NOTE:    1 is returned if I is F-saturated, otherwise 0 is returned.
407@*   we check indeed that Ker(D --F--> D/I) is (0)
408EXAMPLE: example isFsat; shows examples
409"
410{
411  /* checks whether I is F-saturated, that is Ke  (D -F-> D/I) is 0 */
412  /* works in any algebra */
413  /*  for simplicity : later check attrib */
414  /* returns -1 if true */
415  if (attrib(I,"isSB")!=1)
416  {
417    I = groebner(I);
418  }
419  matrix @M = matrix(I);
420  matrix @F[1][1] = F;
421  module S = modulo(@F,@M);
422  S = NF(S,I);
423  S = groebner(S);
424  return( (gkdim(S) == -1) );
425}
426example
427{
428  "EXAMPLE:"; echo = 2;
429  ring r = 0,(x,y),dp;
430  poly G = x*(x-y)*y;
431  def A = annfs(G);
432  setring A;
433  poly F = x3-y2;
434  isFsat(LD,F);
435  ideal J = LD*F;
436  isFsat(J,F);
437}
438
439proc DLoc(ideal I, poly F)
440"USAGE:  DLoc(I, F);  I an ideal, F a poly
441RETURN: nothing (exports objects instead)
442ASSUME: the basering is a Weyl algebra
443PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s
444NOTE:    In the basering, the following objects are exported:
445@* the ideal LD0 (in Groebner basis) is the presentation of the localization
446@* the list BS contains roots with multiplicities of Bernstein polynomial of (D/I)_f.
447DISPLAY: If printlevel=1, progress debug messages will be printed,
448@*       if printlevel>=2, all the debug messages will be printed.
449EXAMPLE: example DLoc; shows examples
450"
451{
452  /* runs SDLoc and DLoc0 */
453  /* assume: run from Weyl algebra */
454  if (dmodappassumeViolation())
455  {
456    ERROR("Basering is inappropriate: characteristic>0 or qring present");
457  }
458  if (!isWeyl())
459  {
460    ERROR("Basering is not a Weyl algebra");
461  }
462  if (defined(LD0) || defined(BS))
463  {
464    ERROR("Reserved names LD0 and/or BS are used. Please rename the objects.");
465  }
466  int old_printlevel = printlevel;
467  printlevel=printlevel+1;
468  def @R = basering;
469  def @R2 = SDLoc(I,F);
470  setring @R2;
471  poly F = imap(@R,F);
472  def @R3 = DLoc0(LD,F);
473  setring @R3;
474  ideal bs = BS[1];
475  intvec m = BS[2];
476  setring @R;
477  ideal LD0 = imap(@R3,LD0);
478  export LD0;
479  ideal bs = imap(@R3,bs);
480  list BS; BS[1] = bs; BS[2] = m;
481  export BS;
482  kill @R3;
483  printlevel = old_printlevel;
484}
485example;
486{
487  "EXAMPLE:"; echo = 2;
488  ring r = 0,(x,y,Dx,Dy),dp;
489  def R = Weyl();    setring R; // Weyl algebra in variables x,y,Dx,Dy
490  poly F = x2-y3;
491  ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F;
492  // I is not holonomic, since its dimension is not 4/2=2
493  gkdim(I);
494  DLoc(I, x2-y3); // exports LD0 and BS
495  LD0; // localized module (R/I)_f is isomorphic to R/LD0
496  BS; // description of b-function for localization
497}
498
499proc DLoc0(ideal I, poly F)
500"USAGE:  DLoc0(I, F);  I an ideal, F a poly
501RETURN:  ring
502PURPOSE: compute the presentation of the localization of D/I w.r.t. f^s,
503@*           where D is a Weyl Algebra, based on the output of procedure SDLoc
504ASSUME: the basering is similar to the output ring of SDLoc procedure
505NOTE:    activate this ring with the @code{setring} command. In this ring,
506@* the ideal LD0 (in Groebner basis) is the presentation of the localization
507@* the list BS contains roots and multiplicities of Bernstein polynomial of (D/I)_f.
508DISPLAY: If printlevel=1, progress debug messages will be printed,
509@*       if printlevel>=2, all the debug messages will be printed.
510EXAMPLE: example DLoc0; shows examples
511"
512{
513  if (dmodappassumeViolation())
514  {
515    ERROR("Basering is inappropriate: characteristic>0 or qring present");
516  }
517  /* assume: to be run in the output ring of SDLoc */
518  /* doing: add F, eliminate vars*Dvars, factorize BS */
519  /* analogue to annfs0 */
520  def @R2 = basering;
521  // we're in D_n[s], where the elim ord for s is set
522  ideal J = NF(I,std(F));
523  // make leadcoeffs positive
524  int i;
525  for (i=1; i<= ncols(J); i++)
526  {
527    if (leadcoef(J[i]) <0 )
528    {
529      J[i] = -J[i];
530    }
531  }
532  J = J,F;
533  ideal M = groebner(J);
534  int Nnew = nvars(@R2);
535  ideal K2 = nselect(M,1..Nnew-1);
536  int ppl = printlevel-voice+2;
537  dbprint(ppl,"// -1-1- _x,_Dx are eliminated in basering");
538  dbprint(ppl-1, K2);
539  // the ring @R3 and the search for minimal negative int s
540  ring @R3 = 0,s,dp;
541  dbprint(ppl,"// -2-1- the ring @R3 = K[s] is ready");
542  ideal K3 = imap(@R2,K2);
543  poly p = K3[1];
544  dbprint(ppl,"// -2-2- attempt the factorization");
545  list PP = factorize(p);          //with constants and multiplicities
546  ideal bs; intvec m;             //the Bernstein polynomial is monic, so we are not interested in constants
547  for (i=2; i<= size(PP[1]); i++)  //we delete P[1][1] and P[2][1]
548  {
549    bs[i-1] = PP[1][i];
550    m[i-1]  = PP[2][i];
551  }
552  ideal bbs; int srat=0; int HasRatRoots = 0;
553  int sP;
554  for (i=1; i<= size(bs); i++)
555  {
556    if (deg(bs[i]) == 1)
557    {
558      bbs = bbs,bs[i];
559    }
560  }
561  if (size(bbs)==0)
562  {
563    dbprint(ppl-1,"// -2-3- factorization: no rational roots");
564    //    HasRatRoots = 0;
565    HasRatRoots = 1; // s0 = -1 then
566    sP = -1;
567    // todo: return ideal with no subst and a b-function unfactorized
568  }
569  else
570  {
571    // exist rational roots
572    dbprint(ppl-1,"// -2-3- factorization: rational roots found");
573    HasRatRoots = 1;
574    //    dbprint(ppl-1,bbs);
575    bbs = bbs[2..ncols(bbs)];
576    ideal P = bbs;
577    dbprint(ppl-1,P);
578    srat = size(bs) - size(bbs);
579    // define minIntRoot on linear factors or find out that it doesn't exist
580    intvec vP;
581    number nP;
582    P = normalize(P); // now leadcoef = 1
583    P = lead(P)-P;
584    sP = size(P);
585    int cnt = 0;
586    for (i=1; i<=sP; i++)
587    {
588      nP = leadcoef(P[i]);
589      if ( (nP - int(nP)) == 0 )
590      {
591        cnt++;
592        vP[cnt] = int(nP);
593      }
594    }
595//     if ( size(vP)>=2 )
596//     {
597//       vP = vP[2..size(vP)];
598//     }
599    if ( size(vP)==0 )
600    {
601      // no roots!
602      dbprint(ppl,"// -2-4- no integer root, setting s0 = -1");
603      sP = -1;
604      //      HasRatRoots = 0; // older stuff, here we do substitution
605      HasRatRoots = 1;
606    }
607    else
608    {
609      HasRatRoots = 1;
610      sP = -Max(-vP);
611      dbprint(ppl,"// -2-4- minimal integer root found");
612      dbprint(ppl-1, sP);
613      //    int sP = minIntRoot(bbs,1);
614//       P =  normalize(P);
615//       bs = -subst(bs,s,0);
616      if (sP >=0)
617      {
618        dbprint(ppl,"// -2-5- nonnegative root, setting s0 = -1");
619        sP = -1;
620      }
621      else
622      {
623        dbprint(ppl,"// -2-5- the root is negative");
624      }
625    }
626  }
627
628  if (HasRatRoots)
629  {
630    setring @R2;
631    K2 = subst(I,s,sP);
632    // IF min int root exists ->
633    // create the ordinary Weyl algebra and put the result into it,
634    // thus creating the ring @R5
635    // ELSE : return the same ring with new objects
636    // keep: N, i,j,s, tmp, RL
637    Nnew = Nnew - 1; // former 2*N;
638    // list RL = ringlist(save);  // is defined earlier
639    //  kill Lord, tmp, iv;
640    list L = 0;
641    list Lord, tmp;
642    intvec iv;
643    list RL = ringlist(basering);
644    L[1] = RL[1];
645    L[4] = RL[4];  //char, minpoly
646    // check whether vars have admissible names -> done earlier
647    // list Name = RL[2]M
648    // DName is defined earlier
649    list NName; // = RL[2]; // skip the last var 's'
650    for (i=1; i<=Nnew; i++)
651    {
652      NName[i] =  RL[2][i];
653    }
654    L[2] = NName;
655    // dp ordering;
656    string s = "iv=";
657    for (i=1; i<=Nnew; i++)
658    {
659      s = s+"1,";
660    }
661    s[size(s)] = ";";
662    execute(s);
663    tmp     = 0;
664    tmp[1]  = "dp";  // string
665    tmp[2]  = iv;  // intvec
666    Lord[1] = tmp;
667    kill s;
668    tmp[1]  = "C";
669    iv = 0;
670    tmp[2]  = iv;
671    Lord[2] = tmp;
672    tmp     = 0;
673    L[3]    = Lord;
674    // we are done with the list
675    // Add: Plural part
676    def @R4@ = ring(L);
677    setring @R4@;
678    int N = Nnew/2;
679    matrix @D[Nnew][Nnew];
680    for (i=1; i<=N; i++)
681    {
682      @D[i,N+i]=1;
683    }
684    def @R4 = nc_algebra(1,@D);
685    setring @R4;
686    kill @R4@;
687    dbprint(ppl,"// -3-1- the ring @R4 is ready");
688    dbprint(ppl-1, @R4);
689    ideal K4 = imap(@R2,K2);
690    intvec vopt = option(get);
691    option(redSB);
692    dbprint(ppl,"// -3-2- the final cosmetic std");
693    K4 = groebner(K4);  // std does the job too
694    option(set,vopt);
695    // total cleanup
696    setring @R2;
697    ideal bs = imap(@R3,bs);
698    bs = -normalize(bs); // "-" for getting correct coeffs!
699    bs = subst(bs,s,0);
700    kill @R3;
701    setring @R4;
702    ideal bs = imap(@R2,bs); // only rationals are the entries
703    list BS; BS[1] = bs; BS[2] = m;
704    export BS;
705    //    list LBS = imap(@R3,LBS);
706    //    list BS; BS[1] = sbs; BS[2] = m;
707    //    BS;
708    //    export BS;
709    ideal LD0 = K4;
710    export LD0;
711    return(@R4);
712  }
713  else
714  {
715    /* SHOULD NEVER GET THERE */
716    /* no rational/integer roots */
717    /* return objects in the copy of current ring */
718    setring @R2;
719    ideal LD0 = I;
720    poly BS = normalize(K2[1]);
721    export LD0;
722    export BS;
723    return(@R2);
724  }
725}
726example;
727{
728  "EXAMPLE:"; echo = 2;
729  ring r = 0,(x,y,Dx,Dy),dp;
730  def R = Weyl();    setring R; // Weyl algebra in variables x,y,Dx,Dy
731  poly F = x2-y3;
732  ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F;
733  // moreover I is not holonomic, since its dimension is not 2 = 4/2
734  gkdim(I); // 3
735  def W = SDLoc(I,F);  setring W; // creates ideal LD in W = R[s]
736  def U = DLoc0(LD, x2-y3);  setring U; // compute in R
737  LD0; // Groebner basis of the presentation of localization
738  BS; // description of b-function for localization
739}
740
741
742proc SDLoc(ideal I, poly F)
743"USAGE:  SDLoc(I, F);  I an ideal, F a poly
744RETURN:  ring
745PURPOSE: compute a generic presentation of the localization of D/I w.r.t. f^s
746ASSUME: the basering D is a Weyl algebra
747NOTE:    activate this ring with the @code{setring} command. In this ring,
748@*  the ideal LD (in Groebner basis) is the presentation of the localization
749DISPLAY: If printlevel=1, progress debug messages will be printed,
750@*       if printlevel>=2, all the debug messages will be printed.
751EXAMPLE: example SDLoc; shows examples
752"
753{
754  /* analogue to Sannfs */
755  /* printlevel >=4 gives debug info */
756  /* assume: we're in the Weyl algebra D  in x1,x2,...,d1,d2,... */
757
758  if (dmodappassumeViolation())
759  {
760    ERROR("Basering is inappropriate: characteristic>0 or qring present");
761  }
762  if (!isWeyl())
763  {
764    ERROR("Basering is not a Weyl algebra");
765  }
766  def save = basering;
767  /* 1. create D <t, dt, s > as in LOT */
768  /* ordering: eliminate t,dt */
769  int ppl = printlevel-voice+2;
770  int N = nvars(save); N = N div 2;
771  int Nnew = 2*N + 3; // t,Dt,s
772  int i,j;
773  string s;
774  list RL = ringlist(save);
775  list L, Lord;
776  list tmp;
777  intvec iv;
778  L[1] = RL[1]; // char
779  L[4] = RL[4]; // char, minpoly
780  // check whether vars have admissible names
781  list Name  = RL[2];
782  list RName;
783  RName[1] = "@t";
784  RName[2] = "@Dt";
785  RName[3] = "@s";
786  for(i=1;i<=N;i++)
787  {
788    for(j=1; j<=size(RName);j++)
789    {
790      if (Name[i] == RName[j])
791      {
792        ERROR("Variable names should not include @t,@Dt,@s");
793      }
794    }
795  }
796  // now, create the names for new vars
797  tmp    =  0;
798  tmp[1] = "@t";
799  tmp[2] = "@Dt";
800  list SName ; SName[1] = "@s";
801  list NName = tmp + Name + SName;
802  L[2]   = NName;
803  tmp    = 0;
804  kill NName;
805  // block ord (a(1,1),dp);
806  tmp[1]  = "a"; // string
807  iv      = 1,1;
808  tmp[2]  = iv; //intvec
809  Lord[1] = tmp;
810  // continue with dp 1,1,1,1...
811  tmp[1]  = "dp"; // string
812  s       = "iv=";
813  for(i=1;i<=Nnew;i++)
814  {
815    s = s+"1,";
816  }
817  s[size(s)]= ";";
818  execute(s);
819  tmp[2]    = iv;
820  Lord[2]   = tmp;
821  tmp[1]    = "C";
822  iv        = 0;
823  tmp[2]    = iv;
824  Lord[3]   = tmp;
825  tmp       = 0;
826  L[3]      = Lord;
827  // we are done with the list
828  def @R@ = ring(L);
829  setring @R@;
830  matrix @D[Nnew][Nnew];
831  @D[1,2]=1;
832  for(i=1; i<=N; i++)
833  {
834    @D[2+i,N+2+i]=1;
835  }
836  // ADD [s,t]=-t, [s,Dt]=Dt
837  @D[1,Nnew] = -var(1);
838  @D[2,Nnew] = var(2);
839  def @R = nc_algebra(1,@D);
840  setring @R;
841  kill @R@;
842  dbprint(ppl,"// -1-1- the ring @R(@t,@Dt,_x,_Dx,@s) is ready");
843  dbprint(ppl-1, @R);
844  poly  F = imap(save,F);
845  ideal I = imap(save,I);
846  dbprint(ppl-1, "the ideal after map:");
847  dbprint(ppl-1, I);
848  poly p = 0;
849  for(i=1; i<=N; i++)
850  {
851    p = diff(F,var(2+i))*@Dt + var(2+N+i);
852    dbprint(ppl-1, p);
853    I = subst(I,var(2+N+i),p);
854    dbprint(ppl-1, var(2+N+i));
855    p = 0;
856  }
857  I = I, @t - F;
858  // t*Dt + s +1 reduced with t-f gives f*Dt + s
859  I = I, F*var(2) + var(Nnew); // @s
860  // -------- the ideal I is ready ----------
861  dbprint(ppl,"// -1-2- starting the elimination of @t,@Dt in @R");
862  dbprint(ppl-1, I);
863  //  ideal J = engine(I,eng);
864  ideal J = groebner(I);
865  dbprint(ppl-1,"// -1-2-1- result of the  elimination of @t,@Dt in @R");
866  dbprint(ppl-1, J);;
867  ideal K = nselect(J,1..2);
868  dbprint(ppl,"// -1-3- @t,@Dt are eliminated");
869  dbprint(ppl-1, K);  // K is without t, Dt
870  K = groebner(K);  // std does the job too
871  // now, we must change the ordering
872  // and create a ring without t, Dt
873  setring save;
874  // ----------- the ring @R3 ------------
875  // _x, _Dx,s;  elim.ord for _x,_Dx.
876  // keep: N, i,j,s, tmp, RL
877  Nnew = 2*N+1;
878  kill Lord, tmp, iv, RName;
879  list Lord, tmp;
880  intvec iv;
881  L[1] = RL[1];
882  L[4] = RL[4];  // char, minpoly
883  // check whether vars hava admissible names -> done earlier
884  // now, create the names for new var
885  tmp[1] = "s";
886  list NName = Name + tmp;
887  L[2] = NName;
888  tmp = 0;
889  // block ord (dp(N),dp);
890  // string s is already defined
891  s = "iv=";
892  for (i=1; i<=Nnew-1; i++)
893  {
894    s = s+"1,";
895  }
896  s[size(s)]=";";
897  execute(s);
898  tmp[1] = "dp";  // string
899  tmp[2] = iv;   // intvec
900  Lord[1] = tmp;
901  // continue with dp 1,1,1,1...
902  tmp[1] = "dp";  // string
903  s[size(s)] = ",";
904  s = s+"1;";
905  execute(s);
906  kill s;
907  kill NName;
908  tmp[2]      = iv;
909  Lord[2]     = tmp;
910  tmp[1]      = "C";  iv  = 0;  tmp[2]=iv;
911  Lord[3]     = tmp;  tmp = 0;
912  L[3]        = Lord;
913  // we are done with the list. Now add a Plural part
914  def @R2@ = ring(L);
915  setring @R2@;
916  matrix @D[Nnew][Nnew];
917  for (i=1; i<=N; i++)
918  {
919    @D[i,N+i]=1;
920  }
921  def @R2 = nc_algebra(1,@D);
922  setring @R2;
923  kill @R2@;
924  dbprint(ppl,"//  -2-1- the ring @R2(_x,_Dx,s) is ready");
925  dbprint(ppl-1, @R2);
926  ideal MM = maxideal(1);
927  MM = 0,s,MM;
928  map R01 = @R, MM;
929  ideal K = R01(K);
930  // total cleanup
931  ideal LD = K;
932  // make leadcoeffs positive
933  for (i=1; i<= ncols(LD); i++)
934  {
935    if (leadcoef(LD[i]) <0 )
936    {
937      LD[i] = -LD[i];
938    }
939  }
940  export LD;
941  kill @R;
942  return(@R2);
943}
944example;
945{
946  "EXAMPLE:"; echo = 2;
947  ring r = 0,(x,y,Dx,Dy),dp;
948  def R = Weyl(); // Weyl algebra on the variables x,y,Dx,Dy
949  setring R;
950  poly F = x2-y3;
951  ideal I = Dx*F, Dy*F;
952  // note, that I is not holonomic, since it's dimension is not 2
953  gkdim(I); // 3, while dim R = 4
954  def W = SDLoc(I,F);
955  setring W; // = R[s], where s is a new variable
956  LD; // Groebner basis of s-parametric presentation
957}
958
959proc annRat(poly g, poly f)
960"USAGE:  annRat(g,f);  f, g polynomials
961RETURN:  ring
962PURPOSE: compute the annihilator of the rational function g/f in the Weyl algebra D
963NOTE: activate the output ring with the @code{setring} command.
964@*      In the output ring, the ideal LD (in Groebner basis) is the annihilator.
965@*      The algorithm uses the computation of ann f^{-1} via D-modules.
966DISPLAY: If printlevel=1, progress debug messages will be printed,
967@*       if printlevel>=2, all the debug messages will be printed.
968SEE ALSO: annPoly
969EXAMPLE: example annRat; shows examples
970"
971{
972
973  if (dmodappassumeViolation())
974  {
975    ERROR("Basering is inappropriate: characteristic>0 or qring present");
976  }
977
978  // assumptions: f is not a constant
979  if (f==0) { ERROR("Denominator cannot be zero"); }
980  if (leadexp(f) == 0)
981  {
982    // f = const, so use annPoly
983    g = g/f;
984    def @R = annPoly(g);
985    return(@R);
986  }
987    // computes the annihilator of g/f
988  def save = basering;
989  int ppl = printlevel-voice+2;
990  dbprint(ppl,"// -1-[annRat] computing the ann f^s");
991  def  @R1 = SannfsBM(f);
992  setring @R1;
993  poly f = imap(save,f);
994  int i,mir;
995  int isr = 0; // checkRoot1(LD,f,1); // roots are negative, have to enter positive int
996  if (!isr)
997  {
998    // -1 is not the root
999    // find the m.i.r iteratively
1000    mir = 0;
1001    for(i=nvars(save)+1; i>=1; i--)
1002    {
1003      isr =  checkRoot1(LD,f,i);
1004      if (isr) { mir =-i; break; }
1005    }
1006    if (mir ==0)
1007    {
1008      "No integer root found! Aborting computations, inform the authors!";
1009      return(0);
1010    }
1011    // now mir == i is m.i.r.
1012  }
1013  else
1014  {
1015    // -1 is the m.i.r
1016    mir = -1;
1017  }
1018  dbprint(ppl,"// -2-[annRat] the minimal integer root is ");
1019  dbprint(ppl-1, mir);
1020  // use annfspecial
1021  dbprint(ppl,"// -3-[annRat] running annfspecial ");
1022  ideal AF = annfspecial(LD,f,mir,-1); // ann f^{-1}
1023  //  LD = subst(LD,s,j);
1024  //  LD = engine(LD,0);
1025  // modify the ring: throw s away
1026  // output ring comes from SannfsBM
1027  list U = ringlist(@R1);
1028  list tmp; // variables
1029  for(i=1; i<=size(U[2])-1; i++)
1030  {
1031    tmp[i] = U[2][i];
1032  }
1033  U[2] = tmp;
1034  tmp = 0;
1035  tmp[1] = U[3][1]; // x,Dx block
1036  tmp[2] = U[3][3]; // module block
1037  U[3] = tmp;
1038  tmp = 0;
1039  tmp = U[1],U[2],U[3],U[4];
1040  def @R2 = ring(tmp);
1041  setring @R2;
1042  // now supply with Weyl algebra relations
1043  int N = nvars(@R2)/2;
1044  matrix @D[2*N][2*N];
1045  for(i=1; i<=N; i++)
1046  {
1047    @D[i,N+i]=1;
1048  }
1049  def @R3 = nc_algebra(1,@D);
1050  setring @R3;
1051  dbprint(ppl,"// - -[annRat] ring without s is ready:");
1052  dbprint(ppl-1,@R3);
1053  poly g = imap(save,g);
1054  matrix G[1][1] = g;
1055  matrix LL = matrix(imap(@R1,AF));
1056  kill @R1;   kill @R2;
1057  dbprint(ppl,"// -4-[annRat] running modulo");
1058  ideal LD  = modulo(G,LL);
1059  dbprint(ppl,"// -4-[annRat] running GB on the final result");
1060  LD  = engine(LD,0);
1061  export LD;
1062  return(@R3);
1063}
1064example
1065{
1066  "EXAMPLE:"; echo = 2;
1067  ring r = 0,(x,y),dp;
1068  poly g = 2*x*y;  poly f = x^2 - y^3;
1069  def B = annRat(g,f);
1070  setring B;
1071  LD;
1072  // Now, compare with the output of Macaulay2:
1073  ideal tst = 3*x*Dx + 2*y*Dy + 1, y^3*Dy^2 - x^2*Dy^2 + 6*y^2*Dy + 6*y,
10749*y^2*Dx^2*Dy-4*y*Dy^3+27*y*Dx^2+2*Dy^2, 9*y^3*Dx^2-4*y^2*Dy^2+10*y*Dy -10;
1075 option(redSB); option(redTail);
1076 LD = groebner(LD);
1077 tst = groebner(tst);
1078 print(matrix(NF(LD,tst)));  print(matrix(NF(tst,LD)));
1079 // So, these two answers are the same
1080}
1081
1082/*
1083//static proc ex_annRat()
1084{
1085  // more complicated example for annRat
1086  ring r = 0,(x,y,z),dp;
1087  poly f = x3+y3+z3; // mir = -2
1088  poly g = x*y*z;
1089  def A = annRat(g,f);
1090  setring A;
1091}
1092*/
1093
1094proc annPoly(poly f)
1095"USAGE:  annPoly(f);  f a poly
1096RETURN:  ring
1097PURPOSE: compute the complete annihilator ideal of f in the Weyl algebra D
1098NOTE:  activate the output ring with the @code{setring} command.
1099@*   In the output ring, the ideal LD (in Groebner basis) is the annihilator.
1100DISPLAY: If printlevel=1, progress debug messages will be printed,
1101@*       if printlevel>=2, all the debug messages will be printed.
1102SEE ALSO: annRat
1103EXAMPLE: example annPoly; shows examples
1104"
1105{
1106  // computes a system of linear PDEs with polynomial coeffs for f
1107  def save = basering;
1108  list L = ringlist(save);
1109  list Name = L[2];
1110  int N = nvars(save);
1111  int i;
1112  for (i=1; i<=N; i++)
1113  {
1114    Name[N+i] = "D"+Name[i]; // concat
1115  }
1116  L[2] = Name;
1117  def @R = ring(L);
1118  setring @R;
1119  def @@R = Weyl();
1120  setring @@R;
1121  kill @R;
1122  matrix M[1][N];
1123  for (i=1; i<=N; i++)
1124  {
1125    M[1,i] = var(N+i);
1126  }
1127  matrix F[1][1] = imap(save,f);
1128  ideal I = modulo(F,M);
1129  ideal LD = groebner(I);
1130  export LD;
1131  return(@@R);
1132}
1133example
1134{
1135  "EXAMPLE:"; echo = 2;
1136  ring r = 0,(x,y,z),dp;
1137  poly f = x^2*z - y^3;
1138  def A = annPoly(f);
1139  setring A; // A is the 3rd Weyl algebra in 6 variables
1140  LD; // the Groebner basis of annihilator
1141  gkdim(LD); // must be 3 = 6/2, since A/LD is holonomic module
1142  NF(Dy^4, LD); // must be 0 since Dy^4 clearly annihilates f
1143}
1144
1145/* DIFFERENT EXAMPLES
1146
1147//static proc exCusp()
1148{
1149  "EXAMPLE:"; echo = 2;
1150  ring r = 0,(x,y,Dx,Dy),dp;
1151  def R = Weyl();   setring R;
1152  poly F = x2-y3;
1153  ideal I = (y^3 - x^2)*Dx - 2*x, (y^3 - x^2)*Dy + 3*y^2; // I = Dx*F, Dy*F;
1154  def W = SDLoc(I,F);
1155  setring W;
1156  LD;
1157  def U = DLoc0(LD,x2-y3);
1158  setring U;
1159  LD0;
1160  BS;
1161  // the same with DLoc:
1162  setring R;
1163  DLoc(I,F);
1164}
1165
1166//static proc exWalther1()
1167{
1168  // p.18 Rem 3.10
1169  ring r = 0,(x,Dx),dp;
1170  def R = nc_algebra(1,1);
1171  setring R;
1172  poly F = x;
1173  ideal I = x*Dx+1;
1174  def W = SDLoc(I,F);
1175  setring W;
1176  LD;
1177  ideal J = LD, x;
1178  eliminate(J,x*Dx); // must be [1]=s // agree!
1179  // the same result with Dloc0:
1180  def U = DLoc0(LD,x);
1181  setring U;
1182  LD0;
1183  BS;
1184}
1185
1186//static proc exWalther2()
1187{
1188  // p.19 Rem 3.10 cont'd
1189  ring r = 0,(x,Dx),dp;
1190  def R = nc_algebra(1,1);
1191  setring R;
1192  poly F = x;
1193  ideal I = (x*Dx)^2+1;
1194  def W = SDLoc(I,F);
1195  setring W;
1196  LD;
1197  ideal J = LD, x;
1198  eliminate(J,x*Dx); // must be [1]=s^2+2*s+2 // agree!
1199  // the same result with Dloc0:
1200  def U = DLoc0(LD,x);
1201  setring U;
1202  LD0;
1203  BS;
1204  // almost the same with DLoc
1205  setring R;
1206  DLoc(I,F);
1207  LD0;  BS;
1208}
1209
1210//static proc exWalther3()
1211{
1212  // can check with annFs too :-)
1213  // p.21 Ex 3.15
1214  LIB "nctools.lib";
1215  ring r = 0,(x,y,z,w,Dx,Dy,Dz,Dw),dp;
1216  def R = Weyl();
1217  setring R;
1218  poly F = x2+y2+z2+w2;
1219  ideal I = Dx,Dy,Dz,Dw;
1220  def W = SDLoc(I,F);
1221  setring W;
1222  LD;
1223  ideal J = LD, x2+y2+z2+w2;
1224  eliminate(J,x*y*z*w*Dx*Dy*Dz*Dw); // must be [1]=s^2+3*s+2 // agree
1225  ring r2 =  0,(x,y,z,w),dp;
1226  poly F = x2+y2+z2+w2;
1227  def Z = annfs(F);
1228  setring Z;
1229  LD;
1230  BS;
1231  // the same result with Dloc0:
1232  setring W;
1233  def U = DLoc0(LD,x2+y2+z2+w2);
1234  setring U;
1235  LD0;  BS;
1236  // the same result with DLoc:
1237  setring R;
1238  DLoc(I,F);
1239  LD0;  BS;
1240}
1241
1242*/
1243
1244proc engine(def I, int i)
1245"USAGE:  engine(I,i);  I  ideal/module/matrix, i an int
1246RETURN:  the same type as I
1247PURPOSE: compute the Groebner basis of I with the algorithm, chosen via i
1248NOTE: By default and if i=0, slimgb is used; otherwise std does the job.
1249EXAMPLE: example engine; shows examples
1250"
1251{
1252  /* std - slimgb mix */
1253  def J;
1254  //  ideal J;
1255  if (i==0)
1256  {
1257    J = slimgb(I);
1258  }
1259  else
1260  {
1261    // without options -> strange! (ringlist?)
1262    intvec v = option(get);
1263    option(redSB);
1264    option(redTail);
1265    J = std(I);
1266    option(set, v);
1267  }
1268  return(J);
1269}
1270example
1271{
1272  "EXAMPLE:"; echo = 2;
1273  ring r = 0,(x,y),Dp;
1274  ideal I  = y*(x3-y2),x*(x3-y2);
1275  engine(I,0); // uses slimgb
1276  engine(I,1); // uses std
1277}
1278
1279proc insertGenerator (list #)
1280"USAGE:  insertGenerator(id,p[,k]);  id an ideal/module, p a poly/vector, k an optional int
1281RETURN:  same as id
1282PURPOSE: inserts p into the first argument at k-th index position and returns the enlarged object
1283NOTE:    If k is given, p is inserted at position k, otherwise (and by default),
1284@*       p is inserted at the beginning.
1285EXAMPLE: example insertGenerator; shows examples
1286"
1287{
1288  if (size(#) < 2)
1289  {
1290    ERROR("insertGenerator has to be called with at least 2 arguments (ideal/module,poly/vector)");
1291  }
1292  string inp1 = typeof(#[1]);
1293  if (inp1 == "ideal" || inp1 == "module")
1294  {
1295    if (inp1 == "ideal") { ideal id = #[1]; }
1296    else { module id = #[1]; }
1297  }
1298  else { ERROR("first argument has to be of type ideal or module"); }
1299  string inp2 = typeof(#[2]);
1300  if (inp2 == "poly" || inp2 == "vector")
1301  {
1302    if (inp2 =="poly") { poly f = #[2]; }
1303    else
1304    {
1305      if (inp1 == "ideal")
1306      {
1307        ERROR("second argument has to be a polynomial if first argument is an ideal");
1308      }
1309      else { vector f = #[2]; }
1310    }
1311  }
1312  else { ERROR("second argument has to be of type poly/vector"); }
1313  int n = ncols(id);
1314  int k = 1; // default
1315  if (size(#)>=3)
1316  {
1317    if (typeof(#[3]) == "int")
1318    {
1319      k = #[3];
1320      if (k<=0)
1321      {
1322        ERROR("third argument has to be positive");
1323      }
1324    }
1325    else { ERROR("third argument has to be of type int"); }
1326  }
1327  execute(inp1 +" J;");
1328  if (k == 1) { J = f,id; }
1329  else
1330  {
1331    if (k>n)
1332    {
1333      J = id;
1334      J[k] = f;
1335    }
1336    else // 1<k<=n
1337    {
1338      J[1..k-1] = id[1..k-1];
1339      J[k] = f;
1340      J[k+1..n+1] = id[k..n];
1341    }
1342  }
1343  return(J);
1344}
1345example
1346{
1347  "EXAMPLE:"; echo = 2;
1348  ring r = 0,(x,y,z),dp;
1349  ideal I = x^2,z^4;
1350  insertGenerator(I,y^3);
1351  insertGenerator(I,y^3,2);
1352  module M = I;
1353  insertGenerator(M,[x^3,y^2,z],2);
1354}
1355
1356proc deleteGenerator (list #)
1357"USAGE:  deleteGenerator(id,k);  id an ideal/module, k an int
1358RETURN:  same as id
1359PURPOSE: deletes the k-th generator from the first argument and returns the altered object
1360EXAMPLE: example insertGenerator; shows examples
1361"
1362{
1363  if (size(#) < 2)
1364  {
1365    ERROR("deleteGenerator has to be called with 2 arguments (ideal/module,int)");
1366  }
1367  string inp1 = typeof(#[1]);
1368  if (inp1 == "ideal" || inp1 == "module")
1369  {
1370    if (inp1 == "ideal") { ideal id = #[1]; }
1371    else { module id = #[1]; }
1372  }
1373  else { ERROR("first argument has to be of type ideal or module"); }
1374  string inp2 = typeof(#[2]);
1375  if (inp2 == "int" || inp2 == "number") { int k = int(#[2]); }
1376  else { ERROR("second argument has to be of type int"); }
1377  int n = ncols(id);
1378  if (n == 1) { ERROR(inp1+" must have more than one generator"); }
1379  if (k<=0 || k>n) { ERROR("second argument has to be in the range 1,...,"+string(n)); }
1380  execute(inp1 +" J;");
1381  if (k == 1) { J = id[2..n]; }
1382  else
1383  {
1384    if (k == n) { J = id[1..n-1]; }
1385    else
1386    {
1387      J[1..k-1] = id[1..k-1];
1388      J[k..n-1] = id[k+1..n];
1389    }
1390  }
1391  return(J);
1392}
1393example
1394{
1395  "EXAMPLE:"; echo = 2;
1396  ring r = 0,(x,y,z),dp;
1397  ideal I = x^2,y^3,z^4;
1398  deleteGenerator(I,2);
1399  module M = [x,y,z],[x2,y2,z2],[x3,y3,z3];;
1400  deleteGenerator(M,2);
1401}
1402
1403proc fl2poly(list L, string s)
1404"USAGE:  fl2poly(L,s); L a list, s a string
1405RETURN:  poly
1406PURPOSE: reconstruct a monic polynomial in one variable from its factorization
1407ASSUME:  s is a string with the name of some variable and
1408@*         L is supposed to consist of two entries:
1409@*         L[1] of the type ideal with the roots of a polynomial
1410@*         L[2] of the type intvec with the multiplicities of corr. roots
1411EXAMPLE: example fl2poly; shows examples
1412"
1413{
1414  if (varNum(s)==0)
1415  {
1416    ERROR("no such variable found in the basering"); return(0);
1417  }
1418  poly x = var(varNum(s));
1419  poly P = 1;
1420  int sl = size(L[1]);
1421  ideal RR = L[1];
1422  intvec IV = L[2];
1423  for(int i=1; i<= sl; i++)
1424  {
1425    if (IV[i] > 0)
1426    {
1427      P = P*((x-RR[i])^IV[i]);
1428    }
1429    else
1430    {
1431      printf("Ignored the root with incorrect multiplicity %s",string(IV[i]));
1432    }
1433  }
1434  return(P);
1435}
1436example
1437{
1438  "EXAMPLE:"; echo = 2;
1439  ring r = 0,(x,y,z,s),Dp;
1440  ideal I = -1,-4/3,-5/3,-2;
1441  intvec mI = 2,1,1,1;
1442  list BS = I,mI;
1443  poly p = fl2poly(BS,"s");
1444  p;
1445  factorize(p,2);
1446}
1447
1448static proc safeVarName (string s, string cv)
1449{
1450  string S;
1451  if (cv == "v")  { S = "," + "," + varstr(basering)  + ","; }
1452  if (cv == "c")  { S = "," + "," + charstr(basering) + ","; }
1453  if (cv == "cv") { S = "," + charstr(basering) + "," + varstr(basering) + ","; }
1454  s = "," + s + ",";
1455  while (find(S,s) <> 0)
1456  {
1457    s[1] = "@";
1458    s = "," + s;
1459  }
1460  s = s[2..size(s)-1];
1461  return(s)
1462}
1463
1464proc initialIdealW (ideal I, intvec u, intvec v, list #)
1465"USAGE:  initialIdealW(I,u,v [,s,t,w]);  I ideal, u,v intvecs, s,t optional ints,
1466@*       w an optional intvec
1467RETURN:  ideal, GB of initial ideal of the input ideal w.r.t. the weights u and v
1468ASSUME:  The basering is the n-th Weyl algebra in characteristic 0 and  for all
1469@*       1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the
1470@*       sequence of variables is given by x(1),...,x(n),D(1),...,D(n),
1471@*       where D(i) is the differential operator belonging to x(i).
1472PURPOSE: computes the initial ideal with respect to given weights.
1473NOTE:    u and v are understood as weight vectors for x(1..n) and D(1..n)
1474@*       respectively.
1475@*       If s<>0, @code{std} is used for Groebner basis computations,
1476@*       otherwise, and by default, @code{slimgb} is used.
1477@*       If t<>0, a matrix ordering is used for Groebner basis computations,
1478@*       otherwise, and by default, a block ordering is used.
1479@*       If w consist of 2n strictly positive entries, w is used for weighted
1480@*       homogenization, otherwise, and by default, no weights are used.
1481DISPLAY: If printlevel=1, progress debug messages will be printed,
1482@*       if printlevel>=2, all the debug messages will be printed.
1483EXAMPLE: example initialIdealW; shows examples
1484"
1485{
1486  if (dmodappassumeViolation())
1487  {
1488    ERROR("Basering is inappropriate: characteristic>0 or qring present");
1489  }
1490  if (!isWeyl())
1491  {
1492    ERROR("Basering is not a Weyl algebra.");
1493  }
1494  int ppl = printlevel - voice +2;
1495  def save = basering;
1496  int n = nvars(save)/2;
1497  int N = 2*n+1;
1498  list RL = ringlist(save);
1499  int i;
1500  int whichengine = 0;           // default
1501  int methodord   = 0;           // default
1502  intvec homogweights = 1:(2*n); // default
1503  if (size(#)>0)
1504  {
1505    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1506    {
1507      whichengine = int(#[1]);
1508    }
1509    if (size(#)>1)
1510    {
1511      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1512      {
1513        methodord = int(#[2]);
1514      }
1515      if (size(#)>2)
1516      {
1517        if (typeof(#[3])=="intvec")
1518        {
1519          if (size(#[3])==2*n && allPositive(#[3])==1)
1520          {
1521            homogweights = #[3];
1522          }
1523          else
1524          {
1525            print("// Homogenization weight vector must consist of positive entries and be");
1526            print("// of size " + string(n) + ". Using no weights.");
1527          }
1528        }
1529      }
1530    }
1531  }
1532  for (i=1; i<=n; i++)
1533  {
1534    if (bracket(var(i+n),var(i))<>1)
1535    {
1536      ERROR(string(var(i+n)) + " is not a differential operator for " + string(var(i)));
1537    }
1538  } 
1539  // 1. create  homogenized Weyl algebra
1540  // 1.1 create ordering
1541  intvec uv = u,v,0;
1542  homogweights = homogweights,1;
1543  list Lord = list(list("a",homogweights));
1544  list C0 = list("C",intvec(0));
1545  if (methodord == 0) // default: blockordering
1546  {
1547    Lord[2] = list("a",uv);
1548    Lord[3] = list("dp",intvec(1:(N-1)));
1549    Lord[4] = list("lp",intvec(1));
1550    Lord[5] = C0;
1551  }
1552  else // M() ordering
1553  {
1554    intmat @Ord[N][N];
1555    @Ord[1,1..N] = uv; @Ord[2,1..N] = 1:(N-1);
1556    for (i=1; i<=N-2; i++) { @Ord[2+i,N - i] = -1; }
1557    dbprint(ppl-1,"// the ordering matrix:",@Ord);
1558    Lord[2] = list("M",intvec(@Ord));
1559    Lord[3] = C0;
1560  }
1561  // 1.2 the new var
1562  list Lvar = RL[2]; Lvar[N] = safeVarName("h","cv");
1563  // 1.3 create commutative ring
1564  list L@@Dh = RL; L@@Dh = L@@Dh[1..4];
1565  L@@Dh[2] = Lvar; L@@Dh[3] = Lord;
1566  def @Dh = ring(L@@Dh); kill L@@Dh;
1567  setring @Dh;
1568  dbprint(ppl-1,"// the ring @Dh:",@Dh);
1569  // 1.4 create non-commutative relations
1570  matrix @relD[N][N];
1571  for (i=1; i<=n; i++) { @relD[i,n+i] = var(N)^(homogweights[i]+homogweights[n+i]); }
1572  def Dh = nc_algebra(1,@relD);
1573  setring Dh; kill @Dh;   
1574  dbprint(ppl-1,"// computing in ring",Dh);
1575  // 2. Compute the initial ideal
1576  ideal I = imap(save,I);
1577  I = homog(I,h);
1578  // 2.1 the hard part: Groebner basis computation
1579  dbprint(ppl, "// starting Groebner basis computation with engine: "+string(whichengine));
1580  I = engine(I, whichengine);
1581  dbprint(ppl, "// finished Groebner basis computation");
1582  dbprint(ppl-1, "// I before dehomogenization is " +string(I));
1583  I = subst(I,var(N),1); // dehomogenization
1584  dbprint(ppl-1, "I after dehomogenization is " +string(I));
1585  // 2.2 the initial form
1586  I = inForm(I,uv);
1587  setring save;
1588  I = imap(Dh,I); kill Dh;
1589  // 2.3 the final GB
1590  dbprint(ppl, "// starting cosmetic Groebner basis computation with engine: "+string(whichengine));
1591  I = engine(I, whichengine);
1592  dbprint(ppl,"// finished cosmetic Groebner basis computation");
1593  return(I);
1594}
1595example
1596{
1597  "EXAMPLE:"; echo = 2;
1598  ring @D = 0,(x,Dx),dp;
1599  def D = Weyl();
1600  setring D;
1601  intvec u = -1; intvec v = 2;
1602  ideal I = x^2*Dx^2,x*Dx^4;
1603  ideal J = initialIdealW(I,u,v); J;
1604}
1605
1606proc initialMalgrange (poly f,list #)
1607"USAGE:  initialMalgrange(f,[,a,b,v]); f poly, a,b optional ints, v opt. intvec
1608RETURN:  ring, Weyl algebra induced by basering, extended by two new vars t,Dt
1609PURPOSE: computes the initial Malgrange ideal of a given polynomial w.r.t. the weight
1610@*       vector (-1,0...,0,1,0,...,0) such that the weight of t is -1 and the
1611@*       weight of Dt is 1.
1612ASSUME:  The basering is commutative and of characteristic 0.
1613NOTE:    Activate the output ring with the @code{setring} command.
1614@*       The returned ring contains the ideal \"inF\", being the initial ideal
1615@*       of the Malgrange ideal of f.
1616@*       Varnames of the basering should not include t and Dt.
1617@*       If a<>0, @code{std} is used for Groebner basis computations,
1618@*       otherwise, and by default, @code{slimgb} is used.
1619@*       If b<>0, a matrix ordering is used for Groebner basis computations,
1620@*       otherwise, and by default, a block ordering is used.
1621@*       If a positive weight vector v is given, the weight
1622@*       (d,v[1],...,v[n],1,d+1-v[1],...,d+1-v[n]) is used for homogenization
1623@*       computations, where d denotes the weighted degree of f.
1624@*       Otherwise and by default, v is set to (1,...,1). See Noro, 2002.
1625DISPLAY: If printlevel=1, progress debug messages will be printed,
1626@*       if printlevel>=2, all the debug messages will be printed.
1627EXAMPLE: example initialMalgrange; shows examples
1628"
1629{
1630  if (dmodappassumeViolation())
1631  {
1632    ERROR("Basering is inappropriate: characteristic>0 or qring present");
1633  }
1634  if (!isCommutative())
1635  {
1636    ERROR("Basering must be commutative.");
1637  }
1638  int ppl = printlevel - voice +2;
1639  def save = basering;
1640  int n = nvars(save);
1641  int i;
1642  int whichengine = 0; // default
1643  int methodord   = 0; // default
1644  intvec u0 = 1:n;     // default
1645  if (size(#)>0)
1646  {
1647    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1648    {
1649      whichengine = int(#[1]);
1650    }
1651    if (size(#)>1)
1652    {
1653      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1654      {
1655        methodord = int(#[2]);
1656      }
1657      if (size(#)>2)
1658      {
1659        if (typeof(#[3])=="intvec" && size(#[3])==n && allPositive(#[3])==1)
1660        {
1661          u0 = #[3];
1662        }
1663      }
1664    }
1665  }
1666  // creating the homogenized extended Weyl algebra
1667  list RL = ringlist(save);
1668  RL = RL[1..4]; // if basering is commutative nc_algebra
1669  list C0 = list("C",intvec(0));
1670  // 1. create homogenization weights
1671  // 1.1. get the weighted degree of f
1672  list Lord = list(list("wp",u0),C0);
1673  list L@r = RL;
1674  L@r[3] = Lord;
1675  def r = ring(L@r); kill L@r,Lord;
1676  setring r;
1677  poly f = imap(save,f);
1678  int d = deg(f);
1679  setring save; kill r;
1680  // 1.2 the homogenization weights
1681  intvec homogweights = d;
1682  homogweights[n+2] = 1;
1683  for (i=1; i<=n; i++)
1684  {
1685    homogweights[i+1]   = u0[i];
1686    homogweights[n+2+i] = d+1-u0[i];
1687  }
1688  // 2. create extended Weyl algebra
1689  int N = 2*n+2;
1690  // 2.1 create names for vars
1691  string vart  = safeVarName("t","cv");
1692  string varDt = safeVarName("D"+vart,"cv");
1693  while (varDt <> "D"+vart)
1694  {
1695    vart  = safeVarName("@"+vart,"cv");
1696    varDt = safeVarName("D"+vart,"cv");
1697  }
1698  list Lvar;
1699  Lvar[1] = vart; Lvar[n+2] = varDt;
1700  for (i=1; i<=n; i++)
1701  {
1702    Lvar[i+1]   = string(var(i));
1703    Lvar[i+n+2] = safeVarName("D" + string(var(i)),"cv");
1704  }
1705  //  2.2 create ordering
1706  list Lord = list(list("dp",intvec(1:N)),C0);
1707  // 2.3 create the (n+1)-th Weyl algebra
1708  list L@D = RL; L@D[2] = Lvar; L@D[3] = Lord;
1709  def @D = ring(L@D); kill L@D;
1710  setring @D;
1711  def D = Weyl();
1712  setring D; kill @D;
1713  dbprint(ppl,"// the (n+1)-th Weyl algebra :" ,D);
1714  // 3. compute the initial ideal
1715  // 3.1 create the Malgrange ideal
1716  poly f = imap(save,f);
1717  ideal I = var(1)-f;
1718  for (i=1; i<=n; i++)
1719  {
1720    I = I, var(n+2+i)+diff(f,var(i+1))*var(n+2);
1721  }
1722  // I = engine(I,whichengine); // todo is it efficient to compute a GB first wrt dp first?
1723  // 3.2 computie the initial ideal
1724  intvec w = 1,0:n;
1725  I = initialIdealW(I,-w,w,whichengine,methodord,homogweights);
1726  ideal inF = I; attrib(inF,"isSB",1);
1727  export(inF);
1728  setring save;
1729  return(D);
1730}
1731example
1732{
1733  "EXAMPLE:"; echo = 2;
1734  ring r = 0,(x,y),dp;
1735  poly f = x^2+y^3+x*y^2;
1736  def D = initialMalgrange(f);
1737  setring D;
1738  inF;
1739  setring r;
1740  intvec v = 3,2;
1741  def D2 = initialMalgrange(f,1,1,v);
1742  setring D2;
1743  inF;
1744}
1745
1746proc dmodappassumeViolation()
1747{
1748  // returns Boolean : yes/no [for assume violation]
1749  // char K = 0
1750  // no qring
1751  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
1752  {
1753    //    "ERROR: no qring is allowed";
1754    return(1);
1755  }
1756  return(0);
1757}
1758
1759proc bFactor (poly F)
1760"USAGE:  bFactor(f);  f poly
1761RETURN:  list
1762PURPOSE: tries to compute the roots of a univariate poly f
1763NOTE:    The output list consists of two or three entries:
1764@*       roots of f as an ideal, their multiplicities as intvec, and,
1765@*       if present, a third one being the product of all irreducible factors
1766@*       of degree greater than one, given as string.
1767DISPLAY: If printlevel=1, progress debug messages will be printed,
1768@*       if printlevel>=2, all the debug messages will be printed.
1769EXAMPLE: example bFactor; shows examples
1770"
1771{
1772  int ppl = printlevel - voice +2;
1773  def save = basering;
1774  ideal LI = variables(F);
1775  list L;
1776  int i = size(LI);
1777  if (i>1) { ERROR("poly has to be univariate")}
1778  if (i == 0)
1779  {
1780    dbprint(ppl,"// poly is constant");
1781    L = list(ideal(0),intvec(0),string(F));
1782    return(L);
1783  }
1784  poly v = LI[1];
1785  L = ringlist(save); L = L[1..4];
1786  L[2] = list(string(v));
1787  L[3] = list(list("dp",intvec(1)),list("C",intvec(0)));
1788  def @S = ring(L);
1789  setring @S;
1790  poly ir = 1; poly F = imap(save,F);
1791  list L = factorize(F);
1792  ideal I = L[1];
1793  intvec m = L[2];
1794  ideal II; intvec mm;
1795  for (i=2; i<=ncols(I); i++)
1796  {
1797    if (deg(I[i]) > 1)
1798    {
1799      ir = ir * I[i]^m[i];
1800    }
1801    else
1802    {
1803      II[size(II)+1] = I[i];
1804      mm[size(II)] = m[i];
1805    }
1806  }
1807  II = normalize(II);
1808  II = subst(II,var(1),0);
1809  II = -II;
1810  if (size(II)>0)
1811  {
1812    dbprint(ppl,"// found roots");
1813    dbprint(ppl-1,string(II));
1814  }
1815  else
1816  {
1817    dbprint(ppl,"// found no roots");
1818  } 
1819  L = list(II,mm);
1820  if (ir <> 1)
1821  {
1822    dbprint(ppl,"// found irreducible factors");
1823    ir = cleardenom(ir);
1824    dbprint(ppl-1,string(ir));
1825    L = L + list(string(ir));
1826  }
1827  else
1828  {
1829    dbprint(ppl,"// no irreducible factors found");
1830  } 
1831  setring save;
1832  L = imap(@S,L);
1833  return(L);
1834}
1835example
1836{
1837  "EXAMPLE:"; echo = 2;
1838  ring r = 0,(x,y),dp;
1839  bFactor((x^2-1)^2);
1840  bFactor((x^2+1)^2); 
1841  bFactor((y^2+1/2)*(y+9)*(y-7));
1842}
Note: See TracBrowser for help on using the repository browser.