declare verbose FindpCubed, 5;
SetVerbose("FindpCubed", 0);


/* Computes the idempotent corresponding to the subgroup H of G.*/
intrinsic SubgroupIdempotent(QG :: AlgGrp, H :: GrpPerm) -> AlgGrpElt
{}
    local i;

    G := Group(QG);
    Q := CoefficientRing(QG);
    Gelts := [g : g in G];
    v := [Q!0 : g in Gelts];
    for i:=1 to #Gelts do
        if Gelts[i] in H then
	    v[i] := 1/#H;
	end if;
    end for;
    return QG ! v;
end intrinsic;
	    

/******************************** Frobenius ***************************************
 *
 * OK ring of ints in K, h map provided by AutomorphismGroup or ArtinRepresentations
 * h:G -> Aut(K/Q), p is a rational prime
 *
 * Computes  a lift of Frobenius of \frp / p in the Galois group G, where \frp is a chosen prime ideal
 * in K above p.
 *
 *******************************************************************************************/
intrinsic Frobenius(OK :: RngOrd, h :: Map, p :: RngIntElt) -> GrpPermElt
{}
    local G, b, is_frob, vals, i;

    G := Domain(h);

    P := Factorization(p*OK)[1,1];
    cnt := 1;
    for g in G do
        // print "cnt = ", cnt;
        vals := [Valuation( h(g)(b) - b^p, P ) : b in Basis(OK)];
        is_frob := true;
        i := 2;
        while i le #vals and is_frob do
            if vals[i] lt 1 then
                is_frob := false;
            end if;
            i :=i+1;
        end while;
        if is_frob then
            return(g);
        end if;
        cnt := cnt+1;
    end for;
end intrinsic;

/*
intrinsic Frobenius(OK :: RngOrd, h :: Map, p :: RngIntElt, I :: Any) -> GrpPermElt
{}
    local G, b, is_frob, vals, i;

    G := Domain(h);
    Q, mQ := quo< G | I >;

    P := Factorization(p*OK)[1,1];
    cnt := 1;
    for q in Q do
        g := (mQ^-1)(q);
        print "cnt = ", cnt;
        vals := [Valuation( h(g)(b) - b^p, P ) : b in Basis(OK)];
        is_frob := true;
        i := 2;
        while i le #vals and is_frob do
            if vals[i] lt 1 then
                is_frob := false;
            end if;
            i :=i+1;
        end while;
        if is_frob then
            return(g);
        end if;
        cnt := cnt+1;
    end for;
end intrinsic;
*/

/******************************** DecompositionGroup ***************************************
 *
 * OK ring of ints in K, h map provided by AutomorphismGroup or ArtinRepresentations, 
 * p a rational prime
 *
 * Computes the decomposition group of a chosen prime \frp above p. 
 *
 *******************************************************************************************/
intrinsic DecompositionGroup(OK :: RngOrd, h :: Map, p :: RngIntElt) -> GrpPerm, Map
{}
    local G, P, D, g, vals, is_in_D, i;

    G := Domain(h);

    P := Factorization(p*OK)[1,1];
    D := [];
    cnt := 1;
    for g in G do
        // print "cnt = ", cnt;
        Q := ideal<OK | [h(g)(b) : b in Basis(P)]>;
        if Q eq P then
            Append(~D, g);
        end if;
        cnt := cnt+1;
    end for;
    return sub<G | D>;
end intrinsic;

/* Here U is a subgroup of a ray class group rcgp which lives in E/Q, h the usual map from G=Gal(E/Q) to the automorphisms and 
   m the map from rcgp to fractional ideals of E. We check if the subgroup U is invariant under the action of G.
*/
intrinsic IsInvariant(U :: GrpAb, h :: Map, m :: Map) -> BoolElt
{}

    local c, I, g, J, z;

    for c in Generators(U) do
        I := m(c);
        for g in Generators(Domain(h)) do
            J := ideal< Order(I) | [h(g)(z) : z in Generators(I)]>;
	    if not (m^-1)(J) in U then
	        return false;
	    end if;
        end for;
    end for;
    return true;
