Package $contrib/invariants -- Invariant Theory

/*-----[    suggestions for use:        ]-----*\

     Run CoCoA and type:

     Alias Inv := $contrib/invariants;
     Inv.About();
     Inv.Man();

\*--------------------------------------------*/

Alias Inv := $contrib/invariants;

Define About()
  PrintLn "
  Topic    : Invariant Theory
  KeyWords : Invariant Theory, Linearly Reductive Groups
  Author   : A.Del Padrone
  E-mail   : delpadrone@dima.unige.it
  Version  : CoCoA 4.2
  Date     : 23 July 2002
  ";
EndDefine; -- About
  
------[ Manual ]--------

Define Man()
  PrintLn "
Suggested alias for this package:
    Alias Inv := $contrib/invariants;

SYNTAX

    Inv.IsInvAlgGroup(F: POLY, EquationsG: LIST, MatActionG: MAT ): BOOL
    Inv.GensInvIdealLinRedGroup(EquationsG: LIST, MatActionG: MAT): LIST
    Inv.GensInvAlgLinRedGroup(EquationsG: LIST, MatActionG: MAT): LIST
    Inv.GensInvAlgLinRedGroup(EquationsG: LIST, MatActionG: MAT,
                              PrintInfo: BOOL): LIST

DESCRIPTION

Let K be the field Q of the rational numbers, or the field Z/(p) with
p elements, where p is a prime integer.  Given a rational action of a
linear algebraic (K-affine) group G (embedded in some affine space K^L
as an affine algebraic set) on a finite dimensional vector space V,
over K, consider the induced (rational) action of G on the K-algebra
K[V] of the polynomial functions over V, if N is the dimension of V
then K[V]=K[x[1..N]] is the ring of N-variate polynomials over K. We
are interested in those polynomials of K[V] which are invariant under
the action of G.

Here a Linear Reductive Algebraic Group G will always be represented
by the radical ideal I(G)=(H_1,..,H_T) of K[z[1..L]] such that G is
the set of zeros of I(G) in K^L.

The (linear) action of G on a N-dimensional K-module V will be encoded
via a rational representation R=R(G,V):G --> GL(V)=GL_N(K) < Mat_N(K).
Then R is represented by a polynomial square matrix A=A(G,V) in
Mat_N(K[z[1..L]]), where N is the dimension of V; if A=(A_[H,K]), with
A_[H,K] in K[z[1..L]], then the action of G on V is given by the
(modulo I(G)) K-automorphisms Phi_E of K[V]=K[x[1..N]] (E is any
element of G) defined by
 Phi_E: x[I] --> Sum(K=1 To L; A_[K,I](E) * x[K]) 
Henceforth a polynomial F in K[X[1..N]] is G-invariant iff 'F=F(A*x)
modulo I(G)'  i.e. NF(F-F(A*x), I(G))=0.
 
The set K[V]^G of G-invariant polynomials of K[V] is a K-subalgebra of
K[V], moreover it is a graded subalgebra w.r.t. the standard
graduation on K[V] because of we consider linear actions.
It is known that K[V]^G is always finitely generated if the group G is
REDUCTIVE, otherwise K[V]^G can be not Noetherian (both facts are due
to NAGATA).

This package deals with two problems:
1) decide whether a polynomial F in K[V] is G-invariant,
2) calculate a minimal system of homogeneous generators for K[V]^G 
   when G is a LINEARLY REDUCTIVE. The algorithm used here is based 
   on a result of HARM DERKSEN.

LITERATURE

All can be found in:
 Derksen, Harm and Kemper, Gregor;
 Computational Invariant Theory;
 Encyclopaedia of Mathematical Sciences 130,
 Springer-Verlag, Berlin, Heidelberg, New York 2002.   
see also the reference given there.

MAIN FUNCTIONS

The function: 

   Inv.GensInvAlgLinRedGroup(EquationsG, MatActionG)
   Inv.GensInvAlgLinRedGroup(EquationsG, MatActionG, PrintInfo)

computes a minimal system of homogeneous generators for the K-algebra
K[V]^G, when G is a linearly reductive algebraic group embedded in
some affine space, K^l, as an algebraic set and:

- EquationsG: a list of polynomials in K[z[1..l]] which generate
  the vanishing ideal of G as an algebraic set; 
- MatActionG: a N><N polynomial matrix in Mat_N( K[z[1..l]])
  defining the action of G on a K-vector space V s.t. V has dimension
  N over K.
- PrintInfo: if TRUE progress information is displayed

The function:

   Inv.IsInvAlgGroup(F, EquationsG, MatActionG)

tests if the polynomial F in K[V] is G-invariant.

