///////////////////////////////////////////////////////////////////////////////
version="$Id: nctools.lib,v 1.10 2005-02-23 18:10:46 levandov Exp $";
category="Noncommutative";
info="
LIBRARY: nctools.lib     General tools for noncommutative algebras
AUTHORS:   Levandovskyy V.,     levandov@mathematik.uni-kl.de,
@*         Lobillo, F.J.,       jlobillo@ugr.es,
@*	   Rabelo, C.,          crabelo@ugr.es.

SUPPORT: DFG (Deutsche Forschungsgesellschaft) and Metodos algebraicos y efectivos en grupos cuanticos, BFM2001-3141, MCYT, Jose Gomez-Torrecillas (Main researcher).

PROCEDURES:
Gweights(r);		computes weights for a compatible ordering in a G-algebra,
weightedRing(r);	changes the ordering of a ring to a weighted one,
ndc([v]);               the ideal of non-degeneracy conditions in G-algebra,
RootOfUnity(n);         the minimal polynomial for the n-th primitive root of unity,
Weyl([p]);              creates Weyl algebra structure in a basering (two different realizations),
CreateWeyl(n, [p]);     returns n-th Weyl algebra in (x(i),D(i)) presentation; in char p, if an integer p is given,
Heisenberg(N, [p,d]);   returns N-th  Heisenberg algebra in (x(i),y(i),h) realization,
Exterior();             returns qring, the exterior algebra of a basering,
Fin_dim_algebra(matrix M, list #); creates finite dimensional algebra structure from the basering and the multiplication matrix M,

AUXILIARY PROCEDURES:
Newton(f);  		Newton diagram of a polynomial f,
NCRelations(r);		recovers the non-commutative relations of a G-algebra,
IsCentral(p,[v]);       check for the commutativity of polynomial p with the G-algebra,
Is_NC();                checks whether basering is noncommutative,
UpOneMatrix(N);         returns NxN matrix with 1's in the whole upper triagle,

ALIAS PROCEDURES:
wRing(r);		alias to weightedRing,
";

///////////////////////////////////////////////////////////////////////////////

// This procedure computes a weights vector for a G-Algebra
// r must be a G-algebra

proc Gweights(def r)
"USAGE:   Gweights(r); r a ring or a square matrix
RETURN:   a weight vector for the G-algebra:
          r itself, if it is of the type ring,
          or for a G-algebra, defined by the square matrix r
NOTE:    with Gweights you obtain a vector, which must be used to redefine the G-Algebra. If the input is a matrix and the output is the zero vector then there is not a G-algebra structure associated to these relations with respect to the given variables. Another possibility is to use wRing or weightedRing to obtain directly the G-Algebra.
EXAMPLE: example Gweights; shows examples
SEE ALSO: wRing, weightedRing
"{
  matrix tails;
  int novalid=0;
  if (typeof(r)=="ring") //a ring is admissible as input
  {
    setring r;
    def l = NCRelations(r);
    tails = l[2]; // l=C,D we need D, the tails of the relations
  }
  else
  {
    if ( (typeof(r)=="matrix") || (typeof(r)=="intmat") )
    {
      if ( nrows(r)==ncols(r) ) //the input is a square matrix
      {
        tails = matrix(r);
      }
      else
      {
        novalid = 1;
      }
    }
    else
    {
      novalid=1;
    }
  }
  if (novalid==0)
  {
    intmat IM = SimplMat(tails);
    if ( size(IM)>1 )
    {
      int n  = ncols(tails);
      int m  = nrows(IM)-1;
      int m1 = 0;
      int m2 = m;
      int m3 = 0;
      ring simplexring=(real,10),(x),lp;// The simplex procedure requires a basering of this type
      matrix M = IM;
      list sol = simplex (M,m,n,m1,m2,m3);
      return(weightvector(sol));
    }
    else
    {
      "Invalid input"; //usually because the input is a one variable ring
      return();
    }
  }
  else
  {
    "The input must be a ring or a square matrix";
    return();
  }
}
example
{
  "EXAMPLE:";echo=2;
  ring r = (0,q),(a,b,c,d),lp;
  matrix C[4][4]; 
  C[1,2]=q; C[1,3]=q; C[1,4]=1; C[2,3]=1; C[2,4]=q; C[3,4]=q;
  matrix D[4][4]; 
  D[1,4]=(q-1/q)*b*c;
  ncalgebra(C,D);
  setring r; r;
  Gweights(r);
  Gweights(D);
}

