source: git/dmodloc.lib @ 3dcd82

//
version="$Id$";
category="Noncommutative";
info="
LIBRARY: dmodloc.lib     Localization of algebraic D-modules and applications
AUTHORS: Daniel Andres,  daniel.andres@math.rwth-aachen.de

Support: DFG Graduiertenkolleg 1632 'Experimentelle und konstruktive Algebra'


OVERVIEW:
Let I be a left ideal I in the n-th polynomial Weyl algebra D=K[x]<d> and
let f be a polynomial in K[x].

If D/I is a holonomic module over D, it is known that the localization of D/I
at f is also holonomic. The procedure Dlocalization computes an ideal J in D
such that this localization is isomorphic to D/J as D-modules.

If one regards I as an ideal in the rational Weyl algebra as above, K(x)<d>*I,
and intersects with K[x]<d>, the result is called the Weyl closure of I.
The procedures WeylClosure (if I has finite holonomic rank) and WeylClosure1
(if I is in the first Weyl algebra) can be used for computations.

As an application of the Weyl closure, the procedure annRatSyz computes a
holonomic part of the annihilator of a rational function by computing certain
syzygies. The full annihilator can be obtained by taking the Weyl closure of
the result.

If one regards the left ideal I as system of linear PDEs, one can find its
polynomial solutions with polSol (if I is holonomic) or polSolFiniteRank
(if I is of finite holonomic rank). Rational solutions can be obtained
with ratSol.

The procedure bfctBound computes a possible multiple of the b-function for
f^s*u at a generic root of f. Here, u stands for [1] in D/I.

This library also offers the procedures holonomicRank and DsingularLocus to
compute the holonomic rank and the singular locus of the D-module D/I.


References:
@* (OT)  Oaku, Takayama: 'Algorithms for D-modules',
         Journal of Pure and Applied Algebra, 1998
@* (OTT) Oaku, Takayama, Tsai: 'Polynomial and rational solutions of holonomic 
         systems', Journal of Pure and Applied Algebra, 2001
@* (OTW) Oaku, Takayama, Walther: 'A Localization Algorithm for D-modules',
         Journal of Symbolic Computation, 2000
@* (Tsa) Tsai: 'Algorithms for algebraic analysis', PhD thesis, 2000


PROCEDURES:

holonomicRank(I);  computes the holonomic rank of I 
DsingularLocus(I); computes the singular locus of a D-module
Dlocalization(I,f[,k,e]);  computes the localization of a D-module
WeylClosure1(L);   computes the Weyl closure of operator in first Weyl algebra
WeylClosure(I);    computes the Weyl closure of an ideal in the Weyl algebra
polSol(I[,w,m]);   computes basis of polynomial solutions to the given system
polSolFiniteRank(I[,w]); computes basis of polynomial solutions to given system
ratSol(I);         computes basis of rational solutions to the given system
bfctBound(I,f[,primdec]); computes multiple of b-function for f^s*u
annRatSyz(f,g[,db,eng]);  computes part of annihilator of rational function g/f

dmodGeneralAssumptionCheck();   check general assumptions
safeVarName(s);    finds an untaken name to use for a new variable
extendWeyl(S);     extends basering (Weyl algebra) by given vars
polyVars(f,v);     checks whether f contains only variables indexed by v
monomialInIdeal(I);    computes all monomials appearing in generators of ideal
vars2pars(v);      converts variables specified by v into parameters
killTerms(p,k);    delete all terms of p containing var(k)
minIntRoot2(L);    finds minimal integer root in a list of roots
maxIntRoot(L);     finds maximal integer root in a list of roots
dmodAction(id,f[,v]);  computes the natural action of a D-module on K[x]
dmodActionRat(id,w);   computes the natural action of a D-module on K(x)
simplifyRat(v);    simplifies rational function
addRat(v,w);       adds rational functions
multRat(v,w);      multiplies rational functions
diffRat(v,j);      derives rational function
commRing();        deletes non-commutative relations from ring
rightNFWeyl(id,k); computes right NF wrt. right ideal (var(k)) in Weyl algebra


KEYWORDS: 
D-module; 
holonomic rank; singular locus of D-module; 
D-localization; localization of D-module; characteristic variety;
Weyl closure;
polynomial solutions; rational solutions;
annihilator of rational function;

SEE ALSO: bfun_lib, dmod_lib, dmodapp_lib, dmodvar_lib, gmssing_lib
";


/*
CHANGELOG:
12.11.12: bugfixes, updated docu
*/


LIB "bfun.lib";    // for pIntersect etc
LIB "dmodapp.lib"; // for GBWeight, charVariety etc
LIB "nctools.lib"; // for Weyl, isWeyl etc
// TODO uncomment this once chern.lib is ready
// LIB "chern.lib";   // for orderedPartition


// testing for consistency of the library /////////////////////////////////////
proc testdmodloc()
{
  example dmodGeneralAssumptionCheck;
  example safeVarName;
  example extendWeyl;
  example polyVars;
  example monomialInIdeal;
  example vas2pars;
  example killTerms;
  example minIntRoot2;
  example maxIntRoot;
  example dmodAction;
  example dmodActionRat;
  example simplifyRat;
  example addRat;
  example multRat;
  example diffRat;
  example commRing;
  example holonomicRank;
  example DsingularLocus;
  example rightNFWeyl;
  example Dlocalization;
  example WeylClosure1;
  example WeylClosure;
  example polSol;
  example polSolFiniteRank;
  example ratSol;
  example bfctBound;
  example annRatSyz;
}


// tools //////////////////////////////////////////////////////////////////////

proc dmodGeneralAssumptionCheck ()
"
USAGE:    dmodGeneralAssumptionCheck();
RETURN:   nothing, but checks general assumptions on the basering
NOTE:     This procedure checks the following conditions on the basering R
          and prints an error message if any of them is violated:
@*         - R is the n-th Weyl algebra over a field of characteristic 0,
@*         - R is not a qring,
@*         - for all 1<=i<=n the identity var(i+n)*var(i)=var(i)*var(i+1)+1
             holds, i.e. the sequence of variables is given by
             x(1),...,x(n),D(1),...,D(n), where D(i) is the differential 
             operator belonging to x(i).
EXAMPLE:  example dmodGeneralAssumptionCheck; shows examples
"
{
  // char K <> 0, qring
  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
  {
    ERROR("Basering is inappropriate: characteristic>0 or qring present");
  }
  // no Weyl algebra
  if (isWeyl() == 0)
  {
    ERROR("Basering is not a Weyl algebra");
  }
  // wrong sequence of vars
  int i,n;
  n = nvars(basering) div 2;
  for (i=1; i<=n; i++)
  {
    if (bracket(var(i+n),var(i))<>1)
    {
      ERROR(string(var(i+n))+" is not a differential operator for " +string(var(i)));
    }
  }
  return();
}
example
{
  "EXAMPLE"; echo=2;
  ring r = 0,(x,D),dp;
  dmodGeneralAssumptionCheck(); // prints error message
  def W = Weyl();
  setring W;
  dmodGeneralAssumptionCheck(); // returns nothing
}


proc safeVarName (string s)
"
USAGE:    safeVarName(s);  s string
RETURN:   string, returns s if s is not the name of a par/var of basering
          and '@' + s otherwise
EXAMPLE:  example safeVarName; shows examples
"
{
  string S = "," + charstr(basering) + "," + varstr(basering) + ",";
  s = "," + s + ",";
  while (find(S,s) <> 0)
  {
    s[1] = "@";
    s = "," + s;
  }
  s = s[2..size(s)-1];
  return(s);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = (0,a),(w,@w,x,y),dp;
  safeVarName("a");
  safeVarName("x");
  safeVarName("z");
  safeVarName("w");
}


proc extendWeyl (def newVars)
"
USAGE:    extendWeyl(S);  S string or list of strings
ASSUME:   The basering is the n-th Weyl algebra over a field of
          characteristic 0 and for all 1<=i<=n the identity
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
          is the differential operator belonging to x(i).
RETURN:   ring, Weyl algebra extended by vars given by S
EXAMPLE:  example extendWeyl; shows examples
"
{
  dmodGeneralAssumptionCheck();
  int i,s;
  string inpt = typeof(newVars);
  list L;
  if (inpt=="string")
  {
    s = 1;
    L = newVars;
  }
  else
  {
    if (inpt=="list")
    {
      s = size(newVars);
      if (s<1)
      {
        ERROR("No new variables specified.");
      }
      for (i=1; i<=s; i++)
      {
        if (typeof(newVars[i]) <> "string")
        {
          ERROR("Entries of input list must be of type string.");
        }
      }
      L = newVars;
    }
    else
    {
      ERROR("Expected string or list of strings as input.");
    }
  }
  def save = basering;
  int n = nvars(save) div 2;
  list RL = ringlist(save);
  RL = RL[1..4];
  list Ltemp = L;
  for (i=s; i>0; i--)
  {
    Ltemp[n+s+i] = "D" + newVars[i];
  }
  for (i=n; i>0; i--)
  {
    Ltemp[s+i]     = RL[2][i];
    Ltemp[n+2*s+i] = RL[2][n+i];
  }
  RL[2] = Ltemp;
  Ltemp = list();
  Ltemp[1] = list("dp",intvec(1:(2*n+2*s)));
  Ltemp[2] = list("C",intvec(0));
  RL[3] = Ltemp;
  kill Ltemp;
  def @Dv = ring(RL);
  setring @Dv;
  def Dv = Weyl();
  setring save;
  return(Dv);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring @D2 = 0,(x,y,Dx,Dy),dp;
  def D2 = Weyl();
  setring D2;
  def D3 = extendWeyl("t");
  setring D3; D3;
  list L = "u","v";
  def D5 = extendWeyl(L);
  setring D5;
  D5;
}


