# source:git/Singular/LIB/dmodapp.lib@238c959

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