//////////////////////////////////////////////////////////////////////////////
version="version perron.lib 4.0.0.0 Jun_2013 "; // $Id$
category="Noncommutative";
info="
LIBRARY:  perron.lib         computation of algebraic dependences
AUTHOR:   Oleksandr Motsak   U@D, where U={motsak}, D={mathematik.uni-kl.de}

PROCEDURES:
perron(L[, D]);  relations between pairwise commuting polynomials

KEYWORDS:  algebraic dependence; relations
";

LIB "central.lib";

//////////////////////////////////////////////////////////////////////////////
proc perron( ideal L, list # )
"USAGE:  perron( L [, D] )
RETURN:  commutative ring with ideal `Relations`
PURPOSE: computes polynomial relations ('Relations') between pairwise
         commuting polynomials of L [, up to a given degree bound D]
NOTE:    the implementation was partially inspired by the Perron's theorem.
EXAMPLE: example perron; shows an example
"
{
  int N, D, i;
  N = size(L);
  if( N == 0 )
  {
    ERROR( "Input ideal must be non-zero!" );
  }
  intvec W; // weights
  for ( i = N; i > 0; i-- )
  {
    W[i] = deg(L[i]);
  }
  ////////////////////////////////////////////////////////////////////////
  D = -1;
  // Check whether the degree bound 'D' is given:
  if( size(#)>0 )
  {
    if ( typeof(#[1]) == typeof(D) )
    {
      D = #[1];
      if( D <= 0 )
      {
        ERROR( "An optional parameter D must be positive!" );
      }
    }
  }

  // , otherwise we try to estimate it according to Perron's Th:
  if( D < 0 )
  {
    D = 1;
    int d;
    int min = -1;

    for ( i = size(L); i > 0 ; i-- )
    {
      d = W[i];
      D = D * d;
      if( min == -1)
      {
        min = d;
      }
      else
      {
        if( min > d )
        {
          min = d;
        }
      }
    }
    if( (D == 0) or (min <= 0) )
    {
      ERROR( "Wrong set of polynomials!" );
    }
    D = D div min;
  }
  ////////////////////////////////////////////////////////////////////////
  def NCRING = basering;
  def CurrentField = ringlist( NCRING )[1];

  // We are going to construct a commutative ring in N variables F(i),
  // with the field specified by 'CurrentField':

  ring TEMPRING = 0, ( F(1..N) ), dp;
  list RingList = ringlist( TEMPRING );
  setring NCRING;

  if( !defined(RingList) )
  {
    list RingList = imap( TEMPRING, RingList );
  }
  RingList[1] = CurrentField;

  // New Commutative Ring with correct field!
  def COMMUTATIVERING = ring( RingList );

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

  setring COMMUTATIVERING; // we are in COMMUTATIVERING now

  ideal PBWBasis = PBW_maxDeg( D ); // All monomials of degree(!) <= D.

  // TODO: it would be better to compute weighted monomials of weight
  // <= W[1] \cdots W[N].

  setring NCRING; // and back to NCRING

  map Psi = COMMUTATIVERING, L; // F(i) \mapsto L[i]

  ideal Images = Psi( PBWBasis ); // Corresponding products of polynomials
                                  // from L

  // ::MAIN STEP:: // Compute relations in NC ring:
  def T = linearMapKernel( Images );

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

  // check the output of 'linearMapKernel':
  int t = 0;

  if( (typeof(T) != "module") and (typeof(T) != "int" ) )
  {
    ERROR( "Wrong output from function 'linearMapKernel'!" );
  }

  if( typeof(T) == "int" )
  {
    t = 1;
    if( T != 0 )
    {
      ERROR( "Wrong output from function 'linearMapKernel'!" );
    }
  }

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

  // Go back to commutative case in both cases:
  setring COMMUTATIVERING;

  ideal Relations; // And generate Relations:

  if( t == 0 ) // T is a module
  {
    module KER = imap( NCRING, T );
    Relations = linearCombinations( PBWBasis, KER );
  }
  else
  { // T == int(0) => all images are zero =>
    Relations = PBWBasis;
  }

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

  // we compute an std basis of the relations:

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

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

  // reduce everything in as far as possible
  Relations = simplify( groebner( Relations ), 1 + 2 + 8 );

  // restore options
  option( set, v );

  //    Relations;
  export Relations;

  return( COMMUTATIVERING );
}
example
{ "EXAMPLE:"; echo = 2;
  int p = 3;
  ring AA = p,(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 I = x^p, y^p, z^p-z, 4*x*y+z^2-2*z; // the center
  def RA = perron( I, p );
  setring RA;
  RA;
  Relations; // it was exported from perron to be in the returned ring.

  // perron can be also used in a commutative case, for example:
  ring B = 0,(x,y,z),dp;
  ideal J = xy+z2, z2+y2, x2y2-2xy3+y4;
  def RB = perron(J);
  setring RB;
  Relations;
  // one more test:
  setring A;
  map T=RA,I;
  T(Relations);  // should be zero
}