source: git/Singular/LIB/central.lib @ 66d68c

///////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Noncommutative";
info="
LIBRARY:  central.lib      Computation of central elements of GR-algebras

AUTHOR:  Oleksandr Motsak,    U@D, where U={motsak}, D={mathematik.uni-kl.de}

OVERVIEW: A library for computing elements of the center and centralizers of sets of elements in GR-algebras.

KEYWORDS:  center; centralizer; reduce; centralize; PBW

PROCEDURES:
  centralizeSet(F, V)          v.s. basis of the centralizer of F within V
  centralizerVS(F, D)          v.s. basis of the centralizer of F
  centralizerRed(F, D[, N])    reduced basis of the centralizer of F
  centerVS(D)                  v.s. basis of the center
  centerRed(D[, k])            reduced basis of the center
  center(D[, k])               reduced basis of the center
  centralizer(F, D[, k])       reduced bais of the centralizer of F
  sa_reduce(V)              's.a. reduction' of pairwise commuting elements
  sa_poly_reduce(p, V)      's.a. reduction' of p by pairwise commuting elements
  inCenter(T)                  checks the centrality of list/ideal/poly T
  inCentralizer(T, S)          checks whether list/ideal/poly T commute with S
  isCartan(p)                  checks whether polynomial p is a Cartan element
  applyAdF(Basis, f)           images of elements under the k-linear map Ad_f
  linearMapKernel(Images)      kernel of a linear map given by images
  linearCombinations(Basis, C)  k-linear combinations of elements
  variablesStandard()          set of algebra generators in their natural order
  variablesSorted()            heuristically sorted set of algebra generators
  PBW_eqDeg(Deg)               PBW monomials of given degree
  PBW_maxDeg(MaxDeg)           PBW monomials up to given degree
  PBW_maxMonom(MaxMonom)       PBW monomials up to given maximal monomial
";

LIB "general.lib" // for "sort"
LIB "poly.lib" // for "maxdeg"


/******************************************************/
// ::DefaultStuff:: Shortcuts to useful short functions. Just to avoid if( if( if( ... ))).
/******************************************************/


/******************************************************/
static proc DefaultValue ( def s, list # ) // Process general variable parameters list
  "
RETURN: s or (typeof(s))(#)
"
{
  def @p = s;
  if ( size(#) > 0 )
  {
    if ( typeof(#[1]) == typeof(s) )
      {
        @p = #[1];
      }
  }
  return( @p );
}

/******************************************************/
static proc DefaultInt( list # ) // Process variable parameters list with 'int' default value
  "
RETURN: 0 or int(#)
"
{
  int @p = 0;
  return( DefaultValue( @p, # ) );
}

/******************************************************/
static proc DefaultIdeal ( list # ) // Process variable parameters list with 'ideal' default value
  "
RETURN: 0 or ideal(#)
"
{
  ideal @p = 0;
  return( DefaultValue( @p, # ) );
}



/******************************************************/
// ::Debug:: Shortcuts to used debugging functions.
/******************************************************/


/******************************************************/
static proc toprint( int pl ) // To print or not to print?!? The answer is here!
  "
RETURN: 1 means to print, otherwise 0.
"
{
  return( printlevel >= ( 3 -  pl) ); // voice + ?
}

/******************************************************/
static proc DBPrint( int pl, list # ) // My 'dbprint' which uses toprint(i).
  "
    USAGE:
"
{
  if( toprint(pl) )
    {
      dbprint(1, #);
    }
}

/******************************************************/
static proc BCall( string Name, list # ) // This function must be called at the beginning of every 'mathematical' function.
  "
USAGE: Name is a name of a mathematical function to trace. # means parameters into the function.
"
{
  if( toprint(0) )
    {
      "CALL: ", Name, "( ";
      dbprint(1, #);
      "     )";
    }
}

/******************************************************/
static proc ECall(string Name, list #) // This function must be called at the end of every 'mathematical' function.
  "
USAGE: Name is a name of a mathematical function to trace. # means result of the function.
"
{
  if( toprint(0) )
    {
      "RET : ", Name, " => Result: {";
      dbprint(1, #);
      "    }";
    }
}



/******************************************************/
// ::Helpers:: Small functions used in this library.
/******************************************************/

/******************************************************/
static proc makeNice( def l )
  "
RETURN: the same as input
PURPOSE: make 'nice' polynomials, kill scalars
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "makeNice", l ); }; /*4DEBUG*/

  poly p;

  if( typeof(l) == "poly" )
    {
      // "normal" polynomial form == no denominators, gcd of coeffs is a unit
      l = cleardenom( l );
      if ( maxdegInt(l) > 0 )
        {
          l = cleardenom( l / leadcoef(l) );
        }
    } else
      {
        if( typeof(l) == "ideal" )
          {
            for( int i = 1; i <= size(l); i++ )
              {
                p = l[i];
                   p = cleardenom( p );

                // Now make polynomials look nice:
                if ( maxdegInt(p) > 0 ) // throw away scalars!
                  {
                        // "normal" polynomial form == no denominators, gcd of coeffs is a unit
                        p = cleardenom( p / leadcoef(p) );
                  } else
                    {
                      p = 0; // no constants
                    }
                l[i] = p;

              }

            l = simplify(l, 2 + 8);
          }
      }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "makeNice", l ); }; /*4DEBUG*/
  return( l );
}



/******************************************************/
static proc monomialForm( def p )
  "
: p is either poly or ideal
RETURN: polynomial with all monomials from p but without coefficients.
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "monomialForm", p ); }; /*4DEBUG*/

  poly result = 0; int k, j; poly m;

  if( typeof(p) == "ideal" ) //
    {
      if( ncols(p) > 0 )
        {
          result = uni_poly( p[1] );

          for ( k = ncols(p); k > 1; k -- )
            {
              for( j = size(p[k]); j > 0; j-- )
                {
                  m = leadmonom( p[k][j] );

                  if( size(result + m) > size(result) ) // trick!
                    {
                      result = result + m;
                    }
                }

            }
        }
    }
  else
    {
      if( typeof(p) == "poly" ) //
        {
          result = uni_poly(p);
        } else
          {
            ERROR( "ERROR: Wrong input! Expected polynomial or ideal!" );
          }
    }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "monomialForm", result ); }; /*4DEBUG*/
  return( result );
}

/******************************************************/
static proc uni_poly( poly p )
  "
    returns polynomial with the same monomials but without coefficients.
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "uni_poly", p ); }; /*4DEBUG*/

  poly result = 0;

  for ( int k = size(p); k > 0; k -- )
    {
      result = result + leadmonom( p[k] );
    }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "uni_poly", result ); }; /*4DEBUG*/
  return( result );
}





/******************************************************/
static proc smoothQideal( ideal I, list #)
  "
PURPOSE: smooths the ideal in a current QUOTIENT(!) ring.
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "smoothQideal", I ); }; /*4DEBUG*/

  ideal A = matrix(I) - matrix(NF( I, twostd(DefaultIdeal(#)), 1 )); // component wise

  if( size(A) > 0 ) // Were there any changes (any non-zero component)?
    {
      ideal T;

      int j = 1;

      for( int i = 1; i <= ncols(I); i++ )
        {
          if( size(A[i]) == 0 ) // keep only unchanged elements
            {
              T[ j ] = I[i]; j++;
            }
        }
      I = T;
    }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "smoothQideal", I ); }; /*4DEBUG*/

  return( I );
}




/******************************************************/
// ::PBW:: PBW basis construction.
/******************************************************/




/******************************************************/
proc PBW_maxDeg( int MaxDeg )
"
USAGE: PBW_maxDeg(MaxDeg); MaxDeg int
PURPOSE: Compute PBW elements up to a given maximal degree.
RETURN: ideal consisting of found elements.
NOTE: unit is omitted. Weights are ignored!
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "PBW_maxDeg", MaxDeg ); }; /*4DEBUG*/

  ideal Basis = 0;

  for (int k = 1; k <= MaxDeg; k ++ )
    {
      Basis = Basis + maxideal(k);
    }

  Basis = smoothQideal( Basis );

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "PBW_maxDeg", Basis ); }; /*4DEBUG*/
  return( Basis );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  def A = nc_algebra(1,D); // this algebra is U(sl_2)
  setring A;

  // PBW Basis of A_2 - monomials of degree <= 2, without unit:
  PBW_maxDeg( 2 );
}