end intrinsic;


intrinsic IsTrivialAction(Delta :: GrpPerm, cl :: GrpAb, U :: GrpAb, h :: Map, m :: Map) -> BoolElt
{}
    local i;

    H, mq := quo<cl | U>;
    for c in Generators(H) do
        I := m((mq^-1)(c));
	for g in Delta do
	    J := ideal< Order(I) | [h(g)(z) : z in Generators(I)]>;
	    if not mq((m^-1)(J)) eq c then
	        return false;
	    end if;
	end for;
     end for;
     return true;
end intrinsic;



/* Compute the i-th ramification subgroup in lower numbering.*/
intrinsic RamificationGroup(P_F :: RngOrdIdl, i :: RngIntElt, hG :: Map) -> GrpPerm
{}
    local G, I, B, g, is_in_Gi, j, a;

    G := Domain(hG);
    Igens := [One(G)];
    I := sub<G | Igens >;
    B := Basis(Order(P_F));
    cnt := 1;
    for g in G do
        if not g in I then
            is_in_Gi := true;
	    j := 1;
	    while is_in_Gi and j le #B do
	        a := B[j];
	        if Valuation(hG(g)(a) - a, P_F) le i then
	            is_in_Gi := false;
	        end if;
	        j := j+1; 
	    end while;
	    if is_in_Gi then
                Append(~Igens, g);
                I := sub<G | Igens >;
	    end if;
            cnt := cnt+1;
        end if;
    end for;
    return I;
end intrinsic;



intrinsic ComputeDQ2(OE :: RngOrd, P :: RngOrdIdl, h :: Map) -> GrpAb, map
{}
    local R, mR, U, mU, gens, iota, u, Q2, DQ2, mDQ2, j;

    R, mR := quo< OE | P^2 >;
    U, mU := UnitGroup(R);
    gens :=[];
    iota := (mU*(mR^-1));
    for u in U do
        if Valuation(iota(u) - 1, P) ge 1 then
            Append(~gens, u);
        end if;
    end for;
    Q2 := sub<U | gens >;
    DQ2 := [(iota^-1)( h(gamma)(iota(q)) ) - q : q in Q2, gamma in Domain(h)];
    DQ2, mDQ2 := quo<Q2 | DQ2 >;
    j := (mDQ2^-1)*iota;
    return DQ2, j;
end intrinsic;


