source: git/Singular/LIB/center.lib @ ec9b96

///////////////////////////////////////////////////////////////////////////////
version="$Id: center.lib,v 1.20 2006-05-11 18:53:25 motsak Exp $"
category="Noncommutative"
info="
LIBRARY:  center.lib      computation of central elements of GR-algebras
AUTHOR:  Oleksandr Motsak,        motsak@mathematik.uni-kl.de.
OVERVIEW: Computing elements of center and centralizers of sets of elements
in non-commutative algebras.

MAIN 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

AUXILIARY PROCEDURES:
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

KEYWORDS:  center; centralizer; cartan; reduce; centralize; PBW
"

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 ( deg(p) > 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 ( deg(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 = I - 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); int MaxDeg
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 A = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  ncalgebra(1,D); // this algebra is U(sl_2)

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


/******************************************************/
proc PBW_eqDeg( int Deg )
"USAGE: PBW_eqDeg(Deg); int Deg
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 A = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  ncalgebra(1,D); // 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); poly m
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 A = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  ncalgebra(1,D); // this algebra is U(sl_2)

  // 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( Basis, f); ideal Basis, poly f
PURPOSE: Apply Ad_{f} to every element of Basis
RETURN: ideal, Ad_{f}(Basis)
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 A = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  ncalgebra(1,D); // 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:
  module C = linearMapKernel( Image ); C;

  // Now we can compute the kernel of Ad_e by means 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 ); ideal Images
PURPOSE: Computes the kernel of a linear map: e_i \mapsto Images[i], 
@* where e_i are certain basis vectors
RETURN: syzygy module, or 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( std( 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 A = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  ncalgebra(1,D); // 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:
  module C = linearMapKernel( Image ); C;

  // Now we can compute the kernel of Ad_e by means 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 ); ideal Basis, module C
PURPOSE: computes linear combinations of Basis vectors with coefficients from 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 A = 0,(e,f,h),dp;
  matrix D[3][3]=0;
  D[1,2]=-h;  D[1,3]=2*e;  D[2,3]=-2*f;
  ncalgebra(1,D); // 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:
  module C = linearMapKernel( Image ); C;

  // Now we can compute the kernel of Ad_e by means 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") and (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
/******************************************************/



/******************************************************/
proc variablesStandard() // Returns an ideal of variables in a current base ring.
"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 A = 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 algebra is U(sl_2)
  // Variables in their natural order:
  variablesStandard(); 
}

/******************************************************/
proc variablesSorted() // Sorts variables into an ideal. This is a kind of heuristics!
"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 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 A = 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 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.
/******************************************************/
/******************************************************/





/******************************************************/
proc centralizeSet( ideal F, ideal V ) // HL 'core' function
"USAGE:      centralizeSet( F, V ); ideal F, ideal V
INPUT:      a finite set of elements F, vector space basis V
RETURN:     ideal, generated by base elements
PURPOSE:    computes a v.s. basis of the centralizer of F in v.s. V
NOTE:       Cen(F,V) is computed by the formula Cen(F[N],..,Cen(F[1],V)..)
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_4_1 = 0,(e(1..4)),dp;
  matrix D[4][4]=0;
  D[2,4] = -e(1);
  D[3,4] = -e(2);
  // This is $A_{41}$ - the first real Lie algebra of dimension $4$.
  ncalgebra(1,D); 
  
  ideal F = variablesSorted(); F;
 
  // the center of $A_{41}$ 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 PBW basis consisting of 
  //    monomials of exponent <= (1,2,1,0)
  ideal V = PBW_maxMonom( e(1) * e(2)^ 2 * e(3) ); 
 
  // 2. Compute the centralizer of F within vector space 
  //    spanned by these monomials:
  ideal C = centralizeSet( F, V ); C;
 
  inCenter(C);
}



/******************************************************/
proc centralizerVS( ideal F, int d )
  "USAGE:      centralizerVS( F, D ); ideal F, int D
RETURN:     ideal, generated by elements of degree <= D
PURPOSE:    computes a v.s. basis of the centralizer of 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 A = 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 algebra is U(sl_2)
  ideal F = x, y;
  // find all elements commuting with x and y of degree <= 4:
  ideal C = centralizerVS(F, 4); 
  C; 
  inCentralizer(C, F);
}




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




/******************************************************/
proc centerVS( int D )
"USAGE:      centerVS( D ); int D
RETURN:     ideal, generated by elements of degree <= D
PURPOSE:    computes a v.s. 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 A = 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 algebra is U(sl_2)
  // find all central elements of degree <= 4
  ideal Z = centerVS(4); 
  Z; 
  // note that the second element is the square of the first 
  // plus the multiple of the first:
  Z[2] - Z[1]^2 + 8*Z[1];
  inCenter(Z);
}