/******************************************************/
proc PBW_eqDeg( int Deg )
"
USAGE: PBW_eqDeg(Deg); Deg int
PURPOSE: Compute PBW elements of a given degree.
RETURN: ideal consisting of found elements.
NOTE: Unit is omitted. Weights are ignored!
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "PBW_eqDeg", Deg ); }; /*4DEBUG*/

  ideal Basis = smoothQideal( maxideal( Deg ) );

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "PBW_eqDeg", Basis ); }; /*4DEBUG*/
  return( Basis );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)

  // PBW Basis of A_2 \ A_1 - monomials of degree == 2:
  PBW_eqDeg( 2 );
}


/******************************************************/
proc PBW_maxMonom( poly MaxMonom )
"
USAGE: PBW_maxMonom(m); m poly
PURPOSE: Compute PBW elements up to a given maximal one.
RETURN: ideal consisting of found elements.
NOTE: Unit is omitted. Weights are ignored!
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "PBW_maxMonom", MaxMonom ); }; /*4DEBUG*/

  ideal K = 0;

  intvec exp = leadexp( MaxMonom );

  for ( int k = 1; k <= size(exp); k ++ )
    {
      K[ 1 + size(K) ] = var(k)^( 1 + exp[k] );
    }

  attrib(K, "isSB", 1);

  K = kbase(K);

  K = K[ (ncols(K)-1)..1]; // reverse, kill last 1

  K = smoothQideal( K );

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "PBW_maxMonom", K ); }; /*4DEBUG*/

  return( K );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  def A = nc_algebra(1,D); // this algebra is U(sl_2)
  setring A;

  // At most 1st degree in e, h and at most 2nd degree in f, unit is omitted:
  PBW_maxMonom( e*(f^2)* h );
}




/******************************************************/
// ::CORE:: Core procedures...
/******************************************************/



/******************************************************/
proc applyAdF( ideal I, poly p )
  "
USAGE: applyAdF(B, f); B ideal, f poly
PURPOSE: Apply Ad_f to every element of B
RETURN: ideal, generated by Ad_f(B[i]), 1<=i<=size(B)
NOTE:  Ad_f(v) := [f, v] = f*v - v*f
SEE ALSO:   linearMapKernel; linearMapKernel
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "applyAdF", I, p ); }; /*4DEBUG*/

  poly t; ideal II = 0;

  for ( int k = ncols(I); k > 0; k --)
    {
      t = I[k];
      II[k] = p * t - t * p; // we have to reduce smooth images in Qrings...
    }

  II = NF( II, twostd(0) ); // ... now!

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "applyAdF", II ); }; /*4DEBUG*/
  return( II );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)

  // Let us consider the linear map Ad_{e} from A_2 into A.

  // Compute the PBW basis of A_2:
  ideal Basis = PBW_maxDeg( 2 ); Basis;

  // Compute images of basis elements under the linear map Ad_e:
  ideal Image = applyAdF( Basis, e ); Image;

  // Now we have a linear map given by: Basis_i --> Image_i
  // Let's compute its kernel K:

  // 1. compute syzygy module C:
  module C = linearMapKernel( Image ); C;

  // 2. compute corresponding combinations of basis vectors:
  ideal K = linearCombinations(Basis, C); K;

  // Let's check that Ad_e(K) is zero:
  applyAdF( K, e );
}



/******************************************************/
proc linearMapKernel( ideal Images )
"
USAGE: linearMapKernel( Images ); Images ideal
PURPOSE: Computes the syzygy module of the linear map given by Images.
RETURN: syzygy module, or int(0) if all images are zeroes
SEE ALSO:   applyAdF; linearMapKernel
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "linearMapKernel", Images ); }; /*4DEBUG*/

  // This must be a list of monomials in a form of polynomial (sum with coeffs == 1)
  poly Monomials = monomialForm( Images );

  int N = size( Monomials ); // number of different monomials

  if ( N == 0 ) // & ncols( Images ) > 0 => all Images == 0
    {
      int result = 0;

      /*4DEBUG*/        if( defined( @@@DEBUG ) ){ ECall( "linearMapKernel", result ); }; /*4DEBUG*/
      return( result );
    }

  // Compute matrix MD
  module MD; // zero

  int x, y;

  vector w;

  poly p, m;

  int V = ncols(Images); // must be equal to ncols(Basis) and size(Basis)!

  // We take monomials as vector space basis of <Image>_k...

  // TODO: Is there any other way to compute a basis of it and represent images as
  // linear combination of them???

  // Maybe some 'free resolution' stuff??? But we need linear maps only!!!

  for ( x = 1; x <= V; x++ ) // For every Image vector
    {
      w = 0;

      p = Images[x];

      y = 1; // from 1st monomial in Monomials...

      while( size(p) > 0 )
        {
          m = leadmonom(p);

          // y < N!
          while( Monomials[y] != m )
            // There MUST be this monomial because of the construction of Monomials polynomial!
            {
              y++; // to size()
            }

          // found monomial m at position y.

          w = w + leadcoef(p) * gen(y); // leadcoef(p) MUST be nonzero!
          p = p - lead(p); // kill lead term
        }

      MD[x] = w;
    }

  /*******************************************/

  // save options
  intvec v = option( get );

  // set right options
  option( redSB );
  option( redTail );

  // compute everything in a right form
  MD = simplify( groebner( syz(MD) ), 1 + 2 + 8 );
  // note that MD is a matrix of numbers - no polynomials...

  // restore options
  option( set, v );

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "linearMapKernel", MD ); }; /*4DEBUG*/

  return( MD );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  def A = nc_algebra(1,D); // this algebra is U(sl_2)
  setring A;

  // Let us consider the linear map Ad_{e} from A_2 into A.

  // Compute the PBW basis of A_2:
  ideal Basis = PBW_maxDeg( 2 ); Basis;

  // Compute images of basis elements under the linear map Ad_e:
  ideal Image = applyAdF( Basis, e ); Image;

  // Now we have a linear map given by: Basis_i --> Image_i
  // Let's compute its kernel K:

  // 1. compute syzygy module C:
  module C = linearMapKernel( Image ); C;

  // 2. compute corresponding combinations of basis vectors:
  ideal K = linearCombinations(Basis, C); K;

  // Let's check that Ad_e(K) is zero:
  ideal Z = applyAdF( K, e ); Z;

  // Now linearMapKernel will return a single integer 0:
  def CC  = linearMapKernel(Z); typeof(CC); CC;
}


