// version string automatically expanded by CVS
version="$Id: powers.lib,v 1.3 2008-06-26 13:17:43 motsak Exp $";
category="Linear algebra";

info="
LIBRARY:   powers.lib  Symmetric and Exterior k-powers
AUTHOR:    Oleksandr Motsak, email: {U@D}, where U={motsak}, D={mathematik.uni-kl.de}.

KEYWORDS: symmetric power, exterior power

PROCEDURES:
    symmetricBasis(int, int)       return a basis of a k-th symmetric power
    symmetricPower(module,int)     return k-th symmetric power of a matrix

    exteriorBasis(int, int)        return a basis of a k-th exterior power
    exteriorPower(module,int)      return k-th exterior power of a matrix
       
";
    

LIB "general.lib"; // for sort
LIB "nctools.lib"; // for SuperCommutative

// thanks to Oleksandr Iena for his persistence ;-)

proc symmetricBasis(int n, int k)
"USAGE:    symmetricBasis(n, k); n int, k int
RETURN:   ring: the basis of the k-th symmetric power of a n-dim vector space
NOTE:     The output consists of a polynomial ring in n variables, together with the ideal called I.
          The ideal I is the basis of the k-th symmetric power of a n-dim vector space (ordered lexicographically).
          The char. of the returned ring is the same as of current basis ring or zero.
SEE ALSO: exteriorBasis
KEYWORDS: symmetric basis
EXAMPLE:  example symmetricBasis; shows an example"
{
  int p = 0;
  if( nameof(basering) != "basering" )
  {
    p = char(basering);
  }

  ring T = (p), (e(1..n)), dp;
  ideal I = simplify( NF(maxideal(k), std(0)), 1 + 2 + 8 );
  int N = ncols(I);
  I = sort(I)[1]; // lex
  I = I[N..1];

  export I;
  return(T);
}
example
{ "EXAMPLE:"; echo = 2;

def r = symmetricBasis(2, 3); setring r; r;
I; // basis of 3rd sym. power of a 2-dim v.s.
kill r;

ring r = (0),(a, b, c, d), dp; r;
def g = symmetricBasis(3, 2); setring g; g; I;

kill g, r;

ring r = (32003),(a, b, c, d), dp; r;
def g = symmetricBasis(4, 2); setring g; g; I;
}

proc exteriorBasis(int n, int k)
"USAGE:    exteriorBasis(n, k); n int, k int
RETURN:   ring: the basis of the k-th exterior power of a n-dim vector space
NOTE:     The output consists of a polynomial ring in n variables, together with the ideal called I.
          The ideal I is the basis of the k-th exterior power of a n-dim vector space (ordered lexicographically).
          The char. of the returned ring is the same as of current basis ring or zero.

SEE ALSO: symmetricBasis
KEYWORDS: exterior basis
EXAMPLE:  example exteriorBasis; shows an example"
{
  int p = 0;
  if( nameof(basering) != "basering" )
  {
    p = char(basering);
  }
    
  ring S = (p), (e(1..n)), dp;
  def T = SuperCommutative(); setring T;
  ideal I = simplify( NF(maxideal(k), std(0)), 1 + 2 + 8 );
  int N = ncols(I);
  I = sort(I)[1];
  I = I[N..1 ];

  export I;
  return(T);
}
example
{ "EXAMPLE:"; echo = 2;

def r = exteriorBasis(2, 3); setring r; r;
"Basis: ", I; // basis of 3rd sym. power of a 2-dim v.s.
kill r;

ring r = (0),(a, b, c, d), dp; r;
def g = exteriorBasis(3, 2); setring g; g; "Basis: ", I;

kill g, r;

ring r = (32003),(a, b, c, d), dp; r;
def g = exteriorBasis(4, 2); setring g; g; "Basis: ", I;
}



proc symmetricPower(module A, int k)
"USAGE:    symmetricPower(A, k); A module, k int
RETURN:   module: the k-th symmetric power of A
NOTE:     the chosen bases (ordered lexicographically) and
          temporary data  may be shown if printlevel is big enough