proc polyVars (poly f, intvec v)
"
USAGE:    polyVars(f,v);  f poly, v intvec
RETURN:   int, 1 if f contains only variables indexed by v, 0 otherwise
EXAMPLE:  example polyVars; shows examples
"
{
  ideal varsf = variables(f); // vars contained in f
  ideal V;
  int i;
  int n = nvars(basering);
  for (i=1; i<=nrows(v); i++)
  {
    if ( (v[i]<0) || (v[i]>n) )
    {
      ERROR("var(" + string(v[i]) + ") out of range");
    }
    V[i] = var(v[i]);
  }
  attrib(V,"isSB",1);
  ideal N = NF(varsf,V);
  N = simplify(N,2);
  if (N[1]==0)
  {
    return(1);
  }
  else
  {
    return(0);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  poly f = y^2+zy;
  intvec v = 1,2;
  polyVars(f,v); // does f depend only on x,y?
  v = 2,3;
  polyVars(f,v); // does f depend only on y,z?
}


proc monomialInIdeal (ideal I)
"
USAGE:    monomialInIdeal(I);  I ideal
RETURN:   ideal consisting of all monomials appearing in generators of ideal
EXAMLPE:  example monomialInIdeal; shows examples
"
{
  // returns ideal consisting of all monomials appearing in generators of ideal
  I = simplify(I,2+8);
  int i;
  poly p;
  ideal M;
  for (i=1; i<=size(I); i++)
  {
    p = I[i];
    while (p<>0)
    {
      M[size(M)+1] = leadmonom(p);
      p = p - lead(p);
    }
  }
  M = simplify(M,4+2);
  return(M);
}
example
{
  "EXAMPLE"; echo=2;
  ring r = 0,(x,y),dp;
  ideal I = x2+5x3y7, x-x2-6xy;
  monomialInIdeal(I);
}


proc vars2pars (intvec v)
"
USAGE:    vars2pars(v);  v intvec
ASSUME:   The basering is commutative.
RETURN:   ring with variables specified by v converted into parameters
EXAMPLE:  example vars2pars; shows examples
"
{
  if (isCommutative() == 0)
  {
    ERROR("The basering must be commutative.");
  }
  v = sortIntvec(v)[1];
  int sv = size(v);
  if ( (v[1]<1) || (v[sv]<1) )
  {
    ERROR("Expected entries of intvec in the range 1.."+string(n));
  }
  def save = basering;
  int i,j,n;
  n = nvars(save);
  list RL = ringlist(save);
  list Lp,Lv,L1;
  if (typeof(RL[1]) == "list")
  {
    L1 = RL[1];
    Lp = L1[2];
  }
  else
  {
    L1[1] = RL[1];
    L1[4] = ideal(0);
  }
  j = sv;
  for (i=1; i<=n; i++)
  {
    if (j>0)
    {
      if (v[j]==i)
      {
        Lp[size(Lp)+1] = string(var(i));
        j--;
      }
      else
      {
        Lv[size(Lv)+1] = string(var(i));
      }
    }
    else
    {
      Lv[size(Lv)+1] = string(var(i));
    }
  }
  RL[2] = Lv;
  L1[2] = Lp;
  L1[3] = list(list("lp",intvec(1:size(Lp))));
  RL[1] = L1;
  L1 = list();
  L1[1] = list("dp",intvec(1:sv));
  L1[2] = list("C",intvec(0));
  RL[3] = L1;
//   RL;
  def R = ring(RL);
  return(R);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z,a,b,c),dp;
  intvec v = 4,5,6;
  def R = vars2pars(v);
  setring R;
  R;
  v = 1,2;
  def RR = vars2pars(v);
  setring RR;
  RR;
}


proc killTerms (poly p, int k)
"
USAGE:    killTerms(p,k);  p poly, k int
RETURN:   poly, p minus all terms containing the k-th variable
NOTE:     No Groebner basis computations are used.
EXAMPLE:  example killTerms; shows examples
"
{
  if ( (k<0) || (k>nvars(basering)) )
  {
    ERROR("Given int is out of range");
  }
  int j;
  poly q;
  intvec v;
  for (j=1; j<=size(p); j++)
  {
    v = leadexp(p[j]);
    if (v[k]==0)
    {
      q = q + p[j];
    }
  }
  return(q);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,z),dp;
  poly p = 8xyz + 7xy + 6xz + 5yz + 4x + 3y + 2z + 1;
  killTerms(p,1);
  killTerms(p,2);
  killTerms(p,3);
}


static proc minMaxIntRoot (list L, string minmax)
{
  int win;
  if (size(L)>1)
  {
    if ( (typeof(L[1])<>"ideal") || (typeof(L[2])<>"intvec") )
    {
      win = 1;
    }
  }
  else
  {
    win = 1;
  }
  if (win)
  {
    ERROR("Expected list in the format of bFactor");
  }
  if (size(L)>2)
  {
    if ( (L[3]=="1") || (L[3]=="0") )
    {
      print("// Warning: Constant poly. Returning 0.");
      return(int(0));
    }
  }
  ideal I = L[1];
  int i,k,b;
  if (minmax=="min")
  {
    i = ncols(I);
    k = -1;
    b = 0;
  }
  else // minmax=="max"
  {
    i = 1;
    k = 1;
    b = ncols(I);
  }
  for (; k*i<k*b; i=i+k)
  {
    if (isInt(leadcoef(I[i])))
    {
      return(int(leadcoef(I[i])));
    }
  }
  print("// Warning: No integer root found. Returning 0.");
  return(int(0));
}


//TODO rename? minIntRoot is name of proc in dmod.lib
proc minIntRoot2 (list L)
"
USAGE:    minIntRoot2(L);  L list
ASSUME:   L is the output of bFactor.
RETURN:   int, the minimal integer root in a list of roots
SEE ALSO: minIntRoot, maxIntRoot, bFactor
EXAMPLE:  example minIntRoot2; shows examples
"
{
  return(minMaxIntRoot(L,"min"));
}
example
{
  "EXAMPLE"; echo=2;
  ring r = 0,x,dp;
  poly f = x*(x+1)*(x-2)*(x-5/2)*(x+5/2);
  list L = bFactor(f);
  minIntRoot2(L);
}


proc maxIntRoot (list L)
"
USAGE:    maxIntRoot(L);  L list
ASSUME:   L is the output of bFactor.
RETURN:   int, the maximal integer root in a list of roots
SEE ALSO: minIntRoot2, bFactor
EXAMPLE:  example maxIntRoot; shows examples
"
{
  return(minMaxIntRoot(L,"max"));
}
example
{
  "EXAMPLE"; echo=2;
  ring r = 0,x,dp;
  poly f = x*(x+1)*(x-2)*(x-5/2)*(x+5/2);
  list L = bFactor(f);
  maxIntRoot(L);
}