/******************************************************/
proc linearCombinations( ideal Basis, module KER )
  "
USAGE:  linearCombinations( Basis, C ); Basis ideal, C module
PURPOSE: forms linear combinations of elements from Basis by replacing gen(i) by Basis[i] in C
RETURN: ideal generated by computed linear combinations
SEE ALSO:   linearMapKernel; applyAdF
"
{

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "linearCombinations", Basis, KER ); }; /*4DEBUG*/


  number c;

  int x, y;

  vector w;

  poly p;

  ideal result = 0;

  // Kernel' basis translation
  for ( x = 1; x <= ncols(KER); x++ )
    {
      p = 0;
      w = KER[x];

      for ( y = 1; y <= nrows(w); y++ )
        {
          c = leadcoef( w[y] );

          if ( c != 0 )
            {
              p = p + c * Basis[y]; // linear combination of base vectors { Basis[y] }
            }
        }
      result[ x ]  = p;
    }


  // no reduction in quotient algebras is needed. No multiplications were done!


  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "linearCombinations", result ); }; /*4DEBUG*/

  return( result );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)

  // Let us consider the linear map Ad_{e} from A_2 into A.

  // Compute the PBW basis of A_2:
  ideal Basis = PBW_maxDeg( 2 ); Basis;

  // Compute images of basis elements under the linear map Ad_e:
  ideal Image = applyAdF( Basis, e ); Image;

  // Now we have a linear map given by: Basis_i --> Image_i
  // Let's compute its kernel K:

  // 1. compute syzygy module C:
  module C = linearMapKernel( Image ); C;

  // 2. compute corresponding combinations of basis vectors:
  ideal K = linearCombinations(Basis, C); K;

  // Let's check that Ad_e(K) is zero:
  applyAdF( K, e );
}



/******************************************************/
static proc LINEAR_MAP_KERNEL(ideal Basis, ideal Images ) // Ker of the linear map given by its values on basis vectors
  "
PURPOSE: Computation of the kernel basis of the linear map given by the list @given
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "LINEAR_MAP_KERNEL", Basis, Images ); }; /*4DEBUG*/

  ideal result = 0;

  if ( size( Basis ) == 0 )
    {
      /*4DEBUG*/        if( defined( @@@DEBUG ) ){ ECall( "LINEAR_MAP_KERNEL", result ); }; /*4DEBUG*/
      return( result );
    }

  // compute fundamental solutions system
  def T = linearMapKernel( Images );


  // check result of linearMapKernel
  if( typeof(T) == "int" )
  {
    if( T == 0 )
    {
      // All zeroes! Return Basis:
      /*4DEBUG*/        if( defined( @@@DEBUG ) ){ ECall( "LINEAR_MAP_KERNEL", Basis ); }; /*4DEBUG*/
      return( Basis );
    }
  }
  else
  {
    if( typeof(T) != "module" )
    {
       ERROR( "Wrong output from the 'linearMapKernel' function!" );
    }
  }

  result = linearCombinations( Basis, T );

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "LINEAR_MAP_KERNEL", result ); }; /*4DEBUG*/
  return( result );
}




/******************************************************/
static proc ZeroKer( ideal Basis, ideal Images ) // VS Basis of a Kernel of the linear map AD_h, h is a Cartan element
"
PURPOSE: Computes VS Basis of a Kernel of the linear map AD_h, when h is a Cartan element
NOTE: the result is a set of all basis vectors having a zero image
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "ZeroKer", Basis, Images ); }; /*4DEBUG*/

  ideal result = 0;

  for( int i = 1; i <= ncols( Basis ); i++ )
    {
      if( size( Images[i] ) == 0 ) // zero image?
        {
          result[ 1 + size(result) ] = Basis[i]; // take this basis vector!
        }
    }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "ZeroKer", result ); }; /*4DEBUG*/
  return( result );
}




/******************************************************/
// ::Variables:: Computes a set of variables
/******************************************************/



/******************************************************/
// Returns an ideal of variables in a current base ring.
proc variablesStandard()
"
USAGE:      variablesStandard();
RETURN:     ideal, generated by algebra variables
PURPOSE:    computes the set of algebra variables taken in their natural order
SEE ALSO:   variablesSorted
EXAMPLE:    example variablesStandard; shows an example
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "variablesStandard" ); }; /*4DEBUG*/

  ideal result = maxideal(1);

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "variablesStandard", result ); }; /*4DEBUG*/
  return( result );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=-z;  D[1,3]=2*x;  D[2,3]=-2*y;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)
  // Variables in their natural order:
  variablesStandard();
}

/******************************************************/
// Sorts variables into an ideal. This is a kind of heuristics!
proc variablesSorted()
"
USAGE:      variablesSorted();
RETURN:     ideal, generated by sorted algebra variables
PURPOSE:    computes the set of algebra variables sorted so that
@* Cartan variables go first
NOTE:       This is a heuristics for the computation of the center:
@* it is better to compute centralizers of Cartan variables first since in this
@* case we can omit solving the system of equations.
SEE ALSO:   variablesStandard
EXAMPLE:    example variablesSorted; shows an example
"{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "variablesSorted" ); }; /*4DEBUG*/

  ideal V   = variablesStandard();
  int  N    = size( V ); // == nvars( basering )

  ideal result = 0;

  int  r_begin = 1;
  int  r_end   = N;

  poly v;

  for( int k = 1; k <= N; k++ )
    {
      v = V[k];

      if( isCartan(v) == 1 ) // Cartan elements go 1st
        {
          result[r_begin] = v;
          r_begin++;
        } else // Other - in the end...
          {
            result[r_end] = v;
            r_end--;
          }
    }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "variablesSorted", result ); }; /*4DEBUG*/
  return( result );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=-z;  D[1,3]=2*x;  D[2,3]=-2*y;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)
  // There is only one Cartan variable - z in U(sl_2),
  // it must go 1st:
  variablesSorted();
}





