##############################################################################
##
##  Fortran friendly group representations in GAP
##
##  by Jack Schmidt (c) 2012 (CC BY-SA 3.0)
##
##
##  Requires GAP 4.5.x, repsn package, and Frank's futil package.
##
##  Download http://www.math.rwth-aachen.de/~Frank.Luebeck/gap/FUtil/
##  and change "Exp" to "ExpFutil" in two places in gap/Decimal.gi.
##
##  Functions:
##    RealCharacters,                - irred characters that have a real rep
##    RealRepresentation,            - will double dimension if needed
##    UnitaryRepresentation,         - unitary rep
##    RealOrthogonalRepresentation,  - will double dimension if needed
##    FortranStringOfMatrix,         - 2 line string rep for easy READ
##
##  Example usage: ExampleUsage( SymmetricGroup( 4 ) );
##  (see source at bottom to customize it; it outputs a bunch of text files
##  in the current directory, one per matrix)
##
##  Warning: Unitary representations are somewhat foreign to GAP (which instead
##  prefers to find an inner product for which an already given representation
##  becomes unitary).  Don't use them for much other than exporting to fortran.

##############################################################################
##
##  Return a list of characters that have real representations without 
##  dimension doubling (so not complex, and not quaternionic).
##
##  Input: group, a finite group, preferably permutation or pc-group.
##  Output: chis, the list of irreducible real characters.
##

RealCharacters := group ->
  Irr( CharacterTable( group ) ){ 
    Positions( Indicator( CharacterTable( group ), 2 ), 1 ) };

IsIrreducibleRealCharacter := chi ->
  IsClassFunction( chi ) and 
  IsIrreducibleCharacter( chi ) and
  Indicator( UnderlyingCharacterTable( chi ), [ chi ], 2 )[1] = 1;

##############################################################################
##
##  RealRepresentation - convert real irred character into real rep
##
##  Input: chi, a real irreducible character, such as returned 
##    by RealCharacters
##  Output: rep, a homomorphism from the group to a matrix group over real
##    subfield of the Cyclotomic field
##

LoadPackage("repsn");

RealRepresentation := function( chi )
  local rep, src, gen, fld, sub, bas, img, ran, image, gens, imgs, dim,
    vectorSpace, sim, seed, basis;
  Assert( 0, IsClassFunction( chi ) and IsCharacter( chi ) );
  if not IsIrreducibleCharacter( chi ) then
    # split into irreducible components, and then direct sum
    Error("Not done, but easy");
  elif not IsIrreducibleRealCharacter( chi ) then
    # impossible to write over real field, so double dimension
    rep := IrreducibleAffordingRepresentation( chi );
    src := Source( rep );
    gen := GeneratorsOfGroup( src );
    fld := FieldOfMatrixGroup( Image( rep ) );
    sub := NF( Conductor( fld ), Concatenation( GaloisStabilizer( fld ), [ -1 ] ) );
    bas := Basis( AsField( sub, fld ) );
    img := List( gen, x -> BlownUpMat( bas, Image( rep, x ) ) );
    ran := Group( img );
    SetSize( ran, Size( Image( rep ) ) );
    return GroupHomomorphismByImagesNC( src, ran, gen, img );
  else # try random conjugations to get inside the real field
    rep := IrreducibleAffordingRepresentation( chi );
    image := Image( rep );
    gens := GeneratorsOfGroup( Source( rep ) );
    imgs := List( gens, x -> x^rep );
    dim := DimensionOfMatrixGroup( image );
    vectorSpace := FieldOfMatrixGroup( image ) ^ dim;
    sim := One( imgs[1] );
    repeat
      seed := Random( vectorSpace );
      basis := List( [1..dim], i -> seed^Random( image ) );
      if RankMat( basis ) < dim then continue; fi;
      sim := basis^-1;
    until ForAll( Flat( List( imgs, mat -> mat^sim ) ),
      entry -> entry = ComplexConjugate( entry ) );
    image := image^sim;
    return GroupHomomorphismByImagesNC(
      Source( rep ), image^sim,
      gens, List( imgs, x -> x^sim ) );
  fi;
end;

##############################################################################
##
##  CholeskyFactors - indefinite Cholesky factorization
## 
##  Input: A, a Hermitian matrix
##  Output: [ L, D ], such that L is lower triangular, D is diagonal,
##    and A = L D L^H
##

CholeskyFactors := function( A )
  local L, D, i, j, n;
  n := Size( A );
  Assert( 0, IsMatrix( A ) and DimensionsMat( A ) = [ n, n ] and
    A = TransposedMat( List( A, row -> List( row, ent -> ComplexConjugate( ent ) ) ) ) );
  D := [];
  L := [];
  for i in [ 1 .. n ] do
    L[i] := [];
    for j in [ 1 .. i-1 ] do
      L[i][j] := ( A[i][j] - Sum( [1..j-1], k -> L[i][k]*ComplexConjugate(L[j][k])*D[k] ) )/D[j];
    od;
    L[i][i]:=1;
    for j in [ i+1 .. n ] do
      L[i][j] := 0;
    od;
    D[i] := A[i][i] - Sum( [1..i-1],k->L[i][k]*ComplexConjugate(L[i][k])*D[k] );
  od;
  return [L,D];
end;

##############################################################################
##
##  CholeskyFactor - definite Cholesky factorization
##
##  Input: A, a Hermitian positive definite matrix over the cyclotomics
##  Output: L, a lower triangular matrix over the cyclotomics such that
##    A = L L^H
##

LoadPackage("futil");

SquareRoot := function( R )
  return SqrtDecimalApproximation( DecimalApproximation( R, 51 ), 25 );