proc dmodAction (def id, poly f, list #)
"
USAGE:    dmodAction(id,f[,v]);  id ideal or poly, f poly, v optional intvec
ASSUME:   If v is not given, the basering is the n-th Weyl algebra W over a
          field of characteristic 0 and for all 1<=i<=n the identity 
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of 
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
          differential operator belonging to x(i).
          Otherwise, v is assumed to specify positions of variables, which form
          a Weyl algebra as a subalgebra of the basering:
          If size(v) equals 2*n, then bracket(var(v[i]),var(v[j])) must equal
          1 if and only if j equals i+n, and 0 otherwise,  for all 1<=i,j<=n.
@*        Further, assume that f does not contain any D(i).
RETURN:   same type as id, the result of the natural D-module action of id on f
NOTE:     The assumptions made are not checked.
EXAMPLE:  example dmodAction; shows examples
"
{
  string inp1 = typeof(id);
  if ((inp1<>"poly") && (inp1<>"ideal"))
  {
    ERROR("Expected first argument to be poly or ideal but received "+inp1);
  }
  intvec posXD = 1..nvars(basering);
  if (size(#)>0)
  {
    if (typeof(#[1])=="intvec")
    {
      posXD = #[1];
    }
  }
  if ((size(posXD) mod 2)<>0)
  {
    ERROR("Even number of variables expected.")
  }
  int n = (size(posXD)) div 2;
  int i,j,k,l;
  ideal resI = id;
  int sid = ncols(resI);
  intvec v;
  poly P,h;
  for (l=1; l<=sid; l++)
  {
    P = resI[l];
    resI[l] = 0;
    for (i=1; i<=size(P); i++)
    {
      v = leadexp(P[i]);
      h = f;
      for (j=1; j<=n; j++)
      {
        for (k=1; k<=v[posXD[j+n]]; k++) // action of Dx
        {
          h = diff(h,var(posXD[j]));
        }
        h = h*var(posXD[j])^v[posXD[j]]; // action of x
      }
      h = leadcoef(P[i])*h;
      resI[l] = resI[l] + h;
    }
  }
  if (inp1 == "ideal")
  {
    return(resI);
  }
  else
  {
    return(resI[1]);
  }
}
example
{
  ring r = 0,(x,y,z),dp;
  poly f = x^2*z - y^3;
  def A = annPoly(f);
  setring A;
  poly f = imap(r,f);
  dmodAction(LD,f);
  poly P = y*Dy+3*z*Dz-3;
  dmodAction(P,f);
  dmodAction(P[1],f);
}


static proc checkRatInput (vector I)
{
  // check for valid input
  int wrginpt;
  if (nrows(I)<>2)
  {
    wrginpt = 1;
  }
  else
  {
    if (I[2] == 0)
    {
      wrginpt = 1;
    }
  }
  if (wrginpt)
  {
    ERROR("Vector must consist of exactly two components, second one not 0");
  }
  return();
}


proc dmodActionRat(def id, vector w)
"
USAGE:    dmodActionRat(id,w);  id ideal or poly, f vector
ASSUME:   The basering is the n-th Weyl algebra W over a field of
          characteristic 0 and for all 1<=i<=n the identity 
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of 
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
          differential operator belonging to x(i).
@*        Further, assume that w has exactly two components, second one not 0,
          and that w does not contain any D(i).
RETURN:   same type as id, the result of the natural D-module action of id on
          the rational function w[1]/w[2]
EXAMPLE:  example dmodActionRat; shows examples
"
{
  string inp1 = typeof(id);
  if ( (inp1<>"poly") && (inp1<>"ideal") )
  {
    ERROR("Expected first argument to be poly or ideal but received " + inp1);
  }
  checkRatInput(w);
  poly f = w[1];
  finKx(f);
  f = w[2];
  finKx(f);
  def save = basering;
  def r = commRing();
  setring r;
  ideal I = imap(save,id);
  vector w = imap(save,w);
  int i,j,k,l;
  int n = nvars(basering) div 2;
  int sid = ncols(I);
  intvec v;
  poly P;
  vector h,resT;
  module resL;
  for (l=1; l<=sid; l++)
  {
    P = I[l];
    resT = [0,1];
    for (i=1; i<=size(P); i++)
    {
      v = leadexp(P[i]);
      h = w;
      for (j=1; j<=n; j++)
      {
        for (k=1; k<=v[j+n]; k++) // action of Dx
        {
          h = diffRat(h,j);
        }
        h = h + h[1]*(var(j)^v[j]-1)*gen(1);  // action of x
      }
      h = h + (leadcoef(P[i])-1)*h[1]*gen(1);
      resT = addRat(resT,h);
    }
    resL[l] = resT;
  }
  setring save;
  module resL = imap(r,resL);
  return(resL);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  poly f = 2*x;  poly g = y;
  def A = annRat(f,g); setring A;
  poly f = imap(r,f); poly g = imap(r,g);
  vector v = [f,g]; // represents f/g
  // x and y act by multiplication
  dmodActionRat(x,v);
  dmodActionRat(y,v);
  // Dx and Dy act by partial derivation
  dmodActionRat(Dx,v);
  dmodActionRat(Dy,v);
  dmodActionRat(x*Dx+y*Dy,v);
  setring r;
  f = 2*x*y; g = x^2 - y^3;
  def B = annRat(f,g); setring B;
  poly f = imap(r,f); poly g = imap(r,g);
  vector v = [f,g];
  dmodActionRat(LD,v); // hence LD is indeed the annihilator of f/g
}


static proc arithmeticRat (vector I, vector J, string op, list #)
{
  // op = "+": return I+J
  // op = "*": return I*J
  // op = "s": return simplified I
  // op = "d": return diff(I,var(#[1]))
  int isComm = isCommutative();
  if (!isComm)
  {
    def save = basering;
    def r = commRing();
    setring r;
    ideal m = maxideal(1);
    map f = save,m;
    vector I = f(I);
    vector J = f(J);
  }
  vector K;
  poly p;
  if (op == "s")
  {
    p = gcd(I[1],I[2]);
    K = (I[1]/p)*gen(1) + (I[2]/p)*gen(2);
  }
  else
  {
    if (op == "+")
    {
      I = arithmeticRat(I,vector(0),"s");
      J = arithmeticRat(J,vector(0),"s");
      p = lcm(I[2],J[2]);
      K = (I[1]*p/I[2] + J[1]*p/J[2])*gen(1) + p*gen(2);
    }
    else
    {
      if (op == "*")
      {
        K = (I[1]*J[1])*gen(1) + (I[2]*J[2])*gen(2);
      }
      else
      {
        if (op == "d")
        {
          int j = #[1];
          K = (diff(I[1],var(j))*I[2] - I[1]*diff(I[2],var(j)))*gen(1)+ (I[2]^2)*gen(2);
        }
      }
    }
    K = arithmeticRat(K,vector(0),"s");
  }
  if (!isComm)
  {
    setring save;
    vector K = imap(r,K);
  }
  return(K);
}


proc simplifyRat (vector J)
"
USAGE:    simplifyRat(v);  v vector
ASSUME:   Assume that v has exactly two components, second one not 0.
RETURN:   vector, representing simplified rational function v[1]/v[2]
NOTE:     Possibly present non-commutative relations of the basering are
          ignored.
EXAMPLE:  example simplifyRat; shows examples
"
{
  checkRatInput(J);
  return(arithmeticRat(J,vector(0),"s"));
}
example
{
  "EXAMPLE:"; echo = 2;
   ring r = 0,(x,y),dp;
  vector v = [x2-1,x+1];
  simplifyRat(v);
  simplifyRat(v) - [x-1,1];
}


proc addRat (vector I, vector J)
"
USAGE:    addRat(v,w);  v,w vectors
ASSUME:   Assume that v,w have exactly two components, second ones not 0.
RETURN:   vector, representing rational function (v[1]/v[2])+(w[1]/w[2])
NOTE:     Possibly present non-commutative relations of the basering are
          ignored.
EXAMPLE:  example addRat; shows examples
"
{
  checkRatInput(I);
  checkRatInput(J);
  return(arithmeticRat(I,J,"+"));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  vector v = [x,y];
  vector w = [y,x];
  addRat(v,w);
  addRat(v,w) - [x2+y2,xy];
}


proc multRat (vector I, vector J)
"
USAGE:    multRat(v,w);  v,w vectors
ASSUME:   Assume that v,w have exactly two components, second ones not 0.
RETURN:   vector, representing rational function (v[1]/v[2])*(w[1]/w[2])
NOTE:     Possibly present non-commutative relations of the basering are
          ignored.
EXAMPLE:  example multRat; shows examples
"
{
  checkRatInput(I);
  checkRatInput(J);
  return(arithmeticRat(I,J,"*"));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  vector v = [x,y];
  vector w = [y,x];
  multRat(v,w);
  multRat(v,w) - [1,1];
}


proc diffRat (vector I, int j)
"
USAGE:    diffRat(v,j);  v vector, j int
ASSUME:   Assume that v has exactly two components, second one not 0.
RETURN:   vector, representing rational function derivative of rational
          function (v[1]/v[2]) w.r.t. var(j)
NOTE:     Possibly present non-commutative relations of the basering are
          ignored.
EXAMPLE:  example diffRat; shows examples
"
{
  checkRatInput(I);
  if ( (j<1) || (j>nvars(basering)) )
  {
    ERROR("Second argument must be in the range 1.."+string(nvars(basering)));
  }
  return(arithmeticRat(I,vector(0),"d",j));
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y),dp;
  vector v = [x,y];
  diffRat(v,1);
  diffRat(v,1) - [1,y];
  diffRat(v,2);
  diffRat(v,2) - [-x,y2];
}


proc commRing ()
"
USAGE:    commRing();
RETURN:   ring, basering without non-commutative relations
EXAMPLE:  example commRing; shows examples
"
{
  list RL = ringlist(basering);
  if (size(RL)<=4)
  {
    return(basering);
  }
  RL = RL[1..4];
  def r = ring(RL);
  return(r)
}
example
{
  "EXAMPLE:"; echo = 2;
  def W = makeWeyl(3);
  setring W; W;
  def W2 = commRing();
  setring W2; W2;
  ring r = 0,(x,y),dp;
  def r2 = commRing(); // same as r
  setring r2; r2;
}


// TODO remove this proc once chern.lib is ready
static proc orderedPartition(int n, list #)
"
USUAGE:  orderedPartition(n,a); n,a positive ints
         orderedPartition(n,w); n positive int, w positive intvec
RETURN:  list of intvecs
PURPOSE: Computes all partitions of n of length a, if the second
         argument is an int, or computes all weighted partitions 
         w.r.t. w of n of length size(w) if the second argument
         is an intvec.
         In both cases, zero parts are included.
EXAMPLE: example orderedPartition; shows an example
"
{
  int a,wrongInpt,intInpt;
  intvec w = 1;
  if (size(#)>0)
  {
    if (typeof(#[1]) == "int")
    {
      a = #[1];
      intInpt = 1;
    }
    else
    {
      if (typeof(#[1]) == "intvec")
      {
        w = #[1];
        a = size(w);
      }
      else
      {
        wrongInpt = 1;
      }
    }
  }
  else
  {
    wrongInpt = 1;
  }
  if (wrongInpt)
  {
    ERROR("Expected second argument of type int or intvec.");
  }
  kill wrongInpt;
  if (n==0 && a>0)
  {
    return(list(0:a));
  }
  if (n<=0 || a<=0 || allPositive(w)==0)
  {
    ERROR("Positive arguments expected.");
  }
  int baseringdef;
  if (defined(basering)) // if a basering is defined, it should be saved for later use
  {
    def save = basering;
    baseringdef = 1;
  }
  ring r = 0,(x(1..a)),dp; // all variables for partition of length a
  ideal M;
  if (intInpt)
  {
    M = maxideal(n);   // all monomials of total degree n
  }
  else
  {
    M = weightKB(ideal(0),n,w); // all monomials of total weighted degree n
  }  
  list L; 
  int i;
  for (i = 1; i <= ncols(M); i++) {L = insert(L,leadexp(M[i]));} 
  // the leadexp corresponds to a partition
  if (baseringdef) // sets the old ring as basering again
  {
    setring save;
  }
  return(L); //returns the list of partitions
}
example
{
 "EXAMPLE"; echo = 2;
 orderedPartition(4,2);
 orderedPartition(5,3);
 orderedPartition(2,4);
 orderedPartition(8,intvec(2,3));
 orderedPartition(7,intvec(2,2)); // no such partition
}


// applications of characteristic variety /////////////////////////////////////

proc holonomicRank (ideal I, list #)
"
USAGE:    holonomicRank(I[,e]);   I ideal, e optional int
ASSUME:   The basering is the n-th Weyl algebra over a field of
          characteristic 0 and for all 1<=i<=n the identity
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
          is the differential operator belonging to x(i).
RETURN:   int, the holonomic rank of I
REMARKS:  The holonomic rank of I is defined to be the K(x(1..n))-dimension of 
          I regarded as ideal in the rational Weyl algebra K(x(1..n))<D(1..n)>.
          If this dimension is not finite, -1 is returned.
NOTE:     If e<>0, @code{std} is used for Groebner basis computations,
          otherwise (and by default) @code{slimgb} is used.
@*        If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
EXAMPLE:  example holonomicRank; shows examples
"
{
  // assumption check is done by charVariety
  int ppl = printlevel - voice + 2;
  int eng;
  if (size(#)>0)
  {
    if(typeof(#[1])=="int")
    {
      eng = #[1];
    }
  }
  def save = basering;
  dbprint(ppl  ,"// Computing characteristic variety...");
  def grD = charVariety(I);
  setring grD; // commutative ring
  dbprint(ppl  ,"// ...done.");
  dbprint(ppl-1,"// " + string(charVar));
  int n = nvars(save) div 2;
  intvec v = 1..n;
  def R = vars2pars(v);
  setring R;
  ideal J = imap(grD,charVar);
  dbprint(ppl  ,"// Starting GB computation...");
  J = engine(J,0); // use slimgb
  dbprint(ppl  ,"// ...done.");
  dbprint(ppl-1,"// " + string(J));
  int d = vdim(J);
  setring save;
  return(d);
}
example
{
  "EXAMPLE:"; echo = 2;
  // (OTW), Example 8
  ring r3 = 0,(x,y,z,Dx,Dy,Dz),dp;
  def D3 = Weyl();
  setring D3;
  poly f = x^3-y^2*z^2;
  ideal I = f^2*Dx+3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z;
  // I annihilates exp(1/f)
  holonomicRank(I);
}


proc DsingularLocus (ideal I)
"
USAGE:    DsingularLocus(I);  I ideal
ASSUME:   The basering is the n-th Weyl algebra over a field of
          characteristic 0 and for all 1<=i<=n the identity
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
          is the differential operator belonging to x(i).
RETURN:   ideal, the singular locus of the D-module D/I
NOTE:     If printlevel>=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed
EXAMPLE:  example DsingularLocus; shows examples
"
{
  // assumption check is done by charVariety
  int ppl = printlevel - voice + 2;
  def save = basering;
  dbprint(ppl  ,"// Computing characteristic variety...");
  def grD = charVariety(I);
  setring grD;
  dbprint(ppl  ,"// ...done");
  dbprint(ppl-1,"// " + string(charVar));
  poly pDD = 1;
  ideal IDD;
  int i;
  int n = nvars(basering) div 2;
  for (i=1; i<=n; i++)
  {
    pDD = pDD*var(i+n);
    IDD[i] = var(i+n);
  }
  dbprint(ppl  ,"// Computing saturation...");
  ideal S = sat(charVar,IDD)[1];
  dbprint(ppl  ,"// ...done");
  dbprint(ppl-1,"// " + string(S));
  dbprint(ppl  ,"// Computing elimination...");
  S = eliminate(S,pDD);
  dbprint(ppl  ,"// ...done");
  dbprint(ppl-1,"// " + string(S));
  dbprint(ppl  ,"// Computing radical...");
  S = radical(S);
  dbprint(ppl  ,"// ...done");
  dbprint(ppl-1,"// " + string(S));
  setring save;
  ideal S = imap(grD,S);
  return(S);
}
example
{
  "EXAMPLE:"; echo = 2;
  // (OTW), Example 8
  ring @D3 = 0,(x,y,z,Dx,Dy,Dz),dp;
  def D3 = Weyl();
  setring D3;
  poly f = x^3-y^2*z^2;
  ideal I = f^2*Dx + 3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z;
  // I annihilates exp(1/f)
  DsingularLocus(I);
}


// localization ///////////////////////////////////////////////////////////////

static proc finKx(poly f)
{
  int n = nvars(basering) div 2;
  intvec iv = 1..n;
  if (polyVars(f,iv) == 0)
  {
    ERROR("Given poly may not contain differential operators.");
  }
  return();
}


proc rightNFWeyl (def id, int k)
"
USAGE:    rightNFWeyl(id,k);  id ideal or poly, k int
ASSUME:   The basering is the n-th Weyl algebra over a field of
          characteristic 0 and for all 1<=i<=n the identity
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
          is the differential operator belonging to x(i).
RETURN:   same type as id, the right normal form of id with respect to the
          principal right ideal generated by the k-th variable
NOTE:     No Groebner basis computation is used.
EXAMPLE:  example rightNFWeyl; shows examples.
"
{
  dmodGeneralAssumptionCheck();
  string inpt = typeof(id);
  if (inpt=="ideal" || inpt=="poly")
  {
    ideal I = id;
  }
  else
  {
    ERROR("Expected first input to be of type ideal or poly.");
  }
  def save = basering;
  int n = nvars(save) div 2;
  if (0>k || k>2*n)
  {
    ERROR("Expected second input to be in the range 1.."+string(2*n)+".");
  }
  int i,j;
  if (k>n) // var(k) = Dx(k-n)
  {
    // switch var(k),var(k-n)
    list RL = ringlist(save);
    matrix rel = RL[6];
    rel[k-n,k] = -1;
    RL = RL[1..4];
    list L = RL[2];
    string str = L[k-n];
    L[k-n] = L[k];
    L[k] = str;
    RL[2] = L;
    def @W = ring(RL);
    kill L,RL,str;
    ideal @mm = maxideal(1);
    setring @W;
    matrix rel = imap(save,rel);
    def W = nc_algebra(1,rel);
    setring W;
    ideal @mm = imap(save,@mm);
    map mm = save,@mm;
    ideal I = mm(I);
    i = k-n;
  }
  else  // var(k) = x(k)
  {
    def W = save;
    i = k;
  }
  for (j=1; j<=ncols(I); j++)
  {
    I[j] = killTerms(I[j],i);
  }
  setring save;
  I = imap(W,I);
  if (inpt=="poly")
  {
    return(I[1]);
  }
  else
  {
    return(I);
  }
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,y,Dx,Dy),dp;
  def W = Weyl();
  setring W;
  ideal I = x^3*Dx^3, y^2*Dy^2, x*Dy, y*Dx;
  rightNFWeyl(I,1); // right NF wrt principal right ideal x*W
  rightNFWeyl(I,3); // right NF wrt principal right ideal Dx*W
  rightNFWeyl(I,2); // right NF wrt principal right ideal y*W
  rightNFWeyl(I,4); // right NF wrt principal right ideal Dy*W
  poly p = x*Dx+1;
  rightNFWeyl(p,1); // right NF wrt principal right ideal x*W
}


// TODO check OTW for assumptions on holonomicity
proc Dlocalization (ideal J, poly f, list #)
"
USAGE:    Dlocalization(I,f[,k,e]);  I ideal, f poly, k,e optional ints
ASSUME:   The basering is the n-th Weyl algebra over a field of
          characteristic 0 and for all 1<=i<=n the identity
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i)
          is the differential operator belonging to x(i).
@*        Further, assume that f does not contain any D(i) and that I is
          holonomic on K^n\V(f).
RETURN:   ideal or list, computes an ideal J such that D/J is isomorphic
          to D/I localized at f as D-modules.
          If k<>0, a list consisting of J and an integer m is returned,
          such that f^m represents the natural map from D/I to D/J.
          Otherwise (and by default), only the ideal J is returned.
REMARKS:  Reference: (OTW)
NOTE:     If e<>0, @code{std} is used for Groebner basis computations,
          otherwise (and by default) @code{slimgb} is used.
@*        If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
SEE ALSO: DLoc, SDLoc, DLoc0
EXAMPLE: example Dlocalization; shows examples
"
{
  dmodGeneralAssumptionCheck();
  finKx(f);
  int ppl = printlevel - voice + 2;
  int outList,eng;
  if (size(#)>0)
  {
    if (typeof(#[1])=="int" || typeof(#[1])=="number")
    {
      outList = int(#[1]);
    }
    if (size(#)>1)
    {
      if (typeof(#[2])=="int" || typeof(#[2])=="number")
      {
        eng = int(#[2]);
      }
    }
  }
  int i,j;
  def save = basering;
  int n = nvars(save) div 2;
  def Dv = extendWeyl(safeVarName("v"));
  setring Dv;
  poly f = imap(save,f);
  ideal phiI;
  for (i=n; i>0; i--)
  {
    phiI[i+n] = var(i+n+2)-var(1)^2*bracket(var(i+n+2),f)*var(n+2);
    phiI[i]   = var(i+1);
  }
  map phi = save,phiI;
  ideal J = phi(J);
  J = J, 1-f*var(1);
  // TODO original J has to be holonomic only on K^n\V(f), not on all of K^n
  // does is suffice to show that new J is holonomic on Dv??
  if (isHolonomic(J) == 0)
  {
    ERROR("Module is not holonomic.");
  }
  intvec w = 1; w[n+1]=0;
  ideal G = GBWeight(J,w,-w,eng);
  dbprint(ppl  ,"// found GB wrt weight " +string(-w));
  dbprint(ppl-1,"// " + string(G));
  intvec ww = w,-w;
  ideal inG = inForm(G,ww);
  inG = engine(inG,eng);
  poly s = var(1)*var(n+2); // s=v*Dv
  vector intersecvec = pIntersect(s,inG);
  s = vec2poly(intersecvec);
  s = subst(s,var(1),-var(1)-1);
  list L = bFactor(s);
  dbprint(ppl  ,"// found b-function");
  dbprint(ppl-1,"// roots: "+string(L[1]));
  dbprint(ppl-1,"// multiplicities: "+string(L[2]));
  kill inG,intersecvec,s;
  // TODO: use maxIntRoot
  L = intRoots(L);           // integral roots of b-function
  if (L[2]==0:size(L[2]))    // no integral roots
  {
    setring save;
    return(ideal(1));
  }
  intvec iv;
  for (i=1; i<=ncols(L[1]); i++)
  {
    iv[i] = int(L[1][i]);
  }
  int l0 = Max(iv);
  dbprint(ppl,"// maximal integral root is " +string(l0));
  kill L,iv;
  intvec degG;
  ideal Gk;
  for (j=1; j<=ncols(G); j++)
  {
    degG[j] = deg(G[j],ww);
    for (i=0; i<=l0-degG[j]; i++)
    {
      Gk[ncols(Gk)+1] = var(1)^i*G[j];
    }
  }
  Gk = rightNFWeyl(Gk,n+2);
  dbprint(ppl,"// found right normalforms");
  module M = coeffs(Gk,var(1));
  setring save;
  def mer = makeModElimRing(save);
  setring mer;
  module M = imap(Dv,M);
  kill Dv;
  M = engine(M,eng);
  dbprint(ppl  ,"// found GB of free module of rank " + string(l0+1));
  dbprint(ppl-1,"// " + string(M));
  M = prune(M);
  setring save;
  matrix M = imap(mer,M);
  kill mer;
  int ro = nrows(M);
  int co = ncols(M);
  ideal I;
  if (ro == 1) // nothing to do
  {
    I = M;
  }
  else
  {
    matrix zm[ro-1][1]; // zero matrix
    matrix v[ro-1][1];
    for (i=1; i<=co; i++)
    {
      v = M[1..ro-1,i];
      if (v == zm)
      {
        I[size(I)+1] = M[ro,i];
      }
    }
  }
  if (outList<>0)
  {
    return(list(I,l0+2));
  }
  else
  {
    return(I);
  }  
}
example
{
  "EXAMPLE:"; echo = 2;
  // (OTW), Example 8
  ring r = 0,(x,y,z,Dx,Dy,Dz),dp;
  def W = Weyl();
  setring W;
  poly f = x^3-y^2*z^2;
  ideal I = f^2*Dx+3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z;
  // I annihilates exp(1/f)
  ideal J = Dlocalization(I,f);
  J;
  Dlocalization(I,f,1); // The natural map D/I -> D/J is given by 1/f^2
}



// Weyl closure ///////////////////////////////////////////////////////////////

static proc orderFiltrationD1 (poly f)
{
  // returns list of ideal and intvec
  // ideal contains x-parts, intvec corresponding degree in Dx
  poly g,h;
  g = f;
  ideal I;
  intvec v,w,u;
  w = 0,1;
  int i,j;
  i = 1;
  while (g<>0)
  {
    h = inForm(g,w);
    I[i] = 0;
    for (j=1; j<=size(h); j++)
    {
      v = leadexp(h[j]);
      u[i] = v[2];
      v[2] = 0;
      I[i] = I[i] + leadcoef(h[j])*monomial(v);
    }
    g = g-h;
    i++;
  }
  return(list(I,u));
}


static proc kerLinMapD1 (ideal W, poly L, poly p)
{
  // computes kernel of right multiplication with L viewed
  // as homomorphism of K-vector spaces span(W) -> D1/p*D1
  // assume p in K[x], basering is K<x,Dx>
  ideal G,K;
  G = std(p);
  list l;
  int i,j;
  // first, compute the image of span(W)
  if (simplify(W,2)[1] == 0)
  {
    return(K); // = 0
  }
  for (i=1; i<=size(W); i++)
  {
    l = orderFiltrationD1(W[i]*L);
    K[i] = 0;
    for (j=1; j<=size(l[1]); j++)
    {
      K[i] = K[i] + NF(l[1][j],G)*var(2)^(l[2][j]);
    }
  }
  // now, we get the kernel by linear algebra
  l = linReduceIdeal(K,1);
  i = ncols(l[1]) - size(l[1]);
  if (i<>0)
  {
    K = module(W)*l[2];
    K = K[1..i];
  }
  else
  {
    K = 0;
  }
  return(K);
}


static proc leftDivisionKxD1 (poly p, poly L)
{
  // basering is D1 = K<x,Dx>
  // p in K[x]
  // compute p^(-1)*L if p is a left divisor of L
//   if (rightNF(L,ideal(p))<>0)
//   {
//     ERROR("First poly is not a right factor of second poly");
//   }
  def save = basering;
  list l = orderFiltrationD1(L);
  ideal l1 = l[1];
  ring r = 0,x,dp;
  ideal l1 = fetch(save,l1);
  poly p = fetch(save,p);
  int i;
  for (i=1; i<=ncols(l1); i++)
  {
    l1[i] =  division(l1[i],p)[1][1,1];
  }
  setring save;
  ideal I = fetch(r,l1);
  poly f;
  for (i=1; i<=ncols(I); i++)
  {
    f = f + I[i]*var(2)^(l[2][i]);
  }
  return(f);
}


proc WeylClosure1 (poly L)
"
USAGE:    WeylClosure1(L);  L a poly
ASSUME:   The basering is the first Weyl algebra D=K<x,d|dx=xd+1> over a field
          K of characteristic 0.
RETURN:   ideal, the Weyl closure of the principal left ideal generated by L
REMARKS:  The Weyl closure of a left ideal I in the Weyl algebra D is defined
          to be the intersection of I regarded as left ideal in the rational
          Weyl algebra K(x)<d> with D.
@*        Reference: (Tsa), Algorithm 1.2.4
NOTE:     If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
SEE ALSO: WeylClosure
EXAMPLE:  example WeylClosure1; shows examples
"
{    
  dmodGeneralAssumptionCheck(); // assumption check
  int ppl = printlevel - voice + 2;
  def save = basering;
  intvec w = 0,1; // for order filtration
  poly p = inForm(L,w);
  ring @R = 0,var(1),dp;
  ideal mm = var(1),1;
  map m = save,mm;
  ideal @p = m(p);
  poly p = @p[1];
  poly g = gcd(p,diff(p,var(1)));
  if (g == 1)
  {
    g = p;
  }
  ideal facp = factorize(g,1); // g is squarefree, constants aren't interesting
  dbprint(ppl-1,
          "// squarefree part of highest coefficient w.r.t. order filtration:");
  dbprint(ppl-1, "// " + string(facp));
  setring save;
  p = imap(@R,p);
  // 1-1 extend basering by parameter and introduce new var t=x*d
  list RL = ringlist(save);
  RL = RL[1..4];
  list l;
  l[1] = int(0);
  l[2] = list(safeVarName("a"));
  l[3] = list(list("lp",intvec(1)));
  l[4] = ideal(0);
  RL[1] = l;
  l = RL[2] + list(safeVarName("t"));
  RL[2] = l;
  l = list();
  l[1] = list("dp",intvec(1,1));
  l[2] = list("dp",intvec(1));
  l[3] = list("C",intvec(0));
  RL[3] = l;
  def @Wat = ring(RL);
  kill RL,l;
  setring @Wat;
  matrix relD[3][3];
  relD[1,2] = 1;
  relD[1,3] = var(1);
  relD[2,3] = -var(2);
  def Wat = nc_algebra(1,relD);
  setring Wat;
  kill @Wat;
  // 1-2 rewrite L using Euler operators
  ideal mm = var(1)+par(1),var(2);
  map m = save,mm;
  poly L = m(L);
  w = -1,1,0; // for Bernstein filtration
  int i = 1;
  ideal Q;
  poly p = L;
  intvec d;
  while (p<>0)
  {
    Q[i] = inForm(p,w);
    p = p - Q[i];
    d[i] = -deg(Q[i],w);
    i++;
  }
  ideal S = std(var(1)*var(2)-var(3));
  Q = NF(Q,S);
  dbprint(ppl,  "// found Euler representation of operator");
  dbprint(ppl-1,"// " + string(Q));
  Q = subst(Q,var(1),1);
  Q = subst(Q,var(2),1);
  // 1-3 prepare for algebraic extensions with minpoly = facp[i]
  list RL = ringlist(Wat);
  RL = RL[1..4];
  list l;
  l = string(var(3));
  RL[2] = l;
  l = list();
  l[1] = list("dp",intvec(1));
  l[2] = list("C",intvec(0));
  RL[3] = l;
  mm = par(1);
  m = @R,par(1);
  ideal facp = m(facp);
  kill @R,m,mm,l,S;
  intvec maxroots,testroots;
  int sq = size(Q);
  string strQ = "ideal Q = " + string(Q) + ";";
  // TODO do it without string workaround when issue with maps from 
  //   transcendental to algebraic extension fields is fixed
  int j,maxr;
  // 2-1 get max int root of lowest nonzero entry of Q in algebraic extension
  for (i=1; i<=size(facp); i++)
  {
    testroots = 0;
    def Ra = ring(RL);
    setring Ra;
    ideal mm = 1,1,var(1);
    map m = Wat,mm;
    ideal facp = m(facp);
    minpoly = leadcoef(facp[i]);
    execute(strQ);
    if (simplify(Q,2)[1] == poly(0))
    {
      break;
    }
    j = 1;
    while (j<sq)
    {
      if (Q[j]==0)
      {
        j++;
      }
      else
      {
        break;
      }
    }
    maxroots[i] = d[j]; // d[j] = r_k 
    list LR = bFactor(Q[j]);
    LR = intRoots(LR);
    if (LR[2]<>0:size(LR[2])) // there are integral roots
    {
      for (j=1; j<=ncols(LR[1]); j++)
      {
        testroots[j] = int(LR[1][j]);
      }
      maxr = Max(testroots);
      if(maxr<0)
      {
        maxr = 0;
      }
      maxroots[i] = maxroots[i] + maxr;
    }
    kill LR;
    setring Wat;
    kill Ra;
  }
  maxr = Max(maxroots);
  // 3-1 build basis of vectorspace
  setring save;
  ideal KB;
  for (i=0; i<deg(p); i++)  // it's really <, not <=
  {
    for (j=0; j<=maxr; j++) // it's really <=, not <
    {
      KB[size(KB)+1] = monomial(intvec(i,j));
    }
  }
  dbprint(ppl,"// got vector space basis");
  dbprint(ppl-1, "// " + string(KB));
  // 3-2 get kernel of *L: span(KB)->D/pD
  KB = kerLinMapD1(KB,L,p);
  dbprint(ppl,"// got kernel");
  dbprint(ppl-1, "// " + string(KB));
  // 4-1 get (1/p)*f*L where f in KB
  for (i=1; i<=ncols(KB); i++)
  {
    KB[i] = leftDivisionKxD1(p,KB[i]*L);
  }
  KB = L,KB;
  // 4-2 done
  return(KB);
}
example
{
  "EXAMPLE:"; echo = 2;
  ring r = 0,(x,Dx),dp;
  def W = Weyl();
  setring W;
  poly L = (x^3+2)*Dx-3*x^2;
  WeylClosure1(L);
  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);
  WeylClosure1(L);
}


proc WeylClosure (ideal I)
"
USAGE:    WeylClosure(I);  I an ideal
ASSUME:   The basering is the n-th Weyl algebra W over a field of 
          characteristic 0 and for all 1<=i<=n the identity 
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of 
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
          differential operator belonging to x(i).
@*        Moreover, assume that the holonomic rank of W/I is finite.
RETURN:   ideal, the Weyl closure of I
REMARKS:  The Weyl closure of a left ideal I in the Weyl algebra W is defined to
          be the intersection of I regarded as left ideal in the rational Weyl
          algebra K(x(1..n))<D(1..n)> with W.
@*        Reference: (Tsa), Algorithm 2.2.4
NOTE:     If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
SEE ALSO: WeylClosure1
EXAMPLE:  example WeylClosure; shows examples
"
{
  // assumption check
  dmodGeneralAssumptionCheck();
  if (holonomicRank(I)==-1)
  {
    ERROR("Input is not of finite holonomic rank.");
  }
  int ppl = printlevel - voice + 2;
  int eng = 0; // engine
  def save = basering;
  dbprint(ppl  ,"// Starting to compute singular locus...");
  ideal sl = DsingularLocus(I);
  sl = simplify(sl,2);
  dbprint(ppl  ,"// ...done.");
  dbprint(ppl-1,"// " + string(sl));
  if (sl[1] == 0) // can never get here
  {
    ERROR("Can't find polynomial in K[x] vanishing on singular locus.");
  }
  poly f = sl[1];
  dbprint(ppl  ,"// Found poly vanishing on singular locus: " + string(f));
  dbprint(ppl  ,"// Starting to compute localization...");
  list L = Dlocalization(I,f,1);
  ideal G = L[1];
  dbprint(ppl  ,"// ...done.");
  dbprint(ppl-1,"// " + string(G));
  dbprint(ppl  ,"// Starting to compute kernel of localization map...");
  if (eng == 0)
  {
    G = moduloSlim(f^L[2],G);
  }
  else
  {
    G = modulo(f^L[2],G);
  }
  dbprint(ppl  ,"// ...done.");
  return(G);
}
example
{
  "EXAMPLE:"; echo = 2;
  // (OTW), Example 8
  ring r = 0,(x,y,z,Dx,Dy,Dz),dp;
  def D3 = Weyl();
  setring D3;
  poly f = x^3-y^2*z^2;
  ideal I = f^2*Dx + 3*x^2, f^2*Dy-2*y*z^2, f^2*Dz-2*y^2*z;
  // I annihilates exp(1/f)
  WeylClosure(I);
}



// solutions to systems of PDEs ///////////////////////////////////////////////

proc polSol (ideal I, list #)
"
USAGE:    polSol(I[,w,m]);  I ideal, w optional intvec, m optional int
ASSUME:   The basering is the n-th Weyl algebra W over a field of 
          characteristic 0 and for all 1<=i<=n the identity 
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of 
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
          differential operator belonging to x(i).
@*        Moreover, assume that I is holonomic.
RETURN:   ideal, a basis of the polynomial solutions to the given system
REMARKS:  If w is given, w should consist of n strictly negative entries.
          Otherwise and by default, w is set to -1:n.
          In this case, w is used as weight vector for the computation of a
          b-function.
@*        If m is given, m is assumed to be the minimal integer root of the
          b-function of I w.r.t. w. Note that this assumption is not checked.
@*        Reference: (OTT), Algorithm 2.4
NOTE:     If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
SEE ALSO: polSolFiniteRank, ratSol
EXAMPLE:  example polSol; shows examples
"
{
  dmodGeneralAssumptionCheck();
  int ppl = printlevel - voice + 2;
  int mr,mrgiven;
  def save = basering;
  int n = nvars(save);
  intvec w = -1:(n div 2);
  if (size(#)>0)
  {
    if (typeof(#[1])=="intvec")
    {
      if (allPositive(-#[1]))
      {
        w = #[1];
      }
    }
    if (size(#)>1)
    {
      if (typeof(#[2])=="int")
      {
        mr = #[2];
        mrgiven = 1;
      }
    }
  }
  // Step 1: the b-function
  list L;
  if (!mrgiven)
  {
    if (!isHolonomic(I))
    {
      ERROR("Ideal is not holonomic. Try polSolFiniteRank.");
    }
    dbprint(ppl,"// Computing b-function...");
    L = bfctIdeal(I,w);
    dbprint(ppl,"// ...done.");
    dbprint(ppl-1,"//   Roots: " + string(L[1]));
    dbprint(ppl-1,"//   Multiplicities: " + string(L[2]));
    mr = minIntRoot2(L);
    dbprint(ppl,"// Minimal integer root is " + string(mr) + ".");
  }
  if (mr>0)
  {
    return(ideal(0));
  }
  // Step 2: get the form of a solution f
  int i;
  L = list();
  for (i=0; i<=-mr; i++)
  {
    L = L + orderedPartition(i,-w);
  }
  ideal mons;
  for (i=1; i<=size(L); i++)
  {
    mons[i] = monomial(L[i]);
  }
  kill L;
  mons = simplify(mons,2+4); // L might contain lots of 0s by construction
  ring @C = (0,@c(1..size(mons))),dummyvar,dp;
  def WC = save + @C;
  setring WC;
  ideal mons = imap(save,mons);
  poly f;
  for (i=1; i<=size(mons); i++)
  {
    f = f + par(i)*mons[i];
  }
  // Step 3: determine values of @c(i) by equating coefficients
  ideal I = imap(save,I);
  I = dmodAction(I,f,1..n);
  ideal M = monomialInIdeal(I);
  matrix CC = coeffs(I,M);
  int j;
  ideal C;
  for (i=1; i<=nrows(CC); i++)
  {
    f = 0;
    for (j=1; j<=ncols(CC); j++)
    {
      f = f + CC[i,j];
    }
    C[size(C)+1] = f;
  }
  // Step 3.1: solve a linear system
  ring rC = 0,(@c(1..size(mons))),dp;
  ideal C = imap(WC,C);
  matrix M = coeffs(C,maxideal(1));
  module MM = leftKernel(M);
  setring WC;
  module MM = imap(rC,MM);
  // Step 3.2: return the solution
  ideal F = ideal(MM*transpose(mons));
  setring save;
  ideal F = imap(WC,F);
  return(F);
}
example
{
  "EXAMPLE:"; echo=2;
  ring r = 0,(x,y,Dx,Dy),dp;
  def W = Weyl();
  setring W;
  poly tx,ty = x*Dx, y*Dy;
  ideal I =      // Appel F1 with parameters (2,-3,-2,5)
    tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3),
    ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2),
    (x-y)*Dx*Dy+2*Dx-3*Dy;
  intvec w = -1,-1;
  polSol(I,w);
}


static proc ex_polSol()
{  ring r = 0,(x,y,Dx,Dy),dp;
  def W = Weyl();
  setring W;
  poly tx,ty = x*Dx, y*Dy;
  ideal I =      // Appel F1 with parameters (2,-3,-2,5)
    tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3),
    ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2),
    (x-y)*Dx*Dy+2*Dx-3*Dy;
  intvec w = -5,-7;
  // the following gives a bug
  polSol(I,w);
  // this is due to a bug in weightKB, see ticket #339
  // http://www.singular.uni-kl.de:8002/trac/ticket/339
}


proc polSolFiniteRank (ideal I, list #)
"
USAGE:    polSolFiniteRank(I[,w]);  I ideal, w optional intvec
ASSUME:   The basering is the n-th Weyl algebra W over a field of 
          characteristic 0 and for all 1<=i<=n the identity 
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of 
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
          differential operator belonging to x(i).
@*        Moreover, assume that I is of finite holonomic rank.
RETURN:   ideal, a basis of the polynomial solutions to the given system
REMARKS:  If w is given, w should consist of n strictly negative entries.
          Otherwise and by default, w is set to -1:n.
          In this case, w is used as weight vector for the computation of a
          b-function.
@*        Reference: (OTT), Algorithm 2.6
NOTE:     If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
SEE ALSO: polSol, ratSol
EXAMPLE:  example polSolFiniteRank; shows examples
"
{
  dmodGeneralAssumptionCheck();
  if (holonomicRank(I)==-1)
  {
    ERROR("Ideal is not of finite holonomic rank.");
  }
  int ppl = printlevel - voice + 2;
  int n = nvars(basering) div 2;
  int eng;
  intvec w = -1:(n div 2);
  if (size(#)>0)
  {
    if (typeof(#[1])=="intvec")
    {
      if (allPositive(-#[1]))
      {
        w = #[1];
      }
    }
  }
  dbprint(ppl,"// Computing initial ideal...");
  ideal J = initialIdealW(I,-w,w);
  dbprint(ppl,"// ...done.");
  dbprint(ppl,"// Computing Weyl closure...");
  J = WeylClosure(J);
  J = engine(J,eng);
  dbprint(ppl,"// ...done.");
  poly s;
  int i;
  for (i=1; i<=n; i++)
  {
    s = s + w[i]*var(i)*var(i+n);
  }
  dbprint(ppl,"// Computing intersection...");
  vector v = pIntersect(s,J);
  list L = bFactor(vec2poly(v));
  dbprint(ppl-1,"//   roots: " + string(L[1]));
  dbprint(ppl-1,"//   multiplicities: " + string(L[2]));
  dbprint(ppl,"// ...done.");
  int mr =  minIntRoot2(L);
  int pl = printlevel;
  printlevel = printlevel + 1;
  ideal K = polSol(I,w,mr);
  printlevel = printlevel - 1;
  return(K);
}
example
{
  "EXAMPLE:"; echo=2;
  ring r = 0,(x,y,Dx,Dy),dp;
  def W = Weyl();
  setring W;
  poly tx,ty = x*Dx, y*Dy;
  ideal I =      // Appel F1 with parameters (2,-3,-2,5)
    tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3),
    ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2),
    (x-y)*Dx*Dy+2*Dx-3*Dy;
  intvec w = -1,-1;
  polSolFiniteRank(I,w);
}


static proc twistedIdeal(ideal I, poly f, intvec k, ideal F)
{
  // I subset D_n, f in K[x], F = factorize(f,1), size(k) = size(F), k[i]>0
  def save = basering;
  int n = nvars(save) div 2;
  int i,j;
  intvec a,v,w;
  w = (0:n),(1:n);
  for (i=1; i<=size(I); i++)
  {
    a[i] = deg(I[i],w);
  }
  ring FD = 0,(fd(1..n)),dp;
  def @@WFD = save + FD;
  setring @@WFD;
  poly f = imap(save,f);
  list RL = ringlist(basering);
  RL = RL[1..4];
  list L = RL[3];
  v = (1:(2*n)),((deg(f)+1):n);
  L = insert(L,list("a",v));
  RL[3] = L;
  def @WFD = ring(RL);
  setring @WFD;
  poly f = imap(save,f);
  matrix Drel[3*n][3*n];
  for (i=1; i<=n; i++)
  {
    Drel[i,i+n] = 1;     // [D,x]
    Drel[i,i+2*n] = f;   // [fD,x]
    for (j=1; j<=n; j++)
    {
      Drel[i+n,j+2*n] = -diff(f,var(i))*var(j+n);  // [fD,D]
      Drel[j+2*n,i+2*n] = diff(f,var(i))*var(j+2*n) - diff(f,var(j))*var(i+2*n);
      // [fD,fD]
    }
  }
  def WFD = nc_algebra(1,Drel);
  setring WFD;
  kill @WFD,@@WFD;
  ideal I = imap(save,I);
  poly f = imap(save,f);
  for (i=1; i<=size(I); i++)
  {
    I[i] = f^(a[i])*I[i];
  }
  ideal S;
  for (i=1; i<=n; i++)
  {
    S[size(S)+1] = var(i+2*n) - f*var(i+n);
  }
  S = slimgb(S);
  I = NF(I,S);
  if (select1(I,intvec((n+1)..2*n))[1] <> 0)
  {
    // should never get here
    ERROR("Something's wrong. Please inform the author.");
  }
  setring save;
  ideal mm = maxideal(1);
  poly s;
  for (i=1; i<=n; i++)
  {
    s = f*var(i+n);
    for (j=1; j<=size(F); j++)
    {
      s = s + k[j]*(f/F[j])*bracket(var(i+n),F[j]);
    }
    mm[i+2*n] = s;
  }
  map m = WFD,mm;
  ideal J = m(I);
  return(J);
}
example
{
  "EXAMPLE"; echo=2;
  ring r = 0,(x,y,Dx,Dy),dp;
  def W = Weyl();
  setring W;
  poly tx,ty = x*Dx, y*Dy;
  ideal I =      // Appel F1 with parameters (3,-1,1,1) is a solution
    tx*(tx+ty)-x*(tx+ty+3)*(tx-1),
    ty*(tx+ty)-y*(tx+ty+3)*(ty+1);
  kill tx,ty;
  poly f = x^3*y^2-x^2*y^3-x^3*y+x*y^3+x^2*y-x*y^2;
  ideal F = x-1,x,-x+y,y-1,y;
  intvec k = -1,-1,-1,-3,-1;
  ideal T = twistedIdeal(I,f,k,F);
  // TODO change the ordering of WFD
  // introduce new var f??
  //paper:
  poly fx = diff(f,x);
  poly fy = diff(f,y);
  poly fDx = f*Dx;
  poly fDy = f*Dy;
  poly fd(1) = fDx;
  poly fd(2) = fDy;
  ideal K=
    (x^2-x^3)*(fDx)^2+x*((1-3*x)*f-(1-x)*y*fy-(1-x)*x*fx)*(fDx)
    +x*(1-x)*y*(fDy)*(fDx)+x*y*f*(fDy)+3*x*f^2,
    (y^2-y^3)*(fDy)^2+y*((1-5*y)*f-(1-y)*x*fx-(1-y)*y*fy)*(fDy)
    +y*(1-y)*x*(fDx)*(fDy)-y*x*f*(fDx)-3*y*f^2;
}


proc ratSol (ideal I)
"
USAGE:    ratSol(I);  I ideal
ASSUME:   The basering is the n-th Weyl algebra W over a field of 
          characteristic 0 and for all 1<=i<=n the identity 
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of 
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
          differential operator belonging to x(i).
@*        Moreover, assume that I is holonomic.
RETURN:   module, a basis of the rational solutions to the given system.
          Note that each entry has two components, the first one standing for
          the enumerator, the second one for the denominator.
REMARKS:  Reference: (OTT), Algorithm 3.10
NOTE:     If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
SEE ALSO: polSol, polSolFiniteRank
EXAMPLE:  example ratSol; shows examples
"
{
  dmodGeneralAssumptionCheck();
  if (!isHolonomic(I))
  {
    ERROR("Ideal is not holonomic.");
  }
  int ppl = printlevel - voice + 2;
  def save = basering;
  dbprint(ppl,"// computing singular locus...");
  ideal S = DsingularLocus(I);
  dbprint(ppl,"// ...done.");
  poly f = S[1];
  dbprint(ppl,"// considering poly " + string(f));
  int n = nvars(save) div 2;
  list RL = ringlist(save);
  RL = RL[1..4];
  list L = RL[2];
  L = list(L[1..n]);
  RL[2] = L;
  L = list();
  L[1] = list("dp",intvec(1:n));
  L[2] = list("C",intvec(0));
  RL[3] = L;
  def rr = ring(RL);
  setring rr;
  poly f = imap(save,f);
  ideal F = factorize(f,1); // not interested in multiplicities
  dbprint(ppl,"// with irreducible factors " + string(F));
  setring save;
  ideal F = imap(rr,F);
  kill rr,RL;
  int i;
  intvec k;
  ideal FF = 1,1;
  dbprint(ppl,"// computing b-functions of irreducible factors...");
  for (i=1; i<=size(F); i++)
  {
    dbprint(ppl,"//   considering " + string(F[i]) + "...");
    L = bfctBound(I,F[i]);
    if (size(L) == 3) // bfct is constant
    {
      dbprint(ppl,"//   ...got " + string(L[3]));
      if (L[3] == "1")
      {
        return(0); // TODO type // no rational solutions
      }
      else // should never get here
      {
        ERROR("Oops, something went wrong. Please inform the author.");
      }
    }
    else
    {
      dbprint(ppl,"//   ...got roots " + string(L[1]));
      dbprint(ppl,"//      with multiplicities " + string(L[2]));
      k[i] = -maxIntRoot(L)-1;
      if (k[i] < 0)
      {
        FF[2] = FF[2]*F[i]^(-k[i]);
      }
      else
      {
        FF[1] = FF[1]*F[i]^(k[i]);
      }      
    }
  }
  vector v = FF[1]*gen(1) + FF[2]*gen(2);
  kill FF;
  dbprint(ppl,"// ...done");
  ideal twI = twistedIdeal(I,f,k,F);
  intvec w = -1:n;
  dbprint(ppl,"// computing polynomial solutions of twisted system...");
  if (isHolonomic(twI))
  {
    ideal P = polSol(twI,w);
  }
  else
  {
    ideal P = polSolFiniteRank(twI,w);
  }
  module M;
  vector vv;
  for (i=1; i<=ncols(P); i++)
  {
    vv = P[i]*gen(1) + 1*gen(2);
    M[i] = multRat(v,vv);
  }
  dbprint(ppl,"// ...done");
  return (M);
}
example
{
  "EXAMPLE"; echo=2;
  ring r = 0,(x,y,Dx,Dy),dp;
  def W = Weyl();
  setring W;
  poly tx,ty = x*Dx, y*Dy;
  ideal I =      // Appel F1 with parameters (3,-1,1,1) is a solution
    tx*(tx+ty)-x*(tx+ty+3)*(tx-1),
    ty*(tx+ty)-y*(tx+ty+3)*(ty+1);
  module M = ratSol(I);
  // We obtain a basis of the rational solutions to I represented by a
  // module / matrix with two rows.
  // Each column of the matrix represents a rational function, where
  // the first row correspond to the enumerator and the second row to
  // the denominator.
  print(M); 
}


proc bfctBound (ideal I, poly f, list #)
"
USAGE:    bfctBound (I,f[,primdec]);  I ideal, f poly, primdec optional string
ASSUME:   The basering is the n-th Weyl algebra W over a field of 
          characteristic 0 and for all 1<=i<=n the identity 
          var(i+n)*var(i)=var(i)*var(i+1)+1 holds, i.e. the sequence of 
          variables is given by x(1),...,x(n),D(1),...,D(n), where D(i) is the
          differential operator belonging to x(i).
@*        Moreover, assume that I is holonomic.
RETURN:   list of roots (of type ideal) and multiplicities (of type intvec) of
          a multiple of the b-function for f^s*u at a generic root of f. 
          Here, u stands for [1] in D/I.
REMARKS:  Reference: (OTT), Algorithm 3.4
NOTE:     This procedure requires to compute a primary decomposition in a 
          commutative ring. The optional string primdec can be used to specify
          the algorithm to do so. It may either be 'GTZ' (Gianni, Trager,
          Zacharias) or 'SY' (Shimoyama, Yokoyama). By default, 'GTZ' is used.
@*        If printlevel=1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
SEE ALSO: bernstein, bfct, bfctAnn
EXAMPLE:  example bfctBound; shows examples
"
{
  dmodGeneralAssumptionCheck();
  finKx(f);
  if (!isHolonomic(I))
  {
    ERROR("Ideal is not holonomic.");
  }
  int ppl = printlevel - voice + 2;
  string primdec = "GTZ";
  if (size(#)>1)
  {
    if (typeof(#[1])=="string")
    {
      if ( (#[1]=="SY") || (#[1]=="sy") || (#[1]=="Sy") )
      {
        primdec = "SY";
      }
      else
      {
        if ( (#[1]<>"GTZ") && (#[1]<>"gtz") && (#[1]<>"Gtz") )
        {
          print("// Warning: optional string may either be 'GTZ' or 'SY',");
          print("//          proceeding with 'GTZ'.");
          primdec = "GTZ";
        }
      }
    }
  }
  def save = basering;
  int n = nvars(save) div 2;
  // step 1
  ideal mm = maxideal(1);
  def Wt = extendWeyl(safeVarName("t"));
  setring Wt;
  poly f = imap(save,f);
  ideal mm = imap(save,mm);
  int i;
  for (i=1; i<=n; i++)
  {
    mm[i+n] = var(i+n+2) + bracket(var(i+n+2),f)*var(n+2);
  }
  map m = save,mm;
  ideal I = m(I);
  I = I, var(1)-f;
  // step 2
  intvec w = 1,(0:n);
  dbprint(ppl  ,"// Computing initial ideal...");
  I = initialIdealW(I,-w,w);
  dbprint(ppl  ,"// ...done.");
  dbprint(ppl-1,"// " + string(I));
  // step 3: rewrite I using Euler operator t*Dt
  list RL = ringlist(Wt);
  RL = RL[1..4];
  list L = RL[2] + list(safeVarName("s")); // s=t*Dt
  RL[2] = L;
  L = list();
  L[1] = list("dp",intvec(1:(2*n+2)));
  L[2] = list("dp",intvec(1));
  L[3] = list("C",intvec(0));
  RL[3] = L;
  def @Wts = ring(RL);
  kill L,RL;
  setring @Wts;
  matrix relD[2*n+3][2*n+3];
  relD[1,2*n+3] = var(1);
  relD[n+2,2*n+3] = -var(n+2);
  for (i=1; i<=n+1; i++)
  {
    relD[i,n+i+1] = 1;
  }
  def Wts = nc_algebra(1,relD);
  setring Wts;
  ideal I = imap(Wt,I);
  kill Wt,@Wts;
  ideal S = var(1)*var(n+2)-var(2*n+3);
  attrib(S,"isSB",1);
  dbprint(ppl  ,"// Computing Euler representation...");
  // I = NF(I,S);
  int d;
  intvec ww = 0:(2*2+2); ww[1] = -1; ww[n+2] = 1;
  for (i=1; i<=size(I); i++)
  {
    d = deg(I[i],ww);
    if (d>0)
    {
      I[i] = var(1)^d*I[i];
    }
    if (d<0)
    {
      d = -d;
      I[i] = var(n+2)^d*I[i];
    }
  }
  I = NF(I,S); // now there are no t,Dt in I, only s
  dbprint(ppl  ,"// ...done.");
  I = subst(I,var(2*n+3),-var(2*n+3)-1);
  ring Ks = 0,s,dp;
  def Ws = save + Ks;
  setring Ws;
  ideal I = imap(Wts,I);
  kill Wts;
  poly DD = 1;
  for (i=1; i<=n; i++)
  {
    DD = DD * var(n+i);
  }
  dbprint(ppl  ,"// Eliminating differential operators...");
  ideal J = eliminate(I,DD); // J subset K[x,s]
  dbprint(ppl  ,"// ...done.");
  dbprint(ppl-1,"// " + string(J));
  list RL = ringlist(Ws);
  RL = RL[1..4];
  list L = RL[2];
  L = list(L[1..n]) + list(L[2*n+1]);
  RL[2] = L;
  L = list();
  L[1] = list("dp",intvec(1:(n+1)));
  L[2] = list("C",intvec(0));
  RL[3] = L;
  def Kxs = ring(RL);
  setring Kxs;
  ideal J = imap(Ws,J);
  dbprint(ppl  ,"// Computing primary decomposition with engine " 
          + primdec + "...");
  if (primdec == "GTZ")
  {
    list P = primdecGTZ(J);
  }
  else
  {
    list P = primdecSY(J);
  }
  dbprint(ppl  ,"// ...done.");
  dbprint(ppl-1,"// " + string(P));
  ideal GP,Qix,rad,B;
  poly f = imap(save,f);
  vector v;
  for (i=1; i<=size(P); i++)
  {
    dbprint(ppl  ,"// Considering primary component " + string(i) 
            + " of " + string(size(P)) + "...");
    dbprint(ppl  ,"//   Intersecting with K[x] and computing radical...");
    GP = std(P[i][1]);
    Qix = eliminate(GP,var(n+1)); // subset K[x]
    rad = radical(Qix);
    rad = std(rad);
    dbprint(ppl  ,"//   ...done.");
    dbprint(ppl-1,"// " + string(rad));
    if (rad[1]==0 || NF(f,rad)==0)
    {
      dbprint(ppl  ,"//   Intersecting with K[s]...");
      v = pIntersect(var(n+1),GP);
      B[size(B)+1] = vec2poly(v,n+1);
      dbprint(ppl  ,"//   ...done.");
      dbprint(ppl-1,"// " + string(B[size(B)]));
    }
    dbprint(ppl  ,"// ...done.");
  }
  f = lcm(B); // =lcm(B[1],...,B[size(B)])
  list bb = bFactor(f);
  setring save;
  list bb = imap(Kxs,bb);
  return(bb);
}
example
{
  "EXAMPLE"; echo=2;
  ring r = 0,(x,y,Dx,Dy),dp;
  def W = Weyl();
  setring W;
  poly tx,ty = x*Dx, y*Dy;
  ideal I =      // Appel F1 with parameters (2,-3,-2,5)
    tx*(tx+ty+4)-x*(tx+ty+2)*(tx-3),
    ty*(tx+ty+4)-y*(tx+ty+2)*(ty-2),
    (x-y)*Dx*Dy+2*Dx-3*Dy;
  kill tx,ty;
  poly f = x-1;
  bfctBound(I,f);
}


//TODO check f/g or g/f, check Weyl closure of result
proc annRatSyz (poly f, poly g, list #)
"
USAGE:    annRatSyz(f,g[,db,eng]);  f, g polynomials, db,eng optional integers
ASSUME:   The Basering is commutative and over a field of characteristic 0.
RETURN:   ring (a Weyl algebra) containing an ideal 'LD', which is (part of)
          the annihilator of the rational function g/f in the corresponding
          Weyl algebra
REMARKS:  This procedure uses the computation of certain syzygies.
          One can obtain the full annihilator by computing the Weyl closure of
          the ideal LD.
NOTE:     Activate the output ring with the @code{setring} command.
          In the output ring, the ideal 'LD' (in Groebner basis) is (part of)
          the annihilator of g/f.
@*        If db>0 is given, operators of order up to db are considered,
          otherwise, and by default, a minimal holonomic solution is computed.
@*        If eng<>0, @code{std} is used for Groebner basis computations,
          otherwise, and by default, @code{slimgb} is used.
@*        If printlevel =1, progress debug messages will be printed,
          if printlevel>=2, all the debug messages will be printed.
SEE ALSO: annRat, annPoly
EXAMPLE:  example annRatSyz; shows examples
"
{
  // check assumptions
  if (!isCommutative())
  {
    ERROR("Basering must be commutative.");
  }
  if (  (size(ideal(basering)) >0) || (char(basering) >0) )
  {
    ERROR("Basering is inappropriate: characteristic>0 or qring present.");
  }
  if (g == 0)
  {
    ERROR("Second polynomial must not be zero.");
  }
  int db,eng;
  if (size(#)>0)
  {
    if (typeof(#[1]) == "int")
    {
      db = int(#[1]);
    }
    if (size(#)>1)
    {
      if (typeof(#[2]) == "int")
      {
        eng = int(#[1]);
      }
    }
  }
  int ppl = printlevel - voice + 2;
  vector I = f*gen(1)+g*gen(2);
  checkRatInput(I);
  int i,j;
  def R = basering;
  int n = nvars(R);
  list RL = ringlist(R);
  RL = RL[1..4];
  list Ltmp = RL[2];
  for (i=1; i<=n; i++)
  {
    Ltmp[i+n] = safeVarName("D" + Ltmp[i]);
  }
  RL[2] = Ltmp;
  Ltmp = list();
  Ltmp[1] = list("dp",intvec(1:2*n));
  Ltmp[2] = list("C",intvec(0));
  RL[3] = Ltmp;
  kill Ltmp;
  def @D = ring(RL);
  setring @D;
  def D = Weyl();
  setring D;
  ideal DD = 1;
  ideal Dcd,Dnd,LD,tmp;
  Dnd = 1;
  module DS;
  poly DJ;
  kill @D;
  setring R;
  module Rnd,Rcd;
  Rnd[1] = I;
  vector RJ;
  ideal L = I[1];
  module RS;
  poly p,pnew;
  pnew = I[2];
  int k,c;
  while(1)
  {
    k++;
    setring R;
    dbprint(ppl,"// Testing order: " + string(k));
    Rcd = Rnd;
    Rnd = 0;
    setring D;
    Dcd = Dnd;
    Dnd = 0;
    dbprint(ppl-1,"// Current members of the annihilator: " + string(LD));
    setring R;
    c = size(Rcd);
    p = pnew;
    for (i=1; i<=c; i++)
    {
      for (j=1; j<=n; j++)
      {
        RJ = diffRat(Rcd[i],j);
        setring D;
        DJ = Dcd[i]*var(n+j);
        tmp = Dnd,DJ;
        if (size(Dnd) <> size(simplify(tmp,4))) // new element
        {
          Dnd[size(Dnd)+1] = DJ;
          setring R;
          Rnd[size(Rnd)+1] = RJ;
          pnew = lcm(pnew,RJ[2]);
        }
        else // already have DJ in Dnd
        {
          setring R;
        }
      }
    }
    p = pnew/p;
    for (i=1; i<=size(L); i++)
    {
      L[i] = p*L[i];
    }
    for (i=1; i<=size(Rnd); i++)
    {
      L[size(L)+1] = pnew/Rnd[i][2]*Rnd[i][1];
    }
    RS = syz(L);
    setring D;
    DD = DD,Dnd;
    setring R;
    if (RS <> 0)
    {
      setring D;
      DS = imap(R,RS);
      LD = ideal(transpose(DS)*transpose(DD));
    }
    else
    {
      setring D;
    }
    LD = engine(LD,eng);
    // test if we're done
    if (db<=0)
    {
      if (isHolonomic(LD)) { break; }
    }
    else
    {
      if (k==db) { break; }
    }
  }
  export(LD);
  setring R;
  return(D);
}
example
{
  "EXAMPLE:"; echo = 2;
  // printlevel = 3;
  ring r = 0,(x,y),dp;
  poly f = 2*x*y; poly g = x^2 - y^3;
  def A = annRatSyz(f,g);   // compute a holonomic solution
  setring A; A;
  LD;
  setring r;
  def B = annRatSyz(f,g,5); // compute a solution up to degree 5
  setring B;
  LD; // this is the full annihilator as we will check below
  setring r;
  def C = annRat(f,g); setring C;
  LD; // the full annihilator
  ideal BLD = imap(B,LD);
  NF(LD,std(BLD));
}