/******************************************************/
/******************************************************/
// ::BasicCentralizerComputation:: Basic functions for centralize' computation.
/******************************************************/
/******************************************************/





/******************************************************/
// HL 'core' function
proc centralizeSet( ideal F, ideal V )
"
USAGE:      centralizeSet( F, V ); F, V ideals
INPUT:      F, V finite sets of elements of the base algebra
RETURN:     ideal, generated by computed elements
PURPOSE:    computes a vector space basis of the centralizer of the set F in the vector space generated by V over the ground field
SEE ALSO:   centralizerVS; centralizer; inCentralizer
EXAMPLE:    example centralizeSet; shows an example
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "centralizeSet", F, V ); }; /*4DEBUG*/

  int  N = size(F);

  if( N == 0)
    {
      ERROR( "F MUST be non empty!!!" );
    }

  DBPrint(1, "BasisSize: " + string(size(V)) );

  ideal Images;

  for( int v = 1; (v <= N) and (size(V) > 0); v++ )
    {
      DBPrint(1, "Centralizing " + string(F[v]) );

      Images = applyAdF( V, F[v] );

      if( (isCartan(F[v]) == 1) or (size(V) == 1) )
        {
          V = ZeroKer( V, Images );
        } else
          {
            V = LINEAR_MAP_KERNEL( V, Images );
          }

      // Printing...
      DBPrint(1, "Progress: [ " + string(v) + " / " + string(N) + " ]"+
              " => BasisSize: " + string(size(V)) );
    }

  V = makeNice(V);

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "centralizeSet", V ); }; /*4DEBUG*/

  return( V );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring A = 0,(e(1..4)),dp;
  matrix D[4][4]=0;
  D[2,4] = -e(1);
  D[3,4] = -e(2);
  // This is A_4_1 - the first real Lie algebra of dimension 4.
  def A_4_1 = nc_algebra(1,D); setring A_4_1;

  ideal F = variablesSorted(); F;

  // the center of A_4_1 is generated by
  // e(1) and -1/2*e(2)^2+e(1)*e(3)
  // therefore one may consider computing it in the following way:

  // 1. Compute a PBW basis consisting of
  //    monomials with exponent <= (1,2,1,0)
  ideal V = PBW_maxMonom( e(1) * e(2)^ 2 * e(3) );

  // 2. Compute the centralizer of F within the vector space
  //    spanned by these monomials:
  ideal C = centralizeSet( F, V ); C;

  inCenter(C); // check the result
}



/******************************************************/
proc centralizerVS( ideal F, int d )
  "
USAGE:      centralizerVS( F, D ); F ideal, D int
RETURN:     ideal, generated by computed elements
PURPOSE:    computes a vector space basis of centralizer(F) up to degree D
NOTE:       D must be non-negative
SEE ALSO:   centerVS; centralizer; inCentralizer
EXAMPLE:    example centralizerVS; shows an example
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "centralizerVS", F, d ); }; /*4DEBUG*/

  if( size(F) == 0)
    {
      ERROR( "F MUST be non-empty!!!" );
    }

  ideal V = centralizeSet( F, PBW_maxDeg( d ) ); // basis of the Centralizer of S in PBW basis

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "centralizerVS", V ); }; /*4DEBUG*/

  return( V );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=-z;  D[1,3]=2*x;  D[2,3]=-2*y;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)
  ideal F = x, y;
  // find generators of the vector space of elements
  // of degree <= 4 commuting with x and y:
  ideal C = centralizerVS(F, 4);
  C;
  inCentralizer(C, F); // check the result
}




/******************************************************/
// ::CenterAliases:: Basic functions/aliases for center' computation.
/******************************************************/




/******************************************************/
proc centerVS( int D )
"
USAGE:      centerVS( D ); D int
RETURN:     ideal, generated by computed elements
PURPOSE:    computes a vector space basis of the center up to degree D
NOTE:       D must be non-negative
SEE ALSO:   centralizerVS; center; inCenter
EXAMPLE:    example centerVS; shows an example
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "centerVS", D ); }; /*4DEBUG*/


  if( nameof( basering ) == "basering" )
    {
      //        ERROR( "No current ring!" );
    }

  if( D < 0 )
    {
      ERROR( "Degree D must be non-negative!" );
    }

  ideal result = centralizerVS( variablesSorted(), D );

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "centerVS", result ); }; /*4DEBUG*/

  return( result );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=-z;  D[1,3]=2*x;  D[2,3]=-2*y;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)
  // find a basis of the vector space of all
  // central elements of degree <= 4:
  ideal Z = centerVS(4);
  Z;
  // note that the second element is the square of the first
  // plus a multiple of the first:
  Z[2] - Z[1]^2 + 8*Z[1];
  inCenter(Z); // check the result
}