SEE ALSO: exteriorPower
KEYWORDS: symmetric power
EXAMPLE:  example symmetricPower; shows an example"
{
  int p = printlevel - voice + 2;
  def save = basering;

  int M = nrows(A);
  int N = ncols(A);

  int i, j;

  /////////////////////////////////////////////////////////////////
  def Tn = symmetricBasis(N, k); setring Tn;
  ideal Ink = I;
  dbprint(p-3, "Temporary Source Ring", basering);
  dbprint(p, "S^k(Source Basis)", Ink);

  /////////////////////////////////////////////////////////////////
  def Tm = symmetricBasis(M, k); setring Tm;
  ideal Imk = I;
  dbprint(p-3, "Temporary Image Ring", basering);
  dbprint(p, "S^k(Image Basis)", Imk);

  /////////////////////////////////////////////////////////////////
  def Rm = save + Tm;
  setring Rm;
  dbprint(p-2, "Temporary Working Ring", Rm);

  module A = imap(save, A);

  ideal B; poly t;

  for( i = N; i > 0; i-- )
  {
    t = 0;
    for( j = M; j > 0; j-- )
    {
      t = t + A[i][j] * var( nvars(save) + j); // tensor product!!!
    }

    B[i] = t;
  }

  dbprint(p-1, "Matrix of simgle images", B);

  map TMap = Tn, B; ideal C = TMap(Ink); // apply S^k A to Source basis vectors... (Ink)
  C = NF(C, std(0));

  dbprint(p-1, "Image Matrix: ", C);

  // and write it in Image basis (Imk)

  ideal Imk = imap(Tm, Imk);

  module D; poly lm; vector tt;

  for( i = ncols(C); i > 0; i-- )
  {
    t = C[i];
    tt = 0;

    while( t != 0 )
    {
      lm = leadmonom(t);
      //    lm;
      for( j = ncols(Imk); j > 0; j-- )
      {
        if( lm / Imk[j] != 0 )
        {
          tt = tt + (lead(t) / Imk[j]) * gen(j);
          break;
        }
      }
      t = t - lead(t);
    }

    D[i] = tt;
//    tt;
  }

  // D; print(D);


  setring save;

  // basering;

  module AA = imap(Rm, D);

//  nrows(AA) = Nk; // ????
//  ncols(AA) = Mk;

//  Nk, Mk;

//  AA[Mk] = AA[Mk] * 1;

  return(AA);
}
example
{ "EXAMPLE:"; echo = 2;
int save = printlevel; printlevel = 1;

ring r = (0),(a, b, c, d), dp; r;

module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1);

print(symmetricPower(A, 2));
print(symmetricPower(A, 3));

module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2);

print(symmetricPower(B, 2));
print(symmetricPower(B, 3));

kill r;
ring r = (0),(a, b, c, d), dp;def g = SuperCommutative();setring g;

module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1);

print(symmetricPower(A, 2));
print(symmetricPower(A, 3));

module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2);

print(symmetricPower(B, 2));
print(symmetricPower(B, 3));

kill r;

printlevel = save;
}



proc exteriorPower(module A, int k)
"USAGE:    exteriorPower(A, k); A module, k int
RETURN:   module: the k-th exterior power of A
NOTE:     the chosen bases and temporary data
          may be shown if printlevel is big enough.
          Last rows may be invisible if zero.
SEE ALSO: symmetricPower
KEYWORDS: exterior power
EXAMPLE:  example exteriorPower; shows an example"
{
  int p = printlevel - voice + 2;
  def save = basering;

  int M = nrows(A);
  int N = ncols(A);

  int i, j;

  /////////////////////////////////////////////////////////////////
  def Tn = exteriorBasis(N, k); setring Tn;
  ideal Ink = I;
  dbprint(p-3, "Temporary Source Ring", basering);
  dbprint(p, "E^k(Source Basis)", Ink);

  /////////////////////////////////////////////////////////////////
  def Tm = exteriorBasis(M, k); setring Tm;
  ideal Imk = I;
  dbprint(p-3, "Temporary Image Ring", basering);
  dbprint(p, "E^k(Image Basis)", Imk);


  def Rm = save + Tm;
  setring Rm;
  dbprint(p-2, "Temporary Working Ring", Rm);

  module A = imap(save, A);

  ideal B; poly t;

  for( i = N; i > 0; i-- )
  {
    t = 0;
    for( j = M; j > 0; j-- )
    {
      t = t + A[i][j] * var( nvars(save) + j); // tensor product!!!
    }

    B[i] = t;
  }

  dbprint(p-1, "Matrix of simgle images", B);

  map TMap = Tn, B; ideal C = TMap(Ink); // apply S^k A to Source basis vectors... (Ink)
  C = NF(C, std(0));

  dbprint(p-1, "Image Matrix: ", C);

  // and write it in Image basis (Imk)

  ideal Imk = imap(Tm, Imk);

  module D; poly lm; vector tt;

  for( i = ncols(C); i > 0; i-- )
  {
    t = C[i];
    tt = 0;

    while( t != 0 )
    {
      lm = leadmonom(t);
      //    lm;
      for( j = ncols(Imk); j > 0; j-- )
      {
        if( lm / Imk[j] != 0 )
        {
          tt = tt + (lead(t) / Imk[j]) * gen(j);
          break;
        }
      }
      t = t - lead(t);
    }

    D[i] = tt;
//    tt;
  }

  // D; print(D);


  setring save;

  // basering;

  module AA = imap(Rm, D);

//  nrows(AA) = Nk; // ????
//  ncols(AA) = Mk;

//  Nk, Mk;

//  AA[Mk] = AA[Mk] * 1;

  return(AA);
}
example
{ "EXAMPLE:"; echo = 2;
  int save = printlevel; printlevel = 1;

  ring r = (0),(a, b, c, d), dp; r;

  module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1);

  print(exteriorPower(A, 2));
  print(exteriorPower(A, 3));

  module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2);

  print(exteriorPower(B, 2));
  print(exteriorPower(B, 3));

  kill r;
  ring r = (0),(a, b, c, d), dp;def g = SuperCommutative();setring g;

  module A = a*gen(1), b * gen(1), c*gen(1), d * gen(1);

  print(exteriorPower(A, 2));
  print(exteriorPower(A, 3));

  module B = a*gen(1) + c* gen(2), b * gen(1) + d * gen(2);

  print(exteriorPower(B, 2));
  print(exteriorPower(B, 3));

  kill r;

  printlevel = save;
}