source: git/Singular/LIB/dmodloc.lib @ 3686937

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