/******************************************************/
proc centralizerRed( ideal F, int D, list # )
"
USAGE:      centralizerRed( F, D[, N] ); F ideal, D int, N optional int
RETURN:     ideal, generated by computed elements
PURPOSE:    computes subalgebra generators of centralizer(F) up to degree D.
NOTE:       In general, one cannot compute the whole centralizer(F).
@* Hence, one has to specify a termination condition via arguments D and/or N.
@* If D is positive, only centralizing elements up to degree D are computed.
@* If D is negative, the termination is determined by N only.
@* If N is given, the computation stops if at least N elements have been found.
@* Warning: if N is given and bigger than the actual number of generators,
@* the procedure may not terminate.
@* Current ordering must be a degree compatible well-ordering.
SEE ALSO:   centralizerVS; centerRed; centralizer; inCentralizer
EXAMPLE:    example centralizerRed; shows an example
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "centralizerRed", F, D, # ); }; /*4DEBUG*/

  if( nameof( basering ) == "basering" )
    {
      //        ERROR( "No current ring!" );
    }

  if( size(F) == 0)
    {
      ERROR( "F MUST be non-empty!!!" );
    }

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

  int i, j, l, d;

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

  int k = DefaultInt(#);

  int m = (k > 0);

  int @MinDeg = 6; // starting guess for Maximal Bounding Degree, 6
  int @Delta  = 4; // increment of it, 4

  if( m and (D <= 0) )
    {
      // minimal guess
      D = @MinDeg;
    }

  if( !m and D < 0)
    {
      ERROR("Wrong bounding condition!");
    }

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

  def NCRING = basering; // Non-commutative ring
  list L = ringlist( NCRING );
  def L1, L2, L3, L4 = L[1..4]; // General components

  def COMMRING = ring( list( L1, L2, L3, L4 ) ); // Underlying commutative ring


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

  // we keep the list of found leading monomials in the commutative ring!
  setring COMMRING;

  // Init
  list FOUND_LEADING_MONOMIALS = list();

  for( i = 1; i <= D; i++ )
    {
      FOUND_LEADING_MONOMIALS[i] = ideal();
    }

  ideal FLM, NEW, T; // in COMMRING

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

  setring NCRING;

  ideal result, FLM, PBW, NEW, T, P; // in NCRING

  // Main loop:
  i = 1;

  while( i <= D )
    {
      DBPrint( 1, "Current degree is " + string(i) );

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

      // Compute current "reduced" PBW basis...

      // Prepare current found leading monomials
      setring COMMRING;
      FLM = FOUND_LEADING_MONOMIALS[i];

      // And back to NCRing
      setring NCRING;

      FLM = imap(COMMRING, FLM); // We cannot write imap(COMMRING, FOUND_LEADING_MONOMIALS[i]) :(((

      attrib(FLM, "isSB", 1); // just to avoid "no standard basis" warning.

      // degrees should not change,
      // no monomials should be multiplied here
      T = reduce( PBW_eqDeg( i ), FLM, 1 );

      // we simply kill in T monomials occurring in FOUND_LEADING_MONOMIALS[i]
      P = PBW + T; // + simplifies

      // Compute current centralizer
      NEW = centralizeSet( F, P );

      if( size(NEW) > 0 )
        {
          // In order to speedup multiplications we are going into a commutative ring:
          setring COMMRING;

          // we can perform commutative interreduction
          // since no monomials should be multiplied!
          // degrees should not change
          NEW = interred( imap( NCRING, NEW ) );

          // Go back!
          setring NCRING;

          NEW = imap( COMMRING, NEW );

          DBPrint( 1, "Found: ", NEW );

          // Add them to result...
          result = result + NEW;
        }

      // Did we find needed number of generators? Or reached the bound?
      if( (m and (size(result) >= k)) or (!m and (i == D)) )
        {
          break; // Get out of here!!!
        }

      // otherwise we must update FOUND_LEADING_MONOMIALS
      if( size(NEW) > 0 )
        {
          setring COMMRING;

          FLM = 0;

          // We must update FOUND_LEADING_MONOMIALS!!!
          for( j = 1; j <= size(NEW); j++ )
            {
              FLM[j] = leadmonom( NEW[j] ); // we are interested in leading monomials only!
            }

          FOUND_LEADING_MONOMIALS[i] = FOUND_LEADING_MONOMIALS[i] + FLM;

          for( j = 1; j <= D; j = j + i ) // For every degree (j*i) of LNEW, do:
            {
              for( l = j; (l+i) <= D; l++ )
                {
                  FOUND_LEADING_MONOMIALS[l+i] =
                    FOUND_LEADING_MONOMIALS[l+i] + FOUND_LEADING_MONOMIALS[l] * FLM;
                }
            }

          // Return to NCRING
          setring NCRING;

          FLM = imap(COMMRING, FLM);
          attrib(FLM, "isSB", 1);// just to avoid "no standard basis" warning.

          // we simply kill in T monomials occurring in FOUND_LEADING_MONOMIALS[i]
          T = reduce( T, FLM, 1 );

          PBW = PBW + T;
        } else
          {
            PBW = P;
          }


      if( m and (i == D) ) // Was the previous estimation too small???
        {
          // We must update FOUND_LEADING_MONOMIALS in their Commutative world:
          setring COMMRING;

          // Init new grades:
          for( j = D + 1; j <= (D + @Delta); j++ )
            {
              FOUND_LEADING_MONOMIALS[j] = ideal();
            }

          FLM = 0;

          // All previously computed elements in their order!
          NEW = imap( NCRING, result );

          for( j = 1; j <= size(NEW); j++ )
            {
              FLM[j] = leadmonom( NEW[j] ); // we are interested in leading monomials only!
            }

          while( size(FLM) > 0 )
            {
              // minimal degree:
              d = mindegInt(FLM);  /// ### ///

              // take all of minimal degree:
              T = jet( FLM, d );

              // there are size(T) elements of smallest degree (deg(FLM[1])) in FLM!

              // Add them in the same way:
              for( j = 1; j <= (D + @Delta); j = j + d ) // For every degree (j*d) of T, do:
                {
                  for( l = j; (l + d) <= (D + @Delta); l++ )
                    {
                      if( (l + d) > D ) // Only new should be updated!
                        {
                          FOUND_LEADING_MONOMIALS[l+d] =
                            FOUND_LEADING_MONOMIALS[l+d] + FOUND_LEADING_MONOMIALS[l] * T;
                        }
                    }
                }

              // Kill them from FLM:
              if( size(T) < size(FLM) )
                {
                  FLM = FLM[ (size(T)+1) .. size(FLM) ];
                } else
                  {
                    FLM = 0; // break;
                  }

            }

          // Go back...
          setring NCRING;

/*
    if(toprint())
    {
      typeof(@Delta); @Delta;
      typeof(D); D;
    };
*/
          // And set new Bound
          D = D + @Delta;
        }

      i++;
    }

  result = makeNice(result);

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "centralizerRed", result ); }; /*4DEBUG*/

  return( result );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=-z;  D[1,3]=2*x;  D[2,3]=-2*y;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)
  ideal F = x, y;
  // find subalgebra generators of degree <= 4 of the subalgebra of
  // all elements commuting with x and y:
  ideal C = centralizerRed(F, 4);
  C;
  inCentralizer(C, F); // check the result
}


/******************************************************/
proc centerRed( int D, list # )
"
USAGE:      centerRed( D[, N] ); D int, N optional int
RETURN:     ideal, generated by computed elements
PURPOSE:    computes subalgebra generators of the center up to degree D
NOTE:       In general, one cannot compute the whole center.
@* Hence, one has to specify a termination condition via arguments D and/or N.
@* If D is positive, only central elements up to degree D will be found.
@* If D is negative, the termination is determined by N only.
@* If N is given, the computation stops if at least N elements have been found.
@* Warning: if N is given and bigger than the actual number of generators,
@* the procedure may not terminate.
@* Current ordering must be a degree compatible well-ordering.
SEE ALSO:   centralizerRed; centerVS; center; inCenter
EXAMPLE:    example centerRed; shows an example
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "centerRed", D, # ); }; /*4DEBUG*/

  if( nameof( basering ) == "basering" )
    {
      //        ERROR( "No current ring!" );
    }

  ideal result = centralizerRed( variablesSorted(), D, # );

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "centerRed", result ); }; /*4DEBUG*/

  return( result );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=z;
  def A = nc_algebra(1,D); setring A; // it is a Heisenberg algebra
  // find a basis of the vector space of
  // central elements of degree <= 3:
  ideal VSZ = centerVS(3);
  // There should be 3 degrees of z.
  VSZ;
  inCenter(VSZ); // check the result
  // find "minimal" central elements of degree <= 3
  ideal SAZ = centerRed(3);
  // Only 'z' must be computed
  SAZ;
  inCenter(SAZ); // check the result
}