The function:

   Inv.GensInvIdealLinRedGroup(EquationsG, MatActionG)

computes a minimal system of homogeneous generators of the ideal
generated in K[V] by all the homogeneous invariants of positive degree,
the zero set in V of this ideal is the HILBERT NULLCONE.

>EXAMPLE<

Alias Inv := $contrib/invariants;  -- Suggested alias

Use GRing ::= Q[z];

-- G = C_2 the cyclic group of order two, V=Q^2;
EquationsG := [z^2-1];   /* G is an affine group embedded in Q^1 */

MatActionG := Mat([[(z+1)/2, (1-z)/2],
                   [(1-z)/2, (z+1)/2]]);
/* it defines the rational action of G on Q[x[1..2]] which exchanges
x[1] with x[2]*/

Inv.MakeInvRing(Len(MatActionG));
/* This is to create the ambient ring for the invariants */

F := InvRing :: x[1]^2-2x[2];
Inv.IsInvAlgGroup(F, EquationsG, MatActionG);

Inv.GensInvAlgLinRedGroup(EquationsG, MatActionG);
Inv.GensInvAlgLinRedGroup(EquationsG, MatActionG, TRUE);
";
EndDefine;-- Man

------[ Functions ]--------
  
Define MakeInvRing(DimRep)
  InvRing ::= CoeffRing[x[1..DimRep]];
  PrintLn '-- Created: InvRing ::= ', InvRing, ';';
EndDefine;--MakeInvRing(DimRep)


Define DegSort(A,B)
  Return Deg(A)<=Deg(B);
EndDefine;--DegSort(A,B)

  
/*Auxiliary Function:
Compute a vector space basis of Id in degree D
If Id is an homogeneous ideal w.r.t. the standard graduation*/

Define HomogIdealDegreeVectorBasis(Id,D)
  G := SortedBy(Gens(Minimalized(Id)), Function(Inv.PkgName(),'DegSort'));
  InDegId := Deg(G[1]);
  If D < InDegId Then Return [0] End;
  If D = InDegId Then Return [F In G| Deg(F)=InDegId] End;
  GUpD := [F In G | Deg(F)<=D];
  RG :=  ConcatLists([F*Support(DensePoly(D-Deg(F)))| F In GUpD]);
  U := Ideal(RG);
  Minimalize(U);
  Return Gens(U);
EndDefine;--HomogIdealDegreeVectorBasis(Id,D)


/*Auxiliary Function:
Exhaustive search of natural solutions of a linear equations with natural
coefficients: inductive version*/

Define NatSolNatLinEq(Eq,SC,EC,TV)
  /*Eq=[[C_1,..,C_K],D], where Eq[1]=[C_1,..,C_K] is the coefficient
    list of the equation (we suppose 1<=C_1<=..<=C_K), Eq[2]=D is the
    known term of the eq.;SC=first coefficient to consider; EC=K;TV=test
     "vector" in N^K;*/
  SS := [];/*SS=solution set*/
  /*These are trivial cases*/
  If SC>EC Then
    Return SS
  EndIf;/*an exit condition*/
  If Eq[2]=0 Then
    For I := SC To EC Do
      TV[I] := 0;
    End;
    Append(SS, TV);
    Return SS
  End;/*SS=[only the zero "vector"]*/
  If Comp(Eq[1],SC)>Eq[2] Then
    Return SS
  End;/*impossible eq.*/

  N := Max([I In SC..EC | Comp(Eq[1],I)<=Eq[2]]);
  If Mod(Eq[2],GCD([Comp(Eq[1],I)|I In SC..N]))<>0 Then
    Return SS
  End;/*impossible eq.*/

  If N<EC Then
    For I := (N+1) To EC Do
      TV[I] := 0;
    EndFor;
  EndIf; /*if I>N then the I-component of every solution is necessarily
           zero*/

  If SC=N Then
    If Mod(Eq[2],Comp(Eq[1],SC))=0 Then
      TV[SC] := Div(Eq[2],Comp(Eq[1],SC));
      Append(SS,TV)
    EndIf;
    Return SS
  End;/*this is the non-trivial base case: Eq[1][SC]X[SC]=Eq[2]*/

  ToTry := [];
  /*The possible values to test for X[SC]: we have necessarily
     X[SC] In 0..Div(Eq[2],Comp(Eq[1],SC)), but not all these values
     are possible...*/
  For I := 0 To (Div(Eq[2],Comp(Eq[1],SC))-1) Do
    If Comp(Eq[1],SC+1)<=(Eq[2]-Comp(Eq[1],SC)*I) Then
      K := Max([J In (SC+1)..N | Comp(Eq[1],J)<=(Eq[2]-Comp(Eq[1],SC)*I)]);
      If Mod((Eq[2]-Comp(Eq[1],SC)*I),GCD([Comp(Eq[1],H)|H In (SC+1)..K]))=0
  Then Append(ToTry, I)
      EndIf;
    EndIf;
  EndFor;
  If Mod(Eq[2],Comp(Eq[1],SC))=0 Then
    Append(ToTry,Div(Eq[2],Comp(Eq[1],SC)))
  EndIf;
  /*Now the recursive call...*/
  Foreach I In ToTry Do
    TV[SC] := I;
    SS := Concat(SS,
                 $.NatSolNatLinEq([Eq[1],(Eq[2]-Comp(Eq[1],SC)*I)],SC+1,EC,TV));
  EndForeach;
  Return SS;