/******************************************************/
proc centralizerRed( ideal F, int D, list # )
"USAGE:      centralizerRed( F, D[, N] ); ideal F, int D[, int N]
RETURN:     ideal, generated by computed generators
PURPOSE:    computes a 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 has 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 = 1; // starting guess for Maximal Bounding Degree, 6
  int @Delta  = 1; // 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 A = 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 algebra is U(sl_2)
  ideal F = x, y;
  // find subalgebra generators degree <= 4 of an algebra of
  // all elements commuting with x and y:
  ideal C = centralizerRed(F, 4); 
  C; 
  inCentralizer(C, F);
}


/******************************************************/
proc centerRed( int D, list # )
"USAGE:      centerRed( D[, N] ); int D[, int N]
RETURN:     ideal, generated by computed generators
PURPOSE:    computes a 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 has 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 A = 0,(x,y,z),dp;
  matrix D[3][3]=0;
  D[1,2]=z;
  ncalgebra(1,D); // it is a Heisenberg algebra
  // find vector space basis of center of degree <= 3
  ideal VSZ = centerVS(3); 
  // There should be 3 degrees of z.
  VSZ;
  inCenter(VSZ);
  // find "minimal" central elements of degree <= 3
  ideal SAZ = centerRed(3); 
  // Only 'z' must be computed
  SAZ; 
  inCenter(SAZ);
}


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

/******************************************************/
static proc INTERRED( ideal S )
  "USAGE:      INTERRED( S ); ideal S
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 wrt 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 );    
}

/******************************************************/
proc sa_reduce( ideal V ) // 'subalgebra basis' computation
"USAGE:     sa_reduce(V); ideal V
RETURN:     ideal, generated by found elements
PURPOSE:    compute a s.a. basis of an algebra generated by V
NOTE:       May produce wrong result in quotient 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 )
    {
      // 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); 
        
      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 );            
        } 
        
      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 A = 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 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
}



/******************************************************/
proc sa_poly_reduce( poly p, ideal V ) // subalgebra reduction of a polynomial
"USAGE:      sa_poly_reduce(p, V); poly p, ideal V
RETURN:     polynomial, a reduction of p wrt V
PURPOSE:    computes a reduction of polynomial p wrt V
NOTE:       May produce wrong result in quotient 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 A = 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 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 poly {" + 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);
}


/******************************************************************************/
proc inCentralizer( def a, ideal S ) // Checks the commutativity of polynomials of a with the polynomials in S
"USAGE:   inCentralizer(a, S); a poly/list/ideal, S poly/ideal
RETURN:  integer, 1 if a in the centralizer(S), 0 otherwise
PURPOSE: check whether a is centralizing with respect to elements of 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;
  ncalgebra(1,D); // the Heisenberg algebra
  poly f = x^2;
  poly a = z; // 'z' is central => it lies in any 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
}

/******************************************************/
proc inCenter( def a ) // Checks the centrality of a
  "USAGE:   inCenter(a); a poly/list/ideal
RETURN:  integer, 1 if a in the center, 0 otherwise
PURPOSE: check whether a given element is 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;
  ncalgebra(1,D); // 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);
}


/******************************************************/
proc isCartan( poly f ) // Checks whether f is a Cartan element.
"USAGE:       isCartan(f); poly f
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 
@* iff for all g in A there exists C in K such that [f, g] = C * g
@* iff 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;
  ncalgebra(1,D); // 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.
/******************************************************/
/******************************************************/




/******************************************************/
proc center( int D, list # ) // Computes the generators of the center of a basering
"USAGE:      center(D[, N]); int D, int N
RETURN:     ideal, generated by elements of degree at most D
PURPOSE:    computes a 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 has 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 A = 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;
  ncalgebra(1,D); // this algebra is U(sl_2) tensored with K[t]
  ideal Z = center(3); // find all central elements of degree <= 3
  Z;
  inCenter(Z);
  ideal ZZ = center(-1, 1); // find one central element of the lowest degree
  ZZ;
  inCenter(ZZ);
}

/******************************************************/
proc centralizer( ideal S, int D, list # ) // Computes the generators of the centralizer of S in a basering
"USAGE:      centralizer(F, D[, N]); poly/ideal F, int D[, int N]
RETURN:     ideal, generated by computed generators
PURPOSE:    computes a 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 has 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 A = 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 algebra is U(sl_2)
  poly f = 4*x*y+z^2-2*z; // a central polynomial
  f;
  ideal c = centralizer(f, 2); // find all elements of the centralizer of f
  // of degree <= 2
  c;  // since f is central, the answer consists of generators of A
  inCentralizer(c, f);
  ideal cc = centralizer(f,-1,2); // find at least two elements of the centralizer of f
  cc;
  inCentralizer(cc, f);
  poly g = z^2-2*z; // some non-central polynomial
  c = centralizer(g, 2); // find all elements of the centralizer of g
  // of degree <= 2
  c;
  inCentralizer(c, g);
  centralizer(g,-1,1); // find the element of the lowest degree in the centralizer
  cc = centralizer(g,-1,2); // find at least two elements of the centralizer of g
  cc;
  inCentralizer(cc, g);
}


/*******************************************************
 // 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 "center.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;
*******************************************************/