/******************************************************/
/******************************************************/
// ::SubAlgebraReduction:: A kind of subalgebra reduction...
/******************************************************/
/******************************************************/

/******************************************************/
static proc INTERRED( ideal S )
  "
USAGE:      INTERRED( S ); S ideal
RETURN:      ideal, interreduced S
PURPOSE:     interreduction without monomial multiplication,
    just make every leading monomial occur in a single polynomial
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "INTERRED", S ); }; /*4DEBUG*/

  ideal result = S;

  int flag = 1;

  int i, j, N;

  poly p, lm;

  while( flag == 1 )
    {
      flag = 0;

      result = sort( simplify( result, 1 + 2 + 8) )[1];
      // sorting w.r.t. actual monomial ordering
      // generators with SMALLER(!) leading term come FIRST

      N = size(result);

      // kill leading monomials:

      i = 1;
      while( i < N )
        {
          p = result[i];
          lm = leadmonom(p);

          j = i + 1;
          while( leadmonom(result[j]) == lm )
            {
              result[j] = result[j] - p; // leadcoefs are 1 because of simplify.
              flag = 1; // we have changed something => we do still need to care about it...
              j++;

              if( j > N )
                {
                  break;
                }
            }

          i = j;
        }
    }

  // We are done! No common leading monomials!
  // The result is sorted

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "INTERRED", result ); }; /*4DEBUG*/

  return( result );
}


/******************************************************/
static proc SANF( poly p, list FOUND_LEADING_MONOMIALS )
  "
    reduce p w.r.t. found multiples without ANY polynomial multiplications!
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "SANF", p, FOUND_LEADING_MONOMIALS); }; /*4DEBUG*/

  poly q = p;
  poly head = 0;

  int d; int N = size(FOUND_LEADING_MONOMIALS);

  while( size(q) > 0 )
    {
      d = maxdegInt(p); /// ### ///

      if( (0 < d) and (d <= N) )
        {
          if( size(FOUND_LEADING_MONOMIALS[d]) > 0 )
            {
              attrib( FOUND_LEADING_MONOMIALS[d], "isSB", 1);
              q = reduce( p, FOUND_LEADING_MONOMIALS[d] );
            }

          DBPrint(1, string(p) + " --> " + string(q) );
        }

      if( q == p )
        {
          p = lead(q);

          if( d > 0 )
            {
              // No scalars!
              head = head + p;
            }

          q = q - p;
        }

      p = q;
    }



  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "SANF", head ); }; /*4DEBUG*/

  return( head );
}


/******************************************************/
static proc maxdegInt( ideal I )
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "maxdegInt", I ); }; /*4DEBUG*/

  intmat D = maxdeg(I);

  int max = D[1, 1]; int m;

  for( int c = 2; c <= ncols(D); c++ )
    {
      m = D[1, c];

      if( m > max )
        {
          max = m;
        }
    }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "maxdegInt", max ); }; /*4DEBUG*/

  return( max );
}


/******************************************************/
static proc mindegInt( ideal I )
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "mindegInt", I ); }; /*4DEBUG*/

  intmat D = mindeg(I);

  int min = D[1, 1]; int m;

  for( int c = 2; c <= ncols(D); c++ )
    {
      m = D[1, c];

      if( m < min )
        {
          min = m;
        }
    }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "mindegInt", min ); }; /*4DEBUG*/

  return( min );
}

/******************************************************/
// 'subalgebra basis' computation
proc sa_reduce( ideal V )
"
USAGE:     sa_reduce(V); V ideal
RETURN:     ideal, generated by computed elements
PURPOSE:    compute a subalgebra basis of an algebra generated by the elements of V
NOTE:       At the moment the usage of this procedure is limited to G-algebras
SEE ALSO:   sa_poly_reduce
EXAMPLE:    example sa_reduce; shows an example
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "sa_reduce", V ); }; /*4DEBUG*/

  ideal result = ideal();

  ideal FLM = INTERRED( V ); // The output is sorted "[1]<[2]<[3]<..."

  // We are bounded by maximal degree!!!
  int D = maxdegInt( FLM );

  // Init
  list FOUND_LEADING_MONOMIALS = list();

  int i;

  for( i = 1; i <= D; i++ )
    {
      FOUND_LEADING_MONOMIALS[i] = ideal();
    }

  int d, j, l;

  poly p, q; ideal T;


  int c = 1;  // polynomials in FLM commute pairwise

  for( j = 1; (j < size(FLM)) and (c == 1); j++ )
    {
      p = FLM[j];

      for( l = j+1; (l <= size(FLM)) and (c == 1); l++ )
        {
          q = FLM[l];

          if( NF(p*q - q*p, twostd(0)) != 0  )
            {
              c = 0; // There exists non-commuting pair
            }
        }
    }

  while( size(FLM) > 0 )
    {
//    FLM;

      // Take the 1st element of FLM...
      p = FLM[1]; // SANF( FLM[1], FOUND_LEADING_MONOMIALS );

      FLM[1] = 0; // ...and kill it from FLM

      d = maxdegInt( p );
      T = ideal(p);

//      d; size(FOUND_LEADING_MONOMIALS);

    if( d > 0 )
    {

      FOUND_LEADING_MONOMIALS[d] = FOUND_LEADING_MONOMIALS[d] + T;

      for( j = 1; j <= D; j = j + d ) // For every degree (j*d) of T, do:
        {
          for( l = j; (l + d) <= D; l++ )
            {
              FOUND_LEADING_MONOMIALS[l+d] =
                FOUND_LEADING_MONOMIALS[l+d] + FOUND_LEADING_MONOMIALS[l] * T;

              if( c != 1 )
                {
                  FOUND_LEADING_MONOMIALS[l+d] =
                    FOUND_LEADING_MONOMIALS[l+d] + T * FOUND_LEADING_MONOMIALS[l];
                }
            }
        }
    }

      if( size(FLM) > 0 )
        {
          for( i = 2; i <= ncols(FLM); i++ )
            {
              FLM[i] = SANF( FLM[i], FOUND_LEADING_MONOMIALS );
            }
          FLM = INTERRED( FLM );
        }

        if( size(T) > 0 )
        {
          DBPrint(1, "Found: " + string(T) );
          result = result + T;
        }

    }

  result = makeNice(result);

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "sa_reduce", result ); }; /*4DEBUG*/

  return( result );
}
example
{ "EXAMPLE:"; echo = 2;
 ring AA = 0,(x,y,z),dp;
 matrix D[3][3]=0;
 D[1,2]=-z; D[1,3]=2*x; D[2,3]=-2*y;
 def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)
 poly f = 4*x*y+z^2-2*z; // a central polynomial
 ideal I = f, f*f, f*f*f - 10*f*f, f+3*z^3; I;
 sa_reduce(I); // should be just f and z^3
}



