# source:git/Singular/LIB/dmodapp.lib@8c66fb

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