///////////////////////////////////////////////////////////////////////////////

// This procedure take a ring r, call to Gweights(r) and use the output 
// of Gweights(r) to make a change of order in r
// The output is a new ring, equal to r but the order
// r must be a G-algebra

proc weightedRing(def r)
"USAGE:   weightedRing(r); r a ring
RETURN:  ring with the relations of r with order changed to comply with
the ordering condition for G-algebras
NOTE:    you have to activate this ring with the "setring" command
EXAMPLE: example weightedRing; shows examples
SEE ALSO: wRing, Gweights
"{
  def wv=Gweights(r);
  if (typeof(wv)=="intvec")
  {
    setring r;
    int n=nvars(r);
    // Generating an nxn-intmat order
    intmat m[n][n];
    m[1,1]=wv[1];
    int i;
    for (i=2; i<=n; i++)
    {
      m[1,i]=wv[i];
      m[i,n+2-i]=1;
    }
    // End of generation.
    def lr=NCRelations(r);
    string newringstring="ring newring=("+charstr(r)+"),("+varstr(r)+"),M("+string(m)+")";
    execute (newringstring);
    def lnewring=imap(r,lr);
    ncalgebra(lnewring[1],lnewring[2]);
    return(newring);
  }
  else
  {
    "Invalid input.";//usually because the input is a one variable ring
    return();
  }
}
example
{
  "EXAMPLE:";echo=2;
  ring r = (0,q),(a,b,c,d),lp;
  matrix C[4][4]; 
  C[1,2]=q; C[1,3]=q; C[1,4]=1; C[2,3]=1; C[2,4]=q; C[3,4]=q;
  matrix D[4][4]; 
  D[1,4]=(q-1/q)*b*c;
  ncalgebra(C,D);
  setring r; r;
  def t=weightedRing(r);
  setring t; t;
}

///////////////////////////////////////////////////////////////////////////////

proc wRing(def r)
"USAGE:   see weightedRing
"{
  def a=weightedRing(r);
  if (typeof(a)=="ring") { return(a); }
  "Error";
  return();
}
example
{
  "EXAMPLE:";echo=2;
  LIB "qmatrix.lib";
  def r=quant(2); // generate quant(2) and store it in r
  setring r; // set the ring r the active ring
  r;
  def s=wRing(r);
  setring s;
  s;
}

///////////////////////////////////////////////////////////////////////////////

// This procedure computes ei+ej-f with f running in Newton(pij) and deletes the zero rows

static proc Cij(intmat M, int i,j)
{
  M=(-1)*M;
  int nc=ncols(M);
  intvec N;
  int k;
  for (k=1; k<=nrows(M); k++)
  {
    M[k,i]=M[k,i]+1;
    M[k,j]=M[k,j]+1;
    if (intvec(M[k,1..nc])!=0)
    {
      N=N,intvec(M[k,1..nc]);
    } // we only want non-zero rows
  }
  if (size(N)>1)
  {
    N=N[2..size(N)]; // Deleting the zero added in the definition of N
    M=intmat(N,size(N)/nc,nc); // Conversion from vector to matrix
  }
  else
  {
    intmat M[1][1]=0;
  }
  return (M);
}

///////////////////////////////////////////////////////////////////////////////

// This procedure run over the matrix of pij calculating Cij