EndDefine;--NatSolNatLinEq(Eq,SC,EC,TV)

  /*Auxiliary Function:
     We want to find a basis for Id_D modulo the vector space spanned by the
     products, with degree D, of those minimal generators of the homogeneous ideal
     Id in lower degrees which are invariant: i.e. If Id is minimally generated
     by {F1,F2,...,Ft}, with invariants G1,...,Gs in {F1,F2,...,Ft} of degrees
     D1,...,Ds I want to find a basis of the vector space
     Id_D/<{G1^X1***Gs^Xs | X1D1+...+XsDs=D}>*/

Define QuotientIdealBasis(Id,ToElim,D)
  L := SortedBy(ToElim, Function(Inv.PkgName(),'DegSort'));
  NewToElim := [F In L | Deg(F)<=D];
  EC := Len(NewToElim);
  If EC>=1 Then
    Eq := [[Deg(NewToElim[I])| I In 1..EC ],D];
    SC := 1;
    TV := NewList(EC,0);
    Exp := Inv.NatSolNatLinEq(Eq,SC,EC,TV);
    GensW := 
    [Product([ NewToElim[I]^S[I] | I In 1..EC])| S In Exp];
    RedundantBasis := Inv.HomogIdealDegreeVectorBasis(Id,D);
    M := Mat([GenRepr(F,Ideal(RedundantBasis))| F In GensW]);
    U := LinKer(M);
    B := [ScalarProduct(V,RedundantBasis)| V In U];
    Return B
  Else
    Return Inv.HomogIdealDegreeVectorBasis(Id,D)
  EndIf;
EndDefine;--QuotientIdealBasis(Id,ToElim,D)

  /*Auxiliary Function:
     Remember that a polynomial
     F in K[x[1..N]] is G-invariant iff "F=F(A*x) modulo I(G)"
     i.e. NF(F-F(A*x), I(G))=0 */

Define  Extended(F,X,A,V)
 Return F-Subst(F,[[X[I],Comp(Flatten(List(A*V)),I)]|I In 1..Len(A)])
EndDefine;--Extended(F,X,A,V)

  /*Auxiliary Function: Test of invariance w.r.t. an algebraic group*/

Define IsInvAlgGroup(F,EquationsG, MatActionG)
  EmbDimG := Len(Indets());
  VDim := Len(MatActionG);
  OverRingTest ::= CoeffRing[z[1..EmbDimG],x[1..VDim]];
  Using OverRingTest Do
    W := RMap(z);
    H := Image(EquationsG,W);
    A := Image(MatActionG,W);
    WW := RMap(x);
    F := Image(F,WW);
    X := x;
    V := Transposed(Mat(Vector(X)));
    B := (NF(Inv.Extended(F,X,A,V),Ideal(H))=0);
  EndUsing;
  Destroy OverRingTest;
  Return B;
EndDefine;--IsInvAlgGroup(F,EquationsG, MatActionG)

  /*Auxiliary Function*/

Define GrafIdeal(H,Y,A,V)
  Return Ideal(Concat(H,Y-Flatten(List(A*V))));
EndDefine;--GrafIdeal(H,Y,A,V)

  /*Auxiliary Function:
     Compute a minimal system of generators of the ideal generated by
     all the invariants of positive degree of a linearly reductive
     algebraic group: i.e. 
     I is the ideal of K[V]=K[x[1],...,x[DimV]] generated by
     (K[V]^G)_+={F in K[V]^G | Deg(F)>0} */