end;

CholeskyFactor := function( A )
  local LD, L, D, Li, Di, i, j, n;
  LD := CholeskyFactors( A );
  L := LD[1];
  D := LD[2];
  Li := L^-1;
  Di := List(D,INV);
  D := List( D, SquareRoot );
  Di := List( Di, SquareRoot );
  n := Size( A );
  for i in [1..n] do
    for j in [1..n] do
      L[i][j] := DecimalApproximation(L[i][j],Precision(D[j])+2)*D[j];
      Li[i][j] := Di[i]*DecimalApproximation(Li[i][j],Precision(Di[i])+2);
    od;
  od;
  return [ L, Li ];
end;

InstallMethod( ComplexConjugate, [ IsDecimalApproximation ], x -> x );

##############################################################################
##
##  UnitaryForm - the Gram matrix of the G-invariant inner product
##
##  Input: rep, a homomorphism into a matrix group of cyclotomics
##  Output: mat, a Gramm matrix of the G-invariant inner product
##

UnitaryForm := function( rep )
  local img, dim, one;
  img := Image( rep );
  dim := DimensionOfMatrixGroup( img );
  one := One( img );
  return List( [1..dim], i -> List( [1..dim], j ->
    Sum( img, mat -> ScalarProduct( one[i]*mat, ComplexConjugate( one[j]*mat ) ) ) ) );
end;

##############################################################################
##
##  UnitaryRepresentation - rewrite a rep as a unitary rep
##

UnitaryRepresentation := function( rep )
  local sim, inv;
  sim := CholeskyFactor( UnitaryForm( rep ) );
  inv := sim[2];
  sim := sim[1];
  return GroupHomomorphismByFunction( Source(rep), Range(rep),
    elt -> inv *
    List(Image(rep,elt),row->List(row,ent->
    DecimalApproximation(ent,Precision(sim[1][1])+2))) * sim );
end;

##############################################################################
##
##  IsRealOrthogonalRep
##

UnitaryRepError := rep ->
  Sum( Source(rep), function( elt )
    local mat, err;
    mat := Image(rep,elt);
    err := mat*TransposedMat( List( mat, row -> List( row, ComplexConjugate ) ) ) - One(mat);
    return SqrtDecimalApproximation( Sum( err, row -> Sum( row, ent ->
      ent*ComplexConjugate(ent) ) ), Int(Precision(mat[1][1])/2) ); end ) / Size( Source(rep) );

IsUnitaryRep := function( rep )
  if IsDecimalApproximation( Image( rep, One( Source( rep ) ) )[1][1] )
  then return TruelyLessDecimal( UnitaryRepError( rep ), DecimalApproximation( 1/10^5, 7 ) );
  else return IsGroupHomomorphism(rep) and
    IsFinite( Image( rep ) ) and
    IsMatrixGroup( Image( rep ) ) and
    Characteristic( FieldOfMatrixGroup( Image( rep ) ) ) = 0 and
    ForAll( GeneratorsOfGroup( Image( rep ) ), img ->
      IsOne( img * TransposedMat( List( img, row -> 
      List( row, ComplexConjugate ) ) ) ) );
  fi;
end;

IsRealRep := function( rep )
  if IsDecimalApproximation( Image( rep, One( Source( rep ) ) )[1][1] )
  then return true;
  else return
  IsGroupHomomorphism(rep) and
    IsFinite( Image( rep ) ) and
    IsMatrixGroup( Image( rep ) ) and
    Characteristic( FieldOfMatrixGroup( Image( rep ) ) ) = 0 and
    ForAll( GeneratorsOfGroup( Image( rep ) ), img ->
      ForAll( img, row -> ForAll( row, ent -> ent = ComplexConjugate( ent ) ) ) );
  fi;
end;

IsRealOrthogonalRep := rep -> IsRealRep( rep ) and IsUnitaryRep( rep );

##############################################################################
##
##  FortranStringOfMatrix - string form of a matrix for Fortran
##
##  Input - a matrix, probably square
##  Output - a string, first line the column length, second line is
##    all of the entries as flaoting point numbers in column major
##    order.
##
FortranStringOfMatrix := function( mat )
  return Concatenation( String( Size( mat[1] ) ), "\n",
    JoinStringsWithSeparator( List( Flat( TransposedMat( mat ) ),
      e -> DecimalApproximation( e, 25 ) ), " " ), "\n" );
end;

##############################################################################
##
##  ExampleUsage - create one file per rep per group element
##

ExampleUsage := function( grp )
  local chis, reps, elts, filename, file, repnr, eltnr;

  # Find the characters we can handle
  chis := RealCharacters(grp);;

  # Convert them to representations
  reps := List( chis, RealRepresentation );;

  # Convert them to unitary even?
  reps := List( reps, UnitaryRepresentation );

  # Which group elements do we want? Just generators?
  elts := GeneratorsOfGroup( grp );
  # or maybe all elements?
  elts := AsSet( grp );

  # Now write them to files:
  for repnr in [ 1 .. Size( reps ) ] do
    for eltnr in [ 1.. Size( elts ) ] do

      # name the file based on the rep and the elt
      filename := Concatenation( "rep", String(repnr), "elt", String(eltnr), ".txt" );

      # open a file for writing, no append
      file := OutputTextFile( filename, false );

      # turn off insanity
      SetPrintFormattingStatus( file, false );

      # print the matrix on lines 1 and 2, and the element on line 3
      PrintTo( file, FortranStringOfMatrix( Image( reps[repnr], elts[eltnr] ) ),
        "# ", String( elts[eltnr] ), "\n" );

      # close the file
      CloseStream( file );

    od;
  od;
end;