//Determines constituents of the permutation character of a finite group G
//wrt a subgroup H

//load "D:\\magmadoc\\PermutationConstituents.mag";

KroDel:=function(i,j)
        if i eq j then
            return 1;
        else
            return 0;
        end if;
end function;

SmallerPermutationRepresentation:=function(grp,depth)
    local methods,natdeg,degs,subgr,orbgrp,orbpos,fthorb,min;
    methods:=["nat"];
    natdeg:=Degree(grp);
    degs:=[natdeg];

    orbgrp:=Orbits(grp);
    orbgrp:=[x:x in orbgrp];
    minorbreps:=[Minimum(o):o in orbgrp];

    fthorb:=[x:x in [1..#orbgrp]|#Core(grp,Stabilizer(grp,minorbreps[x])) eq 1];
    if fthorb ne [] then
        Append(~methods,"orb");
        min,orbpos:=Minimum([#orbgrp[x]:x in fthorb]);
        Append(~degs,min);
    end if;

    pigrp:=PrimeDivisors(#grp);
    fthpri:=[p:p in pigrp|#pCore(grp,p) eq 1];
    if fthpri ne [] then
        Append(~methods,"syl");
        min,sylpos:=Minimum([Index(grp,SylowSubgroup(grp,p)):p in fthpri]);
        Append(~degs,min);
    end if;

    min,j:=Minimum(degs);
    method:=methods[j];
    case method:
    when "nat":
        print "Try with maximal subgroups";
        subgr2:=grp;
        llact:=<>;
    when "orb":
        print "Starting orbit method";
        stab:=Stabilizer(grp,minorbreps[fthorb[orbpos]]);
        act2,subgr2:=CosetAction(grp,stab);
        llact:=<act2>;
    when "syl":
        print "Starting Sylow method";
        sylp:=SylowSubgroup(grp,fthpri[sylpos]);
        act2,subgr2:=CosetAction(grp,sylp);
        llact:=<act2>;
    end case;

    print "Starting maximal subgroups method";
    dpt:=1;
    max2:=MaximalSubgroups(subgr2);
    max2:=[x`subgroup:x in max2];
    lmax2:=[1..#max2];
    crf:=[x:x in lmax2|#Core(subgr2,max2[x]) eq 1];
    check:=-1;
    if crf ne [] then
        min,p:=Minimum([Index(subgr2,max2[x]):x in crf]);
        act,subgr:=CosetAction(subgr2,max2[crf[p]]);
        check:=1;
    else

        dpt:=1;
        repeat
            print "No Core-free " cat IntegerToString(dpt) cat "-maximal subgroups available";
            max3lis:={};
            for m in max2 do
                max3:=MaximalSubgroups(m);
                max3:=[x`subgroup:x in max3];
                for x in max3 do
                    Include(~max3lis,x);
                end for;
                lmax3:=[1..#max3];
                crf:=[x:x in lmax3|#Core(subgr2,max3[x]) eq 1];
                if crf ne [] then
                    min,p:=Minimum([Index(subgr2,max3[x]):x in crf]);
                    act,subgr:=CosetAction(subgr2,max3[crf[p]]);
                    check:=1;
                    break;
                end if;
            end for;
            dpt +:=1;
            max2:=max3lis;
        until check eq 1 or dpt gt depth;
    end if;


    if check eq 1 then
        Append(~llact,act);
    else
        subgr:=subgr2;
    end if;

    if #llact ne 0 then
        phi:=&*[x:x in llact];
        return phi,subgr;
    else
        print "no smaller permutation found";
        phi:=hom<grp->grp|[<grp.i,grp.i>:i in [1..Ngens(grp)]]>;
        return phi,subgr2;
    end if;
end function;


PermutationConstituents:=function(grp1,subg1:SMALL:=false)
    local K,qu,grp2,spr,grp,subg,ZZ,ccgrp1,nccg1,ccr1,ccgrp,nccg,ccr,ccgord,pm,cmgrp,ca,G1,base,Gf,fG,
    H1,orbs,minreps,dcr,Hf,fH,tr,dcd,ndcr,pp,F,R,V,b,Dlis,D,VD,bd,pos,XV,SE,KSE,dimz,ZD,bz,RZz,zet,regens,
    den,zz,esp,lowbound,p,E,PE,w,red,RE,redzz,egz,egzv,lagpol,z0,zz0,PCI,mult,Am,EGOM,goc,chilis,mat,mat1;

    ///////////////////////////////////////////////////////////////////////////////////////////
    //We implicitly assume that grp is given as a permutation group of small degree
    //and we try get a hold on subgroups whose monomial representation has many constituents.
    //So we get conjugacy classes in the group grp.
    ///////////////////////////////////////////////////////////////////////////////////////////

    if SMALL then
        K:=Core(grp1,subg1);
        qu:=CosetAction(grp1,K);
        grp2:=qu(grp1);
        spr:=SmallerPermutationRepresentation(grp2,3);
        grp:=spr(grp2);
        subg:=spr(qu(subg1));
        ccgrp1:=ConjugacyClasses(grp1);
        nccg1:=#ccgrp1;
        ccr1:=[x[3]:x in ccgrp1];
        cmgrp:=ClassMap(grp);
    else
        grp:=grp1;
        subg:=subg1;
    end if;

    ccgrp:=ConjugacyClasses(grp);
    nccg:=#ccgrp;
    ccr:=[x[3]:x in ccgrp];
    ccgord:=[x[1]:x in ccgrp];
    pm:=PowerMap(grp);


    ///////////////////////////////////////////////////////////////////////////////////////////
    //We build a coset action of grp on subg in order to obtain double coset representatives.
    //We use a strong base in order to avoid storage problems.
    ///////////////////////////////////////////////////////////////////////////////////////////

    ca:=CosetAction(grp,subg);
    G1:=ca(grp);

    base := Base( G1 );
    if base[1] ne 1 then
        ChangeBase( ~G1,[1] );
        base := Base( G1 );
    end if;
    Gf,fG:=FPGroupStrong(G1);

    H1:=ca(subg);
    orbset:=Orbits(H1);
    orbs:=[Sort([x:x in o]):o in orbset];
    minreps:=[Min(x):x in orbs];

    dcr:=[];
    for m in minreps do
        b,x:=IsConjugate(G1,1,m);
        Append(~dcr,x@@fG);
    end for;
    ndcr:=#dcr;

    dcd := [ 0 : i in [1..Degree(G1)] ];
    for dcn in [1..#orbs] do
        for p in orbs[dcn] do
            dcd[p] := dcn;
        end for;
    end for;

    ////////////////////////////////////////////////////////////////////////////////
    //We get a regular representation D for C=End_G(V) using the pp coefficients.
    //Use it to get Z(D) and its regular representation.
    //Finally we build a cyclic generator for Z.
    ////////////////////////////////////////////////////////////////////////////////

    invo:=Sym(#orbs)![dcd[1^fG(dcr[j]^-1)]:j in [1..#dcr]];
    pp:=function(i,j,k)
            return #(orbset[i] meet orbset[j^invo]^fG(dcr[k]));
    end function;

    F:=Rationals();
    R:=MatrixAlgebra(F,ndcr);
    GenIntersectionMatrix:=function(j)
        return R![pp(i,j,k):i,k in [1..ndcr]];
    end function;

    V:=RSpace(F,ndcr);
    b:=Basis(V);
    Dlis:=[R!1];
    D:=sub<R|Dlis>;
    VD:=sub<V|V.1>;
    poslis:=[2..ndcr];
    repeat
        bvd:=Basis(VD);
        poslis:=[i:i in poslis|not b[i] in VD];
        minpos:=Minimum(poslis);
        Dlis:=Append(Dlis,R!GenIntersectionMatrix(minpos));
        D:=sub<R|Dlis>;
        VD:=sub<V|[b*d:b in bvd, d in Dlis]>;
        dim:=Dimension(VD);
        dim;
    until dim eq ndcr;

    bd:=Basis(D);
    pos:=[Position(bd,D.i):i in [2..Ngens(D)]];
    XV:=KMatrixSpace(F,#pos*dim,dim);
    SE:=XV![bd[j,i,k]-bd[i,j,k]:i in [1..dim],j in pos,k in [1..dim]];
    KSE:=Kernel(Transpose(SE));
    dimz:=Dimension(KSE);
    ZD:=sub<D|[&+[(KSE.t)[i]*bd[i]:i in [1..dim]]:t in [1..dimz]]>;
    bz:=[ZD.i:i in [1..dimz]];
    RZz:=MatrixAlgebra(F,dimz);
    zet:=function(z)
        return RZz!&cat[Coordinates(ZD,b*z):b in bz];
    end function;
    regens:=[zet(b):b in bz];

    ///////////////////////////////////////////////////////////////////////////////
    //Here comes one hard part, namely finding the eigenvalues for z.
    //Even Dixon's modular approach does not seem to be feasible.
    //Maybe we need to combine with it factorization techniques.
    ///////////////////////////////////////////////////////////////////////////////


    esp:=Exponent(grp);
    lowbound:=Max(Min(Floor(2*Sqrt(#grp)),2*Index(grp,subg)),dimz);

    MinPrimeCongToOne:=function(d)
        p:=lowbound +(d+1-lowbound mod d);
        if not IsPrime(p) then
            repeat
                p +:=d;
            until IsPrime(p);
        end if;
        return p;
    end function;

    pri:=MinPrimeCongToOne(esp);
    E:=GF(pri);
    PE<x>:=PolynomialRing(E);
    w:=PrimitiveElement(E);
    ZZ:=Integers();
    red:=hom<ZZ->E|x:->x mod pri>;
    RE:=MatrixAlgebra(E,dimz);

    repeat
        coe:=[Random(-2,2):i in [1..dimz]];
        zz:=&+[coe[i]*regens[i]:i in [1..dimz]];
        den:=LCM([Denominator(zz[i,j]):i,j in [1..dimz]]);
        zz:=den*zz;
        redzz:=RE![red(zz[i,j]):i,j in [1..dimz]];
        pol:=MinimalPolynomial(redzz);
    until Degree(pol) eq dimz and Evaluate(pol,0) ne 0;

    egz:=Eigenvalues(redzz);
    egzv:=[e[1]:e in egz];
    Sort(~egzv);
    lagpol:=[Interpolation(egzv,[E!KroDel(i,j):i in [1..dimz]]):j in [1..dimz]];
    z0:=&+[coe[i]*ZD.i:i in [1..dimz]];
    z0:=den*z0;
    zz0:=MatrixAlgebra(E,ndcr)![red(ZZ!z0[i,j]):i,j in [1..ndcr]];
    PCI:=[Evaluate(l,zz0):l in lagpol];
    mult:=[];
    for p in PCI do
        b,a:=IsSquare(Rank(p));
        mult:=Append(mult,a);
    end for;
    Am:=[[p[1,i]:i in [1..ndcr]]:p in PCI];

    ///////////////////////////////////////////////////////////////////////////////
    //This seems to be the most time consuming routine, namely determining
    //the Gollan--Ostermann matrix
    ////////////////////////////////////////////////////////////////////////////////

    EGOM:=[[0:j in [1..nccg]]:i in [1..ndcr]];
    for p in [1..Degree(G1)] do
       inv := SVPermutation(G1,1,p);
       for s in [1..nccg] do
          img := p^(ca(ccr[s]^-1)*inv);
             EGOM[dcd[img],s] +:= 1;
       end for;
    end for;

    goc:=[[red(EGOM[j,v]):v in [1..nccg]]:j in [1..ndcr]];
    chilis:=[[&+[Am[i,j]*goc[j,v]:j in [1..ndcr]]/mult[i]:v in [1..nccg]]:i in [1..dimz]];

    E1:=CyclotomicField(esp:Sparse);
    mat:=<[E1!0:j in [1..nccg]]:i in [1..dimz]>;
    for j in [1..nccg] do
        ord:=ccgord[j];
        z2:=w^((pri-1) div ord);
        F1<q>:=CyclotomicField(ord);
        for i in [1..dimz] do
            val:=&+[ZZ!(&+[chilis[i,pm(j,n)]*z2^(-s*n):n in [0..ord-1]]/ord)*q^s:s in [0..ord-1]];
            val:=Minimise(val);
            mat[i,j]:=val;
        end for;
    end for;

    if SMALL then
    mat1:=<[E1!0:j in [1..nccg1]]:i in [1..dimz]>;
    for k in [1..nccg1] do
        j:=cmgrp(spr(qu(ccr1[k])));
        for i in [1..dimz] do
            mat1[i,k]:=mat[i,j];
        end for;
    end for;
    mat:=mat1;
    end if;
    return mult,mat;
end function;


PartCharTable2Gap:=procedure(grp,grpstr,lischar)
    ccgrp:=Classes(grp);
    ccgs:=[x[2]:x in ccgrp];
    ccgord:=[x[1]:x in ccgrp];
    pm:=PowerMap(grp);
    maxord:=Max(ccgord);
    nccg:=#ccgrp;
    pmout:=[[pm(i,n):i in [1..nccg]]:n in [1..maxord]];
    filename:="D:\\gapdoc\\CharTab" cat grpstr cat ".gap";
    names:=<"tbl.Size:=","tbl.UnderlyingCharacteristic:=","tbl.SizesConjugacyClasses:=","tbl.OrdersClassRepresentatives:=",
        "tbl.ComputedPowerMaps:=","tbl.Identifier:=","lischar:=">;
    values:=<#grp,0,ccgs,ccgord,pmout,grpstr,lischar>;
    strtmp:="Read(\"D:\\gapdoc\\CharTab" cat grpstr cat ".gap\");";
    PrintFile(filename,strtmp:Overwrite:=true);
    PrintFile(filename,"tbl:= rec();;");
    for x in [1..#names] do
        PrintFile(filename,names[x]);
        PrintFile(filename,values[x]);
        PrintFile(filename,";;");
    end for;
    PrintFile(filename,"ConvertToLibraryCharacterTableNC( tbl );;");
    PrintFile(filename,"lll:=LLL( tbl, lischars );;");
    PrintFile(filename,"irreducibles:=lll.vecs;;");
end procedure;