Define GensInvIdealLinRedGroup(EquationsG, MatActionG)
  EmbDimG := Len(Indets());-- I suppose that the current ring is the ring
  --of the affine space in which G is embedded
  VDim := Len(MatActionG);--this gives the dimension of the vector space V
  --over which G acts
  GrafOverRing ::= CoeffRing[z[1..EmbDimG],x[1..VDim],y[1..VDim]];
  /*this is the ring of the affine space in which the grafic, T, of
   the morphisms (g,v)-->(v,g*v) is embedded N.B. we cosider T in
           G >< V >< V instead of T in (G >< V) >< (V >< V)*/
  Using GrafOverRing Do
    W := RMap(z);
    H := Image(EquationsG,W);
    A := Image(MatActionG,W);
    V := Transposed(Mat(Vector(x)));
    Y := y;
    GrafIdeal := Inv.GrafIdeal(H,Y,A,V);
    GB := Gens(Elim(z,GrafIdeal));
  EndUsing;
  Using InvRing Do
    Target := Concat(NewList(EmbDimG,0), x, NewList(VDim,0));
    WW := RMap(Target);
    EvalGB := Image(GB,WW);
    GensMin := Gens(Minimalized(Ideal(EvalGB)));
    GensInvIdeal := SortedBy([P | P In GensMin And P<>0], 
			     Function(Inv.PkgName(),'DegSort'));
  EndUsing;
  Destroy GrafOverRing;
  Return GensInvIdeal;
EndDefine;--GensInvIdealLinRedGroup(EquationsG,MatActionG)

  /*Main Function:
     Compute a minimal set of homogeneous generators for the algebra of
     invariants of a linearly reductive algebraic group*/

--Define GensInvAlgLinRedGroup(EquationsG, MatActionG,<PrintInfo>)
Define GensInvAlgLinRedGroup(...)
  EquationsG := ARGV[1];
  MatActionG := ARGV[2];
  If Len(ARGV)=2 Then PrintInfo := FALSE; Else PrintInfo := ARGV[3]; EndIf;
  EmbDimG := Len(Indets());
  VDim := Len(MatActionG);
  OverRing ::= CoeffRing[z[1..EmbDimG],x[1..VDim]];
  GensInvIdeal := Inv.GensInvIdealLinRedGroup(EquationsG,MatActionG);
  If PrintInfo Then PrintLn 'Hilbert NullCone=', GensInvIdeal; EndIf;
  InvIn := [F In GensInvIdeal |
    Inv.IsInvAlgGroup(F,EquationsG, MatActionG)=TRUE];
  If PrintInfo Then 
    PrintLn 'Not Invariants=', Diff(GensInvIdeal,InvIn) 
  EndIf;
  Using InvRing Do
    If Len(InvIn)<Len(GensInvIdeal) Then
      ToReplace := Set([Deg(F)| F In Diff(GensInvIdeal,InvIn)]);
      N := Len(ToReplace);
      Id := Ideal(GensInvIdeal);
      For I := 1 To N Do
	If PrintInfo Then 
	  PrintLn 'Searching Invariants in Degre d=', ToReplace[I] 
	EndIf;
	InvUpD := [F In InvIn | Deg(F)<=ToReplace[I]];
	VectorBasis := Inv.QuotientIdealBasis(Id,InvUpD,ToReplace[I]);
	Using OverRing Do
	  W := RMap(z);
	  H := Image(EquationsG,W);
	  A := Image(MatActionG,W);
	  V := Transposed(Mat(Vector(x)));
	  VectorBasis := BringIn(VectorBasis);
	  InvIn := BringIn(InvIn);
	  InvVB := [F In VectorBasis | NF(Inv.Extended(F,x,A,V),Ideal(H))=0];
	  Missing := Diff(VectorBasis,InvVB);
	  If Len(Missing)<>0 Then
	    P := [NF(Inv.Extended(F,x,A,V),Ideal(H)) | F In Missing];
	    T := Set(ConcatLists([Support(F)|F In P]));
	    If PrintInfo Then 
	      PrintLn 'Dimensions of the Linear System:' 
	    EndIf;
	    UnKnowns := Len(VectorBasis);
	    Equations := Len(T);
	    If PrintInfo Then PrintLn 'UnKnowns=', UnKnowns; EndIf;
	    If PrintInfo Then PrintLn 'Equations=', Equations; EndIf;
	    M := Mat([ [ CoeffOfTerm(TT, X) | X In P ] | TT In T]);
	    M := LinKer(M);
	    If PrintInfo Then PrintLn 'Basis of solutions=', M; EndIf;
	    InvIn := Concat(InvIn,
			    [ScalarProduct(S,VectorBasis)| S In M]);
	  EndIf;
	EndUsing/*OverRing*/;
	InvIn := BringIn(InvIn);
	Destroy OverRing;
      EndFor;
      GensInv := InvIn;
    Else
      GensInv := GensInvIdeal;
    EndIf;
  EndUsing/*InvRing*/;
  If PrintInfo Then PrintLn 'Invariant Algebra=', GensInv; EndIf;
  Return GensInv; 
EndDefine;--GensInvAlgLinRedGroup(...)

EndPackage;