/******************************************************/
// subalgebra reduction of a polynomial
proc sa_poly_reduce( poly p, ideal V )
"
USAGE:      sa_poly_reduce(p, V); p poly, V ideal
RETURN:     polynomial, a reduction of p w.r.t. V
PURPOSE:    computes a reduction of the polynomial p w.r.t. the subalgebra generated by elements of V
NOTE:       At the moment the usage of this procedure is limited to G-algebras
SEE ALSO:   sa_reduce
EXAMPLE:    example sa_poly_reduce; shows an example
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "sa_poly_reduce", p, V ); }; /*4DEBUG*/
  // As previous...

  ideal FLM = INTERRED( V ); // The output is sorted "[1]<[2]<[3]<..."

  // We are bounded by maximal degree!!!
  int D = maxdegInt( FLM + ideal(p)  );

  // Init
  list FOUND_LEADING_MONOMIALS = list();

  int i;

  for( i = 1; i <= D; i++ )
    {
      FOUND_LEADING_MONOMIALS[i] = ideal();
    }

  int d, j, l;

  poly f, q; ideal T;


  int c = 1;  // polynomials in FLM commute pairwise

  for( j = 1; (j < size(FLM)) and (c == 1); j++ )
    {
      f = FLM[j];

      for( l = j+1; (l <= size(FLM)) and (c == 1); l++ )
        {
          q = FLM[l];

          if( NF(f*q - q*f, twostd(0)) != 0 )
            {
              c = 0;
            }
        }
    }


  while( size(FLM) > 0 )
    {
      // Take the 1st element of FLM...
      q = SANF( FLM[1], FOUND_LEADING_MONOMIALS );

      FLM[1] = 0; // ...and kill it from FLM

      d = maxdegInt(q);
      T = ideal(q);

      FOUND_LEADING_MONOMIALS[d] = FOUND_LEADING_MONOMIALS[d] + T;

      for( j = 1; j <= D; j = j + d ) // For every degree (j*d) of T, do:
        {
          for( l = j; (l + d) <= D; l++ )
            {
              FOUND_LEADING_MONOMIALS[l+d] =
                FOUND_LEADING_MONOMIALS[l+d] + FOUND_LEADING_MONOMIALS[l] * T;

              if( c != 1 )
                {
                  FOUND_LEADING_MONOMIALS[l+d] =
                    FOUND_LEADING_MONOMIALS[l+d] + T * FOUND_LEADING_MONOMIALS[l];
                }
            }
        }

      if( size(FLM) > 0 )
        {
          for( i = 2; i <= ncols(FLM); i++ )
            {
              FLM[i] = SANF( FLM[i], FOUND_LEADING_MONOMIALS );
            }
          FLM = INTERRED( FLM );
        }
    }

  poly result = SANF(p, FOUND_LEADING_MONOMIALS);

  result = makeNice( result );


  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "sa_poly_reduce", result ); }; /*4DEBUG*/

  return( result );
}
example
{ "EXAMPLE:"; echo = 2;
 ring AA = 0,(x,y,z),dp;
 matrix D[3][3]=0;
 D[1,2]=-z; D[1,3]=2*x; D[2,3]=-2*y;
 def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)
 poly f = 4*x*y+z^2-2*z; // a central polynomial
 sa_poly_reduce(f + 3*f*f + x, ideal(f) ); // should be just 'x'
}







/******************************************************/
// ::inStuff:: inCentralizer, inCenter, isCartan helpers
/******************************************************/


/******************************************************/
static proc inCentralizer_poly( poly p, ideal S )
  "
    if p in centralizer(S) => return 1, otherwise return 0
"
{
  poly f;

  for( int k = 1; k <= size(S); k++ )
    {
      f = S[k];

      if( NF( f * p - p * f, twostd(0) ) != 0 )
        {
          DBPrint( 1, "POLY: " + string (p) +
                   " is NOT in the centralizer of polynomial {" + string(f) + "}" );
          return (0);
        }
    }

  return( 1 );
}

/******************************************************/
static proc inCentralizer_list( def l, ideal S )
{
  for( int @i = 1; @i <= size(l); @i++ )
    {
      if( (typeof(l[@i])=="poly") or (typeof(l[@i]) == "int") or (typeof(l[@i]) == "number") )
        {
          if(! inCentralizer_poly(l[@i], S) )
            {
              return(0);
            }

        } else
          {
            if( (typeof(l[@i])=="list") or (typeof(l[@i])=="ideal") )
              {
                if(! inCentralizer_list(l[@i], S) )
                  {
                    return(0);
                  }
              }
          }
    }
  return(1);
}


/******************************************************************************/
// Checks the commutativity of polynomials of a with the polynomials in S
proc inCentralizer( def a, ideal S )
"
USAGE:   inCentralizer(E, S); E poly/list/ideal, S poly/ideal
RETURN:  integer, 1 if E is in the centralizer(S), 0 otherwise
PURPOSE: check whether the elements of E are in the centralizer(S)
EXAMPLE: example inCentralizer; shows examples
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "inCentralizer", a, S ); }; /*4DEBUG*/

  if( nameof( basering ) == "basering" )
    {
      //        ERROR( "No current ring!" );
    }


  int res;

  if( (typeof(a) == "poly") or (typeof(a) == "int") or (typeof(a) == "number") )
    {
      res = inCentralizer_poly(a, S);
    } else
      {
        if( (typeof(a)=="list") or (typeof(a)=="ideal") )
          {
            res = inCentralizer_list(a, S);
          } else
            {
              res = -1;
            }
      }

  if( res == -1 )
    {
      ERROR( "Wrong argument!" );
    }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "inCentralizer", res ); }; /*4DEBUG*/

  return (res);
}
example
{
  "EXAMPLE:";echo=2;
  ring R = 0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=-z;
  def r = nc_algebra(1,D); setring r; // the Heisenberg algebra
  poly f = x^2;
  poly a = z; // 'z' is central => it lies in every centralizer!
  poly b = y^2;
  inCentralizer(a, f);
  inCentralizer(b, f);
  list  l = list(1, a);
  inCentralizer(l, f);
  ideal I = a, b;
  inCentralizer(I, f);
  printlevel = 2;
  inCentralizer(a, f); // yes
  inCentralizer(b, f); // no
}

/******************************************************/
// Checks the centrality of a
proc inCenter( def a )
  "
USAGE:   inCenter(E); E poly/list/ideal
RETURN:  integer, 1 if E is in the center, 0 otherwise
PURPOSE: check whether the elements of E are central
EXAMPLE: example inCenter; shows examples
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "inCenter", a ); }; /*4DEBUG*/

  if( nameof( basering ) == "basering" )
    {
      //        ERROR( "No current ring!" );
    }

  int result = inCentralizer( a, variablesStandard() );

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "inCenter", result ); }; /*4DEBUG*/

  return( result );
}
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;
  def r = nc_algebra(1,D); setring r; // this is U(sl_2)
  poly p=4*x*y+z^2-2*z;
  inCenter(p);
  poly f=4*x*y;
  inCenter(f);
  list l= list( 1, p, p^2, p^3);
  inCenter(l);
  ideal I= p, f;
  inCenter(I);
}