static proc Ct(matrix P)
{
  int    k = ncols(P);
  intvec T = 0;
  int    i,j;
//  int notails=1;
  def S;
  for (j=2; j<=k; j++)
  {
    for (i=1; i<j; i++)
    {
      if ( P[i,j] != 0 )
      {
//        notails=0;
        S = Newton(P[i,j]);
        S = Cij(S,i,j);
        if ( size(S)>1 )
        {
          T = T,S;
        }
      }
    }
  }
  if ( size(T)==1 )
  {
    intmat C[1][1] = 0;
  }
  else
  {
    T=T[2..size(T)]; // Deleting the zero added in the definition of T
    intmat C = intmat(T,size(T)/k,k); // Conversion from vector to matrix
  }
  return (C);
}

///////////////////////////////////////////////////////////////////////////////

// The purpose of this procedure is to produce the input matrix required by simplex procedure

static proc SimplMat(matrix P)
{
  intmat C=Ct(P);
  if (size(C)>1)
  {
    int r = nrows(C);
    int n = ncols(C);
    int f = 1+n+r;
    intmat M[f][n+1]=0;
    int i;
    for (i=2; i<=(n+1); i++) 
    {
      M[1,i]=-1; // (0,-1,-1,-1,...) objective function in the first row
    }
    for (i=2; i<=f; i++) {M[i,1]=1;} // All the independent terms are 1
    for (i=2; i<=(n+1); i++) {M[i,i]=-1;} // wi>=1 is an identity matrix
    M[(n+2)..f,2..(n+1)]=(-1)*intvec(C); // <wi,a> >= 1, a in C ...
  }
  else
  {
    int n = ncols(P);
    int f = 1+n;
    intmat M[f][n+1]=0;
    int i;
    for (i=2; i<=(n+1); i++) {M[1,i]=-1;} // (0,-1,-1,-1,...) objective function in the first row
    for (i=2; i<=f; i++) {M[i,1]=1;} // All the independent terms are 1
    for (i=2; i<=(n+1); i++) {M[i,i]=-1;} // wi>=1 is an identity matrix
  }
  return (M);
}

///////////////////////////////////////////////////////////////////////////////

// This procedure generates a nice output of the simplex method consisting of a vector
// with the solutions. The vector is ordered.

static proc weightvector(list l)
"ASSUME:  l is the output of simplex.
RETURNS: if there is a solution, an intvec with the solution."
{
  matrix m=l[1];
  intvec nv=l[3];
  int sol=l[2];
  int rows=nrows(m);
  int N=l[6];
  intmat wv[1][N]=0;
  int i;
  if (sol)
  {
    "no solution satisfies the given constraints";
  }
  else
  {
    for ( i = 2; i <= rows; i++ )
    {
      if ( nv[i-1] <= N )
      {
	wv[1,nv[i-1]]=int(m[i,1]);
      }
    }
  }
  return (intvec(wv));
}

///////////////////////////////////////////////////////////////////////////////

// This procedure calculates the Newton diagram of the polynomial f
// The output is a intmat M, each row of M is the exp of a monomial in f

proc Newton(poly f)
"USAGE:   Newton(f); f a poly
RETURN:  intmat, representing the Newton diagram of f
NOTE:    each row is the exponent of a monomial of f
EXAMPLE: example Newton; shows examples
"{
  int n=nvars(basering);
  intvec N=0;
  if ( f != 0 )
  {
    while ( f != 0 )
    {
      N = N, leadexp(f);
      f = f-lead(f);
    }
  }
  else
  {
    N=N, leadexp(f);
  }
  N = N[2..size(N)]; // Deletes the zero added in the definition of T
  intmat M=intmat(N,(size(N)/n),n); // Conversion from vector to matrix
  return (M);
}
example
{
  "EXAMPLE:";echo=2;
  ring r = 0,(x,y,z),lp;
  poly f=x2y+3xz-5y+3;
  Newton(f);
}

///////////////////////////////////////////////////////////////////////////////

// This procedure recover the non-conmutative relations (matrices C and D)

