// 2008-03-19

intersectionMatrices:=function(G)
   H:=Stabilizer(G,1);
   orbs:=Orbits(H);
   d:=#orbs;
   reps:=[ Random(Set(o)) : o in orbs ];
   trsv:=[G|];
   for i in [1..d] do
      t,x:=IsConjugate(G,reps[i],1);
      trsv[i]:=x;
   end for;
   pair:=[ Random({ j : j in [1..d] | 1^trsv[i] in orbs[j] })
      : i in [1..d] ];
   return [ Matrix(Integers(),d,d,
      [ [ #(orbs[i] meet orbs[pair[j]]^(trsv[k]^(-1)))
         : k in [1..d] ] : j in [1..d] ] ) : i in [1..d] ];
end function;

eigenMatrixOverF:=function(G,F)
    Bs:=intersectionMatrices(G);
    Bs:=[ MatrixAlgebra(F,#Bs) | B : B in Bs ];
    egs:={ { Eigenspace(Transpose(Bs[i]),p[1]) 
        : p in Eigenvalues(Bs[i]) } : i in [1..#Bs] };
    md:=Maximum({ #d : d in egs });
    d1:=#Bs;
    dec:=Random({ d : d in egs | #d eq md });
    Exclude(~egs,dec);
    while #dec lt d1 do
       newd:=Random(egs);
       Exclude(~egs,newd);
       dec:={ v meet v1 : v in dec, v1 in newd 
          | Dimension(v meet v1) gt 0 };
    end while;
    bss:={ Basis(v)[1] : v in dec };
    bss0:=Random({ b : b in bss | &+Eltseq(b) ne 0 });
    bss1:=[ bss0 ] cat Setseq(Exclude(bss,bss0));
    return Matrix(F,#Bs,#Bs,[ 1/v[1]*v : v in bss1 ]);
end function;    

eigenMatrix:=function(G)
    Bs:=intersectionMatrices(G);
    egs:={ { Eigenspace(Transpose(Bs[i]),p[1]) 
        : p in Eigenvalues(Bs[i]) } : i in [1..#Bs] };
    md:=Maximum({ #d : d in egs });
    d1:=#Bs;
    dec:=Random({ d : d in egs | #d eq md });
    Exclude(~egs,dec);
    while #dec lt d1 do
       newd:=Random(egs);
       Exclude(~egs,newd);
       dec:={ v meet v1 : v in dec, v1 in newd 
          | Dimension(v meet v1) gt 0 };
    end while;
    bss:={ Basis(v)[1] : v in dec };
    bss0:=Random({ b : b in bss | &+Eltseq(b) ne 0 });
    bss1:=[ bss0 ] cat Setseq(Exclude(bss,bss0));
    return Matrix(Rationals(),#Bs,#Bs,[ 1/v[1]*v : v in bss1 ]);
end function;    


/*
eigenMatrix:=function(G)
    Bs:=intersectionMatrices(G);
    egs:=[ [ Eigenspace(Transpose(Bs[i]),p[1]) 
        : p in Eigenvalues(Bs[i]) ] : i in [1..#Bs] ];
    cegs:={ &meet{ r[i] : i in [1..#Bs] } 
        : r in CartesianProduct(egs) };
    bss:={ Basis(e)[1] : e in cegs | Dimension(e) ne 0 };
    bss0:=Random({ b : b in bss | &+Eltseq(b) ne 0 });
    bss1:=[ bss0 ] cat Setseq(Exclude(bss,bss0));
    return Matrix(Rationals(),#Bs,#Bs,[ 1/v[1]*v : v in bss1 ]);
end function;    
*/

kreinMatrices:=function(P)
   d1:=Nrows(P);
   n:=&+Eltseq(P[1]);
   Q:=n*P^(-1);
   return [ Matrix(Parent(P[1][1]),d1,d1,
       [ [ 1/(n*Q[1,k])*&+[ P[1,l]*Q[l,i]*Q[l,j]*Q[l,k] : l in [1..d1] ]
       : k in [1..d1] ] : j in [1..d1] ]
       ) : i in [1..d1] ];
end function;
         
abstractCentralizerAlgebra:=function(G)
   int:=intersectionMatrices(G);
   strc:=[ Eltseq(M) : M in int ];
   return Algebra< Rationals(),#int | strc >;
end function;
        
isMultiplicityFree:=function(G)
   return IsCommutative(abstractCentralizerAlgebra(G));
end function;

adjacencyMatrices:=function(G)
   n:=Degree(G);
   X2:={ [i,j] : i in {1..n}, j in {1..n} };
   orbitals:=Orbits(G,GSet(G,X2));
   adjM:=[];
   for i in [1..#orbitals] do
      A:=Matrix(Rationals(),n,n,[0:i in {1..n^2}]);
      for x in orbitals[i] do
	     A[x[1],x[2]]:=1;
      end for;
	  Append(~adjM,A);
   end for;
   return adjM;
end function;