/******************************************************/
// Checks whether f is a Cartan element.
proc isCartan( poly f )
"
USAGE:       isCartan(f); f poly
PURPOSE:     check whether f is a Cartan element.
RETURN:      integer, 1 if f is a Cartan element and 0 otherwise.
NOTE:        f is a Cartan element of the algebra A
@* if and only if for all g in A there exists C in K such that [f, g] = C * g
@* if and only if for all variables v_i there exist C in K such that [f, v_i] = C * v_i.
"
{
  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ BCall( "isCartan", f ); }; /*4DEBUG*/

  if( nameof( basering ) == "basering" )
    {
      //        ERROR( "No current ring!" );
    }


  ideal V = variablesStandard();

  int r = 1; poly v, g;

  for( int i = size(V); i > 0; i-- )
    {
      v = leadmonom(V[i]); // V[i] must be just a variable, but...

      g = NF( f*v - v*f, twostd(0) ); // [f, V[i]]

      if( size(g) > 0 )
        {
          if( size(g) > 1 ) // it is not just \alpha * v_i.
            {
              r = 0;
              break;
            }

          if( leadmonom(g) != v ) // g = \alpha * v_j, j != i.
            {
              r = 0;
              break;
            }

        } // else \alpha = 0
    }

  /*4DEBUG*/    if( defined( @@@DEBUG ) ){ ECall( "isCartan", r ); }; /*4DEBUG*/
  return( r );
}
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;
  def r = nc_algebra(1,D); setring r; // this is U(sl_2) with cartan - z
  isCartan(z); // yes!
  poly p=4*x*y+z^2-2*z;
  isCartan(p); // central elements are Cartan elements!
  poly f=4*x*y;
  isCartan(f); // no way!
  isCartan( 10 + p + z ); // scalar + central + cartan
}




/******************************************************/
/******************************************************/
// ::MainAliases:: The main non-static functions, visible to user are here. They are wrappers around basic functions.
/******************************************************/
/******************************************************/




/******************************************************/
// Computes the generators of the center of a basering
proc center( int D, list # )
"
USAGE:      center(D[, N]); D int, N optional int
RETURN:     ideal, generated by computed elements
PURPOSE:    computes subalgebra generators of the center up to degree D
NOTE:       In general, one cannot compute the whole center.
@* Hence, one has to specify a termination condition via arguments D and/or N.
@* If D is positive, only central elements up to degree D will be found.
@* If D is negative, the termination is determined by N only.
@* If N is given, the computation stops if at least N elements have been found.
@* Warning: if N is given and bigger than the actual number of generators,
@* the procedure may not terminate.
@* Current ordering must be a degree compatible well-ordering.
SEE ALSO:   centralizer; inCenter
EXAMPLE:    example center; shows an example
"
{
  if( nameof( basering ) == "basering" )
    {
      //        ERROR( "No current ring!" );
    }

  if( DefaultInt( # ) > 0 )
    {
      return( centerRed( D, # ) );
    }

  if( D >= 0 )
    {
      return( sa_reduce( centerVS(D) ) ); // Experimental! May be wrong!!!
    }

  ERROR( "Wrong arguments!" );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(x,y,z,t),dp;
  matrix D[4][4]=0;
  D[1,2]=-z;  D[1,3]=2*x;  D[2,3]=-2*y;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2) tensored with K[t]
  // find generators of the center of degree <= 3:
  ideal Z = center(3);
  Z;
  inCenter(Z); // check the result
  // find at least one generator of the center:
  ideal ZZ = center(-1, 1);
  ZZ;
  inCenter(ZZ); // check the result
}

/******************************************************/
// Computes the generators of the centralizer of S in a basering
proc centralizer( ideal S, int D, list # )
"
USAGE:      centralizer(F, D[, N]); F poly/ideal, D int, N optional int
RETURN:     ideal, generated by computed elements
PURPOSE:    computes subalgebra generators of centralizer(F) up to degree D
NOTE:       In general, one cannot compute the whole centralizer(F).
@* Hence, one has to specify a termination condition via arguments D and/or N.
@* If D is positive, only centralizing elements up to degree D will be found.
@* If D is negative, the termination is determined by N only.
@* If N is given, the computation stops if at least N elements have been found.
@* Warning: if N is given and bigger than the actual number of generators,
@* the procedure may not terminate.
@* Current ordering must be a degree compatible well-ordering.
SEE ALSO:   center; inCentralizer
EXAMPLE:    example centralizer; shows an example
"
{
  if( nameof( basering ) == "basering" )
    {
      //        ERROR( "No current ring!" );
    }

  if( DefaultInt( # ) > 0 )
    {
      return( centralizerRed( S, D, # ) );
    }

  if( D >= 0 )
    {
      return( sa_reduce( centralizerVS(S, D) ) ); // Experimental! May be wrong!!!
    }

  ERROR( "Wrong arguments!" );
}
example
{
  "EXAMPLE:"; echo = 2;
  ring AA = 0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=-z; D[1,3]=2*x; D[2,3]=-2*y;
  def A = nc_algebra(1,D); setring A; // this algebra is U(sl_2)
  poly f = 4*x*y+z^2-2*z; // a central polynomial
  f;
  // find generators of the centralizer of f of degree <= 2:
  ideal c = centralizer(f, 2);
  c;  // since f is central, the answer consists of generators of A
  inCentralizer(c, f); // check the result
  // find at least two generators of the centralizer of f:
  ideal cc = centralizer(f,-1,2);
  cc;
  inCentralizer(cc, f); // check the result
  poly g = z^2-2*z; // some non-central polynomial
  // find generators of the centralizer of g of degree <= 2:
  c = centralizer(g, 2);
  c;
  inCentralizer(c, g); // check the result
  // find at least one generator of the centralizer of g:
  centralizer(g,-1,1);
  // find at least two generators of the centralizer of g:
  cc = centralizer(g,-1,2);
  cc;
  inCentralizer(cc, g); // check the result
}


/*******************************************************
 // normally one should use this library together with ncalg.lib in the following way:

LIB "ncalg.lib";
def Usl3 = makeUsl(3); // U(sl_3)
setring Usl3;

// show current ring:
basering;

LIB "central.lib";

// easy example(few seconds), must compute two polynomials of degrees 2 and 3.
center(3);

kill Usl3;

def Ug2 = makeUg2(); // U(g_2)
setring Ug2;

// show current ring:
basering;

// easy example(few seconds), must compute one polynomial of degree 2.
center(2);

// hard example (~hours), must compute two polynomials of degrees 2 and 6.
center(6);

quit;
*******************************************************/