proc NCRelations(def r)
"USAGE:   NCRelations(r); r a ring
RETURN:  a list with two elements, both elements are of type matrix and
	 represent the matrices C,D defining the non-commutative relations
	 of the G-algebra r
EXAMPLE: example NCRelations; shows examples
"{
  list l;
  if (typeof(r)=="ring")
  {
    int n=nvars(r);
    matrix C[n][n]=0;
    matrix D[n][n]=0;
    poly f; poly g;
    if (n>1)
    {
      int i,j;
      for (i=2; i<=n; i++)
      {
	for (j=1; j<i; j++)
	{
	  f=var(i)*var(j); // yx=c*xy+...
	  g=var(j)*var(i); // xy
	  while (C[j,i]==0)
	  {
	    if (leadmonom(f)==leadmonom(g))
	    {
	      C[j,i]=leadcoef(f);
	      D[j,i]=D[j,i]+f-lead(f);
	    }
	    else
	    {
	      D[j,i]=D[j,i]+lead(f);
	      f=f-lead(f);
	    }
	  }
	}
      }
      l=C,D;
    }
    else { "The ring must have two or more variables"; }
  }
  else { "The input must be of a type ring";}
  return (l);
}
example
{
  "EXAMPLE:";echo=2;
  ring r = 0,(x,y,z),dp;
  matrix C[3][3]=0,1,2,0,0,-1,0,0,0;
  print(C);
  matrix D[3][3]=0,1,2y,0,0,-2x+y+1;
  print(D);
  ncalgebra(C,D);
  r;
  def l=NCRelations(r);
  print (l[1]);
  print (l[2]);
}

///////////////////////////////////////////////////////////////////////////////

