//Determines constituents of the monomial character lam^G of a finite group G
//wrt a subgroup H, lam linear character of H

//load "D:\\magmadoc\\MonomialConstituents.mag";
//load "D:\\magmadoc\\PermutationConstituents.mag";
//load "/mnt/win_d/magmadoc/MonomialConstituents.mag";

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

MonomialConstituents:=function(grp,subg,lam)
    local ccgrp,nccg,ccr,ccgord,pm,ca,G1,base,Gf,fG,
    H1,Hf,fH,K,orbs,minreps,dcr,dcd,ndcr,Hlis,centdcr,ncentdcr,lis,perm,tr,tot,
    R,V,bv,Dlis,D,VD,poslis,bvd,minpos,U, bd,pos,XV,SE,KSE,dimz,ZD,bz,RZz,zet,regens,
    den,coe,zz,pol,esp,lowbound,p,E,PE,w,z1,red,RE,redzz,egz,egzv,lagpol,z0,zz0,PCI,mult,
    b,a,Am,EGOM,ni,u,g,j,h,goc,chilis,ZZ,mat,ord,z2,F1,q,val,E1;


    ///////////////////////////////////////////////////////////////////////////////////////////
    //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.
    ///////////////////////////////////////////////////////////////////////////////////////////

    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);
    Hf,fH:=FPGroupStrong(H1);
    orbs:=Orbits(H1);
    orbs:=[Sort([x:x in o]):o in orbs];
    minreps:=[Min(x):x in orbs];

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

    K:=Kernel(lam);
    eta:=Exponent(subg/K);
    F<z>:=CyclotomicField(eta);

    Hlis:=[subg meet subg^(fG(d)@@ca):d in dcr];
    centdcr:=[d:d in [1..#dcr]|&and[(h,(fG(dcr[d])@@ca)^-1) in K:h in Generators(Hlis[d])]];
    ncentdcr:=#centdcr;
    if ncentdcr eq 1 then
        ccsubg:=Classes(subg);
        ccsubgr:=[x[3]:x in ccsubg];
        ncsubg:=#ccsubg;
        ordsubg:=#subg;
        ordgrp:=#grp;
        ccsubgsize:=[x[2]:x in ccsubg];
        ccsubgcent:=[ordsubg div x:x in ccsubgsize];
        ccgsize:=[x[2]:x in ccgrp];
        ccgcent:=[ordgrp div x:x in ccgsize];
        cmgrp:=ClassMap(grp);
        fus:=[cmgrp(x):x in ccsubgr];
        rhu:=[F!0:j in [1..nccg]];
        for i in [1..ncsubg] do
            rhu[fus[i]] +:=lam(ccsubgr[i])/ccsubgcent[i];
        end for;
        rho:=[rhu[j]*ccgcent[j]:j in [1..nccg]];
        return [1],[rho];
    else
        noncentdcr:=[x:x in [1..#dcr]|x notin centdcr];
        lis:=centdcr cat noncentdcr;
        perm:=Sym(#dcr)!lis;
        orbs:=[orbs[i^perm]:i in [1..#dcr]];
        minreps:=[minreps[i^perm]:i in [1..#dcr]];
        dcr:=[dcr[i^perm]:i in [1..#dcr]];

        tr:=[];
        for i in [1..#orbs] do
            o:=orbs[i];
            m:=minreps[i];
            tot:=[];
            for j in o do
                b,h:=IsConjugate(H1,m,j);
                tot:=Append(tot,h@@fH);
            end for;
            tr:=Append(tr,tot);
        end for;

        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).
        //We exploit the shape of the first row to speed up this calculation avoiding the calculation of
        //some generalized intersection matrices.
        //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]];
        R:=MatrixAlgebra(F,ncentdcr);
        GenIntersectionMatrix:=function(j)
            local mat,xj,mj,j1,xi,rightmul,b,h,lis;
            mat:=R![0:i,j in [1..ncentdcr]];
            xj:=fG(dcr[j]);
            mj:=minreps[j];
            mat[1,j]:=1;
            j1:=j^invo;
            b,h:=IsConjugate(H1,mj,1^(fG(dcr[j1])^-1));
            mat[j1,1]:=#orbs[j]*lam((fG(dcr[j])*h*fG(dcr[j1])*h^-1)@@ca);
            for i in [2..ncentdcr] do
                xi:=fG(dcr[i]);
                rightmul:=tr[i];
                for r in rightmul do
                    lis:=[k:k in [2..ncentdcr]|minreps[k] in (mj^H1)^(xi*fH(r))];
                    for k in lis do
                        b,h:=IsConjugate(H1,mj,minreps[k]^((xi*fH(r))^-1));
                        mat[i,k] +:=lam((((xj*h*xi)^fH(r))*(h*fG(dcr[k]))^-1)@@ca);
                    end for;
                end for;
            end for;
            return R!mat;
        end function;

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

        ///////////////////////////////////////////////
        //Improving the determination of the center.
        //We obtain as the kernel of the linear map
        //obtained annihilating only the first row of elements in D
        ////////////////////////////////////////////////

        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];

        /////////////////////////////////////////////////////////////////////////////////
        //The eigenvalue problem leads us to consider a modular a la Dixon approach.
        //We simplify considering suitable primes p so that pol factorizes completely
        //over GF(p).
        /////////////////////////////////////////////////////////////////////////////////

        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);
        z1:=w^((pri-1) div eta);
        red:=hom<F->E|z1>;
        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,ncentdcr)![red(z0[i,j]):i,j in [1..ncentdcr]];
        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..ncentdcr]]:p in PCI];

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

        EGOM:=[[F!0:v in [1..nccg]]:j in [1..ncentdcr]];
        EGOM[1,1]:=Degree(G1);
        for i in [1..ndcr] do
            ni:=#orbs[i];
            for s in [1..ni] do
                u:=fG(dcr[i])*fH(tr[i,s]);
                for v in [2..nccg] do
                    g:=ca(ccr[v]);
                    pos:=1^(u*g^-1*u^-1);
                    j:=dcd[pos];
                    if j le ncentdcr then
                        b,h:=IsConjugate(H1,minreps[j],pos);
                        EGOM[j,v] +:=lam((fG(dcr[j])*g^((h*u)^-1))@@ca);
                    end if;
                end for;
            end for;
        end for;


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

        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;
        return mult,mat;
    end if;
end function;

MonomialConstituentsViaPerm:=function(grp,subg,lam:SMALL:=false)
    K:=Kernel(lam);
    eta:=Exponent(subg/K);
    F<z>:=CyclotomicField(eta);
    mult,mat:=PermutationConstituents(grp,K:SMALL);

    ccgrp:=ConjugacyClasses(grp);
    nccg:=#ccgrp;
    ccgsize:=[x[2]:x in ccgrp];
    ordgrp:=#grp;
    ccgcent:=[ordgrp div x:x in ccgsize];
    cmgrp:=ClassMap(grp);

    ccsubg:=Classes(subg);
    ccsubgr:=[x[3]:x in ccsubg];
    ncsubg:=#ccsubg;
    ordsubg:=#subg;

    ccsubgsize:=[x[2]:x in ccsubg];
    ccsubgcent:=[ordsubg div x:x in ccsubgsize];

    fus:=[cmgrp(x):x in ccsubgr];
    rhu:=[F!0:j in [1..nccg]];
    for i in [1..ncsubg] do
        rhu[fus[i]] +:=lam(ccsubgr[i])/ccsubgcent[i];
    end for;
    rho:=[rhu[j]*ccgcent[j]:j in [1..nccg]];

    mut:=[m:m in mat|&+[ccgsize[j]*rho[j]*ComplexConjugate(m[j]):j in [1..nccg]] ne 0];
    return mult,mut;
end function;