source:git/dmodloc.lib@3dcd82

spielwiese
Last change on this file since 3dcd82 was 3dcd82, checked in by Daniel Andres <daniel.andres@…>, 10 years ago
updated docu + bugfixes
• Property mode set to 100644
File size: 65.9 KB
Line
1//
2version="\$Id\$";
3category="Noncommutative";
4info="
5LIBRARY: dmodloc.lib     Localization of algebraic D-modules and applications
6AUTHORS: Daniel Andres,  daniel.andres@math.rwth-aachen.de
7
8Support: DFG Graduiertenkolleg 1632 'Experimentelle und konstruktive Algebra'
9
10
11OVERVIEW:
12Let I be a left ideal I in the n-th polynomial Weyl algebra D=K[x]<d> and
13let f be a polynomial in K[x].
14
15If D/I is a holonomic module over D, it is known that the localization of D/I
16at f is also holonomic. The procedure Dlocalization computes an ideal J in D
17such that this localization is isomorphic to D/J as D-modules.
18
19If one regards I as an ideal in the rational Weyl algebra as above, K(x)<d>*I,
20and intersects with K[x]<d>, the result is called the Weyl closure of I.
21The procedures WeylClosure (if I has finite holonomic rank) and WeylClosure1
22(if I is in the first Weyl algebra) can be used for computations.
23
24As an application of the Weyl closure, the procedure annRatSyz computes a
25holonomic part of the annihilator of a rational function by computing certain
26syzygies. The full annihilator can be obtained by taking the Weyl closure of
27the result.
28
29If one regards the left ideal I as system of linear PDEs, one can find its
30polynomial solutions with polSol (if I is holonomic) or polSolFiniteRank
31(if I is of finite holonomic rank). Rational solutions can be obtained
32with ratSol.
33
34The procedure bfctBound computes a possible multiple of the b-function for
35f^s*u at a generic root of f. Here, u stands for [1] in D/I.
36
37This library also offers the procedures holonomicRank and DsingularLocus to
38compute the holonomic rank and the singular locus of the D-module D/I.
39
40
41References:
42@* (OT)  Oaku, Takayama: 'Algorithms for D-modules',
43         Journal of Pure and Applied Algebra, 1998
44@* (OTT) Oaku, Takayama, Tsai: 'Polynomial and rational solutions of holonomic
45         systems', Journal of Pure and Applied Algebra, 2001
46@* (OTW) Oaku, Takayama, Walther: 'A Localization Algorithm for D-modules',
47         Journal of Symbolic Computation, 2000
48@* (Tsa) Tsai: 'Algorithms for algebraic analysis', PhD thesis, 2000
49
50
51PROCEDURES:
52
53holonomicRank(I);  computes the holonomic rank of I
54DsingularLocus(I); computes the singular locus of a D-module
55Dlocalization(I,f[,k,e]);  computes the localization of a D-module
56WeylClosure1(L);   computes the Weyl closure of operator in first Weyl algebra
57WeylClosure(I);    computes the Weyl closure of an ideal in the Weyl algebra
58polSol(I[,w,m]);   computes basis of polynomial solutions to the given system
59polSolFiniteRank(I[,w]); computes basis of polynomial solutions to given system
60ratSol(I);         computes basis of rational solutions to the given system
61bfctBound(I,f[,primdec]); computes multiple of b-function for f^s*u
62annRatSyz(f,g[,db,eng]);  computes part of annihilator of rational function g/f
63
64dmodGeneralAssumptionCheck();   check general assumptions
65safeVarName(s);    finds an untaken name to use for a new variable
66extendWeyl(S);     extends basering (Weyl algebra) by given vars
67polyVars(f,v);     checks whether f contains only variables indexed by v
68monomialInIdeal(I);    computes all monomials appearing in generators of ideal
69vars2pars(v);      converts variables specified by v into parameters
70killTerms(p,k);    delete all terms of p containing var(k)
71minIntRoot2(L);    finds minimal integer root in a list of roots
72maxIntRoot(L);     finds maximal integer root in a list of roots
73dmodAction(id,f[,v]);  computes the natural action of a D-module on K[x]
74dmodActionRat(id,w);   computes the natural action of a D-module on K(x)
75simplifyRat(v);    simplifies rational function
77multRat(v,w);      multiplies rational functions
78diffRat(v,j);      derives rational function
79commRing();        deletes non-commutative relations from ring
80rightNFWeyl(id,k); computes right NF wrt. right ideal (var(k)) in Weyl algebra
81
82
83KEYWORDS:
84D-module;
85holonomic rank; singular locus of D-module;
86D-localization; localization of D-module; characteristic variety;
87Weyl closure;
88polynomial solutions; rational solutions;
89annihilator of rational function;
90
92";
93
94
95/*
96CHANGELOG:
9712.11.12: bugfixes, updated docu
98*/
99
100
101LIB "bfun.lib";    // for pIntersect etc
102LIB "dmodapp.lib"; // for GBWeight, charVariety etc
103LIB "nctools.lib"; // for Weyl, isWeyl etc
104// TODO uncomment this once chern.lib is ready
105// LIB "chern.lib";   // for orderedPartition
106
107
108// testing for consistency of the library /////////////////////////////////////
109proc testdmodloc()
110{
111  example dmodGeneralAssumptionCheck;
112  example safeVarName;
113  example extendWeyl;
114  example polyVars;
115  example monomialInIdeal;
116  example vas2pars;
117  example killTerms;
118  example minIntRoot2;
119  example maxIntRoot;
120  example dmodAction;
121  example dmodActionRat;
122  example simplifyRat;
124  example multRat;
125  example diffRat;
126  example commRing;
127  example holonomicRank;
128  example DsingularLocus;
129  example rightNFWeyl;
130  example Dlocalization;
131  example WeylClosure1;
132  example WeylClosure;
133  example polSol;
134  example polSolFiniteRank;
135  example ratSol;
136  example bfctBound;
137  example annRatSyz;
138}
139
140
141// tools //////////////////////////////////////////////////////////////////////
142
143proc dmodGeneralAssumptionCheck ()
144"
145USAGE:    dmodGeneralAssumptionCheck();
146RETURN:   nothing, but checks general assumptions on the basering
147NOTE:     This procedure checks the following conditions on the basering R
148          and prints an error message if any of them is violated:
149@*         - R is the n-th Weyl algebra over a field of characteristic 0,
150@*         - R is not a qring,
151@*         - for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
152             holds, i.e. the sequence of variables is given by
153             x(1),...,x(n),D(1),...,D(n), where D(i) is the differential
154             operator belonging to x(i).
155EXAMPLE:  example dmodGeneralAssumptionCheck; shows examples
156"
157{
158  // char K <> 0, qring
159  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
160  {
161    ERROR("Basering is inappropriate: characteristic>0 or qring present");
162  }
163  // no Weyl algebra
164  if (isWeyl() == 0)
165  {
166    ERROR("Basering is not a Weyl algebra");
167  }
168  // wrong sequence of vars
169  int i,n;
170  n = nvars(basering) div 2;
171  for (i=1; i<=n; i++)
172  {
173    if (bracket(var(i+n),var(i))<>1)
174    {
175      ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i)));
176    }
177  }
178  return();
179}
180example
181{
182  "EXAMPLE"; echo=2;
183  ring r = 0,(x,D),dp;
184  dmodGeneralAssumptionCheck(); // prints error message
185  def W = Weyl();
186  setring W;
187  dmodGeneralAssumptionCheck(); // returns nothing
188}
189
190
191proc safeVarName (string s)
192"
193USAGE:    safeVarName(s);  s string
194RETURN:   string, returns s if s is not the name of a par/var of basering
195          and '@' + s otherwise
196EXAMPLE:  example safeVarName; shows examples
197"
198{
199  string S = "," + charstr(basering) + "," + varstr(basering) + ",";
200  s = "," + s + ",";
201  while (find(S,s) <> 0)
202  {
203    s[1] = "@";
204    s = "," + s;
205  }
206  s = s[2..size(s)-1];
207  return(s);
208}
209example
210{
211  "EXAMPLE:"; echo = 2;
212  ring r = (0,a),(w,@w,x,y),dp;
213  safeVarName("a");
214  safeVarName("x");
215  safeVarName("z");
216  safeVarName("w");
217}
218
219
220proc extendWeyl (def newVars)
221"
222USAGE:    extendWeyl(S);  S string or list of strings
223ASSUME:   The basering is the n-th Weyl algebra over a field of
224          characteristic 0 and for all 1<=i<=n the identity
225          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
226          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
227          is the differential operator belonging to x(i).
228RETURN:   ring, Weyl algebra extended by vars given by S
229EXAMPLE:  example extendWeyl; shows examples
230"
231{
232  dmodGeneralAssumptionCheck();
233  int i,s;
234  string inpt = typeof(newVars);
235  list L;
236  if (inpt=="string")
237  {
238    s = 1;
239    L = newVars;
240  }
241  else
242  {
243    if (inpt=="list")
244    {
245      s = size(newVars);
246      if (s<1)
247      {
248        ERROR("No new variables specified.");
249      }
250      for (i=1; i<=s; i++)
251      {
252        if (typeof(newVars[i]) <> "string")
253        {
254          ERROR("Entries of input list must be of type string.");
255        }
256      }
257      L = newVars;
258    }
259    else
260    {
261      ERROR("Expected string or list of strings as input.");
262    }
263  }
264  def save = basering;
265  int n = nvars(save) div 2;
266  list RL = ringlist(save);
267  RL = RL[1..4];
268  list Ltemp = L;
269  for (i=s; i>0; i--)
270  {
271    Ltemp[n+s+i] = "D" + newVars[i];
272  }
273  for (i=n; i>0; i--)
274  {
275    Ltemp[s+i]     = RL[2][i];
276    Ltemp[n+2*s+i] = RL[2][n+i];
277  }
278  RL[2] = Ltemp;
279  Ltemp = list();
280  Ltemp[1] = list("dp",intvec(1:(2*n+2*s)));
281  Ltemp[2] = list("C",intvec(0));
282  RL[3] = Ltemp;
283  kill Ltemp;
284  def @Dv = ring(RL);
285  setring @Dv;
286  def Dv = Weyl();
287  setring save;
288  return(Dv);
289}
290example
291{
292  "EXAMPLE:"; echo = 2;
293  ring @D2 = 0,(x,y,Dx,Dy),dp;
294  def D2 = Weyl();
295  setring D2;
296  def D3 = extendWeyl("t");
297  setring D3; D3;
298  list L = "u","v";
299  def D5 = extendWeyl(L);
300  setring D5;
301  D5;
302}
303
304
305proc polyVars (poly f, intvec v)
306"
307USAGE:    polyVars(f,v);  f poly, v intvec
308RETURN:   int, 1 if f contains only variables indexed by v, 0 otherwise
309EXAMPLE:  example polyVars; shows examples
310"
311{
312  ideal varsf = variables(f); // vars contained in f
313  ideal V;
314  int i;
315  int n = nvars(basering);
316  for (i=1; i<=nrows(v); i++)
317  {
318    if ( (v[i]<0) || (v[i]>n) )
319    {
320      ERROR("var(" + string(v[i]) + ") out of range");
321    }
322    V[i] = var(v[i]);
323  }
324  attrib(V,"isSB",1);
325  ideal N = NF(varsf,V);
326  N = simplify(N,2);
327  if (N[1]==0)
328  {
329    return(1);
330  }
331  else
332  {
333    return(0);
334  }
335}
336example
337{
338  "EXAMPLE:"; echo = 2;
339  ring r = 0,(x,y,z),dp;
340  poly f = y^2+zy;
341  intvec v = 1,2;
342  polyVars(f,v); // does f depend only on x,y?
343  v = 2,3;
344  polyVars(f,v); // does f depend only on y,z?
345}
346
347
348proc monomialInIdeal (ideal I)
349"
350USAGE:    monomialInIdeal(I);  I ideal
351RETURN:   ideal consisting of all monomials appearing in generators of ideal
352EXAMLPE:  example monomialInIdeal; shows examples
353"
354{
355  // returns ideal consisting of all monomials appearing in generators of ideal
356  I = simplify(I,2+8);
357  int i;
358  poly p;
359  ideal M;
360  for (i=1; i<=size(I); i++)
361  {
362    p = I[i];
363    while (p<>0)
364    {
366      p = p - lead(p);
367    }
368  }
369  M = simplify(M,4+2);
370  return(M);
371}
372example
373{
374  "EXAMPLE"; echo=2;
375  ring r = 0,(x,y),dp;
376  ideal I = x2+5x3y7, x-x2-6xy;
377  monomialInIdeal(I);
378}
379
380
381proc vars2pars (intvec v)
382"
383USAGE:    vars2pars(v);  v intvec
384ASSUME:   The basering is commutative.
385RETURN:   ring with variables specified by v converted into parameters
386EXAMPLE:  example vars2pars; shows examples
387"
388{
389  if (isCommutative() == 0)
390  {
391    ERROR("The basering must be commutative.");
392  }
393  v = sortIntvec(v)[1];
394  int sv = size(v);
395  if ( (v[1]<1) || (v[sv]<1) )
396  {
397    ERROR("Expected entries of intvec in the range 1.."+string(n));
398  }
399  def save = basering;
400  int i,j,n;
401  n = nvars(save);
402  list RL = ringlist(save);
403  list Lp,Lv,L1;
404  if (typeof(RL[1]) == "list")
405  {
406    L1 = RL[1];
407    Lp = L1[2];
408  }
409  else
410  {
411    L1[1] = RL[1];
412    L1[4] = ideal(0);
413  }
414  j = sv;
415  for (i=1; i<=n; i++)
416  {
417    if (j>0)
418    {
419      if (v[j]==i)
420      {
421        Lp[size(Lp)+1] = string(var(i));
422        j--;
423      }
424      else
425      {
426        Lv[size(Lv)+1] = string(var(i));
427      }
428    }
429    else
430    {
431      Lv[size(Lv)+1] = string(var(i));
432    }
433  }
434  RL[2] = Lv;
435  L1[2] = Lp;
436  L1[3] = list(list("lp",intvec(1:size(Lp))));
437  RL[1] = L1;
438  L1 = list();
439  L1[1] = list("dp",intvec(1:sv));
440  L1[2] = list("C",intvec(0));
441  RL[3] = L1;
442//   RL;
443  def R = ring(RL);
444  return(R);
445}
446example
447{
448  "EXAMPLE:"; echo = 2;
449  ring r = 0,(x,y,z,a,b,c),dp;
450  intvec v = 4,5,6;
451  def R = vars2pars(v);
452  setring R;
453  R;
454  v = 1,2;
455  def RR = vars2pars(v);
456  setring RR;
457  RR;
458}
459
460
461proc killTerms (poly p, int k)
462"
463USAGE:    killTerms(p,k);  p poly, k int
464RETURN:   poly, p minus all terms containing the k-th variable
465NOTE:     No Groebner basis computations are used.
466EXAMPLE:  example killTerms; shows examples
467"
468{
469  if ( (k<0) || (k>nvars(basering)) )
470  {
471    ERROR("Given int is out of range");
472  }
473  int j;
474  poly q;
475  intvec v;
476  for (j=1; j<=size(p); j++)
477  {
479    if (v[k]==0)
480    {
481      q = q + p[j];
482    }
483  }
484  return(q);
485}
486example
487{
488  "EXAMPLE:"; echo = 2;
489  ring r = 0,(x,y,z),dp;
490  poly p = 8xyz + 7xy + 6xz + 5yz + 4x + 3y + 2z + 1;
491  killTerms(p,1);
492  killTerms(p,2);
493  killTerms(p,3);
494}
495
496
497static proc minMaxIntRoot (list L, string minmax)
498{
499  int win;
500  if (size(L)>1)
501  {
502    if ( (typeof(L[1])<>"ideal") || (typeof(L[2])<>"intvec") )
503    {
504      win = 1;
505    }
506  }
507  else
508  {
509    win = 1;
510  }
511  if (win)
512  {
513    ERROR("Expected list in the format of bFactor");
514  }
515  if (size(L)>2)
516  {
517    if ( (L[3]=="1") || (L[3]=="0") )
518    {
519      print("// Warning: Constant poly. Returning 0.");
520      return(int(0));
521    }
522  }
523  ideal I = L[1];
524  int i,k,b;
525  if (minmax=="min")
526  {
527    i = ncols(I);
528    k = -1;
529    b = 0;
530  }
531  else // minmax=="max"
532  {
533    i = 1;
534    k = 1;
535    b = ncols(I);
536  }
537  for (; k*i<k*b; i=i+k)
538  {
540    {
542    }
543  }
544  print("// Warning: No integer root found. Returning 0.");
545  return(int(0));
546}
547
548
549//TODO rename? minIntRoot is name of proc in dmod.lib
550proc minIntRoot2 (list L)
551"
552USAGE:    minIntRoot2(L);  L list
553ASSUME:   L is the output of bFactor.
554RETURN:   int, the minimal integer root in a list of roots
556EXAMPLE:  example minIntRoot2; shows examples
557"
558{
559  return(minMaxIntRoot(L,"min"));
560}
561example
562{
563  "EXAMPLE"; echo=2;
564  ring r = 0,x,dp;
565  poly f = x*(x+1)*(x-2)*(x-5/2)*(x+5/2);
566  list L = bFactor(f);
567  minIntRoot2(L);
568}
569
570
571proc maxIntRoot (list L)
572"
573USAGE:    maxIntRoot(L);  L list
574ASSUME:   L is the output of bFactor.
575RETURN:   int, the maximal integer root in a list of roots
577EXAMPLE:  example maxIntRoot; shows examples
578"
579{
580  return(minMaxIntRoot(L,"max"));
581}
582example
583{
584  "EXAMPLE"; echo=2;
585  ring r = 0,x,dp;
586  poly f = x*(x+1)*(x-2)*(x-5/2)*(x+5/2);
587  list L = bFactor(f);
588  maxIntRoot(L);
589}
590
591
592proc dmodAction (def id, poly f, list #)
593"
594USAGE:    dmodAction(id,f[,v]);  id ideal or poly, f poly, v optional intvec
595ASSUME:   If v is not given, the basering is the n-th Weyl algebra W over a
596          field of characteristic 0 and for all 1<=i<=n the identity
597          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
598          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
599          differential operator belonging to x(i).
600          Otherwise, v is assumed to specify positions of variables, which form
601          a Weyl algebra as a subalgebra of the basering:
602          If size(v) equals 2*n, then bracket(var(v[i]),var(v[j])) must equal
603          1 if and only if j equals i+n, and 0 otherwise,  for all 1<=i,j<=n.
604@*        Further, assume that f does not contain any D(i).
605RETURN:   same type as id, the result of the natural D-module action of id on f
606NOTE:     The assumptions made are not checked.
607EXAMPLE:  example dmodAction; shows examples
608"
609{
610  string inp1 = typeof(id);
611  if ((inp1<>"poly") && (inp1<>"ideal"))
612  {
613    ERROR("Expected first argument to be poly or ideal but received "+inp1);
614  }
615  intvec posXD = 1..nvars(basering);
616  if (size(#)>0)
617  {
618    if (typeof(#[1])=="intvec")
619    {
620      posXD = #[1];
621    }
622  }
623  if ((size(posXD) mod 2)<>0)
624  {
625    ERROR("Even number of variables expected.")
626  }
627  int n = (size(posXD)) div 2;
628  int i,j,k,l;
629  ideal resI = id;
630  int sid = ncols(resI);
631  intvec v;
632  poly P,h;
633  for (l=1; l<=sid; l++)
634  {
635    P = resI[l];
636    resI[l] = 0;
637    for (i=1; i<=size(P); i++)
638    {
640      h = f;
641      for (j=1; j<=n; j++)
642      {
643        for (k=1; k<=v[posXD[j+n]]; k++) // action of Dx
644        {
645          h = diff(h,var(posXD[j]));
646        }
647        h = h*var(posXD[j])^v[posXD[j]]; // action of x
648      }
650      resI[l] = resI[l] + h;
651    }
652  }
653  if (inp1 == "ideal")
654  {
655    return(resI);
656  }
657  else
658  {
659    return(resI[1]);
660  }
661}
662example
663{
664  ring r = 0,(x,y,z),dp;
665  poly f = x^2*z - y^3;
666  def A = annPoly(f);
667  setring A;
668  poly f = imap(r,f);
669  dmodAction(LD,f);
670  poly P = y*Dy+3*z*Dz-3;
671  dmodAction(P,f);
672  dmodAction(P[1],f);
673}
674
675
676static proc checkRatInput (vector I)
677{
678  // check for valid input
679  int wrginpt;
680  if (nrows(I)<>2)
681  {
682    wrginpt = 1;
683  }
684  else
685  {
686    if (I[2] == 0)
687    {
688      wrginpt = 1;
689    }
690  }
691  if (wrginpt)
692  {
693    ERROR("Vector must consist of exactly two components, second one not 0");
694  }
695  return();
696}
697
698
699proc dmodActionRat(def id, vector w)
700"
701USAGE:    dmodActionRat(id,w);  id ideal or poly, f vector
702ASSUME:   The basering is the n-th Weyl algebra W over a field of
703          characteristic 0 and for all 1<=i<=n the identity
704          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
705          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
706          differential operator belonging to x(i).
707@*        Further, assume that w has exactly two components, second one not 0,
708          and that w does not contain any D(i).
709RETURN:   same type as id, the result of the natural D-module action of id on
710          the rational function w[1]/w[2]
711EXAMPLE:  example dmodActionRat; shows examples
712"
713{
714  string inp1 = typeof(id);
715  if ( (inp1<>"poly") && (inp1<>"ideal") )
716  {
717    ERROR("Expected first argument to be poly or ideal but received " + inp1);
718  }
719  checkRatInput(w);
720  poly f = w[1];
721  finKx(f);
722  f = w[2];
723  finKx(f);
724  def save = basering;
725  def r = commRing();
726  setring r;
727  ideal I = imap(save,id);
728  vector w = imap(save,w);
729  int i,j,k,l;
730  int n = nvars(basering) div 2;
731  int sid = ncols(I);
732  intvec v;
733  poly P;
734  vector h,resT;
735  module resL;
736  for (l=1; l<=sid; l++)
737  {
738    P = I[l];
739    resT = [0,1];
740    for (i=1; i<=size(P); i++)
741    {
743      h = w;
744      for (j=1; j<=n; j++)
745      {
746        for (k=1; k<=v[j+n]; k++) // action of Dx
747        {
748          h = diffRat(h,j);
749        }
750        h = h + h[1]*(var(j)^v[j]-1)*gen(1);  // action of x
751      }
752      h = h + (leadcoef(P[i])-1)*h[1]*gen(1);
754    }
755    resL[l] = resT;
756  }
757  setring save;
758  module resL = imap(r,resL);
759  return(resL);
760}
761example
762{
763  "EXAMPLE:"; echo = 2;
764  ring r = 0,(x,y),dp;
765  poly f = 2*x;  poly g = y;
766  def A = annRat(f,g); setring A;
767  poly f = imap(r,f); poly g = imap(r,g);
768  vector v = [f,g]; // represents f/g
769  // x and y act by multiplication
770  dmodActionRat(x,v);
771  dmodActionRat(y,v);
772  // Dx and Dy act by partial derivation
773  dmodActionRat(Dx,v);
774  dmodActionRat(Dy,v);
775  dmodActionRat(x*Dx+y*Dy,v);
776  setring r;
777  f = 2*x*y; g = x^2 - y^3;
778  def B = annRat(f,g); setring B;
779  poly f = imap(r,f); poly g = imap(r,g);
780  vector v = [f,g];
781  dmodActionRat(LD,v); // hence LD is indeed the annihilator of f/g
782}
783
784
785static proc arithmeticRat (vector I, vector J, string op, list #)
786{
787  // op = "+": return I+J
788  // op = "*": return I*J
789  // op = "s": return simplified I
790  // op = "d": return diff(I,var(#[1]))
791  int isComm = isCommutative();
792  if (!isComm)
793  {
794    def save = basering;
795    def r = commRing();
796    setring r;
797    ideal m = maxideal(1);
798    map f = save,m;
799    vector I = f(I);
800    vector J = f(J);
801  }
802  vector K;
803  poly p;
804  if (op == "s")
805  {
806    p = gcd(I[1],I[2]);
807    K = (I[1]/p)*gen(1) + (I[2]/p)*gen(2);
808  }
809  else
810  {
811    if (op == "+")
812    {
813      I = arithmeticRat(I,vector(0),"s");
814      J = arithmeticRat(J,vector(0),"s");
815      p = lcm(I[2],J[2]);
816      K = (I[1]*p/I[2] + J[1]*p/J[2])*gen(1) + p*gen(2);
817    }
818    else
819    {
820      if (op == "*")
821      {
822        K = (I[1]*J[1])*gen(1) + (I[2]*J[2])*gen(2);
823      }
824      else
825      {
826        if (op == "d")
827        {
828          int j = #[1];
829          K = (diff(I[1],var(j))*I[2] - I[1]*diff(I[2],var(j)))*gen(1)+ (I[2]^2)*gen(2);
830        }
831      }
832    }
833    K = arithmeticRat(K,vector(0),"s");
834  }
835  if (!isComm)
836  {
837    setring save;
838    vector K = imap(r,K);
839  }
840  return(K);
841}
842
843
844proc simplifyRat (vector J)
845"
846USAGE:    simplifyRat(v);  v vector
847ASSUME:   Assume that v has exactly two components, second one not 0.
848RETURN:   vector, representing simplified rational function v[1]/v[2]
849NOTE:     Possibly present non-commutative relations of the basering are
850          ignored.
851EXAMPLE:  example simplifyRat; shows examples
852"
853{
854  checkRatInput(J);
855  return(arithmeticRat(J,vector(0),"s"));
856}
857example
858{
859  "EXAMPLE:"; echo = 2;
860   ring r = 0,(x,y),dp;
861  vector v = [x2-1,x+1];
862  simplifyRat(v);
863  simplifyRat(v) - [x-1,1];
864}
865
866
867proc addRat (vector I, vector J)
868"
870ASSUME:   Assume that v,w have exactly two components, second ones not 0.
871RETURN:   vector, representing rational function (v[1]/v[2])+(w[1]/w[2])
872NOTE:     Possibly present non-commutative relations of the basering are
873          ignored.
875"
876{
877  checkRatInput(I);
878  checkRatInput(J);
879  return(arithmeticRat(I,J,"+"));
880}
881example
882{
883  "EXAMPLE:"; echo = 2;
884  ring r = 0,(x,y),dp;
885  vector v = [x,y];
886  vector w = [y,x];
889}
890
891
892proc multRat (vector I, vector J)
893"
894USAGE:    multRat(v,w);  v,w vectors
895ASSUME:   Assume that v,w have exactly two components, second ones not 0.
896RETURN:   vector, representing rational function (v[1]/v[2])*(w[1]/w[2])
897NOTE:     Possibly present non-commutative relations of the basering are
898          ignored.
899EXAMPLE:  example multRat; shows examples
900"
901{
902  checkRatInput(I);
903  checkRatInput(J);
904  return(arithmeticRat(I,J,"*"));
905}
906example
907{
908  "EXAMPLE:"; echo = 2;
909  ring r = 0,(x,y),dp;
910  vector v = [x,y];
911  vector w = [y,x];
912  multRat(v,w);
913  multRat(v,w) - [1,1];
914}
915
916
917proc diffRat (vector I, int j)
918"
919USAGE:    diffRat(v,j);  v vector, j int
920ASSUME:   Assume that v has exactly two components, second one not 0.
921RETURN:   vector, representing rational function derivative of rational
922          function (v[1]/v[2]) w.r.t. var(j)
923NOTE:     Possibly present non-commutative relations of the basering are
924          ignored.
925EXAMPLE:  example diffRat; shows examples
926"
927{
928  checkRatInput(I);
929  if ( (j<1) || (j>nvars(basering)) )
930  {
931    ERROR("Second argument must be in the range 1.."+string(nvars(basering)));
932  }
933  return(arithmeticRat(I,vector(0),"d",j));
934}
935example
936{
937  "EXAMPLE:"; echo = 2;
938  ring r = 0,(x,y),dp;
939  vector v = [x,y];
940  diffRat(v,1);
941  diffRat(v,1) - [1,y];
942  diffRat(v,2);
943  diffRat(v,2) - [-x,y2];
944}
945
946
947proc commRing ()
948"
949USAGE:    commRing();
950RETURN:   ring, basering without non-commutative relations
951EXAMPLE:  example commRing; shows examples
952"
953{
954  list RL = ringlist(basering);
955  if (size(RL)<=4)
956  {
957    return(basering);
958  }
959  RL = RL[1..4];
960  def r = ring(RL);
961  return(r)
962}
963example
964{
965  "EXAMPLE:"; echo = 2;
966  def W = makeWeyl(3);
967  setring W; W;
968  def W2 = commRing();
969  setring W2; W2;
970  ring r = 0,(x,y),dp;
971  def r2 = commRing(); // same as r
972  setring r2; r2;
973}
974
975
976// TODO remove this proc once chern.lib is ready
977static proc orderedPartition(int n, list #)
978"
979USUAGE:  orderedPartition(n,a); n,a positive ints
980         orderedPartition(n,w); n positive int, w positive intvec
981RETURN:  list of intvecs
982PURPOSE: Computes all partitions of n of length a, if the second
983         argument is an int, or computes all weighted partitions
984         w.r.t. w of n of length size(w) if the second argument
985         is an intvec.
986         In both cases, zero parts are included.
987EXAMPLE: example orderedPartition; shows an example
988"
989{
990  int a,wrongInpt,intInpt;
991  intvec w = 1;
992  if (size(#)>0)
993  {
994    if (typeof(#[1]) == "int")
995    {
996      a = #[1];
997      intInpt = 1;
998    }
999    else
1000    {
1001      if (typeof(#[1]) == "intvec")
1002      {
1003        w = #[1];
1004        a = size(w);
1005      }
1006      else
1007      {
1008        wrongInpt = 1;
1009      }
1010    }
1011  }
1012  else
1013  {
1014    wrongInpt = 1;
1015  }
1016  if (wrongInpt)
1017  {
1018    ERROR("Expected second argument of type int or intvec.");
1019  }
1020  kill wrongInpt;
1021  if (n==0 && a>0)
1022  {
1023    return(list(0:a));
1024  }
1025  if (n<=0 || a<=0 || allPositive(w)==0)
1026  {
1027    ERROR("Positive arguments expected.");
1028  }
1029  int baseringdef;
1030  if (defined(basering)) // if a basering is defined, it should be saved for later use
1031  {
1032    def save = basering;
1033    baseringdef = 1;
1034  }
1035  ring r = 0,(x(1..a)),dp; // all variables for partition of length a
1036  ideal M;
1037  if (intInpt)
1038  {
1039    M = maxideal(n);   // all monomials of total degree n
1040  }
1041  else
1042  {
1043    M = weightKB(ideal(0),n,w); // all monomials of total weighted degree n
1044  }
1045  list L;
1046  int i;
1047  for (i = 1; i <= ncols(M); i++) {L = insert(L,leadexp(M[i]));}
1048  // the leadexp corresponds to a partition
1049  if (baseringdef) // sets the old ring as basering again
1050  {
1051    setring save;
1052  }
1053  return(L); //returns the list of partitions
1054}
1055example
1056{
1057 "EXAMPLE"; echo = 2;
1058 orderedPartition(4,2);
1059 orderedPartition(5,3);
1060 orderedPartition(2,4);
1061 orderedPartition(8,intvec(2,3));
1062 orderedPartition(7,intvec(2,2)); // no such partition
1063}
1064
1065
1066// applications of characteristic variety /////////////////////////////////////
1067
1068proc holonomicRank (ideal I, list #)
1069"
1070USAGE:    holonomicRank(I[,e]);   I ideal, e optional int
1071ASSUME:   The basering is the n-th Weyl algebra over a field of
1072          characteristic 0 and for all 1<=i<=n the identity
1073          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
1074          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
1075          is the differential operator belonging to x(i).
1076RETURN:   int, the holonomic rank of I
1077REMARKS:  The holonomic rank of I is defined to be the K(x(1..n))-dimension of
1078          I regarded as ideal in the rational Weyl algebra K(x(1..n))<D(1..n)>.
1079          If this dimension is not finite, -1 is returned.
1080NOTE:     If e<>0, @code{std} is used for Groebner basis computations,
1081          otherwise (and by default) @code{slimgb} is used.
1082@*        If printlevel=1, progress debug messages will be printed,
1083          if printlevel>=2, all the debug messages will be printed.
1084EXAMPLE:  example holonomicRank; shows examples
1085"
1086{
1087  // assumption check is done by charVariety
1088  int ppl = printlevel - voice + 2;
1089  int eng;
1090  if (size(#)>0)
1091  {
1092    if(typeof(#[1])=="int")
1093    {
1094      eng = #[1];
1095    }
1096  }
1097  def save = basering;
1098  dbprint(ppl  ,"// Computing characteristic variety...");
1099  def grD = charVariety(I);
1100  setring grD; // commutative ring
1101  dbprint(ppl  ,"// ...done.");
1102  dbprint(ppl-1,"// " + string(charVar));
1103  int n = nvars(save) div 2;
1104  intvec v = 1..n;
1105  def R = vars2pars(v);
1106  setring R;
1107  ideal J = imap(grD,charVar);
1108  dbprint(ppl  ,"// Starting GB computation...");
1109  J = engine(J,0); // use slimgb
1110  dbprint(ppl  ,"// ...done.");
1111  dbprint(ppl-1,"// " + string(J));
1112  int d = vdim(J);
1113  setring save;
1114  return(d);
1115}
1116example
1117{
1118  "EXAMPLE:"; echo = 2;
1119  // (OTW), Example 8
1120  ring r3 = 0,(x,y,z,Dx,Dy,Dz),dp;
1121  def D3 = Weyl();
1122  setring D3;
1123  poly f = x^3-y^2*z^2;
1124  ideal I = f^2*Dx+3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z;
1125  // I annihilates exp(1/f)
1126  holonomicRank(I);
1127}
1128
1129
1130proc DsingularLocus (ideal I)
1131"
1132USAGE:    DsingularLocus(I);  I ideal
1133ASSUME:   The basering is the n-th Weyl algebra over a field of
1134          characteristic 0 and for all 1<=i<=n the identity
1135          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
1136          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
1137          is the differential operator belonging to x(i).
1138RETURN:   ideal, the singular locus of the D-module D/I
1139NOTE:     If printlevel>=1, progress debug messages will be printed,
1140          if printlevel>=2, all the debug messages will be printed
1141EXAMPLE:  example DsingularLocus; shows examples
1142"
1143{
1144  // assumption check is done by charVariety
1145  int ppl = printlevel - voice + 2;
1146  def save = basering;
1147  dbprint(ppl  ,"// Computing characteristic variety...");
1148  def grD = charVariety(I);
1149  setring grD;
1150  dbprint(ppl  ,"// ...done");
1151  dbprint(ppl-1,"// " + string(charVar));
1152  poly pDD = 1;
1153  ideal IDD;
1154  int i;
1155  int n = nvars(basering) div 2;
1156  for (i=1; i<=n; i++)
1157  {
1158    pDD = pDD*var(i+n);
1159    IDD[i] = var(i+n);
1160  }
1161  dbprint(ppl  ,"// Computing saturation...");
1162  ideal S = sat(charVar,IDD)[1];
1163  dbprint(ppl  ,"// ...done");
1164  dbprint(ppl-1,"// " + string(S));
1165  dbprint(ppl  ,"// Computing elimination...");
1166  S = eliminate(S,pDD);
1167  dbprint(ppl  ,"// ...done");
1168  dbprint(ppl-1,"// " + string(S));
1171  dbprint(ppl  ,"// ...done");
1172  dbprint(ppl-1,"// " + string(S));
1173  setring save;
1174  ideal S = imap(grD,S);
1175  return(S);
1176}
1177example
1178{
1179  "EXAMPLE:"; echo = 2;
1180  // (OTW), Example 8
1181  ring @D3 = 0,(x,y,z,Dx,Dy,Dz),dp;
1182  def D3 = Weyl();
1183  setring D3;
1184  poly f = x^3-y^2*z^2;
1185  ideal I = f^2*Dx + 3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z;
1186  // I annihilates exp(1/f)
1187  DsingularLocus(I);
1188}
1189
1190
1191// localization ///////////////////////////////////////////////////////////////
1192
1193static proc finKx(poly f)
1194{
1195  int n = nvars(basering) div 2;
1196  intvec iv = 1..n;
1197  if (polyVars(f,iv) == 0)
1198  {
1199    ERROR("Given poly may not contain differential operators.");
1200  }
1201  return();
1202}
1203
1204
1205proc rightNFWeyl (def id, int k)
1206"
1207USAGE:    rightNFWeyl(id,k);  id ideal or poly, k int
1208ASSUME:   The basering is the n-th Weyl algebra over a field of
1209          characteristic 0 and for all 1<=i<=n the identity
1210          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
1211          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
1212          is the differential operator belonging to x(i).
1213RETURN:   same type as id, the right normal form of id with respect to the
1214          principal right ideal generated by the k-th variable
1215NOTE:     No Groebner basis computation is used.
1216EXAMPLE:  example rightNFWeyl; shows examples.
1217"
1218{
1219  dmodGeneralAssumptionCheck();
1220  string inpt = typeof(id);
1221  if (inpt=="ideal" || inpt=="poly")
1222  {
1223    ideal I = id;
1224  }
1225  else
1226  {
1227    ERROR("Expected first input to be of type ideal or poly.");
1228  }
1229  def save = basering;
1230  int n = nvars(save) div 2;
1231  if (0>k || k>2*n)
1232  {
1233    ERROR("Expected second input to be in the range 1.."+string(2*n)+".");
1234  }
1235  int i,j;
1236  if (k>n) // var(k) = Dx(k-n)
1237  {
1238    // switch var(k),var(k-n)
1239    list RL = ringlist(save);
1240    matrix rel = RL[6];
1241    rel[k-n,k] = -1;
1242    RL = RL[1..4];
1243    list L = RL[2];
1244    string str = L[k-n];
1245    L[k-n] = L[k];
1246    L[k] = str;
1247    RL[2] = L;
1248    def @W = ring(RL);
1249    kill L,RL,str;
1250    ideal @mm = maxideal(1);
1251    setring @W;
1252    matrix rel = imap(save,rel);
1253    def W = nc_algebra(1,rel);
1254    setring W;
1255    ideal @mm = imap(save,@mm);
1256    map mm = save,@mm;
1257    ideal I = mm(I);
1258    i = k-n;
1259  }
1260  else  // var(k) = x(k)
1261  {
1262    def W = save;
1263    i = k;
1264  }
1265  for (j=1; j<=ncols(I); j++)
1266  {
1267    I[j] = killTerms(I[j],i);
1268  }
1269  setring save;
1270  I = imap(W,I);
1271  if (inpt=="poly")
1272  {
1273    return(I[1]);
1274  }
1275  else
1276  {
1277    return(I);
1278  }
1279}
1280example
1281{
1282  "EXAMPLE:"; echo = 2;
1283  ring r = 0,(x,y,Dx,Dy),dp;
1284  def W = Weyl();
1285  setring W;
1286  ideal I = x^3*Dx^3, y^2*Dy^2, x*Dy, y*Dx;
1287  rightNFWeyl(I,1); // right NF wrt principal right ideal x*W
1288  rightNFWeyl(I,3); // right NF wrt principal right ideal Dx*W
1289  rightNFWeyl(I,2); // right NF wrt principal right ideal y*W
1290  rightNFWeyl(I,4); // right NF wrt principal right ideal Dy*W
1291  poly p = x*Dx+1;
1292  rightNFWeyl(p,1); // right NF wrt principal right ideal x*W
1293}
1294
1295
1296// TODO check OTW for assumptions on holonomicity
1297proc Dlocalization (ideal J, poly f, list #)
1298"
1299USAGE:    Dlocalization(I,f[,k,e]);  I ideal, f poly, k,e optional ints
1300ASSUME:   The basering is the n-th Weyl algebra over a field of
1301          characteristic 0 and for all 1<=i<=n the identity
1302          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
1303          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
1304          is the differential operator belonging to x(i).
1305@*        Further, assume that f does not contain any D(i) and that I is
1306          holonomic on K^n\V(f).
1307RETURN:   ideal or list, computes an ideal J such that D/J is isomorphic
1308          to D/I localized at f as D-modules.
1309          If k<>0, a list consisting of J and an integer m is returned,
1310          such that f^m represents the natural map from D/I to D/J.
1311          Otherwise (and by default), only the ideal J is returned.
1312REMARKS:  Reference: (OTW)
1313NOTE:     If e<>0, @code{std} is used for Groebner basis computations,
1314          otherwise (and by default) @code{slimgb} is used.
1315@*        If printlevel=1, progress debug messages will be printed,
1316          if printlevel>=2, all the debug messages will be printed.
1318EXAMPLE: example Dlocalization; shows examples
1319"
1320{
1321  dmodGeneralAssumptionCheck();
1322  finKx(f);
1323  int ppl = printlevel - voice + 2;
1324  int outList,eng;
1325  if (size(#)>0)
1326  {
1327    if (typeof(#[1])=="int" || typeof(#[1])=="number")
1328    {
1329      outList = int(#[1]);
1330    }
1331    if (size(#)>1)
1332    {
1333      if (typeof(#[2])=="int" || typeof(#[2])=="number")
1334      {
1335        eng = int(#[2]);
1336      }
1337    }
1338  }
1339  int i,j;
1340  def save = basering;
1341  int n = nvars(save) div 2;
1342  def Dv = extendWeyl(safeVarName("v"));
1343  setring Dv;
1344  poly f = imap(save,f);
1345  ideal phiI;
1346  for (i=n; i>0; i--)
1347  {
1348    phiI[i+n] = var(i+n+2)-var(1)^2*bracket(var(i+n+2),f)*var(n+2);
1349    phiI[i]   = var(i+1);
1350  }
1351  map phi = save,phiI;
1352  ideal J = phi(J);
1353  J = J, 1-f*var(1);
1354  // TODO original J has to be holonomic only on K^n\V(f), not on all of K^n
1355  // does is suffice to show that new J is holonomic on Dv??
1356  if (isHolonomic(J) == 0)
1357  {
1358    ERROR("Module is not holonomic.");
1359  }
1360  intvec w = 1; w[n+1]=0;
1361  ideal G = GBWeight(J,w,-w,eng);
1362  dbprint(ppl  ,"// found GB wrt weight " +string(-w));
1363  dbprint(ppl-1,"// " + string(G));
1364  intvec ww = w,-w;
1365  ideal inG = inForm(G,ww);
1366  inG = engine(inG,eng);
1367  poly s = var(1)*var(n+2); // s=v*Dv
1368  vector intersecvec = pIntersect(s,inG);
1369  s = vec2poly(intersecvec);
1370  s = subst(s,var(1),-var(1)-1);
1371  list L = bFactor(s);
1372  dbprint(ppl  ,"// found b-function");
1373  dbprint(ppl-1,"// roots: "+string(L[1]));
1374  dbprint(ppl-1,"// multiplicities: "+string(L[2]));
1375  kill inG,intersecvec,s;
1376  // TODO: use maxIntRoot
1377  L = intRoots(L);           // integral roots of b-function
1378  if (L[2]==0:size(L[2]))    // no integral roots
1379  {
1380    setring save;
1381    return(ideal(1));
1382  }
1383  intvec iv;
1384  for (i=1; i<=ncols(L[1]); i++)
1385  {
1386    iv[i] = int(L[1][i]);
1387  }
1388  int l0 = Max(iv);
1389  dbprint(ppl,"// maximal integral root is " +string(l0));
1390  kill L,iv;
1391  intvec degG;
1392  ideal Gk;
1393  for (j=1; j<=ncols(G); j++)
1394  {
1395    degG[j] = deg(G[j],ww);
1396    for (i=0; i<=l0-degG[j]; i++)
1397    {
1398      Gk[ncols(Gk)+1] = var(1)^i*G[j];
1399    }
1400  }
1401  Gk = rightNFWeyl(Gk,n+2);
1402  dbprint(ppl,"// found right normalforms");
1403  module M = coeffs(Gk,var(1));
1404  setring save;
1405  def mer = makeModElimRing(save);
1406  setring mer;
1407  module M = imap(Dv,M);
1408  kill Dv;
1409  M = engine(M,eng);
1410  dbprint(ppl  ,"// found GB of free module of rank " + string(l0+1));
1411  dbprint(ppl-1,"// " + string(M));
1412  M = prune(M);
1413  setring save;
1414  matrix M = imap(mer,M);
1415  kill mer;
1416  int ro = nrows(M);
1417  int co = ncols(M);
1418  ideal I;
1419  if (ro == 1) // nothing to do
1420  {
1421    I = M;
1422  }
1423  else
1424  {
1425    matrix zm[ro-1][1]; // zero matrix
1426    matrix v[ro-1][1];
1427    for (i=1; i<=co; i++)
1428    {
1429      v = M[1..ro-1,i];
1430      if (v == zm)
1431      {
1432        I[size(I)+1] = M[ro,i];
1433      }
1434    }
1435  }
1436  if (outList<>0)
1437  {
1438    return(list(I,l0+2));
1439  }
1440  else
1441  {
1442    return(I);
1443  }
1444}
1445example
1446{
1447  "EXAMPLE:"; echo = 2;
1448  // (OTW), Example 8
1449  ring r = 0,(x,y,z,Dx,Dy,Dz),dp;
1450  def W = Weyl();
1451  setring W;
1452  poly f = x^3-y^2*z^2;
1453  ideal I = f^2*Dx+3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z;
1454  // I annihilates exp(1/f)
1455  ideal J = Dlocalization(I,f);
1456  J;
1457  Dlocalization(I,f,1); // The natural map D/I -> D/J is given by 1/f^2
1458}
1459
1460
1461
1462// Weyl closure ///////////////////////////////////////////////////////////////
1463
1464static proc orderFiltrationD1 (poly f)
1465{
1466  // returns list of ideal and intvec
1467  // ideal contains x-parts, intvec corresponding degree in Dx
1468  poly g,h;
1469  g = f;
1470  ideal I;
1471  intvec v,w,u;
1472  w = 0,1;
1473  int i,j;
1474  i = 1;
1475  while (g<>0)
1476  {
1477    h = inForm(g,w);
1478    I[i] = 0;
1479    for (j=1; j<=size(h); j++)
1480    {
1482      u[i] = v[2];
1483      v[2] = 0;
1484      I[i] = I[i] + leadcoef(h[j])*monomial(v);
1485    }
1486    g = g-h;
1487    i++;
1488  }
1489  return(list(I,u));
1490}
1491
1492
1493static proc kerLinMapD1 (ideal W, poly L, poly p)
1494{
1495  // computes kernel of right multiplication with L viewed
1496  // as homomorphism of K-vector spaces span(W) -> D1/p*D1
1497  // assume p in K[x], basering is K<x,Dx>
1498  ideal G,K;
1499  G = std(p);
1500  list l;
1501  int i,j;
1502  // first, compute the image of span(W)
1503  if (simplify(W,2)[1] == 0)
1504  {
1505    return(K); // = 0
1506  }
1507  for (i=1; i<=size(W); i++)
1508  {
1509    l = orderFiltrationD1(W[i]*L);
1510    K[i] = 0;
1511    for (j=1; j<=size(l[1]); j++)
1512    {
1513      K[i] = K[i] + NF(l[1][j],G)*var(2)^(l[2][j]);
1514    }
1515  }
1516  // now, we get the kernel by linear algebra
1517  l = linReduceIdeal(K,1);
1518  i = ncols(l[1]) - size(l[1]);
1519  if (i<>0)
1520  {
1521    K = module(W)*l[2];
1522    K = K[1..i];
1523  }
1524  else
1525  {
1526    K = 0;
1527  }
1528  return(K);
1529}
1530
1531
1532static proc leftDivisionKxD1 (poly p, poly L)
1533{
1534  // basering is D1 = K<x,Dx>
1535  // p in K[x]
1536  // compute p^(-1)*L if p is a left divisor of L
1537//   if (rightNF(L,ideal(p))<>0)
1538//   {
1539//     ERROR("First poly is not a right factor of second poly");
1540//   }
1541  def save = basering;
1542  list l = orderFiltrationD1(L);
1543  ideal l1 = l[1];
1544  ring r = 0,x,dp;
1545  ideal l1 = fetch(save,l1);
1546  poly p = fetch(save,p);
1547  int i;
1548  for (i=1; i<=ncols(l1); i++)
1549  {
1550    l1[i] =  division(l1[i],p)[1][1,1];
1551  }
1552  setring save;
1553  ideal I = fetch(r,l1);
1554  poly f;
1555  for (i=1; i<=ncols(I); i++)
1556  {
1557    f = f + I[i]*var(2)^(l[2][i]);
1558  }
1559  return(f);
1560}
1561
1562
1563proc WeylClosure1 (poly L)
1564"
1565USAGE:    WeylClosure1(L);  L a poly
1566ASSUME:   The basering is the first Weyl algebra D=K<x,d|dx=xd+1> over a field
1567          K of characteristic 0.
1568RETURN:   ideal, the Weyl closure of the principal left ideal generated by L
1569REMARKS:  The Weyl closure of a left ideal I in the Weyl algebra D is defined
1570          to be the intersection of I regarded as left ideal in the rational
1571          Weyl algebra K(x)<d> with D.
1572@*        Reference: (Tsa), Algorithm 1.2.4
1573NOTE:     If printlevel=1, progress debug messages will be printed,
1574          if printlevel>=2, all the debug messages will be printed.
1576EXAMPLE:  example WeylClosure1; shows examples
1577"
1578{
1579  dmodGeneralAssumptionCheck(); // assumption check
1580  int ppl = printlevel - voice + 2;
1581  def save = basering;
1582  intvec w = 0,1; // for order filtration
1583  poly p = inForm(L,w);
1584  ring @R = 0,var(1),dp;
1585  ideal mm = var(1),1;
1586  map m = save,mm;
1587  ideal @p = m(p);
1588  poly p = @p[1];
1589  poly g = gcd(p,diff(p,var(1)));
1590  if (g == 1)
1591  {
1592    g = p;
1593  }
1594  ideal facp = factorize(g,1); // g is squarefree, constants aren't interesting
1595  dbprint(ppl-1,
1596          "// squarefree part of highest coefficient w.r.t. order filtration:");
1597  dbprint(ppl-1, "// " + string(facp));
1598  setring save;
1599  p = imap(@R,p);
1600  // 1-1 extend basering by parameter and introduce new var t=x*d
1601  list RL = ringlist(save);
1602  RL = RL[1..4];
1603  list l;
1604  l[1] = int(0);
1605  l[2] = list(safeVarName("a"));
1606  l[3] = list(list("lp",intvec(1)));
1607  l[4] = ideal(0);
1608  RL[1] = l;
1609  l = RL[2] + list(safeVarName("t"));
1610  RL[2] = l;
1611  l = list();
1612  l[1] = list("dp",intvec(1,1));
1613  l[2] = list("dp",intvec(1));
1614  l[3] = list("C",intvec(0));
1615  RL[3] = l;
1616  def @Wat = ring(RL);
1617  kill RL,l;
1618  setring @Wat;
1619  matrix relD[3][3];
1620  relD[1,2] = 1;
1621  relD[1,3] = var(1);
1622  relD[2,3] = -var(2);
1623  def Wat = nc_algebra(1,relD);
1624  setring Wat;
1625  kill @Wat;
1626  // 1-2 rewrite L using Euler operators
1627  ideal mm = var(1)+par(1),var(2);
1628  map m = save,mm;
1629  poly L = m(L);
1630  w = -1,1,0; // for Bernstein filtration
1631  int i = 1;
1632  ideal Q;
1633  poly p = L;
1634  intvec d;
1635  while (p<>0)
1636  {
1637    Q[i] = inForm(p,w);
1638    p = p - Q[i];
1639    d[i] = -deg(Q[i],w);
1640    i++;
1641  }
1642  ideal S = std(var(1)*var(2)-var(3));
1643  Q = NF(Q,S);
1644  dbprint(ppl,  "// found Euler representation of operator");
1645  dbprint(ppl-1,"// " + string(Q));
1646  Q = subst(Q,var(1),1);
1647  Q = subst(Q,var(2),1);
1648  // 1-3 prepare for algebraic extensions with minpoly = facp[i]
1649  list RL = ringlist(Wat);
1650  RL = RL[1..4];
1651  list l;
1652  l = string(var(3));
1653  RL[2] = l;
1654  l = list();
1655  l[1] = list("dp",intvec(1));
1656  l[2] = list("C",intvec(0));
1657  RL[3] = l;
1658  mm = par(1);
1659  m = @R,par(1);
1660  ideal facp = m(facp);
1661  kill @R,m,mm,l,S;
1662  intvec maxroots,testroots;
1663  int sq = size(Q);
1664  string strQ = "ideal Q = " + string(Q) + ";";
1665  // TODO do it without string workaround when issue with maps from
1666  //   transcendental to algebraic extension fields is fixed
1667  int j,maxr;
1668  // 2-1 get max int root of lowest nonzero entry of Q in algebraic extension
1669  for (i=1; i<=size(facp); i++)
1670  {
1671    testroots = 0;
1672    def Ra = ring(RL);
1673    setring Ra;
1674    ideal mm = 1,1,var(1);
1675    map m = Wat,mm;
1676    ideal facp = m(facp);
1678    execute(strQ);
1679    if (simplify(Q,2)[1] == poly(0))
1680    {
1681      break;
1682    }
1683    j = 1;
1684    while (j<sq)
1685    {
1686      if (Q[j]==0)
1687      {
1688        j++;
1689      }
1690      else
1691      {
1692        break;
1693      }
1694    }
1695    maxroots[i] = d[j]; // d[j] = r_k
1696    list LR = bFactor(Q[j]);
1697    LR = intRoots(LR);
1698    if (LR[2]<>0:size(LR[2])) // there are integral roots
1699    {
1700      for (j=1; j<=ncols(LR[1]); j++)
1701      {
1702        testroots[j] = int(LR[1][j]);
1703      }
1704      maxr = Max(testroots);
1705      if(maxr<0)
1706      {
1707        maxr = 0;
1708      }
1709      maxroots[i] = maxroots[i] + maxr;
1710    }
1711    kill LR;
1712    setring Wat;
1713    kill Ra;
1714  }
1715  maxr = Max(maxroots);
1716  // 3-1 build basis of vectorspace
1717  setring save;
1718  ideal KB;
1719  for (i=0; i<deg(p); i++)  // it's really <, not <=
1720  {
1721    for (j=0; j<=maxr; j++) // it's really <=, not <
1722    {
1723      KB[size(KB)+1] = monomial(intvec(i,j));
1724    }
1725  }
1726  dbprint(ppl,"// got vector space basis");
1727  dbprint(ppl-1, "// " + string(KB));
1728  // 3-2 get kernel of *L: span(KB)->D/pD
1729  KB = kerLinMapD1(KB,L,p);
1730  dbprint(ppl,"// got kernel");
1731  dbprint(ppl-1, "// " + string(KB));
1732  // 4-1 get (1/p)*f*L where f in KB
1733  for (i=1; i<=ncols(KB); i++)
1734  {
1735    KB[i] = leftDivisionKxD1(p,KB[i]*L);
1736  }
1737  KB = L,KB;
1738  // 4-2 done
1739  return(KB);
1740}
1741example
1742{
1743  "EXAMPLE:"; echo = 2;
1744  ring r = 0,(x,Dx),dp;
1745  def W = Weyl();
1746  setring W;
1747  poly L = (x^3+2)*Dx-3*x^2;
1748  WeylClosure1(L);
1749  L = (x^4-4*x^3+3*x^2)*Dx^2+(-6*x^3+20*x^2-12*x)*Dx+(12*x^2-32*x+12);
1750  WeylClosure1(L);
1751}
1752
1753
1754proc WeylClosure (ideal I)
1755"
1756USAGE:    WeylClosure(I);  I an ideal
1757ASSUME:   The basering is the n-th Weyl algebra W over a field of
1758          characteristic 0 and for all 1<=i<=n the identity
1759          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
1760          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
1761          differential operator belonging to x(i).
1762@*        Moreover, assume that the holonomic rank of W/I is finite.
1763RETURN:   ideal, the Weyl closure of I
1764REMARKS:  The Weyl closure of a left ideal I in the Weyl algebra W is defined to
1765          be the intersection of I regarded as left ideal in the rational Weyl
1766          algebra K(x(1..n))<D(1..n)> with W.
1767@*        Reference: (Tsa), Algorithm 2.2.4
1768NOTE:     If printlevel=1, progress debug messages will be printed,
1769          if printlevel>=2, all the debug messages will be printed.
1771EXAMPLE:  example WeylClosure; shows examples
1772"
1773{
1774  // assumption check
1775  dmodGeneralAssumptionCheck();
1776  if (holonomicRank(I)==-1)
1777  {
1778    ERROR("Input is not of finite holonomic rank.");
1779  }
1780  int ppl = printlevel - voice + 2;
1781  int eng = 0; // engine
1782  def save = basering;
1783  dbprint(ppl  ,"// Starting to compute singular locus...");
1784  ideal sl = DsingularLocus(I);
1785  sl = simplify(sl,2);
1786  dbprint(ppl  ,"// ...done.");
1787  dbprint(ppl-1,"// " + string(sl));
1788  if (sl[1] == 0) // can never get here
1789  {
1790    ERROR("Can't find polynomial in K[x] vanishing on singular locus.");
1791  }
1792  poly f = sl[1];
1793  dbprint(ppl  ,"// Found poly vanishing on singular locus: " + string(f));
1794  dbprint(ppl  ,"// Starting to compute localization...");
1795  list L = Dlocalization(I,f,1);
1796  ideal G = L[1];
1797  dbprint(ppl  ,"// ...done.");
1798  dbprint(ppl-1,"// " + string(G));
1799  dbprint(ppl  ,"// Starting to compute kernel of localization map...");
1800  if (eng == 0)
1801  {
1802    G = moduloSlim(f^L[2],G);
1803  }
1804  else
1805  {
1806    G = modulo(f^L[2],G);
1807  }
1808  dbprint(ppl  ,"// ...done.");
1809  return(G);
1810}
1811example
1812{
1813  "EXAMPLE:"; echo = 2;
1814  // (OTW), Example 8
1815  ring r = 0,(x,y,z,Dx,Dy,Dz),dp;
1816  def D3 = Weyl();
1817  setring D3;
1818  poly f = x^3-y^2*z^2;
1819  ideal I = f^2*Dx + 3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z;
1820  // I annihilates exp(1/f)
1821  WeylClosure(I);
1822}
1823
1824
1825
1826// solutions to systems of PDEs ///////////////////////////////////////////////
1827
1828proc polSol (ideal I, list #)
1829"
1830USAGE:    polSol(I[,w,m]);  I ideal, w optional intvec, m optional int
1831ASSUME:   The basering is the n-th Weyl algebra W over a field of
1832          characteristic 0 and for all 1<=i<=n the identity
1833          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
1834          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
1835          differential operator belonging to x(i).
1836@*        Moreover, assume that I is holonomic.
1837RETURN:   ideal, a basis of the polynomial solutions to the given system
1838REMARKS:  If w is given, w should consist of n strictly negative entries.
1839          Otherwise and by default, w is set to -1:n.
1840          In this case, w is used as weight vector for the computation of a
1841          b-function.
1842@*        If m is given, m is assumed to be the minimal integer root of the
1843          b-function of I w.r.t. w. Note that this assumption is not checked.
1844@*        Reference: (OTT), Algorithm 2.4
1845NOTE:     If printlevel=1, progress debug messages will be printed,
1846          if printlevel>=2, all the debug messages will be printed.
1848EXAMPLE:  example polSol; shows examples
1849"
1850{
1851  dmodGeneralAssumptionCheck();
1852  int ppl = printlevel - voice + 2;
1853  int mr,mrgiven;
1854  def save = basering;
1855  int n = nvars(save);
1856  intvec w = -1:(n div 2);
1857  if (size(#)>0)
1858  {
1859    if (typeof(#[1])=="intvec")
1860    {
1861      if (allPositive(-#[1]))
1862      {
1863        w = #[1];
1864      }
1865    }
1866    if (size(#)>1)
1867    {
1868      if (typeof(#[2])=="int")
1869      {
1870        mr = #[2];
1871        mrgiven = 1;
1872      }
1873    }
1874  }
1875  // Step 1: the b-function
1876  list L;
1877  if (!mrgiven)
1878  {
1879    if (!isHolonomic(I))
1880    {
1881      ERROR("Ideal is not holonomic. Try polSolFiniteRank.");
1882    }
1883    dbprint(ppl,"// Computing b-function...");
1884    L = bfctIdeal(I,w);
1885    dbprint(ppl,"// ...done.");
1886    dbprint(ppl-1,"//   Roots: " + string(L[1]));
1887    dbprint(ppl-1,"//   Multiplicities: " + string(L[2]));
1888    mr = minIntRoot2(L);
1889    dbprint(ppl,"// Minimal integer root is " + string(mr) + ".");
1890  }
1891  if (mr>0)
1892  {
1893    return(ideal(0));
1894  }
1895  // Step 2: get the form of a solution f
1896  int i;
1897  L = list();
1898  for (i=0; i<=-mr; i++)
1899  {
1900    L = L + orderedPartition(i,-w);
1901  }
1902  ideal mons;
1903  for (i=1; i<=size(L); i++)
1904  {
1905    mons[i] = monomial(L[i]);
1906  }
1907  kill L;
1908  mons = simplify(mons,2+4); // L might contain lots of 0s by construction
1909  ring @C = (0,@c(1..size(mons))),dummyvar,dp;
1910  def WC = save + @C;
1911  setring WC;
1912  ideal mons = imap(save,mons);
1913  poly f;
1914  for (i=1; i<=size(mons); i++)
1915  {
1916    f = f + par(i)*mons[i];
1917  }
1918  // Step 3: determine values of @c(i) by equating coefficients
1919  ideal I = imap(save,I);
1920  I = dmodAction(I,f,1..n);
1921  ideal M = monomialInIdeal(I);
1922  matrix CC = coeffs(I,M);
1923  int j;
1924  ideal C;
1925  for (i=1; i<=nrows(CC); i++)
1926  {
1927    f = 0;
1928    for (j=1; j<=ncols(CC); j++)
1929    {
1930      f = f + CC[i,j];
1931    }
1932    C[size(C)+1] = f;
1933  }
1934  // Step 3.1: solve a linear system
1935  ring rC = 0,(@c(1..size(mons))),dp;
1936  ideal C = imap(WC,C);
1937  matrix M = coeffs(C,maxideal(1));
1938  module MM = leftKernel(M);
1939  setring WC;
1940  module MM = imap(rC,MM);
1941  // Step 3.2: return the solution
1942  ideal F = ideal(MM*transpose(mons));
1943  setring save;
1944  ideal F = imap(WC,F);
1945  return(F);
1946}
1947example
1948{
1949  "EXAMPLE:"; echo=2;
1950  ring r = 0,(x,y,Dx,Dy),dp;
1951  def W = Weyl();
1952  setring W;
1953  poly tx,ty = x*Dx, y*Dy;
1954  ideal I =      // Appel F1 with parameters (2,-3,-2,5)
1955    tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3),
1956    ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2),
1957    (x-y)*Dx*Dy+2*Dx-3*Dy;
1958  intvec w = -1,-1;
1959  polSol(I,w);
1960}
1961
1962
1963static proc ex_polSol()
1964{  ring r = 0,(x,y,Dx,Dy),dp;
1965  def W = Weyl();
1966  setring W;
1967  poly tx,ty = x*Dx, y*Dy;
1968  ideal I =      // Appel F1 with parameters (2,-3,-2,5)
1969    tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3),
1970    ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2),
1971    (x-y)*Dx*Dy+2*Dx-3*Dy;
1972  intvec w = -5,-7;
1973  // the following gives a bug
1974  polSol(I,w);
1975  // this is due to a bug in weightKB, see ticket #339
1976  // http://www.singular.uni-kl.de:8002/trac/ticket/339
1977}
1978
1979
1980proc polSolFiniteRank (ideal I, list #)
1981"
1982USAGE:    polSolFiniteRank(I[,w]);  I ideal, w optional intvec
1983ASSUME:   The basering is the n-th Weyl algebra W over a field of
1984          characteristic 0 and for all 1<=i<=n the identity
1985          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
1986          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
1987          differential operator belonging to x(i).
1988@*        Moreover, assume that I is of finite holonomic rank.
1989RETURN:   ideal, a basis of the polynomial solutions to the given system
1990REMARKS:  If w is given, w should consist of n strictly negative entries.
1991          Otherwise and by default, w is set to -1:n.
1992          In this case, w is used as weight vector for the computation of a
1993          b-function.
1994@*        Reference: (OTT), Algorithm 2.6
1995NOTE:     If printlevel=1, progress debug messages will be printed,
1996          if printlevel>=2, all the debug messages will be printed.
1998EXAMPLE:  example polSolFiniteRank; shows examples
1999"
2000{
2001  dmodGeneralAssumptionCheck();
2002  if (holonomicRank(I)==-1)
2003  {
2004    ERROR("Ideal is not of finite holonomic rank.");
2005  }
2006  int ppl = printlevel - voice + 2;
2007  int n = nvars(basering) div 2;
2008  int eng;
2009  intvec w = -1:(n div 2);
2010  if (size(#)>0)
2011  {
2012    if (typeof(#[1])=="intvec")
2013    {
2014      if (allPositive(-#[1]))
2015      {
2016        w = #[1];
2017      }
2018    }
2019  }
2020  dbprint(ppl,"// Computing initial ideal...");
2021  ideal J = initialIdealW(I,-w,w);
2022  dbprint(ppl,"// ...done.");
2023  dbprint(ppl,"// Computing Weyl closure...");
2024  J = WeylClosure(J);
2025  J = engine(J,eng);
2026  dbprint(ppl,"// ...done.");
2027  poly s;
2028  int i;
2029  for (i=1; i<=n; i++)
2030  {
2031    s = s + w[i]*var(i)*var(i+n);
2032  }
2033  dbprint(ppl,"// Computing intersection...");
2034  vector v = pIntersect(s,J);
2035  list L = bFactor(vec2poly(v));
2036  dbprint(ppl-1,"//   roots: " + string(L[1]));
2037  dbprint(ppl-1,"//   multiplicities: " + string(L[2]));
2038  dbprint(ppl,"// ...done.");
2039  int mr =  minIntRoot2(L);
2040  int pl = printlevel;
2041  printlevel = printlevel + 1;
2042  ideal K = polSol(I,w,mr);
2043  printlevel = printlevel - 1;
2044  return(K);
2045}
2046example
2047{
2048  "EXAMPLE:"; echo=2;
2049  ring r = 0,(x,y,Dx,Dy),dp;
2050  def W = Weyl();
2051  setring W;
2052  poly tx,ty = x*Dx, y*Dy;
2053  ideal I =      // Appel F1 with parameters (2,-3,-2,5)
2054    tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3),
2055    ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2),
2056    (x-y)*Dx*Dy+2*Dx-3*Dy;
2057  intvec w = -1,-1;
2058  polSolFiniteRank(I,w);
2059}
2060
2061
2062static proc twistedIdeal(ideal I, poly f, intvec k, ideal F)
2063{
2064  // I subset D_n, f in K[x], F = factorize(f,1), size(k) = size(F), k[i]>0
2065  def save = basering;
2066  int n = nvars(save) div 2;
2067  int i,j;
2068  intvec a,v,w;
2069  w = (0:n),(1:n);
2070  for (i=1; i<=size(I); i++)
2071  {
2072    a[i] = deg(I[i],w);
2073  }
2074  ring FD = 0,(fd(1..n)),dp;
2075  def @@WFD = save + FD;
2076  setring @@WFD;
2077  poly f = imap(save,f);
2078  list RL = ringlist(basering);
2079  RL = RL[1..4];
2080  list L = RL[3];
2081  v = (1:(2*n)),((deg(f)+1):n);
2082  L = insert(L,list("a",v));
2083  RL[3] = L;
2084  def @WFD = ring(RL);
2085  setring @WFD;
2086  poly f = imap(save,f);
2087  matrix Drel[3*n][3*n];
2088  for (i=1; i<=n; i++)
2089  {
2090    Drel[i,i+n] = 1;     // [D,x]
2091    Drel[i,i+2*n] = f;   // [fD,x]
2092    for (j=1; j<=n; j++)
2093    {
2094      Drel[i+n,j+2*n] = -diff(f,var(i))*var(j+n);  // [fD,D]
2095      Drel[j+2*n,i+2*n] = diff(f,var(i))*var(j+2*n) - diff(f,var(j))*var(i+2*n);
2096      // [fD,fD]
2097    }
2098  }
2099  def WFD = nc_algebra(1,Drel);
2100  setring WFD;
2101  kill @WFD,@@WFD;
2102  ideal I = imap(save,I);
2103  poly f = imap(save,f);
2104  for (i=1; i<=size(I); i++)
2105  {
2106    I[i] = f^(a[i])*I[i];
2107  }
2108  ideal S;
2109  for (i=1; i<=n; i++)
2110  {
2111    S[size(S)+1] = var(i+2*n) - f*var(i+n);
2112  }
2113  S = slimgb(S);
2114  I = NF(I,S);
2115  if (select1(I,intvec((n+1)..2*n))[1] <> 0)
2116  {
2117    // should never get here
2118    ERROR("Something's wrong. Please inform the author.");
2119  }
2120  setring save;
2121  ideal mm = maxideal(1);
2122  poly s;
2123  for (i=1; i<=n; i++)
2124  {
2125    s = f*var(i+n);
2126    for (j=1; j<=size(F); j++)
2127    {
2128      s = s + k[j]*(f/F[j])*bracket(var(i+n),F[j]);
2129    }
2130    mm[i+2*n] = s;
2131  }
2132  map m = WFD,mm;
2133  ideal J = m(I);
2134  return(J);
2135}
2136example
2137{
2138  "EXAMPLE"; echo=2;
2139  ring r = 0,(x,y,Dx,Dy),dp;
2140  def W = Weyl();
2141  setring W;
2142  poly tx,ty = x*Dx, y*Dy;
2143  ideal I =      // Appel F1 with parameters (3,-1,1,1) is a solution
2144    tx*(tx+ty)-x*(tx+ty+3)*(tx-1),
2145    ty*(tx+ty)-y*(tx+ty+3)*(ty+1);
2146  kill tx,ty;
2147  poly f = x^3*y^2-x^2*y^3-x^3*y+x*y^3+x^2*y-x*y^2;
2148  ideal F = x-1,x,-x+y,y-1,y;
2149  intvec k = -1,-1,-1,-3,-1;
2150  ideal T = twistedIdeal(I,f,k,F);
2151  // TODO change the ordering of WFD
2152  // introduce new var f??
2153  //paper:
2154  poly fx = diff(f,x);
2155  poly fy = diff(f,y);
2156  poly fDx = f*Dx;
2157  poly fDy = f*Dy;
2158  poly fd(1) = fDx;
2159  poly fd(2) = fDy;
2160  ideal K=
2161    (x^2-x^3)*(fDx)^2+x*((1-3*x)*f-(1-x)*y*fy-(1-x)*x*fx)*(fDx)
2162    +x*(1-x)*y*(fDy)*(fDx)+x*y*f*(fDy)+3*x*f^2,
2163    (y^2-y^3)*(fDy)^2+y*((1-5*y)*f-(1-y)*x*fx-(1-y)*y*fy)*(fDy)
2164    +y*(1-y)*x*(fDx)*(fDy)-y*x*f*(fDx)-3*y*f^2;
2165}
2166
2167
2168proc ratSol (ideal I)
2169"
2170USAGE:    ratSol(I);  I ideal
2171ASSUME:   The basering is the n-th Weyl algebra W over a field of
2172          characteristic 0 and for all 1<=i<=n the identity
2173          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
2174          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
2175          differential operator belonging to x(i).
2176@*        Moreover, assume that I is holonomic.
2177RETURN:   module, a basis of the rational solutions to the given system.
2178          Note that each entry has two components, the first one standing for
2179          the enumerator, the second one for the denominator.
2180REMARKS:  Reference: (OTT), Algorithm 3.10
2181NOTE:     If printlevel=1, progress debug messages will be printed,
2182          if printlevel>=2, all the debug messages will be printed.
2184EXAMPLE:  example ratSol; shows examples
2185"
2186{
2187  dmodGeneralAssumptionCheck();
2188  if (!isHolonomic(I))
2189  {
2190    ERROR("Ideal is not holonomic.");
2191  }
2192  int ppl = printlevel - voice + 2;
2193  def save = basering;
2194  dbprint(ppl,"// computing singular locus...");
2195  ideal S = DsingularLocus(I);
2196  dbprint(ppl,"// ...done.");
2197  poly f = S[1];
2198  dbprint(ppl,"// considering poly " + string(f));
2199  int n = nvars(save) div 2;
2200  list RL = ringlist(save);
2201  RL = RL[1..4];
2202  list L = RL[2];
2203  L = list(L[1..n]);
2204  RL[2] = L;
2205  L = list();
2206  L[1] = list("dp",intvec(1:n));
2207  L[2] = list("C",intvec(0));
2208  RL[3] = L;
2209  def rr = ring(RL);
2210  setring rr;
2211  poly f = imap(save,f);
2212  ideal F = factorize(f,1); // not interested in multiplicities
2213  dbprint(ppl,"// with irreducible factors " + string(F));
2214  setring save;
2215  ideal F = imap(rr,F);
2216  kill rr,RL;
2217  int i;
2218  intvec k;
2219  ideal FF = 1,1;
2220  dbprint(ppl,"// computing b-functions of irreducible factors...");
2221  for (i=1; i<=size(F); i++)
2222  {
2223    dbprint(ppl,"//   considering " + string(F[i]) + "...");
2224    L = bfctBound(I,F[i]);
2225    if (size(L) == 3) // bfct is constant
2226    {
2227      dbprint(ppl,"//   ...got " + string(L[3]));
2228      if (L[3] == "1")
2229      {
2230        return(0); // TODO type // no rational solutions
2231      }
2232      else // should never get here
2233      {
2234        ERROR("Oops, something went wrong. Please inform the author.");
2235      }
2236    }
2237    else
2238    {
2239      dbprint(ppl,"//   ...got roots " + string(L[1]));
2240      dbprint(ppl,"//      with multiplicities " + string(L[2]));
2241      k[i] = -maxIntRoot(L)-1;
2242      if (k[i] < 0)
2243      {
2244        FF[2] = FF[2]*F[i]^(-k[i]);
2245      }
2246      else
2247      {
2248        FF[1] = FF[1]*F[i]^(k[i]);
2249      }
2250    }
2251  }
2252  vector v = FF[1]*gen(1) + FF[2]*gen(2);
2253  kill FF;
2254  dbprint(ppl,"// ...done");
2255  ideal twI = twistedIdeal(I,f,k,F);
2256  intvec w = -1:n;
2257  dbprint(ppl,"// computing polynomial solutions of twisted system...");
2258  if (isHolonomic(twI))
2259  {
2260    ideal P = polSol(twI,w);
2261  }
2262  else
2263  {
2264    ideal P = polSolFiniteRank(twI,w);
2265  }
2266  module M;
2267  vector vv;
2268  for (i=1; i<=ncols(P); i++)
2269  {
2270    vv = P[i]*gen(1) + 1*gen(2);
2271    M[i] = multRat(v,vv);
2272  }
2273  dbprint(ppl,"// ...done");
2274  return (M);
2275}
2276example
2277{
2278  "EXAMPLE"; echo=2;
2279  ring r = 0,(x,y,Dx,Dy),dp;
2280  def W = Weyl();
2281  setring W;
2282  poly tx,ty = x*Dx, y*Dy;
2283  ideal I =      // Appel F1 with parameters (3,-1,1,1) is a solution
2284    tx*(tx+ty)-x*(tx+ty+3)*(tx-1),
2285    ty*(tx+ty)-y*(tx+ty+3)*(ty+1);
2286  module M = ratSol(I);
2287  // We obtain a basis of the rational solutions to I represented by a
2288  // module / matrix with two rows.
2289  // Each column of the matrix represents a rational function, where
2290  // the first row correspond to the enumerator and the second row to
2291  // the denominator.
2292  print(M);
2293}
2294
2295
2296proc bfctBound (ideal I, poly f, list #)
2297"
2298USAGE:    bfctBound (I,f[,primdec]);  I ideal, f poly, primdec optional string
2299ASSUME:   The basering is the n-th Weyl algebra W over a field of
2300          characteristic 0 and for all 1<=i<=n the identity
2301          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
2302          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
2303          differential operator belonging to x(i).
2304@*        Moreover, assume that I is holonomic.
2305RETURN:   list of roots (of type ideal) and multiplicities (of type intvec) of
2306          a multiple of the b-function for f^s*u at a generic root of f.
2307          Here, u stands for [1] in D/I.
2308REMARKS:  Reference: (OTT), Algorithm 3.4
2309NOTE:     This procedure requires to compute a primary decomposition in a
2310          commutative ring. The optional string primdec can be used to specify
2311          the algorithm to do so. It may either be 'GTZ' (Gianni, Trager,
2312          Zacharias) or 'SY' (Shimoyama, Yokoyama). By default, 'GTZ' is used.
2313@*        If printlevel=1, progress debug messages will be printed,
2314          if printlevel>=2, all the debug messages will be printed.
2316EXAMPLE:  example bfctBound; shows examples
2317"
2318{
2319  dmodGeneralAssumptionCheck();
2320  finKx(f);
2321  if (!isHolonomic(I))
2322  {
2323    ERROR("Ideal is not holonomic.");
2324  }
2325  int ppl = printlevel - voice + 2;
2326  string primdec = "GTZ";
2327  if (size(#)>1)
2328  {
2329    if (typeof(#[1])=="string")
2330    {
2331      if ( (#[1]=="SY") || (#[1]=="sy") || (#[1]=="Sy") )
2332      {
2333        primdec = "SY";
2334      }
2335      else
2336      {
2337        if ( (#[1]<>"GTZ") && (#[1]<>"gtz") && (#[1]<>"Gtz") )
2338        {
2339          print("// Warning: optional string may either be 'GTZ' or 'SY',");
2340          print("//          proceeding with 'GTZ'.");
2341          primdec = "GTZ";
2342        }
2343      }
2344    }
2345  }
2346  def save = basering;
2347  int n = nvars(save) div 2;
2348  // step 1
2349  ideal mm = maxideal(1);
2350  def Wt = extendWeyl(safeVarName("t"));
2351  setring Wt;
2352  poly f = imap(save,f);
2353  ideal mm = imap(save,mm);
2354  int i;
2355  for (i=1; i<=n; i++)
2356  {
2357    mm[i+n] = var(i+n+2) + bracket(var(i+n+2),f)*var(n+2);
2358  }
2359  map m = save,mm;
2360  ideal I = m(I);
2361  I = I, var(1)-f;
2362  // step 2
2363  intvec w = 1,(0:n);
2364  dbprint(ppl  ,"// Computing initial ideal...");
2365  I = initialIdealW(I,-w,w);
2366  dbprint(ppl  ,"// ...done.");
2367  dbprint(ppl-1,"// " + string(I));
2368  // step 3: rewrite I using Euler operator t*Dt
2369  list RL = ringlist(Wt);
2370  RL = RL[1..4];
2371  list L = RL[2] + list(safeVarName("s")); // s=t*Dt
2372  RL[2] = L;
2373  L = list();
2374  L[1] = list("dp",intvec(1:(2*n+2)));
2375  L[2] = list("dp",intvec(1));
2376  L[3] = list("C",intvec(0));
2377  RL[3] = L;
2378  def @Wts = ring(RL);
2379  kill L,RL;
2380  setring @Wts;
2381  matrix relD[2*n+3][2*n+3];
2382  relD[1,2*n+3] = var(1);
2383  relD[n+2,2*n+3] = -var(n+2);
2384  for (i=1; i<=n+1; i++)
2385  {
2386    relD[i,n+i+1] = 1;
2387  }
2388  def Wts = nc_algebra(1,relD);
2389  setring Wts;
2390  ideal I = imap(Wt,I);
2391  kill Wt,@Wts;
2392  ideal S = var(1)*var(n+2)-var(2*n+3);
2393  attrib(S,"isSB",1);
2394  dbprint(ppl  ,"// Computing Euler representation...");
2395  // I = NF(I,S);
2396  int d;
2397  intvec ww = 0:(2*2+2); ww[1] = -1; ww[n+2] = 1;
2398  for (i=1; i<=size(I); i++)
2399  {
2400    d = deg(I[i],ww);
2401    if (d>0)
2402    {
2403      I[i] = var(1)^d*I[i];
2404    }
2405    if (d<0)
2406    {
2407      d = -d;
2408      I[i] = var(n+2)^d*I[i];
2409    }
2410  }
2411  I = NF(I,S); // now there are no t,Dt in I, only s
2412  dbprint(ppl  ,"// ...done.");
2413  I = subst(I,var(2*n+3),-var(2*n+3)-1);
2414  ring Ks = 0,s,dp;
2415  def Ws = save + Ks;
2416  setring Ws;
2417  ideal I = imap(Wts,I);
2418  kill Wts;
2419  poly DD = 1;
2420  for (i=1; i<=n; i++)
2421  {
2422    DD = DD * var(n+i);
2423  }
2424  dbprint(ppl  ,"// Eliminating differential operators...");
2425  ideal J = eliminate(I,DD); // J subset K[x,s]
2426  dbprint(ppl  ,"// ...done.");
2427  dbprint(ppl-1,"// " + string(J));
2428  list RL = ringlist(Ws);
2429  RL = RL[1..4];
2430  list L = RL[2];
2431  L = list(L[1..n]) + list(L[2*n+1]);
2432  RL[2] = L;
2433  L = list();
2434  L[1] = list("dp",intvec(1:(n+1)));
2435  L[2] = list("C",intvec(0));
2436  RL[3] = L;
2437  def Kxs = ring(RL);
2438  setring Kxs;
2439  ideal J = imap(Ws,J);
2440  dbprint(ppl  ,"// Computing primary decomposition with engine "
2441          + primdec + "...");
2442  if (primdec == "GTZ")
2443  {
2444    list P = primdecGTZ(J);
2445  }
2446  else
2447  {
2448    list P = primdecSY(J);
2449  }
2450  dbprint(ppl  ,"// ...done.");
2451  dbprint(ppl-1,"// " + string(P));
2453  poly f = imap(save,f);
2454  vector v;
2455  for (i=1; i<=size(P); i++)
2456  {
2457    dbprint(ppl  ,"// Considering primary component " + string(i)
2458            + " of " + string(size(P)) + "...");
2459    dbprint(ppl  ,"//   Intersecting with K[x] and computing radical...");
2460    GP = std(P[i][1]);
2461    Qix = eliminate(GP,var(n+1)); // subset K[x]
2464    dbprint(ppl  ,"//   ...done.");
2467    {
2468      dbprint(ppl  ,"//   Intersecting with K[s]...");
2469      v = pIntersect(var(n+1),GP);
2470      B[size(B)+1] = vec2poly(v,n+1);
2471      dbprint(ppl  ,"//   ...done.");
2472      dbprint(ppl-1,"// " + string(B[size(B)]));
2473    }
2474    dbprint(ppl  ,"// ...done.");
2475  }
2476  f = lcm(B); // =lcm(B[1],...,B[size(B)])
2477  list bb = bFactor(f);
2478  setring save;
2479  list bb = imap(Kxs,bb);
2480  return(bb);
2481}
2482example
2483{
2484  "EXAMPLE"; echo=2;
2485  ring r = 0,(x,y,Dx,Dy),dp;
2486  def W = Weyl();
2487  setring W;
2488  poly tx,ty = x*Dx, y*Dy;
2489  ideal I =      // Appel F1 with parameters (2,-3,-2,5)
2490    tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3),
2491    ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2),
2492    (x-y)*Dx*Dy+2*Dx-3*Dy;
2493  kill tx,ty;
2494  poly f = x-1;
2495  bfctBound(I,f);
2496}
2497
2498
2499//TODO check f/g or g/f, check Weyl closure of result
2500proc annRatSyz (poly f, poly g, list #)
2501"
2502USAGE:    annRatSyz(f,g[,db,eng]);  f, g polynomials, db,eng optional integers
2503ASSUME:   The Basering is commutative and over a field of characteristic 0.
2504RETURN:   ring (a Weyl algebra) containing an ideal 'LD', which is (part of)
2505          the annihilator of the rational function g/f in the corresponding
2506          Weyl algebra
2507REMARKS:  This procedure uses the computation of certain syzygies.
2508          One can obtain the full annihilator by computing the Weyl closure of
2509          the ideal LD.
2510NOTE:     Activate the output ring with the @code{setring} command.
2511          In the output ring, the ideal 'LD' (in Groebner basis) is (part of)
2512          the annihilator of g/f.
2513@*        If db>0 is given, operators of order up to db are considered,
2514          otherwise, and by default, a minimal holonomic solution is computed.
2515@*        If eng<>0, @code{std} is used for Groebner basis computations,
2516          otherwise, and by default, @code{slimgb} is used.
2517@*        If printlevel =1, progress debug messages will be printed,
2518          if printlevel>=2, all the debug messages will be printed.
2520EXAMPLE:  example annRatSyz; shows examples
2521"
2522{
2523  // check assumptions
2524  if (!isCommutative())
2525  {
2526    ERROR("Basering must be commutative.");
2527  }
2528  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
2529  {
2530    ERROR("Basering is inappropriate: characteristic>0 or qring present.");
2531  }
2532  if (g == 0)
2533  {
2534    ERROR("Second polynomial must not be zero.");
2535  }
2536  int db,eng;
2537  if (size(#)>0)
2538  {
2539    if (typeof(#[1]) == "int")
2540    {
2541      db = int(#[1]);
2542    }
2543    if (size(#)>1)
2544    {
2545      if (typeof(#[2]) == "int")
2546      {
2547        eng = int(#[1]);
2548      }
2549    }
2550  }
2551  int ppl = printlevel - voice + 2;
2552  vector I = f*gen(1)+g*gen(2);
2553  checkRatInput(I);
2554  int i,j;
2555  def R = basering;
2556  int n = nvars(R);
2557  list RL = ringlist(R);
2558  RL = RL[1..4];
2559  list Ltmp = RL[2];
2560  for (i=1; i<=n; i++)
2561  {
2562    Ltmp[i+n] = safeVarName("D" + Ltmp[i]);
2563  }
2564  RL[2] = Ltmp;
2565  Ltmp = list();
2566  Ltmp[1] = list("dp",intvec(1:2*n));
2567  Ltmp[2] = list("C",intvec(0));
2568  RL[3] = Ltmp;
2569  kill Ltmp;
2570  def @D = ring(RL);
2571  setring @D;
2572  def D = Weyl();
2573  setring D;
2574  ideal DD = 1;
2575  ideal Dcd,Dnd,LD,tmp;
2576  Dnd = 1;
2577  module DS;
2578  poly DJ;
2579  kill @D;
2580  setring R;
2581  module Rnd,Rcd;
2582  Rnd[1] = I;
2583  vector RJ;
2584  ideal L = I[1];
2585  module RS;
2586  poly p,pnew;
2587  pnew = I[2];
2588  int k,c;
2589  while(1)
2590  {
2591    k++;
2592    setring R;
2593    dbprint(ppl,"// Testing order: " + string(k));
2594    Rcd = Rnd;
2595    Rnd = 0;
2596    setring D;
2597    Dcd = Dnd;
2598    Dnd = 0;
2599    dbprint(ppl-1,"// Current members of the annihilator: " + string(LD));
2600    setring R;
2601    c = size(Rcd);
2602    p = pnew;
2603    for (i=1; i<=c; i++)
2604    {
2605      for (j=1; j<=n; j++)
2606      {
2607        RJ = diffRat(Rcd[i],j);
2608        setring D;
2609        DJ = Dcd[i]*var(n+j);
2610        tmp = Dnd,DJ;
2611        if (size(Dnd) <> size(simplify(tmp,4))) // new element
2612        {
2613          Dnd[size(Dnd)+1] = DJ;
2614          setring R;
2615          Rnd[size(Rnd)+1] = RJ;
2616          pnew = lcm(pnew,RJ[2]);
2617        }
2618        else // already have DJ in Dnd
2619        {
2620          setring R;
2621        }
2622      }
2623    }
2624    p = pnew/p;
2625    for (i=1; i<=size(L); i++)
2626    {
2627      L[i] = p*L[i];
2628    }
2629    for (i=1; i<=size(Rnd); i++)
2630    {
2631      L[size(L)+1] = pnew/Rnd[i][2]*Rnd[i][1];
2632    }
2633    RS = syz(L);
2634    setring D;
2635    DD = DD,Dnd;
2636    setring R;
2637    if (RS <> 0)
2638    {
2639      setring D;
2640      DS = imap(R,RS);
2641      LD = ideal(transpose(DS)*transpose(DD));
2642    }
2643    else
2644    {
2645      setring D;
2646    }
2647    LD = engine(LD,eng);
2648    // test if we're done
2649    if (db<=0)
2650    {
2651      if (isHolonomic(LD)) { break; }
2652    }
2653    else
2654    {
2655      if (k==db) { break; }
2656    }
2657  }
2658  export(LD);
2659  setring R;
2660  return(D);
2661}
2662example
2663{
2664  "EXAMPLE:"; echo = 2;
2665  // printlevel = 3;
2666  ring r = 0,(x,y),dp;
2667  poly f = 2*x*y; poly g = x^2 - y^3;
2668  def A = annRatSyz(f,g);   // compute a holonomic solution
2669  setring A; A;
2670  LD;
2671  setring r;
2672  def B = annRatSyz(f,g,5); // compute a solution up to degree 5
2673  setring B;
2674  LD; // this is the full annihilator as we will check below
2675  setring r;
2676  def C = annRat(f,g); setring C;
2677  LD; // the full annihilator
2678  ideal BLD = imap(B,LD);
2679  NF(LD,std(BLD));
2680}
2681
Note: See TracBrowser for help on using the repository browser.