proc Fin_dim_algebra(matrix M, list #)
"USAGE:   Fin_dim_algebra(M,[r]); M a matrix, r an optional ring.
RETURN:  nothing. Creates finite dimensional algebra structure in a ring r
(if it is given) or in a basering (by default) from the matrix M. Exports the ideal called Quot for further qring definition. 
NOTE: matrix M will be read according to relations x_j*x_i = M[i,j]; Quot is not given in its two-sided Groebner basis.
EXAMPLE: Fin_dim_algebra; shows examples
"
{
  if (size(#) >0)
  {
    if ( typeof(#[1])!="ring" ) { return();}
    else
    {
      def @R1 = #[1];
      setring @R1;
    }
  }
  int i,j;
  int n=nvars(basering);
  poly p;
  ideal I;
  number c;
  matrix C[n][n];
  matrix D[n][n];
  for (i=1; i<=n; i++)
  {
    for (j=i; j<=n; j++)
    {
      p=var(i)*var(j)-M[i,j];
      if ( (size(I)==1) && (I[1]==0) )   { I=p; }
      else { I=I,p; }
      if (j>i)
      {
	if ((M[i,j]!=0) && (M[j,i]!=0))
	{
	  c = leadcoef(M[j,i])/leadcoef(M[i,j]);
	}
	else
	{
	  c = 1;
	}
	C[i,j]=c;
        D[i,j]= - M[j,i] +c*M[i,j];
      }
    }
  }
  ncalgebra(C,D);
  ideal Quot = I;
  export Quot;
}
example
{
  "EXAMPLE:";echo=2;
  ring r=(0,a,b),(x(1..3)),dp;
  matrix S[3][3];
  S[2,3]=a*x(1); S[3,2]=-b*x(1);
  Fin_dim_algebra(S);
  Quot = twostd(Quot);
  qring Qr = Quot;
  Qr;
}

///////////////////////////////////////////////////////////////////////////////

proc IsCentral(poly p, list #)
"USAGE:   IsCentral(p,[v]); p poly, v an integer (with v!=0 procedure will be verbose)
RETURN:  integer (1 if p commutes with all variables, 0 otherwise)
EXAMPLE: example IsCentral; shows examples
"{
  int N = nvars(basering);
  int in;
  int flag = 1;
  poly   q = 0;
  for (in=1; in<=N; in++)
  {
    q = p*var(in)-var(in)*p;
    if (q!=0)
    {
      if (size(#) >0 )
      {
	"Noncentral at:", var(in);
      }
      flag = 0;
    }
  }
  return(flag);
}
example
{
  "EXAMPLE:";echo=2;
  ring r=0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=-z;
  D[1,3]=2*x;
  D[2,3]=-2*y;
  ncalgebra(1,D); // this is U(sl_2)
  poly c = 4*x*y+z^2-2*z;
  IsCentral(c,0);
  poly h = x*c;
  IsCentral(h,1);
}

///////////////////////////////////////////////////////////////////////////////

proc UpOneMatrix(int N)
"USAGE:   UpOneMatrix(N); N an integer, the number of columns
RETURN:  intmat, NxN matrix with 1's in the whole upper triagle
NOTE: helpful for setting noncommutative algebras with complicated
coefficient matrices
EXAMPLE: example UpOneMatrix; shows examples
"{
  int ii,jj;
  intmat U[N][N]=0;
  for (ii=1;ii<N;ii++)
  {
    for (jj=ii+1;jj<=N;jj++)
    {
      U[ii,jj]=1;
    }
  }
  return(U);
}
example
{
  "EXAMPLE:";echo=2;
  ring r=(0,q),(x,y,z),dp;
  matrix C=UpOneMatrix(3);
  C[1,3]=q;
  ncalgebra(C,0);
  r;
}

///////////////////////////////////////////////////////////////////////////////
proc ndc(list #)
"USAGE:   ndc([v]); v an optional integer. If v!=0 procedure will be verbose
RETURN:  ideal of non-degeneracy conditions of the basering
EXAMPLE: example ndc; shows examples
"
{
  // if the second argument is given, produces ndc wrt powers x^N
  int N = 1;
  int Verbose = 0;
  if ( size(#)>=1 ) { Verbose = int(#[1]); }
  if ( size(#)>=2 ) { N = int(#[2]); }
  int cnt = 1;
  int numvars = nvars(basering);
  int a,b,c;
  poly p = 1;
  ideal res = 0;
  for (cnt=1; cnt<=N; cnt++)
  {
    if (Verbose) { "Processing degree :",cnt;}
    for (a=1; a<=numvars-2; a++)
    {
      for (b=a+1; b<=numvars-1; b++)
      {
	for(c=b+1; c<=numvars; c++)
	{
	  p = (var(c)^cnt)*(var(b)^cnt);
	  p = p*(var(a)^cnt);
	  p = p-(var(c)^cnt)*((var(b)^cnt)*(var(a)^cnt));
	  if (Verbose) {a,".",b,".",c,".";}
	  if (p!=0)
	  { 
	    if ( res==0 )
	    {
	      res[1] = p;
	    }
	    else
	    {
	      res = res,p;
	    }
	    if (Verbose) { "failed:",p; }
	  }
	}
      }
    }
    if (Verbose) { "done"; }
  }
  return(res);
}
example
{
  "EXAMPLE:";echo=2;
  ring r = (0,q1,q2),(x,y,z),dp;
  matrix C[3][3];
  C[1,2]=q2; C[1,3]=q1; C[2,3]=1;
  matrix D[3][3];
  D[1,2]=x; D[1,3]=z;
  ncalgebra(C,D);
  r;
  ideal j=ndc(); // the silent version
  j;
  ideal i=ndc(1); // the verbose version
  i;
}

///////////////////////////////////////////////////////////////////////////////
proc RootOfUnity(int n)
"USAGE:   RootOfUnity(n); n an integer
RETURN:  number, the n-th primitive root of unity (for use as minpoly)
NOTE: works only in field extensions by one element
EXAMPLE: example RootOfUnity; shows examples
"
{
  if ( npars(basering) !=1 )
  {
    "the procedure works only with one parameter";
    return(0);
  }
  if (n<1) { return(0); }
  number mp = par(1);
  if (n==1) { return(mp-1); }
  if (n==2) { return(mp+1); }
  def OldRing = basering;
  string CH = charstr(basering);
  string MCH;
  int j=1;
  while ( (CH[j] !=",") && (j<=size(CH)))
  {
    MCH=MCH+CH[j]; j++;
  }
  string SR = "ring @@rR="+MCH+","+parstr(basering)+",dp;";
  execute(SR);
  poly @t=var(1)^n-1; // (x^2i-1)=(x^i-1)(x^i+1)
  list l=factorize(@t);
  ideal @l=l[1];
  list @d;
  int s=size(@l);
  int d=deg(@l[s]);
  int cnt=1;
  poly res;
  for (j=s-1; j>=1; j--)
  {
    if ( deg(@l[j]) > d) { d=deg(@l[j]); }
  }
  for (j=1; j<=s; j++)
  {
    if ( deg(@l[j]) == d) { @d[cnt]=@l[j]; cnt++; }
  }
  if ( size(@d)==1 )
  {
    res = poly(@d[1]);
  }
  else
  {
    j=1;
    while  ( j <= size(@d) )
    {
      res = @d[j]-lead(@d[j]);
      if ( leadcoef(res) >=0 ) { j++; }
      else { break; }
    }
    res = @d[j];
  }
  setring OldRing;
  poly I = imap(@@rR,res);
  mp = leadcoef(I);
  kill @@rR;
  return(mp);
}
example
{
  "EXAMPLE:";echo=2;
  ring r8 = (0,q),(x,y,z),dp;
  minpoly = RootOfUnity(8);
  r8;
  ring r7 = (0,q),(x,y,z),dp;
  minpoly = RootOfUnity(7);
  r7;
  ring r6 = (0,q),(x,y,z),dp;
  minpoly = RootOfUnity(6);
  r6;
}

///////////////////////////////////////////////////////////////////////////////
proc Weyl(list #)
"USAGE:   Weyl([p]); p an optional integer.
RETURN:  nothing. Creates Weyl algebra structure in a basering. By default
mimics (x(1..N),d(1..N)) realization. If p is given and is not zero,
uses (x(1),d(1),x(2),d(2),... ) realization.
EXAMPLE: example Weyl; shows examples
"
{
  string rname=nameof(basering);
  if ( rname == "basering") // i.e. no ring has been set yet
  {
    "You have to call the procedure from the ring";
    return();
  }
  int @chr = 0;
  if ( size(#) > 0 )
  {
    if ( typeof( #[1] ) == "int" )
    {
      @chr = #[1];
    }
  }
  int nv = nvars(basering);
  int N = nv div 2;
  if ((nv % 2) != 0)
  {
    "Cannot create Weyl structure for an odd number of generators";
    return();
  }
  matrix @D[nv][nv];
  int i;
  for ( i=1; i<=N; i++ )
  {
    if ( @chr==0 ) // default
    {
      @D[i,N+i]=1;
    }
    else
    {
      @D[2*i-1,2*i]=1;
    }
  }
  ncalgebra(1,@D);
  return();
}
example
{
  "EXAMPLE:";echo=2;
  ring A1=0,(x(1..2),d(1..2)),dp;
  Weyl();
  A1;
  kill A1;
  ring B1=0,(x1,d1,x2,d2),dp;
  Weyl(1);
  B1;
}

///////////////////////////////////////////////////////////////////////////////
proc Heisenberg(int N, list #)
"USAGE:   Heisenberg(N, [p,d]); N an integer (setting 2*N+1 variables), p an optional integer (field characteristic), d an optional integer (power of h in the commutator)
RETURN:  N-th  Heisenberg algebra in x(i),y(i),h realization
NOTE: you have to activate this ring with the "setring" command
EXAMPLE: example Heisenberg; shows examples
"
{
  int @chr = 0;
  int @deg = 1;
  if ( size(#) > 0 )
  {
    if ( typeof( #[1] ) == "int" )
    {
      @chr = #[1];
    }
  }
  if ( size(#) > 1 )
  {
    if ( typeof( #[2] ) == "int" )
    {
      @deg = #[2];
      if (@deg <1) { @deg = 1; }
    }
  }
  ring @@r=@chr,(x(1..N),y(1..N),h),lp;
  matrix D[2*N+1][2*N+1];
  int i;
  for (i=1;i<=N;i++)
  {
    D[i,N+i]=h^@deg;
  }
  ncalgebra(1,D);
  return(@@r);
}
example
{
  "EXAMPLE:";echo=2;
  def a = Heisenberg(2);
  setring a; 
  a;
}

///////////////////////////////////////////////////////////////////////////////

proc Exterior(list #)
"USAGE:   Exterior();
RETURN:  qring, the exterior algebra of a basering
NOTE: you have to activate this qring with the "setring" command
EXAMPLE: example Exterior; shows examples
"
{
  string rname=nameof(basering);
  if ( rname == "basering") // i.e. no ring has been set yet
  {
    "You have to call the procedure from the ring";
    return();
  }
  int N = nvars(basering);
  string NewRing = "ring @R=("+charstr(basering)+"),("+varstr(basering)+"),("+ordstr(basering)+");";
  execute(NewRing);
  matrix @E = UpOneMatrix(N);
  @E = -1*(@E);
  ncalgebra(@E,0);
  int i;
  ideal Q;
  for ( i=1; i<=N; i++ )
  {
    Q[i] = var(i)^2;
  }
  Q = twostd(Q);
  qring @EA = Q;
  return(@EA);
}
example
{
  "EXAMPLE:";echo=2;
  ring R = 0,(x(1..3)),dp;
  def ER = Exterior();
  setring ER;
  ER;
}

///////////////////////////////////////////////////////////////////////////////

proc CreateWeyl(int n, list #)
"USAGE:   CreateWeyl(n,[p]); n an integer, n>0; p an optional integer (field characteristic)
RETURN:  a ring, describing n-th Weyl algebra
NOTE:    You have to activate this ring with the "setring" command.
         The presentation of n-th Weyl algebra is classical:
	 D(i)x(i)=x(i)D(i)+1
SEE ALSO: Weyl
EXAMPLE: example CreateWeyl; shows examples
"{
  if (n<1)
  {
    Print("Incorrect input");
    return();
  }
  int @p = 0;
  if ( size(#) > 0 )
  {
    if ( typeof( #[1] ) == "int" )
    {
      @p = #[1];
    }
  }
  if (n ==1)
  {
    ring @rr = @p,(x,D),dp;
  }
  else
  {
    ring @rr = @p,(x(1..n),D(1..n)),dp;
  }
  setring @rr;
  Weyl();
  return(@rr);
}
example
{ "EXAMPLE:"; echo = 2;
   def a = CreateWeyl(3);
   setring a;
   a;
}

//////////////////////////////////////////////////////////////////////

proc Is_NC()
"USAGE:   Is_NC();
RETURN:   1, if basering is noncommutative; 0 otherwise.
EXAMPLE: example Is_NC; shows examples
"{
  string rname=nameof(basering);
  if ( rname == "basering") // i.e. no ring has been set yet
  {
    "You have to call the procedure from the ring";
    return();
  }
  int n = nvars(basering);
  int i,j;
  poly p;
  for (i=1; i<n; i++)
  {
    for (j=i+1; j<=n; j++)
    {
      p = var(j)*var(i) - var(i)*var(j);
      if (p!=0) { return(1);}
    }
  }
  return(0);
}
example
{ "EXAMPLE:"; echo = 2;
   def a = CreateWeyl(2);
   setring a;
   Is_NC();
   kill a;
   ring r = 17,(x(1..7)),dp;
   Is_NC();
   kill r;
}

//////////////////////////////////////////////////////////////////////