Index: odloc.lib
===================================================================
--- dmodloc.lib	(revision 6b4d4377c3789d9b0d247c80a84385ab6cc0e892)
+++ 	(revision )
@@ -1,2681 +1,0 @@
-//
-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));
-}
-
Index: preim.lib
===================================================================
--- ncpreim.lib	(revision 6b4d4377c3789d9b0d247c80a84385ab6cc0e892)
+++ 	(revision )
@@ -1,785 +1,0 @@
-//////////////////////////////////////////////////////////////////////
-version="$Id$";
-category="Noncommutative";
-info="
-LIBRARY: ncpreim.lib    Non-commutative elimination and preimage computations
-AUTHORS: Daniel Andres, daniel.andres@math.rwth-aachen.de
-
-Support: DFG Graduiertenkolleg 1632 'Experimentelle und konstruktive Algebra'
-
-
-OVERVIEW:
-In G-algebras, elimination of variables is more involved than in the
-commutative case.
-One, not every subset of variables generates an algebra, which is again a
-G-algebra.
-Two, even if the subset of variables in question generates an admissible
-subalgebra, there might be no admissible elimination ordering, i.e. an
-elimination ordering which also satisfies the ordering condition for
-G-algebras.
-
-The difference between the procedure eliminateNC provided in this library and
-the procedure eliminate (plural) from the kernel is that eliminateNC will
-always find an admissible elimination if such one exists. Moreover, the use
-of slimgb for performing Groebner basis computations is possible.
-
-As an application of the theory of elimination, the procedure preimageNC is
-provided, which computes the preimage of an ideal under a homomorphism
-f: A -> B between G-algebras A and B. In contrast to the kernel procedure
-preimage (plural), the assumption that A is commutative is not required.
-
-
-References:
-@* (BGL) J.L. Bueso, J. Gomez-Torrecillas, F.J. Lobillo:
-         'Re-filtering and exactness of the Gelfand-Kirillov dimension',
-         Bull. Sci. math. 125, 8, 689-715, 2001.
-@* (GML) J.I. Garcia Garcia, J. Garcia Miranda, F.J. Lobillo:
-         'Elimination orderings and localization in PBW algebras',
-         Linear Algebra and its Applications 430(8-9), 2133-2148, 2009.
-@* (Lev) Viktor Levandovskyy: 'Intersection of ideals with non-commutative
-         subalgebras', ISSAC'06, 212-219, ACM, 2006.
-
-
-PROCEDURES:
-admissibleSub(v);          checks whether subalgebra is admissible
-isUpperTriangular(M,k)  ;  checks whether matrix is (strictly) upper triangular
-appendWeight2Ord(w);       appends weight to ordering
-elimWeight(v);             computes elimination weight
-eliminateNC(I,v,eng);      elimination in G-algebras
-extendedTensor(A,I);       tensor product of rings with additional relations
-preimageNC(A,f,J[,P,eng]); preimage of ideals under homomorphisms of G-algebras
-
-
-KEYWORDS: 
-preimage; elimination
-
-SEE ALSO: eliminate_lib, preimage (plural)
-";
-
-
-LIB "elim.lib";    // for nselect
-LIB "nctools.lib"; // for makeWeyl etc.
-LIB "dmodapp.lib"; // for sortIntvec
-LIB "ncalg.lib";   // for makeUgl
-LIB "dmodloc.lib"; // for commRing
-
-
-/*
-CHANGELOG
-
-
-*/
-
-
-// -- Tools ------------------------------------------------
-
-
-proc admissibleSub (intvec v)
-"
-USAGE:    admissibleSub(v);  v intvec
-ASSUME:   The entries of v are in the range 1..nvars(basering).
-RETURN:   int, 1 if the variables indexed by the entries of v form an 
-          admissible subalgebra, 0 otherwise
-EXAMPLE:  example admissibleSub; shows examples
-"
-{
-  v = checkIntvec(v);
-  int i,j;
-  list RL = ringlist(basering);
-  if (size(RL) == 4)
-  {
-    return(int(1));
-  }
-  matrix D = RL[6];
-  ideal I;
-  for (i=1; i<=size(v); i++)
-  {
-    for (j=i+1; j<=size(v); j++)
-    {
-      I[size(I)+1] = D[v[j],v[i]];
-    }
-  }
-  ideal M = maxideal(1);
-  ideal J = M[v];
-  attrib(J,"isSB",1);
-  M = NF(M,J);
-  M = simplify(M,2); // get rid of double entries in v
-  intvec opt = option(get);
-  attrib(M,"isSB",1);
-  option("redSB");
-  J = NF(I,M);
-  option(set,opt);
-  for (i=1; i<=ncols(I); i++)
-  {
-    if (J[i]<>I[i])
-    {
-      return(int(0));
-    }
-  }
-  return(int(1));
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(e,f,h),dp;
-  matrix d[3][3];
-  d[1,2] = -h; d[1,3] = 2*e; d[2,3] = -2*f;
-  def A = nc_algebra(1,d);
-  setring A; A; // A is U(sl_2)
-  // the subalgebra generated by e,f is not admissible since [e,f]=h
-  admissibleSub(1..2);
-  // but the subalgebra generated by f,h is admissible since [f,h]=2f
-  admissibleSub(2..3);
-}
-
-
-proc isUpperTriangular(matrix M, list #)
-"
-USAGE:    isUpperTriangular(M[,k]);  M a matrix, k an optional int
-RETURN:   int, 1 if the given matrix is upper triangular, 
-          0 otherwise.
-NOTE:     If k<>0 is given, it is checked whether M is strictly upper
-          triangular.
-EXAMPLE:  example isUpperTriangular; shows examples
-"
-{
-  int strict;
-  if (size(#)>0)
-  {
-    if ((typeof(#[1])=="int") || (typeof(#[1])=="number"))
-    {
-      strict = (0<>int(#[1]));
-    }
-  }
-  int m = Min(intvec(nrows(M),ncols(M)));
-  int j;
-  ideal I;
-  for (j=1; j<=m; j++)
-  {
-    I = M[j..nrows(M),j];
-    if (!strict)
-    {
-      I[1] = 0;
-    }
-    if (size(I)>0)
-    {
-      return(int(0));
-    }
-  }
-  return(int(1));
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,x,dp;
-  matrix M[2][3] = 
-    0,1,2,
-    0,0,3;
-  isUpperTriangular(M);
-  isUpperTriangular(M,1);
-  M[2,2] = 4;
-  isUpperTriangular(M);
-  isUpperTriangular(M,1);
-}
-
-
-proc appendWeight2Ord (intvec w)
-"
-USAGE:    appendWeight2Ord(w);  w an intvec
-RETURN:   ring, the basering equipped with the ordering (a(w),<), where < is
-          the ordering of the basering.
-EXAMPLE:  example appendWeight2Ord; shows examples
-"
-{
-  list RL = ringlist(basering);
-  RL[3] = insert(RL[3],list("a",w),0);
-  def A = ring(RL);
-  return(A);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  ring r = 0,(a,b,x,d),Dp;
-  intvec w = 1,2,3,4;
-  def r2 = appendWeight2Ord(w); // for a commutative ring
-  r2;
-  matrix D[4][4];
-  D[1,2] = 3*a;  D[1,4] = 3*x^2;  D[2,3] = -x;
-  D[2,4] = d;    D[3,4] = 1;
-  def A = nc_algebra(1,D);
-  setring A; A;
-  w = 2,1,1,1;
-  def B = appendWeight2Ord(w);  // for a non-commutative ring
-  setring B; B;
-}
-
-
-static proc checkIntvec (intvec v)
-"
-USAGE:    checkIntvec(v);  v intvec
-RETURN:   intvec consisting of entries of v in ascending order
-NOTE:     Purpose of this proc: check if all entries of v are in the range
-          1..nvars(basering).
-"
-{
-  if (size(v)>1)
-  {
-    v = sortIntvec(v)[1];
-  }
-  int n = nvars(basering);
-  if ( (v[1]<1) || v[size(v)]>n)
-  {
-    ERROR("Entries of intvec must be in the range 1.." + string(n));
-  }
-  return(v);
-}
-
-
-
-// -- Elimination ------------------------------------------
-
-
-/*
-// this is the same as Gweights@nctools.lib
-// 
-// proc orderingCondition (matrix D)
-// "
-// USAGE:    orderingCondition(D);  D a matrix
-// ASSUME:   The matrix D is a strictly upper triangular square matrix.
-// RETURN:   intvec, say w, such that the ordering (a(w),<), where < is
-//           any global ordering, satisfies the ordering condition for
-//           all G-algebras induced by D.
-// NOTE:     If no such ordering exists, the zero intvec is returned.
-// REMARK:   Reference: (BGL)
-// EXAMPLE:  example orderingCondition; shows examples
-// "
-// {
-//   if (ncols(D) <> nrows(D))
-//   {
-//     ERROR("Expected square matrix.");
-//   }
-//   if (isUpperTriangular(D,1)==0)
-//   {
-//     ERROR("Expected strictly upper triangular matrix.");
-//   }
-//   intvec v = 1..nvars(basering);
-//   intvec w = orderingConditionEngine(D,v,0);
-//   return(w);
-// }
-// example
-// {
-//   "EXAMPLE:"; echo = 2;
-//   // (Lev): Example 2
-//   ring r = 0,(a,b,x,d),dp;
-//   matrix D[4][4];
-//   D[1,2] = 3*a;  D[1,4] = 3*x^2;  D[2,3] = -x;
-//   D[2,4] = d;    D[3,4] = 1;
-//   // To create a G-algebra, the ordering condition implies 
-//   // that x^2<a*d must hold (see D[1,4]), which is not fulfilled:
-//   x^2 < a*d;
-//   // Hence, we look for an appropriate weight vector
-//   intwec w = orderingCondition(D); w;
-//   // and use it accordingly.
-//   ring r2 = 0,(a,b,x,d),(a(w),dp);
-//   x^2 < a*d;
-//   matrix D = imap(r,D);
-//   def A = nc_algebra(1,D);
-//   setring A; A;
-// }
-*/
-
-
-proc elimWeight (intvec v)
-"
-USAGE:    elimWeight(v);  v an intvec
-ASSUME:   The basering is a G-algebra.
-@*        The entries of v are in the range 1..nvars(basering) and the
-          corresponding variables generate an admissible subalgebra.
-RETURN:   intvec, say w, such that the ordering (a(w),<), where < is
-          any admissible global ordering, is an elimination ordering
-          for the subalgebra generated by the variables indexed by the
-          entries of the given intvec.
-NOTE:     If no such ordering exists, the zero intvec is returned.
-REMARK:   Reference: (BGL), (GML)
-EXAMPLE:  example elimWeight; shows examples
-"
-{
-  list RL = ringlist(basering);
-  if (size(RL)==4)
-  {
-    ERROR("Expected non-commutative basering.");
-  }
-  matrix D = RL[6];
-  intvec w = orderingConditionEngine(D,v,1);
-  return(w);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  // (Lev): Example 2
-  ring r = 0,(a,b,x,d),Dp;
-  matrix D[4][4];
-  D[1,2] = 3*a;  D[1,4] = 3*x^2;  D[2,3] = -x;
-  D[2,4] = d;    D[3,4] = 1;
-  def A = nc_algebra(1,D);
-  setring A; A;
-  // Since d*a-a*d = 3*x^2, any admissible ordering has to satisfy
-  // x^2 < a*d, while any elimination ordering for {x,d} additionally
-  // has to fulfil a << x and a << d.
-  // Hence neither a block ordering with weights
-  // (1,1,1,1) nor a weighted ordering with weight (0,0,1,1) will do.
-  intvec v = 3,4;
-  elimWeight(v);
-}
-
-
-static proc orderingConditionEngine (matrix D, intvec v, int elimweight)
-{
-  // algorithm from (BGL) and (GML), respectively
-  // solving an LPP via simplex
-  int ppl = printlevel - voice + 1;
-  def save = basering;
-  int n = nvars(save);
-  ideal EV = maxideal(1);
-  EV = EV[v]; // also assumption check for v
-  attrib(EV,"isSB",1);
-  ideal NEV = maxideal(1);
-  NEV = NF(NEV,EV);
-  intmat V1[n-size(NEV)][n+1];
-  if (elimweight)
-  {
-    intmat V2[size(NEV)][n+1];
-  }
-  int rowV1,rowV2;
-  intmat M[1][n];
-  intmat M2,oldM;
-  int i,j,k;
-  for (i=1; i<=n; i++)
-  {
-    if (elimweight)
-    {
-      if (NEV[i]<>0)
-      {
-        V2[rowV2+1,i+1] = 1; // xj == 0
-        rowV2++;
-      }
-      else
-      {
-        V1[rowV1+1,1] = 1; // 1-xi <= 0
-        V1[rowV1+1,i+1] = -1;
-        rowV1++;
-      }
-    }
-    else
-    {
-      V1[i,1] = 1; // 1-xi <= 0
-      V1[i,i+1] = -1;
-      rowV1++;
-    }
-    for (j=i+1; j<=n; j++)
-    {
-      if (deg(D[i,j])>0)
-      {
-        M2 = newtonDiag(D[i,j]);
-        for (k=1; k<=nrows(M2); k++)
-        {
-          M2[k,i] = M2[k,i] - 1; // <beta,x> >= 0
-          M2[k,j] = M2[k,j] - 1;
-        }
-        oldM = M;
-        M = intmat(M,nrows(M)+nrows(M2),n);
-        M = oldM,M2;
-      }
-    }
-  }
-  intvec eq = 0,(-1:n);
-  ring r = 0,x,dp; // to avoid problems with pars or char>0
-  module MM = module(transpose(matrix(M)));
-  MM = simplify(MM,2+4);
-  matrix A;
-  if (MM[1]<>0)
-  {
-    if (elimweight)
-    {
-      MM = 0,transpose(MM);
-    }
-    else
-    {
-      MM = module(matrix(1:ncols(MM)))[1],transpose(MM);
-    }
-    A = transpose(concat(matrix(eq),transpose(-MM)));
-  }
-  else
-  {
-    A = transpose(eq);
-  }
-  A = transpose(concat(transpose(A),matrix(transpose(V1))));
-  if (elimweight)
-  {
-    A = transpose(concat(transpose(A),matrix(transpose(V2))));
-  }
-  int m = nrows(A)-1;
-  ring realr = (real,10),x,lp;
-  matrix A = imap(r,A);
-  dbprint(ppl,"// Calling simplex...");
-  dbprint(ppl-1,"// with the matrix " + print(A));
-  dbprint(ppl-1,"// and parameters " 
-          + string(intvec(m,n,m-rowV1-rowV2,rowV1,rowV2)));
-  list L = simplex(A,m,n,m-rowV1-rowV2,rowV1,rowV2);
-  int se = L[2];
-  if (se==-2)
-  {
-    ERROR("simplex yielded an error. Please inform the authors.");
-  }
-  intvec w = 0:n;
-  if (se==0)
-  {
-    matrix S = L[1];
-    intvec s = L[3];
-    for (i=2; i<=nrows(S); i++)
-    {
-      if (s[i-1]<=n)
-      {
-        w[s[i-1]] = int(S[i,1]);
-      }
-    }
-  }
-  setring save;
-  return(w);
-}
-
-
-proc eliminateNC (ideal I, intvec v, list #)
-"
-USAGE:    eliminateNC(I,v,eng);  I ideal, v intvec, eng optional int
-RETURN:   ideal, I intersected with the subring defined by the variables not
-          index by the entries of v
-ASSUME:   The entries of v are in the range 1..nvars(basering) and the
-          corresponding variables generate an admissible subalgebra.
-REMARKS:  In order to determine the required elimination ordering, a linear
-          programming problem is solved with the simplex algorithm.
-@*        Reference: (GML)
-@*        Unlike eliminate, this procedure will always find an elimination
-          ordering, if such exists.
-NOTE:     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: eliminate (plural)
-EXAMPLE:  example eliminateNC; shows examples
-"
-{
-  int ppl = printlevel - voice + 2;
-  v = checkIntvec(v);
-  if (!admissibleSub(v))
-  {
-    ERROR("Subalgebra is not admissible: no elimination is possible.");
-  }
-  dbprint(ppl,"// Subalgebra is admissible.");
-  int eng;
-  if (size(#)>0)
-  {
-    if (typeof(#[1])=="int" || typeof(#[1])=="number")
-    {
-      eng = int(#[1]);
-    }
-  }
-  def save = basering;
-  int n = nvars(save);
-  dbprint(ppl,"// Computing elimination weight...");
-  intvec w = elimWeight(v);
-  if (w==(0:n))
-  {
-    ERROR("No elimination ordering exists.");
-  }
-  dbprint(ppl,"// ...done.");
-  dbprint(ppl-1,"// Using elimination weight " + string(w) + ".");
-  def r = appendWeight2Ord(w);
-  setring r;
-  ideal I = imap(save,I);
-  dbprint(ppl,"// Computing Groebner basis with engine " + string(eng)+"...");
-  I = engine(I,eng);
-  dbprint(ppl,"// ...done.");
-  dbprint(ppl-1,string(I));
-  I = nselect(I,v);
-  setring save;
-  I = imap(r,I);
-  return(I);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  // (Lev): Example 2
-  ring r = 0,(a,b,x,d),Dp;
-  matrix D[4][4];
-  D[1,2] = 3*a; D[1,4] = 3*x^2;
-  D[2,3] = -x;  D[2,4] = d;     D[3,4] = 1;
-  def A = nc_algebra(1,D);
-  setring A; A;
-  ideal I = a,x;
-  // Since d*a-a*d = 3*x^2, any admissible ordering has to satisfy
-  // x^2 < a*d, while any elimination ordering for {x,d} additionally
-  // has to fulfil a << x and a << d.
-  // Hence, the weight (0,0,1,1) is not an elimination weight for
-  // (x,d) and the call eliminate(I,x*d); will produce an error.
-  eliminateNC(I,3..4);
-  // This call uses the elimination weight (0,0,1,2), which works.
-}
-
-
-
-// -- Preimages ------------------------------------------------
-
-// TODO A or B commutative
-proc extendedTensor(def A, ideal I)
-"
-USAGE:    extendedTensor(A,I);  A ring, I ideal
-RETURN:   ring, A+B (where B denotes the basering) extended with non-
-          commutative relations between the vars of A and B, which arise from
-          the homomorphism A -> B induced by I in the usual sense, i.e. if the
-          vars of A are named x(i) and the vars of B y(j), then putting
-          q(i)(j) = leadcoef(y(j)*I[i])/leadcoef(I[i]*y(j)) and
-          r(i)(j) = y(j)*I[i] - q(i)(j)*I[i]*y(j) yields the relation
-          y(j)*x(i) = q(i)(j)*x(i)*y(j)+r(i)(j).
-REMARK:   Reference: (Lev)
-EXAMPLE:  example extendedTensor; shows examples
-"
-{
-  def B = basering;
-  setring A;
-  int nA = nvars(A);
-  string varA = "," + charstr(A) + "," + varstr(A) + ",";
-  setring B;
-  int nB = nvars(B);
-  list RL = ringlist(B);
-  list L = RL[2];
-  string vB;
-  int i,j;
-  for (i=1; i<=nB; i++)
-  {
-    vB = "," + L[i] + ",";
-    while (find(varA,vB)<>0)
-    {
-      vB[1] = "@";
-      vB = "," + vB;
-    }
-    vB = vB[2..size(vB)-1];
-    L[i] = vB;
-  }
-  RL[2] = L;
-  def @B = ring(RL);
-  kill L,RL;
-  setring @B;
-  ideal I = fetch(B,I);
-  def E = A+@B;
-  setring E;
-  ideal I = imap(@B,I);
-  matrix C = ringlist(E)[5];
-  matrix D = ringlist(E)[6];
-  poly p,q;
-  for (i=1; i<=nA; i++)
-  {
-    for (j=nA+1; j<=nA+nB; j++)
-    {
-      // upper right block: new relations
-      p = var(j)*I[i];
-      q = I[i]*var(j);
-      C[i,j] = leadcoef(p)/leadcoef(q);
-      D[i,j] = p - C[i,j]*q;
-    }
-  }
-  def @EE = commRing();
-  setring @EE;
-  matrix C = imap(E,C);
-  matrix D = imap(E,D);
-  def EE = nc_algebra(C,D);
-  setring B;
-  return(EE);
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  def A = makeWeyl(2);
-  setring A; A;
-  def B = makeUgl(2);
-  setring B; B;
-  ideal I = var(1)*var(3), var(1)*var(4), var(2)*var(3), var(2)*var(4);
-  I;
-  def C = extendedTensor(A,I);
-  setring C; C;
-}
-
-
-proc preimageNC (list #)
-"
-USAGE:    preimageNC(A,f,J[,P,eng]);  A ring, f map or ideal, J ideal, 
-                                      P optional string, eng optional int
-ASSUME:   f defines a map from A to the basering.
-RETURN:   nothing, instead exports an object 'preim' of type ideal to ring A,
-          being the preimage of J under f.
-NOTE:     If P is given and not equal to the empty string, the preimage is
-          exported to A under the name specified by P.
-          Otherwise (and by default), P is set to 'preim'.
-@*        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.
-REMARK:   Reference: (Lev)
-SEE ALSO: preimage (plural)
-EXAMPLE:  example preimageNC; shows examples
-"
-{
-  int ppl = printlevel - voice + 2;
-  if (size(#) <3)
-  {
-    ERROR("Expected 3 arguments.")
-  }
-  def B = basering;
-  if (typeof(#[1])<>"ring")
-  {
-    ERROR("First argument must be a ring.");
-  }
-  def A = #[1];
-  setring A;
-  ideal mm = maxideal(1);
-  setring B;
-  if (typeof(#[2])=="map")
-  {
-    map phi = #[2];
-  }
-  else
-  {
-    if (typeof(#[2])=="ideal")
-    {
-      map phi = A,#[2];
-    }
-    else
-    {
-      ERROR("Second argument must define a map from the specified ring to the basering.");
-    }
-  }
-  if (typeof(#[3])<>"ideal")
-  {
-    ERROR("Third argument must be an ideal in the specified ring");
-  }
-  ideal J = #[3];
-  string str = "preim";
-  int eng;
-  if (size(#)>3)
-  {
-    if (typeof(#[4])=="string")
-    {
-      if (#[4]<>"")
-      {
-        str = #[4];
-      }
-    }
-    if (size(#)>4)
-    {
-      if (typeof(#[5])=="int")
-      {
-        eng = #[5];
-      }
-    }
-  }
-  setring B;
-  ideal I = phi(mm);
-  def E = extendedTensor(A,I);
-  setring E;
-  dbprint(ppl,"// Computing in ring");
-  dbprint(ppl,E);
-  int nA = nvars(A);
-  int nB = nvars(B);
-  ideal @B2E = maxideal(1);
-  @B2E = @B2E[(nA+1)..(nA+nB)];
-  map B2E = B,@B2E;
-  ideal I = B2E(I);
-  ideal Iphi;
-  int i,j;
-  for (i=1; i<=nA; i++)
-  {
-    Iphi[size(Iphi)+1] = var(i) - I[i];
-  }
-  dbprint(ppl,"// I_{phi} is  " + string(Iphi));
-  ideal J = imap(B,J);
-  J = J + Iphi;
-  intvec v = (nA+1)..(nA+nB);
-  dbprint(ppl,"// Starting elimination...");
-  dbprint(ppl-1,string(J));
-  J = eliminateNC(J,v,eng);
-  dbprint(ppl,"// ...done.");
-  dbprint(ppl-1,string(J));
-  J = nselect(J,v);
-  attrib(J,"isSB",1);
-  setring A;
-  dbprint(ppl,"// Writing output to specified ring under the name '" 
-          + str + "'.");
-  str = "ideal " + str + " = imap(E,J); export(" + str + ");";
-  execute(str);
-  setring B;
-  return();
-}
-example
-{
-  "EXAMPLE:"; echo = 2;
-  def A = makeUgl(3); setring A; A; // universal enveloping algebra of gl_3
-  ring r3 = 0,(x,y,z,Dx,Dy,Dz),dp;
-  def B = Weyl(); setring B; B;     // third Weyl algebra
-  ideal ff = x*Dx,x*Dy,x*Dz,y*Dx,y*Dy,y*Dz,z*Dx,z*Dy,z*Dz;
-  map f = A,ff;                     // f: A -> B, e(i,j) |-> x(i)D(j)
-  ideal J = 0;
-  preimageNC(A,f,J,"K");            // compute K := ker(f)
-  setring A;
-  K;
-}
-
-
-// -- Examples ---------------------------------------------
-
-static proc ex1 ()
-{
-  ring r1 = 0,(a,b),dp;
-  int t = 7;
-  def St = nc_algebra(1,t*a);
-  ring r2 = 0,(x,D),dp;
-  def W = nc_algebra(1,1); // W is the first Weyl algebra
-  setring W;
-  map psit = St, x^t,x*D+t;
-  int p = 3;
-  ideal Ip = x^p, x*D+p;
-  preimageNC(St,psit,Ip);
-  setring St; preim;
-}
-
-
-static proc ex2 ()
-{
-  ring r1 = 0,(e,f,h),dp;
-  matrix D1[3][3]; D1[1,2] = -h; D1[1,3] = 2*e; D1[2,3] = -2*f;
-  def U = nc_algebra(1,D1); // D is U(sl_2)
-  ring r2 = 0,(x,D),dp;
-  def W = nc_algebra(1,1); // W is the first Weyl algebra
-  setring W;
-  ideal tau = x,-x*D^2,2*x*D;
-  def E = extendedTensor(U,tau);
-  setring E; E;
-  elimWeight(4..5); 
-  // zero, since there is no elimination ordering for x,D in E
-}
-
-
-static proc ex3 ()
-{
-  ring r1 = 0,(x,d,s),dp;
-  matrix D1[3][3]; D1[1,2] = 1;
-  def A = nc_algebra(1,D1);
-  ring r2 = 0,(X,DX,T,DT),dp;
-  matrix D2[4][4]; D2[1,2] = 1; D2[3,4] = 1;
-  def B = nc_algebra(1,D2);
-  setring B;
-  map phi = A, X,DX,-DT*T;
-  ideal J = T-X^2, DX+2*X*DT;
-  preimageNC(A,phi,J);
-  setring A;
-  preim;
-}
