singular-data 1:4.1.0-p3+ds-2build1 / usr / share / singular / LIB / ncpreim.lib

/////////////////////////////////////////////////////////////////////
version="version ncpreim.lib 4.0.0.0 Jun_2013 "; // $Id: fcca2c37ec1caedea6c6ad3cd032dcf4bedf44e7 $
category="Noncommutative";
info="
LIBRARY: ncpreim.lib    Non-commutative elimination and preimage computations
AUTHOR:  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 @code{eliminateNC} provided in this
library and the procedure @code{eliminate (plural)} from the kernel is that
eliminateNC will always find an admissible elimination if such one exists.
Moreover, the use of @code{slimgb} for performing Groebner basis computations
is possible.

As an application of the theory of elimination, the procedure @code{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
@code{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) V. Levandovskyy: `Intersection of ideals with non-commutative
         subalgebras', ISSAC'06, 212-219, ACM, 2006.


PROCEDURES:
eliminateNC(I,v,eng);      elimination in G-algebras
preimageNC(A,f,J[,P,eng]); preimage of ideals under homomorphisms of G-algebras
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
extendedTensor(A,I);       tensor product of rings with additional relations


KEYWORDS: preimage; elimination


SEE ALSO: elim_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
11.12.12: docu, typos, fixed variable names in extendedTensor,
 moved commRing to dmodloc.lib
12.12.12: typos
17.12.12: docu
24.09.13: bugfix preimageNC naming conflict if f is map from ring called 'B'
*/


// -- Testing for consistency of the library ---------------

static proc testncpreimlib()
{
  example admissibleSub;
  example isUpperTriangular;
  example appendWeight2Ord;
  example elimWeight;
  example eliminateNC;
  example extendedTensor;
  example preimageNC;
}


// -- 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" || typeof(#[2])=="ideal")
  {
    map phi = A,ideal(#[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;
}