/* This function first computes 3 non-abelian extensions L/Q of degree 27 with full decomposition group
   at p=3 such that the completions at the unique prime above p=3 are a set of representatives for
   the non-abelian weakly ramified extensions K/Qp. The approach use the proof of the relevant theorem in Bley/Burns/Hahn.

   We return two lists, both of length 3. The first is the list of global extensions, the second the list of corresponding
   local extensions.

   The global extensions are constructed as extensions of E/Q where E is the composite of the unique subfield of Q(zeta_{p^2}) of degree p
   and the unique subfield of Q(zeta_q1) of degree p.

   We try to find the extesnions of E of degree 3 as subfields of the ray class fields of conductor q*P^2, where P denotes
   the unique prime above 3 in E. We consider primes q le B.
*/ 
intrinsic Find3CubedExtensions(q1 :: RngIntElt, B :: RngIntElt) -> List, List
{}

    p := 3;
    E1 := [S[1] : S in Subfields(CyclotomicField(p^2), p)][1];
    OE1 := MaximalOrder(E1);
    q := q1;
    E2 := [S[1] : S in Subfields(CyclotomicField(q), p)][1];
    OE2 := MaximalOrder(E2);
    // Factorization(p*OE2);
    // Factorization(q1*OE2);
    assert ClassNumber(OE1) eq 1;

    E := OptimizedRepresentation( MergeFields(E1, E2)[1] );
    OE := MaximalOrder(E);
    P := [P[1] : P in Factorization(p*OE)][1];
    assert ClassNumber(OE) eq 1;

    _, s := IsSubfield(E1, E);
    P_E1 := Factorization(p*OE1)[1,1];
    pi := s( Generators(P_E1)[2] );

    Gamma, Aut, h := AutomorphismGroup(E);
    /* Compute the quotient (U_E^1 / U_E^2)_Gamma. See the proof in bbh. */
    DQ2, j := ComputeDQ2(OE, P, h);
    u := j(DQ2.1);

    /* By the proof of the relevant theorem in bbh the order p subgroups of Q := E^\times / T (notation see bbh) correspond to the
       subgroups generated by pi*u^i, i=0,..,p-1.
    */

    qs := [];
    /* In the follwoing we generate a list of primes q such that the quotient of the ray class group of conductor q*P^2
       by <p*cl, I_Gamma*cl> possibly "contains an abelian extension" which is weakly ramified over E of degree 3.
    */ 
    for q := 5 to B do
        if q mod p ne 0 then
            // print "q = ", q;
            cl, m := RayClassGroup(q*P^2); 
            v := [p*c : c in Generators(cl)];
            for c in Generators(cl) do
                for g in Gamma do
                    I := m(c);
                    J := ideal< OE | [h(g)(z) : z in Generators(I)]>;
                    Append(~v, (m^-1)(J) - c);
                end for;
            end for;
            V := sub< cl | v >;
            H, mq := quo<cl | V >;
            // print "#H = ", #H;
            if #H ge 3 then
                m2 := Inverse(mq)*m;
                K := AbelianExtension(m2);
                if Valuation(Conductor(K), P) eq 2 then
                    // print "q found";
                    Append(~qs, q);
                end if;
                // print "********************************************************";
            end if;
        end if;
    end  for;
    
    vprintf FindpCubed, 1: "List of qs computed --> %o\n", qs;

    /* In a second step we compute a list Fs oo lists of the following form [* q, cl, m, U1, U2, ...*]
       where each Ui is a subgroup of cl of index p, invariant under Gamma and weakly wildly ramified.
    */
    Fs := [];
    for q in qs do
        cl, m := RayClassGroup(q*P^2); 
        list := [* q, cl, m  *];

        C, mC := PermutationGroup(cl);
        Vs := [s`subgroup : s in Subgroups(C : OrderEqual:= Integers()!(#cl / p))];
        Us := [sub<cl | [mC(a) : a in Generators(V)]> : V in Vs];
        // print "Check ", #Us, " subgroups.";
        k := 1;
        for U in Us do
            // print "Subgroup ", k;
            k := k+1;
            if IsInvariant(U, h, m) then
                // print "U is invariant under Gamma.";
                H, mq := quo<cl | U>;
                m2 := Inverse(mq)*m;
                K := AbelianExtension(m2);
                if (Valuation(Conductor(K), P) eq 2)  and IsInvariant(U, h, m) then
                    // print("U found");
                    Append(~list, U);
                end if;
            end if;
        end for;
        Append(~Fs, list);
        // print "**************************** qqqqqqqqqqqqqqqqq  ****************************************";
    end for;
    vprintf FindpCubed, 1: "List of Fs computed.\n";

    /* Now we go through the list Fs and try to find a represenative of the 3 local wildly and weakly ramified non
       abelian extensions of Qp.
    */
    found := [false, false, false];
    LocalpCubedFields := [* 1,2,3 *];
    GlobalpCubedFields := [*1,2,3 *];
    i := 1;
    while (not &and found) and i le #Fs do 
        cl := Fs[i,2];
        m := Fs[i,3];
        q := Fs[i, 1];
        vprintf FindpCubed, 1: "----------->               q = %o\n", q;
        for j:=4 to #Fs[i] do
            // print "(i, j) = ", i,j;
            U := Fs[i,j];
            H, mq := quo<cl | U>;
            m2 := Inverse(mq)*m;
            K := AbelianExtension(m2);
            Cond := Conductor(K);
            // print "val_P(cond) = ", Valuation(Cond, P);
            F := Factorization(q*P^2);
            for b:=0 to 2 do
                /* The following computes the Artin symbol of pi*u^b. By the proof in bbh the field K corresponds to pi*u^b
                   if and only if this Artin symbol is trivial.
                */
                xi := CRT([OE!pi, u^b], [q*OE, P^2]);
                I := xi*OE;
                // print "b = ", b, " -----> ", (m2^-1)(I);
                if (m2^-1)(I) eq Domain(m2)!0 and not found[b+1] then
                    vprintf FindpCubed, 1: "Found a field for b = %o\n", b;
                    found[b+1] := true;
		    O := EquationOrder(K);
		    Oabs := SimpleExtension(O);
		    OK := MaximalOrder(Oabs);
		    L0 := NumberField(OK);
		    OL0 := MaximalOrder(L0);
		    fac := Factorization(3*OL0);
		    P1 := fac[1,1];
		    G1, Aut, h1 := AutomorphismGroup(L0);
		    F0 := Completion(L0, P1 : Precision:=200);
		    GlobalpCubedFields[b+1] := L0;
		    LocalpCubedFields[b+1] := F0;
		end if;
            end for;
        end for;
        i := i+1;
     end while;
    return GlobalpCubedFields, LocalpCubedFields;

end intrinsic;

/* Computes a list of global fields with [L:Q] = l = 3 such that the localizations at at p=7 give the
   3 totally ramified extensions of Qp of degree l.
   We return the two lists, a list of the global fileds and secondly, a list of the corresponding completions.
*/ 
intrinsic FindRamifiedExtensionsOfDegree3(ms :: SeqEnum) -> List, List
{}

    found := [false, false, false];
    LocallFields := [* 1,2,3 *];
    GloballFields := [*1,2,3 *];

    l := 3; p := 7; u := 2;
    for m in ms do
        // print "m = ", m, "  found = ", found;
        F<zeta> := CyclotomicField(p*m);
        Q, mQ := CyclotomicAutomorphismGroup(F);
        C, mC := PermutationGroup(Q);
        Vs := [s`subgroup : s in Subgroups(C : OrderEqual:= Integers()!(#Q / l))];
        Us := [sub<Q | [mC(a) : a in Generators(V)]> : V in Vs];
        // print #Us, " subgroups found.";

        /* For the following see the comment in Find3CubedExtensions */ 
        for U in Us do
            H, mH := quo<Q | U>;
            for b := 0 to 2 do
                xi := CRT([p, u^b], [m, p]);
                tau := [g : g in Q | mQ(g)(zeta) eq zeta^xi][1];
                // print "b = ", b, " ------> ", mH(tau);
                if (mH(tau) eq H!0) and (not found[b+1]) then
                        L0 := FixedField(F, U);
                    OL0 := MaximalOrder(L0);
                    h := ClassNumber(L0);
                    // print "h = ", h, "Discriminant = ", Factorization(Discriminant(OL0));
                    // if h eq 1 and GCD(Discriminant(OL0), p) eq p then
                    if GCD(Discriminant(OL0), p) eq p then
                        found[b+1] := true;
                        L0 := FixedField(F, U);
                        OL0 := MaximalOrder(L0);
                        P7 := Factorization(p*OL0)[1,1];
                        F0 := Completion(L0, P7 : Precision:=200);
                        GloballFields[b+1] := L0;
                        LocallFields[b+1] := F0;
                    end if;
                end if;
            end for;
        end for;
        // print "############################################";
    end for;
    return GloballFields, LocallFields;
end intrinsic;

intrinsic FindRamifiedExtensionsOfDegree9(ms :: SeqEnum) -> List, List
{}

    LocallFields := [* 1,2 *];
    GloballFields := [*1,2 *];

    l := 3; p := 7; 
    for m in ms do
        F<zeta> := CyclotomicField(p*m);
        Q, mQ := CyclotomicAutomorphismGroup(F);
        C, mC := PermutationGroup(Q);
        Vs := [s`subgroup : s in Subgroups(C : OrderEqual:= Integers()!(#Q / l^2))];
        Us := [sub<Q | [mC(a) : a in Generators(V)]> : V in Vs];
        print #Us, " subgroups found.";

        b := 1;
        for U in Us do
            H, mH := quo<Q | U>;
            if not Exponent(H) eq 3 then
                L0 := FixedField(F, U);
                OL0 := MaximalOrder(L0);
                P7 := Factorization(p*OL0)[1,1];
                if RamificationIndex(P7) eq l and InertiaDegree(P7) eq 3 then
                        print "Example found m = ", m;
                        F0 := Completion(L0, P7 : Precision:=200);
                        GloballFields[b] := L0;
                        LocallFields[b] := F0;
                        b := b+1;
                end if;
            end if;
        end for;
    end for;
    return GloballFields, LocallFields;
end intrinsic;

intrinsic WrongFindExtensionsOfDeg63(q :: RngIntElt) -> List
{}

    GloballFields, LocallFields := FindRamifiedExtensionsOfDegree3([1,13]);

    ZX<x> := PolynomialRing(Integers());
    F1 := NumberField(x^3 - x^2 - 4*x - 1);
    OF1 := MaximalOrder(F1);

    G63_extensions := [];

    for i in [1..#GloballFields] do
        F2 := GloballFields[i];
        print "---------------> Compute data for the ", i, "th field of degree 63";
        OF2 := MaximalOrder(F2);
        E := OptimizedRepresentation( MergeFields(F1, F2)[1] );
        OE := MaximalOrder(E);
        print "Class number of E: ", ClassNumber(OE);
        Ps := [P[1] : P in Factorization(7*OE)];
        Q, Aut, h := AutomorphismGroup(E);
        /* Wir suchen nun eine Strahlklassnerweiterung von E, die eine Teilerweiterung 
           vom Grad 7 enthaelt und in der
           das Primideal ueber der 7 voll und schwach erzweigt.
        */
    
        q := 2;
        cl, m := RayClassGroup(q*Ps[1]^2); 
        U := sub<cl | [7*u : u in Generators(cl)]>;
        H, mq := quo<cl | U>;
        m2 := Inverse(mq)*m;
        K := AbelianExtension(m2);
        _ := AutomorphismGroup(K:Over := [h(Q.1)]);
        _ := AutomorphismGroup(K:Over := [h(Q.2)]);
        A, lA, mA := AutomorphismGroup(K:Over := [h(Q.1), h(Q.2)]);
        time O := EquationOrder(K);
        time Oabs := SimpleExtension(O);
        time OK := MaximalOrder(Oabs);
        time Ka := NumberField(OK);
        for i in lA do InternalAutomorphismAdd( Ka!i(Ka.1)); end for;
        time Z, _, mZ := AutomorphismGroup(Ka);
	Append(~G63_extensions, [* Ka, OK, Z, mZ *]);
    end for;
    return G63_extensions;
end intrinsic;


intrinsic FindExtensionsOfDeg63(low, high :: RngIntElt) -> List
{}

    GloballFields, LocallFields := FindRamifiedExtensionsOfDegree3([1,13]);
    p := 7;

    ZX<x> := PolynomialRing(Integers());
    F1 := NumberField(x^3 - x^2 - 4*x - 1);
    OF1 := MaximalOrder(F1);

    G63_extensions := [];
    F2 := GloballFields[1];
    E := OptimizedRepresentation( MergeFields(F1, F2)[1] );
    OE := MaximalOrder(E);

    P := Factorization(7*OE)[1,1];
    Gamma, Aut, h := AutomorphismGroup(E);

    /* We compute the subgroups Delta of Gamma = Gal(E/Q) corresponding to the extensions which E^Delta which ramify at p=7 */
    S := [s`subgroup : s in Subgroups(Gamma) | s`order eq 3];
    Deltas := [Delta : Delta in S | Factorization(p*MaximalOrder(FixedField(E, Delta)))[1,2] eq 3];

    qs := [q : q in [low..high] | GCD(7,q) eq 1];
    j := 1;
    found := {};

    while j le #qs and #found ne 3 do
        q := qs[j];
        cl, m := RayClassGroup(q*P^2);
        if #cl mod p eq 0 then
            C, mC := PermutationGroup(cl);
            Vs := [s`subgroup : s in Subgroups(C : OrderEqual:= Integers()!(#cl / p))];
            Us := [sub<cl | [mC(a) : a in Generators(V)]> : V in Vs];
            for U in Us do
                // U := sub<cl | [7*u : u in Generators(cl)]>;
                H, mq := quo<cl | U>;
                if IsInvariant(U, h, m) and #H eq 7 then
                    ind := Index([IsTrivialAction(Delta, cl, U, h, m) : Delta in Deltas], true);
                    if not ind in found then
                        print "q = ", q, "ind = ", ind;
                        found := found join {ind};

                m2 := Inverse(mq)*m;
                K := AbelianExtension(m2);
                _ := AutomorphismGroup(K:Over := [h(Gamma.1)]);
                _ := AutomorphismGroup(K:Over := [h(Gamma.2)]);
                A, lA, mA := AutomorphismGroup(K:Over := [h(Gamma.1), h(Gamma.2)]);
                time O := EquationOrder(K);
                time Oabs := SimpleExtension(O);
                time OK := MaximalOrder(Oabs);
                time Ka := NumberField(OK);
                for i in lA do InternalAutomorphismAdd( Ka!i(Ka.1)); end for;
                time Z, _, mZ := AutomorphismGroup(Ka);
                Append(~G63_extensions, [* Ka, OK, Z, mZ *]); 

                    end if;
                end if;
            end for;
        end if;
        j := j+1;
    end while;
    return G63_extensions;
end intrinsic;


intrinsic FindExtensionsOfDeg63(qs:: SeqEnum) -> List
{}

    GloballFields, LocallFields := FindRamifiedExtensionsOfDegree3([1,13]);
    p := 7;

    ZX<x> := PolynomialRing(Integers());
    F1 := NumberField(x^3 - x^2 - 4*x - 1);
    OF1 := MaximalOrder(F1);

    G63_extensions := [];
    F2 := GloballFields[1];
    E := OptimizedRepresentation( MergeFields(F1, F2)[1] );
    OE := MaximalOrder(E);

    P := Factorization(7*OE)[1,1];
    Gamma, Aut, h := AutomorphismGroup(E);

    /* We compute the subgroups Delta of Gamma = Gal(E/Q) corresponding to the extensions which E^Delta which ramify at p=7 */
    S := [s`subgroup : s in Subgroups(Gamma) | s`order eq 3];
    Deltas := [Delta : Delta in S | Factorization(p*MaximalOrder(FixedField(E, Delta)))[1,2] eq 3];

    j := 1;
    found := {};

    while j le #qs and #found ne 3 do
        q := qs[j];
        cl, m := RayClassGroup(q*P^2);
        U := sub<cl | [7*u : u in Generators(cl)]>;
        H, mq := quo<cl | U>;
        if IsInvariant(U, h, m) and #H eq 7 then
            ind := Index([IsTrivialAction(Delta, cl, U, h, m) : Delta in Deltas], true);
            if not ind in found then
                print "q = ", q, "ind = ", ind;
                found := found join {ind};
                m2 := Inverse(mq)*m;
                K := AbelianExtension(m2);
                _ := AutomorphismGroup(K:Over := [h(Gamma.1)]);
                _ := AutomorphismGroup(K:Over := [h(Gamma.2)]);
                A, lA, mA := AutomorphismGroup(K:Over := [h(Gamma.1), h(Gamma.2)]);
                time O := EquationOrder(K);
                time Oabs := SimpleExtension(O);
                time OK := MaximalOrder(Oabs);
                time Ka := NumberField(OK);
                for i in lA do InternalAutomorphismAdd( Ka!i(Ka.1)); end for;
                time Z, _, mZ := AutomorphismGroup(Ka);
                Append(~G63_extensions, [* Ka, OK, Z, mZ *]); 
            end if;
        end if;
        j := j+1;
    end while;
    return G63_extensions;
end intrinsic;

