/*
 * bbhconj.m
 * (c) 2017 Werner Bley
 * 
 * This implementation computes the invariant \fra_{E/F} for the degree p^3 extension of \Qp,
 * where in principle p could be any odd prime (but it only works for p =3), 
 * and for extensions of Qp of degree 63 with Gal(L^C/Qp) \simeq C_3 x C_3. 
 *
 * For p^3 extensions we also compute the right hand side of the local epsilon constant conjecture (as 
 * formulated in Bley/Burns/Hahn). In this way we get a proof of the local epsilon conjecture.
 *
 * Finally we test the equality \fra_{L/\Qp} = \frc_{L/\Qp}. It is an open question whether this equality holds in all cases.
 * The implementation works for all of the p^3 extensions and the first of the degree 63 extensions. In the remaining two cases
 * of degree 63 extensions the computations take too long, the program does not terminate.
 * 
 * The implementation is to a very large extend an adaption of the implementation of Ruben Debeerst who used it to check the local
 * epsilon constant conjecture. Most if not all of his comments are kept. Intrinsics or functions which are special to the
 * p^3 situation usually start with the prefix pCubed.
 * Many of the intrinsics started with the prefix pCubed can also be used for degree 63 extensions. This is bad style and should
 * be changed in case that someone should try to improve the programs in order to deal with the remaining degree 63 extensions.
 * In my opinion it should be possible to verify \fra_{L/\Qp} = \frc_{L/\Qp} for all degree 63 extensions, however, this will require
 * some additional work.
 *
*/ 


/*
 * epsconj.m
 * (c) 2008-2009 Ruben Debeerst
 *
 * Algorithms to verify the local epsilon constant conjecture.
 * This includes methods for relative algebraic K-groups, for
 * Gaussian sums, for number fields and others.
 *
 * Depends On:
 *   artin.m
 *   characters.m
 *   global.m
 *   brauer.m
 *   ao.m
 *   RelAlgKTheory.m
 *
 * Methods provided:
 *   LECverify
 *   LECcreateRec
 *   LECinitializeCompletion
 *   LECpreparations
 *   LECprepareEps
 *   LEClattice
 *   LECcomputeValues
 *   LECdiscriminant
 *   LECcohomologicalTerm
 *   LECunramifiedTerm
 *   LECepsilonConstant
 *   LECcorrectionTerm
 *   LECimagesK0Rel
 *   LECcheck
 *
 */

/**
 * The functions for the Local Epsilon Constant Conjecture all
 * start with the prefix LEC.
 * 
 * The main function is LECverify. It depends on a global field
 * which is not split at a given prime. If a local field is given,
 * a global representation is computed first.
 *
 * The verification of the conjecture works on a special record-format
 * which holds all neccecary information. To experiment with the values
 * in the conjecture, one can simply give the argument by reference,
 * e.g.
 * 
 * setting := LECcreateRec(L, p);
 * LECverify(~setting);
 * 
 * LECverify will call the following functions:
 *   LECprepareEps    (do preparations for Epsilon Constants),
 *   LECpreparations  (do all other preparations),
 *   LECcomputeValues (compute all values for the conjecture),
 *   LECimagesK0Rel   (read them in the same relative K-group),
 *   LECcheck         (check the conjecture).
 *
 * Some parts are further split, i.e. one can compute each part of
 * the conjecture seperately (LECdiscriminant, LECcorrectionTerm,
 * LECunramifiedTerm, LECcohomologicalTerm, LECepsilonConstant) by
 * either passing a setting-record or all nececary parameters.
 *
 * The values in the record-format that already exist will be used by
 * LECverify, as far as it makes sense. The algorithm will then omit
 * the computation of those values. This allows to reuse values which
 * are already computed.
 *
 * For example ... 
 *
 */

import "RelAlgKTheory.m": K0RelElt;

// local functions
// for epsilon constant conjectures
forward pCubed_lec_verify, lec_verify, star, extension_for_unramified_term,
        unramified_cyclotomic_subfield;
// for cohomological term
forward lec_lattice, lattice_generator_theta, splitting_module_LPmul_modX,
        compute_LPmul_modX, proj_resolution_ZGr, zG_generators,
        tensor_Zmod_with_Q, isomorphism_matrix_KQF, make_QG_linear,
        truncated_exp, truncated_exp_precision;
forward compute_LPmul_modX_vals, truncated_exp_precision_vals;
// for epsilon constants
forward pCubed_lec_unram_char, pCubed_lec_epsilon_constant, G63_lec_epsilon_constant, lec_epsilon_constant, comp_precision, nontriv_ramGrps, ramGrp,
        G63_galois_gauss_sum, pCubed_galois_gauss_sum, galois_gauss_sum, representants, standard_add_char,
        root_of_unity;

// for relative K-groups
forward numden, lin_ext_char;
// for number fields
forward complex_compatible_embedding, embedding, inverse_embedding,
        root_with_complex_embedding, ramification_group,
        inertia_group, fixed_field, frobenius_grp_elt,
        is_unramified_cyclotomic, restrict_aut_to_quotient,
        restricted_aut;
// for local fields
forward completion_with_prec, localized_automorphism_group, local_map,
        basis_pad, elem_to_seq;
// various
forward test_G_compatible, test_G_action, test_cocycle ;

// verbose comments
declare verbose epsConj, 5;
SetVerbose("epsConj", 1);



/**************************************************
 *         Epsilon Constant Conjecture
 **************************************************/


declare attributes ModTupRng : moduli;

/* pCubedLECverify
 * Computes \fra, and also the right hand side of the local epsilon constant conjecture.
 * Finally compares \fra to \frc, the local twisted unramified characteristic. Conjecturally we always
 * have the equaility \fra = \frc.
 *
 * The main work is done in the function  pCubed_lec_verify
 *
*/
intrinsic pCubedLECverify(L::FldNum, p::RngIntElt) -> BoolElt
{ Compute \fra for the extension L/Qp of degree 27 and exponent 3.
  In addition, compute the right hand side of the local epsilon constant conjecture as formulated in Bley/Burns/Hahn.
  Finally compare \fra with the unramified characteristic \frc and check the conjecture \fra = \frc.} 

    local bool;

    require #DecompositionType(L,p) eq 1 :
            "The prime p cannot split in L.";

    t := Cputime();
    bool := pCubed_lec_verify(LECcreateRec(L,p));
    vprintf epsConj, 1: "\n== Total Time: %o\n", Cputime(t);
    return bool;
end intrinsic;


intrinsic pCubedLECverify(~setting::Rec) 
{ Compute \fra for the extension L/Qp of degree 27 and exponent 3.
  In addition, compute the right hand side of the local epsilon constant conjecture as formulated in Bley/Burns/Hahn.
  Finally compare \fra with the unramified characteristic \frc and check the conjecture \fra = \frc.}


    local bool, t;
    t := Cputime();
    bool, setting := pCubed_lec_verify(setting);
    vprintf epsConj, 1: "\n== Total Time: %o\n", Cputime(t);
    print bool;
end intrinsic;


/* pCubed_lec_verify(setting :: Rec)
 * Computes \fra, the right hand side of the local epsilon constant conjecture and \frc.
 * It thus verifies the local epsilon constant conjeture and the conjecture that \fra = \frc.
 * This works well for all of the three extensions of Qp of degree 27.
 *
 * The function can also be used to compute \fra  and \frc for the first of the degree 63 extensions.
 * The remaining two extensions are too big (the defining polynomials have very large coefficients), so
 * that the computations does not terminate. In principle, I would hope that with some additional effort
 * these computations could be done.
*/
function pCubed_lec_verify(settingInput)
    local setting, bool, t;

    setting := settingInput;
    t := Cputime();

    vprint epsConj, 1: "== Compute fra_{L/Q}, the rhs of the epsilon constant conjecture and frc_{L/Q} ==";
    IndentPush();
    vprint epsConj, 1: "L=", setting`L;
    vprint epsConj, 1: "p=", setting`p;
    if assigned setting`N then
        vprint epsConj, 1: "N=", setting`N;
    end if;
    IndentPop();

    L := setting`L; p := setting`p; G := setting`G; h := setting`psiL;
    OL := MaximalOrder(setting`L);
    setting`OL := OL;
    P := Factorisation(p*OL)[1,1];
    if Degree(L) eq 63 then
        AL := P^(-13);
    else
        AL := p^-1*P; /* AL is the square root of the inverse different. */
    end if;

    /* Use pCubed_nb  or G63_nb to compute Ath and theta such that Ath*theta = AL. We do not use the weak normal basis element wnb
       computed in pCubed_nb in the subsequent computations.
    */
    if Degree(L) eq 63 then
        X :=G63_nb(L, G, h: WeakINB:=true);
    else
        X :=pCubed_nb(L, G, AL, h: WeakINB:=true);
    end if;
    theta := X[2]; Ath := X[3]; AssOrd := X[4]; wnb := X[5]; QG := X[7];
    setting`Ath := Ath;
    setting`theta := theta;
    setting`prec := 36;  /* I have forgotten where this comes from and where it is needed. */

    /* Compute setting`brauerInd which is needed for the computation of the Gauss sums. */
    if Degree(L) eq 63 then    
        G63LECprepareEps(~setting);
    else
        pCubedLECprepareEps(~setting);
    end if;

    if not assigned setting`Krel then
        setting`Krel := RelativeGroupC(setting`G, setting`p);
    end if;

    if not assigned setting`dLK then
         LECdiscriminant(~setting);
    end if;
    if not assigned setting`uLK then
        pCubedLECunramifiedTerm(~setting);
    end if;

    /* Compute the Gauss sums \tau(psi_2(chi) - chi) for all chi in IrrG.
       Here psi_2 denotes the second Adams operation. */
    if not assigned setting`tLK then
        inclCE := IdentityFieldMorphism(setting`Qmpt);
        if Degree(L) eq 63 then
            gauss_sums := G63_lec_epsilon_constant(setting`psiL, setting`brauerInd, setting`linChrs, setting`p, setting`t, inclCE, setting`IrrG, 1);
        else        
            gauss_sums := pCubed_lec_epsilon_constant(setting`psiL, setting`brauerInd, setting`linChrs, setting`p, setting`t, inclCE, setting`IrrG, 1);
        end if;
        setting`tLK := gauss_sums;
    end if;

    /* Compute the twisted unramified characteristic \frc_{L/Qp}.
       It is defined by the reduced norm of (1-e_G0) + Frob^-1 * e_G0 . 
    */
    unram_char := pCubed_lec_unram_char(setting);
    setting`unram_char := unram_char;

    setting`LN := Compositum(setting`L, setting`N);

    Qmpt := setting`Qmpt;
    N := setting`N;
    LN := setting`LN;  
    LNG := GroupAlgebra(LN, setting`G);
    GG := [g : g in setting`G];

    /* The following lines replace the use of complex_compatible_embedding which takes a lot of time and often fails.
       Instead we just use coercion for L sseq LN and choose the correct embedding sigma_{i0} such that sigma_{i0}(LN!L.1) = L.1
       (L.1 viewed as an element of C via Conjugates(L.1)[1]).
       We then adapt Coercion(N,LN) by an element of Gal(N/Q).
    */

    precision := 100;
    conj := Conjugates(LN ! L.1);
    i0 := [i : i in [1..#conj] | Abs( conj[i] - Conjugates(L.1)[1]  ) lt 10^(-precision+1) ][1];
    G_N,_,h_N := AutomorphismGroup(N);
    g0 := [g : g in G_N | Abs( Conjugates(LN!h_N(g)(N.1))[i0] - Conjugates(N.1)[1] ) lt 10^(-precision+1)][1];
    incl_L_LN := Coercion(L, LN);
    incl_N_LN := h_N(g0)*Coercion(N,LN);


    dLK_LN := &+([ x[i]*y[i] : i in [1..#x] ]) where y is GG
        where x is [LNG!incl_L_LN(c) : c in Coefficients(setting`dLK)];
    uLK_LN_inv := &+([ x[i]*y[i] : i in [1..#x] ]) where y is GG
        where x is [LNG!incl_N_LN(c) : c in Coefficients(setting`uLK^(-1))];

    /* y is the quotient of discriminant and unramified term. It is known that it lives in the same field as the Gauss sums
       which is small, in our case of degree 18 over the rationals.
       Note: It is not true that y has to live in Qmpt because the reduced norm map is an iso from K_1 to the centre of E[G],
       hence one only knows that y has a representative (maybe even some matrix) mod commutators which is in Qmpt[G].
    */
    y := dLK_LN * uLK_LN_inv;

    /* Compute the subfield M of LN which is generated by the coefficients of y. 
       In our case M is of degree 9.
    */
    coeff := Coefficients(y);
    time M, iota_M := sub<LN | coeff[1]>;
    time coeff_M := [(iota_M^-1)(c) : c in coeff];

    vprintf epsConj, 1: "Compute E.\n";
    fac := Factorization(PolynomialRing(M) ! DefiningPolynomial(Qmpt));
    E := NumberField(fac[1,1]);
    iota_M_E := Coercion(M,E);
    c := Conductor(Qmpt);
    i := [i : i in [1..c] | GCD(i,c) eq 1 and Abs( Conjugates(E.1^i)[1] - Conjugates(Qmpt.1)[1] ) lt 10^-20][1];
    zeta := E.1^i;
    iota_Qmpt_E := hom<Qmpt->E | zeta>;
    iota_Qmpt_E_inv := hom<E->Qmpt | t:->RelativeEmbeddingInverse(t, M, E, Qmpt, zeta)>;

    /* Compute Nrd of y = dLK/uLK as an element of Z(EG)^\times. */
    EG := GroupAlgebra(E, setting`G);
    y_E := &+[(EG ! iota_M_E(coeff_M[i])) * GG[i] : i in [1..#coeff_M]];
    Nrd_y_E := NewtonReducedNormEG(y_E, iota_Qmpt_E);

    /* Compute the final element which must be Q[G]. */
    tLK := setting`tLK;
    tLK_E := [iota_Qmpt_E(z) : z in tLK];
    /* Note that the y-term in the definition of \fra is the inverse of the unramified characteristic. */
    unram_char_E := [iota_Qmpt_E(z) : z in unram_char];
    yLK_E := [iota_Qmpt_E(z^-1) : z in unram_char];
    res := [Nrd_y_E[i] / (tLK_E[i]*yLK_E[i]^-1) : i in [1..#Nrd_y_E]];

    /* It finally remains to compute Corr = [Ath, id, ZG], the correction term which comes from our choice of theta. */
    QG := setting`Krel`QG;
    A := ModuleScalarProd(QG, Ath, QG!(1));
    A := ZGModuleInit(A, Ath`phi);
    BasisA := LocalBasis(A, p);
    ZG := ZGModuleInit(AssOrd, Ath`phi);
    BasisZG := LocalBasis(ZG, p);
    phi := IdentityMatrix(Rationals(), #G);
    S :=  QGMatrix(QG, phi, A, BasisA, ZG, BasisZG);
    lambda := S[1,1];
    Nrd_lambda := NewtonReducedNormEG(lambda, iota_Qmpt_E);

    /* Finally compute \fra as an element of the relative group. */
    Krel := setting`Krel;
    IrrG := setting`IrrG;
    embeddings := [* complex_compatible_embedding(Krel`QG`H[i,1], Qmpt, 100, false) : i in [1..#Krel`QG`H] *];
    unit := [Nrd_lambda[i] * res[i] : i in [1..#res]];
    nr := [Index(IrrG,h[2]): h in Krel`QG`H];
    nr := [ unit[i] :  i in nr ];
    nr := [* (embeddings[i]^(-1))((iota_Qmpt_E_inv)(nr[i])) : i in [1..#nr] *];
    omega := mapToRelKGroup(Krel, nr);
    setting`lhs := omega;
    vprintf epsConj, 1: "lhs computed --------- > %o\n", setting`lhs;


    QG := setting`Krel`QG;
    P := Factorization(p*OL)[1,1];
    I := RamificationGroup(P, 0, h);
    eI := SubgroupIdempotent(QG, I);
    eG := SubgroupIdempotent(QG, G);
    Frob := QG ! Frobenius(OL, h, p) ;
    eins := One(QG);
    eta := (eins - eI) + Frob^-1 * eI;
    incl_id := Coercion(setting`Qmpt, setting`Qmpt);
    Nrd_eta := NewtonReducedNormEG(eta, incl_id);
    nr := [Index(IrrG,h[2]): h in Krel`QG`H];
    nr := [ Nrd_eta[i] :  i in nr ];
    nr := [* (embeddings[i]^(-1))(nr[i]) : i in [1..#nr] *];
    setting`twisted_unram_char := mapToRelKGroup(Krel, nr);
    vprintf epsConj, 1: "twisted_unram_char --------- > %o\n", setting`twisted_unram_char;

    if Degree(L) eq 63 then
        return true, setting;
    end if;


    /* We now start the computation of the right hand side of the local epsilon constant conjecture as formulated in
       Bley/Burns/Hahn.
    */

    setting`N := extension_for_unramified_term(setting`L,setting`P, setting`G);

    /* Computation of the correction term. */
    LECcorrectionTerm(~setting);
    incl_id := Coercion(setting`Qmpt, setting`Qmpt);
    Nrd_mLK := NewtonReducedNormEG(setting`mLK, incl_id);

    vprintf epsConj, 1:  "Nrd_mLK computed.\n";

    /* Computation of the Jacobi sums. */
    if not assigned setting`jLK then
        vprintf epsConj, 1: "Compute the Jacobi sum jLK.\n";
        inclCE := IdentityFieldMorphism(setting`Qmpt);
        gauss_sums := pCubed_lec_epsilon_constant(setting`psiL, setting`brauerInd, setting`linChrs, setting`p, setting`t, inclCE, setting`IrrG, 2);
        setting`jLK := gauss_sums;
    end if;

    vprintf epsConj, 1: "jLK computed.\n";

    // Completion with enough precision
    
    /* theta is replaced by p^s*theta with s minimal such that val(p^s*theta) ge p+2 */
    s := 0;
    while Valuation(setting`theta, setting`P) lt setting`p+2 do 
        setting`theta := setting`theta*p; 
        s := s+1; 
    end while;

    /* For the computation of the cohomological term we need the completion of L at the prime P. */
    if not assigned setting`LP or
       not assigned setting`P /* or
       not test_G_compatible(setting`iota, setting`psiL, setting`psiLP, true, 0, setting`prec+2) or
       not test_G_compatible((setting`iota)^(-1), setting`psiLP, setting`psiL, false, setting`P, setting`prec+2) */
       then
       //not Minimum(test_G_compatible(setting`iota, setting`psiL, setting`psiLP, true, 0)
       //    join test_G_compatible((setting`iota)^(-1), setting`psiLP, setting`psiL, false, setting`P))
       //    ge setting`prec+2 then
        vprintf epsConj, 1: "Compute completion (prec=%o). Takes about 10 and a half minutes.", setting`prec+2;
        vtime   epsConj, 1: LECinitializeCompletion(~setting, setting`prec+2);
    end if;
    IndentPop();
    IndentPush();

    vprintf epsConj, 1: "LP computed.\n";

    /* Computation of the cohomological term. */
    if not assigned setting`eLK then
        vprintf epsConj, 1:  "Compute the cohomological term. Takes about 11 minutes\n";
        vtime epsConj, 1: LECcohomologicalTerm(~setting);
    end if;

    vprintf epsConj, 1: "eLK computed.\n";

    /* It finally remains to compute Corr = [p^-s*Ath, id, ZG], the correction term which comes from our choice of theta. */
    A := ModuleScalarProd(QG, Ath, QG!(p^-s));   
    A := ZGModuleInit(A, Ath`phi);
    BasisA := LocalBasis(A, p);
    ZG := ZGModuleInit(AssOrd, Ath`phi);
    BasisZG := LocalBasis(ZG, p);
    phi := IdentityMatrix(Rationals(), #G);
    S :=  QGMatrix(QG, phi, A, BasisA, ZG, BasisZG);
    lambda := S[1,1];
    incl_id := Coercion(setting`Qmpt, setting`Qmpt);
    Nrd_lambda := NewtonReducedNormEG(lambda, incl_id);

    /* So testet Ruben. */
    Krel := setting`Krel;
    IrrG := setting`IrrG;
    embeddings := [* complex_compatible_embedding(Krel`QG`H[i,1], Qmpt, 100, false) : i in [1..#Krel`QG`H] *];
    /* Note that with the new \frc we must have the inverse of \frc in the rhs of the conjecture. */
    unit := [Nrd_lambda[i] * setting`jLK[i]^-1 * Nrd_mLK[i]^-1 * unram_char[i]^-1: i in [1..#res]];
    nr := [Index(IrrG,h[2]): h in Krel`QG`H];
    nr := [ unit[i] :  i in nr ];
    nr := [* (embeddings[i]^(-1))(nr[i]) : i in [1..#nr] *];
    rhs_omega2 := mapToRelKGroup(Krel, nr);
    rhs_omega1 := K0RelLog(Krel, setting`eLK);
    /* It remains to compute eLK + rhs_omega; */
    vprintf epsConj, 1: "rhs -----> %o\n", rhs_omega1 + rhs_omega2;

    setting`rhs := rhs_omega1 + rhs_omega2;

    /* Finally we compute (once again) the twisted unramified characteristic.
       Note that this is already computed and stored in setting`unram_char.
       So we will better delete it.
     */
    QG := setting`Krel`QG;
    P := Factorization(p*OL)[1,1];
    I := RamificationGroup(P, 0, h);
    eI := SubgroupIdempotent(QG, I);
    eG := SubgroupIdempotent(QG, G);
    Frob := QG ! Frobenius(OL, h, p) ;
    eins := One(QG);
    eta := (eins - eI) + Frob^-1 * eI;  
    incl_id := Coercion(setting`Qmpt, setting`Qmpt);
    Nrd_eta := NewtonReducedNormEG(eta, incl_id);
    nr := [Index(IrrG,h[2]): h in Krel`QG`H];
    nr := [ Nrd_eta[i] :  i in nr ];
    nr := [* (embeddings[i]^(-1))(nr[i]) : i in [1..#nr] *];
    setting`twisted_unram_char := mapToRelKGroup(Krel, nr);
    vprintf epsConj, 1: "twisted_unram_char computed\n";

    bool := true;

    return bool, setting;
end function;



EpsConjSetting := recformat<
    // General information about the situation
    L         : FldNum,         // global field
    K         : Any,         // global base field
    p         : RngIntElt,      // prime
    P         : RngOrdIdl,      // prime ideal above p
    pi        : RngOrdElt,      // uniformizing element of P
    LP        : FldPad,         // completion of L at P
    iota      : Map,            // inclusion L->LP
    G         : Grp,            // Galois group
    psiL      : Map,            // automorphism map G->Aut(L)
    psiLP     : Map,            // automorphism map G->Aut(LP)
    
    // requirements for epsilon constants
    // assigned by LECprepareEps
    IrrG      : SeqEnum,        // irreducible characters
    linChrs   : List,           // linear characters with abelian extensions
    brauerInd : List,           // Brauer inductions in degree 0
    t         : RngIntElt,      // computation precision for eps. constants
    Qmpt      : FldCyc,         // cyclic field for eps. constants
    
    // fields for global computations
    // assigned by LECpreparations
    N         : FldNum,         // unramified extension for unramified term
    LN        : FldNum,         // composite field L*N
    E         : FldNum,         // composite field Qmpt*LN
    Krel      : Rec,            // relative K-group
    
    // requirements for cohomological term
    // assigned by LECpreparations
    //calL      : Lat,            // lattice L
    theta     : RngOrdElt,      // normal basis element generating lattice
    prec      : RngIntElt,      // computation precision for cohom. term
    gamma     : Map,            // cocylcle
    //gammaInfo : Rec,
    
    // results of computations
    // assigned by LECcomputeValues
    uLK       : AlgGrpElt,      // unramified term in N[G]
    mLK       : AlgGrpElt,      // correction term in K[G]
    dLK       : AlgGrpElt,      // equiv. discriminant in L[G]
    tLK       : List,           // epsilon constants (galois gaus sums) (K(chi))_chi
    eLK       : Rec,            // cohomological term in K_0(Zp[G], Qp)
    
    // reduced norms
    // assigned by LECcomputeValues
    uLK_Nr    : List,
    mLK_Nr    : List,
    dLK_Nr    : List,
    
    // reduced extensions
    // assigned by LECimagesK0Rel
    E1        : Any,
    PE1       : RngOrdIdl,
    F         : Any,
    Fm        : Any,
    mipos     : SeqEnum,
    tupFm     : List,
    KrelF     : Rec,
    X         : List,
    Y         : List,
    Z         : List,
    
    // assigned by LECcheck
    result    : BoolElt,

    // new, used for pCubed
    OL        : RngOrd,
    M         : Any,
    iota_M    : Map,
    coeff_M   : SeqEnum,
    unram_char: List,
    Ath       : Rec,
    jLK       : List,
    dLKuLKinv : AlgGrpElt,
    dLKuLKinv_Nr : List,
    lhs       : List,
    rhs       : List,
    twisted_unram_char : List
>;


intrinsic LECcreateRec(L::FldNum, p::RngIntElt) -> Rec
{ Create LEC-record for the number field L and prime p.
  Creates record and computes automorphism map. }
    
    local setting, psiL, g, fac;
    
    setting := rec<EpsConjSetting | >; 
    setting`L := L;
    setting`p := p;
    
    // prime ideal
    fac := Factorization(p*RingOfIntegers(L));
    if #fac eq 1 then
        setting`K := Rationals();
        setting`P := fac[1,1];
    else
        print "The decomposed case is not yet implemented!";
        print "TODO:";
        print "  To compute the completion, L must be an extension of Q.";
        print "  For other algorithms, the extension is needed as tower L|K|Q.";
        error "Number field must be undecomposed at the given prime p!";
        setting`K := DecompositionField(setting`P);
        setting`L := RelativeField(setting`K,setting`L);
        fac := Factorization(p*RingOfIntegers(L));
        setting`P := fac[1,1];
    end if;
    setting`pi := UniformizingElement(setting`P);
    
    // Galois group
    _,_, psiL := AutomorphismGroup(L);
    psiL := map<Domain(psiL)->Codomain(psiL) | g:->psiL(g^-1)>;
    setting`G := Domain(psiL);
    setting`psiL := psiL;
    
    return setting;
end intrinsic;

intrinsic LECcreateRec(L::FldNum, p::RngIntElt, N::FldNum) -> Rec
{ Create LEC-record for the number field L and prime p.
  Creates record and computes automorphism map. }
    setting := LECcreateRec(L,p);
    setting`N := N;
    return setting;
end intrinsic;

intrinsic LECinitializeCompletion(~setting::Rec, prec::RngIntElt)
{ Compute the completion for the setting to the given precision.
  Actually, the completion is computed to a higher precision to ensure
  that Galois-conjugates can be computed to the given precision. }
    setting`LP, setting`iota, setting`psiLP := completion_with_prec(setting`L,setting`P, setting`psiL, prec);
    setting`G := Domain(setting`psiLP);
end intrinsic;

intrinsic LECinitializeCompletion(setting::Rec, prec::RngIntElt) -> Rec
{ Compute the completion for the setting to the given precision.
  Actually, the completion is computed to a higher precision to ensure
  that Galois-conjugates can be computed to the given precision. }
    local newSetting;
    newSetting := setting;
    LECinitializeCompletion(~newSetting, prec);
    newSetting`G := Domain(newSetting`psiLP);
    return newSetting;
end intrinsic;

intrinsic LECimagesK0Rel(~setting)
{ Read all the values of the Epsilon Constant Conjecture,
  as computed by LECcomputeValues, in the same relative K-group. }
    
    local R, x, tm, eLK, dLK_Nr, uLK_Nr, mLK_Nr, tLK_Nr, Krel, Qmpt, IrrG, tup,
          nr, E, r, embeddings, i, charFields, X,Y,Z, mipos,
          subfields, ED, Qm, EDm, incl_EDm_E, KrelED, nr2, t;
    
    // read values
    eLK := setting`eLK;
    dLK_Nr := setting`dLK_Nr;
    uLK_Nr := setting`uLK_Nr;
    mLK_Nr := setting`mLK_Nr;
    tLK_Nr := setting`tLK;

    Krel := setting`Krel;
    Qmpt := setting`Qmpt;
    IrrG := CharacterTable(setting`G);
    
    // tupel
    tup := [ tLK_Nr[i]*dLK_Nr[i]*uLK_Nr[i]*mLK_Nr[i] : i in [1..#tLK_Nr] ];
    // select non-conjugate values in tupel
    nr := [Index(IrrG,h[2]): h in Krel`QG`H];
    nr := [ tup[i] :  i in nr ];
    
    // create smaller extensions E', which still includes
    // the values of the tupel
    E := setting`E;
    vprintf epsConj, 2: "Reduce extension ";
    IndentPush();
    t := Cputime();
    vprintf epsConj, 3: "\n";
    vtime   epsConj, 3: mipos := [MinimalPolynomial(x) : x in nr];
    vprintf epsConj, 3: "Minimal polynomials: %o\n", mipos;
    
    if Max([Degree(f) : f in mipos]) le 4 then
        vtime epsConj, 3: E1 := SplittingField(mipos);
        // vprintf epsConj, 3: "E1: %o\n", E1;
    else
        // start with Q(zeta_m); we need it anyway later 
        vprint epsConj, 3: "using Qm as base field";
        Qm := SplittingField(Polynomial([1] cat [0: i in [1..Exponent(setting`G)-1]] cat [-1]));
        E1 := Qm;
        for f in Reverse(Sort(mipos)) do
            E1 := SplittingField(PolynomialRing(E1)!f);
        end for;
    end if;
    if Degree(E1) gt 1 then
        E1 := OptimizedRepresentation(E1);
        f := DefiningPolynomial(E1);
        c := LCM({Denominator(x) : x in Coefficients(f)});
        if c gt 1 then
            f := c^Degree(f) * Evaluate(f, Polynomial([0, 1/c]));
            E1 := NumberField(f);
            vprintf epsConj, 3: "E': %o\n", E1;
        end if;
    end if;
    t := Cputime(t);
    vprintf epsConj, 2: "Time: %o\n", t;
    IndentPop();
    
    
    if E1 eq Rationals() then
        F := Rationals();
        vprintf epsConj, 3: "No Decomposition in E'\n";
        r := 1;
        
    else
        OE1 := RingOfIntegers(E1);
        vprintf epsConj, 3: "Factorization of p in E'. ";
        vtime   epsConj, 3: PE1 := Factorization(setting`p*RingOfIntegers(OE1))[1,1];
        //F := DecompositionField(PE1);
        // -> bug
        vprintf epsConj, 3: "Compute decomposition field. ";
        vtime   epsConj, 3: F := FixedField(E1, [ sigma : sigma in Automorphisms(E1) |  &and([ sigma(x) in PE1   : x in Generators(PE1) ])]);
        r := Degree(F);
        vprintf epsConj, 3: "Decomposition r=%o\n", r;
        
        if r gt 1 then
            factors := Factorization(setting`p*RingOfIntegers(F));
            // embedding of F into E'
            iota := complex_compatible_embedding(F, E1, 100, false);
            // read prime ideals in E' and choose those belonging to PE1
            factorsE1 := [Factorization(ideal<OE1 | { iota(x) :  x in Generators(F[1])}>) : F in factors ];
            // fix those prime ideal
            PF := [ factors[i] : i in [1..#factors] | #factorsE1[i] eq 1 and factorsE1[i,1,1] eq PE1 ];
            assert #PF eq 1;
            PF := PF[1,1];
        end if;
    end if;
    
    
    // create composite field with character fields
    // F*Q(zeta_m)
    //Qm := CyclotomicField(Exponent(setting`G));
    Qm := SplittingField(Polynomial([1] cat [0: i in [1..Exponent(setting`G)-1]] cat [-1]));
    if Degree(Qm) eq 1 then
        Fm := E1;
    else
        Fm := AbsoluteField(ext<E1 | Factorization(DefiningPolynomial(Qm), E1)[1,1]>);
    end if;
    vprintf epsConj, 3: "Fm: %o\n", Fm;
    
    // read tupel of norms in Fm
    if Degree(E1) eq 1 then
        nr := [* Rationals()!nr[i] : i in [1..#nr]  *];
    else
        vprint epsConj, 3: "compute roots";
        nr := [* root_with_complex_embedding(mipos[i], Fm, nr[i], 100) : i in [1..#nr]  *];
    end if;
    
    setting`E1     := E1;
    setting`F      := F;
    setting`Fm     := Fm;
    setting`mipos  := mipos;
    setting`tupFm  := nr;
    
    vprint epsConj, 3: "compute images in K0rel";
    if r eq 1 then
        print "Case r = 1";
        // no other relative K group needed
        setting`KrelF := setting`Krel;
        
        // no decomposition in E1
        // vprintf epsConj, 3: "No decomposition of %o in E1!\n", setting`p;
        // read tupel directly in character fields Q(chi)
        embeddings := [* complex_compatible_embedding(Krel`QG`H[i,1], Fm, 100, false) : i in [1..#Krel`QG`H] *];
        nr := [* (embeddings[i]^(-1))(nr[i]) : i in [1..#nr] *];
        // read in relative K-group
        X := mapToRelKGroup(Krel, nr);
        print "Before KORelLog eLK = ", eLK;
        Y := K0RelLog(Krel, eLK);
        Z := X+Y;
        // each generator of an ideal is unit in ring of integers ?
        //setting`result := &and([ &and([ IsUnit(Krel`QG`H[i,4]!g) : g in Generators(Z[1,i])]) : i in [1..#Z[1]]])
        //    and Z[2] eq Zero(Parent(Z[2]));

    else
       print "Case r ne 1";
        
        // Consider K_0(Z_p[G], E_Q)
        KrelF := extendGlobalField(Krel, PF);
        setting`KrelF := KrelF;
        embeddings := [* complex_compatible_embedding(KrelF`EG`H[i,1], Fm, 100, false) : i in [1..#KrelF`EG`H] *];
        
        // read tupel in character fields F(chi)
        nr := [* (embeddings[i]^(-1))(nr[i]) : i in [1..#nr] *];
        vprintf epsConj, 4: "nr: %o\n", nr;
        // read in relative K-group
        X := mapToRelKGroup(KrelF, nr);
        
        // cohomological term
        // value in Q(chi)
        nr2 := reducedNorm(Krel, eLK);
        // read in F(chi)
        nr2 := [* KrelF`EG`H[i,3](nr2[i])  : i in [1..#nr2] *];
        // and map to K-group
        Y := mapToRelKGroup(KrelF, nr2);
        
        // sum
        Z := X+Y;
        // each generator of an ideal is unit in ring of integers ?
        //setting`result := &and([ &and([ Valuation(g, KrelF`PrimeIdeals[i,1]) eq 0 : g in Generators(Z[1,i])] ) : i in [1..#Z[1]]])
        //    and Z[2] eq Zero(Parent(Z[2]));
    end if;
    
    setting`X := X;
    setting`Y := Y;
    setting`Z := Z;
end intrinsic;


intrinsic pCubedLECimagesK0Rel(~setting)
{ Read all the values of the Epsilon Constant Conjecture,
  as computed by LECcomputeValues, in the same relative K-group. }
    
    local R, x, tm, eLK, dLK_Nr, uLK_Nr, mLK_Nr, tLK_Nr, Krel, Qmpt, IrrG, tup,
          nr, E, r, embeddings, i, charFields, X,Y,Z, mipos,
          subfields, ED, Qm, EDm, incl_EDm_E, KrelED, nr2, t;
    
    // read values
    eLK := setting`eLK;
    dLKuLKinv_Nr := setting`dLKuLKinv_Nr;
    mLK_Nr := setting`mLK_Nr;
    tLK_Nr := setting`tLK;

    Krel := setting`Krel;
    Qmpt := setting`Qmpt;
    IrrG := CharacterTable(setting`G);
    
    // tupel
    tup := [ tLK_Nr[i]*dLKuLKinv_Nr[i]*mLK_Nr[i] : i in [1..#tLK_Nr] ];
    
idem := QmptIdempotents(setting);
Z := &+[tup[i]*idem[i] : i in [1..#tup]];
coeff_Z := [Rationals()!a : a in Coefficients(Z)];
print "Z ist rational.";

    // select non-conjugate values in tupel
    nr := [Index(IrrG,h[2]): h in Krel`QG`H];
    nr := [ tup[i] :  i in nr ];
    
        print "Case r = 1";
        // no other relative K group needed
        setting`KrelF := setting`Krel;
        
        // no decomposition in E1
        // vprintf epsConj, 3: "No decomposition of %o in E1!\n", setting`p;
        // read tupel directly in character fields Q(chi)
        Fm := setting`E;
        embeddings := [* complex_compatible_embedding(Krel`QG`H[i,1], Fm, 100, false) : i in [1..#Krel`QG`H] *];
        nr := [* (embeddings[i]^(-1))(nr[i]) : i in [1..#nr] *];
        // read in relative K-group
        X := mapToRelKGroup(Krel, nr);
        Y := K0RelLog(Krel, eLK);
        Z := X+Y;
        // each generator of an ideal is unit in ring of integers ?
        //setting`result := &and([ &and([ IsUnit(Krel`QG`H[i,4]!g) : g in Generators(Z[1,i])]) : i in [1..#Z[1]]])
        //    and Z[2] eq Zero(Parent(Z[2]));
    
    setting`X := X;
    setting`Y := Y;
    setting`Z := Z;
end intrinsic;


intrinsic LECcheck(~setting)
{ Verify the local epsilon constant conjecture for the given
  setting, where the reduced norms are already computed. }
    
    setting`result := &and([ &and([ Valuation(g, setting`KrelF`PrimeIdeals[i,1]) eq 0 : g in Generators(Z[1,i])] ) : i in [1..#Z[1]]])
            and Z[2] eq Zero(Parent(Z[2]))
            where Z := setting`Z;
end intrinsic;

intrinsic LECcheck(setting) -> BoolElt
{ Verify the local epsilon constant conjecture for the given
  setting, where the reduced norms are already computed. }
    
    local setting1;
    setting1 := setting;
    LECcheck(~setting1);
    return setting1`result;
end intrinsic;




/*******************************
 *  Equivariant Discriminant
 *******************************/

// Berechne Repraesentant fuer [L, rho, H]
// nach [Bley/Breuning, 4.2.5]
intrinsic LECdiscriminant(psi::Map, theta::RngOrdElt) -> AlgGrpElt
{ Compute the equivariant discriminant of a lattice as described in [Bley/Breuning, 4.2.5]. }
    
    local L, LG, G, v, v1, dLK;
    
    G := Domain(psi);
    L := Domain(psi(G.1));
    LG := GroupAlgebra(L, G);
    dLK := &+([psi(g)(theta)*LG!(g^(-1)) : g in G]);
    
    return dLK;
end intrinsic;

intrinsic LECdiscriminant(psi::Map, P::RngOrdIdl) -> AlgGrpElt
{ Compute the equivariant discriminant of a lattice as described in [Bley/Breuning, 4.2.5]. }
    
    return LECdiscriminant(psi, lattice_generator_theta(psi, P, UniformizingElement(P)));
end intrinsic;

intrinsic LECdiscriminant(~setting::Rec)
{ Compute the equivariant discriminant of a lattice as described in [Bley/Breuning, 4.2.5]. }
    
    if (not assigned setting`theta) then
        vprint epsConj, 5: "Compute normal basis element theta";
        setting`theta := lattice_generator_theta(setting`psiL, setting`P, setting`pi);
    else
        vprint epsConj, 5: "Use saved normal basis element theta";
    end if;
    
    vprintf epsConj, 4: "Valuation v(theta)=%o\n", Valuation(setting`theta, setting`P);
    setting`dLK := LECdiscriminant(setting`psiL, setting`theta);
    vprintf epsConj, 3: "dLK=%o\n", setting`dLK;
end intrinsic;

intrinsic LECdiscriminant(setting::Rec) -> AlgGrpElt
{ Compute the equivariant discriminant of a lattice as described in [Bley/Breuning, 4.2.5]. }
    
    local setting1;
    
    if (not assigned setting`theta) then
        print "!! WARNING !!";
        print "Generator of lattice was not set when computing the equivariant discriminant.";
        print "Pass the parameter by reference, if you want to continue verifying the local epsilon constant conjecture!";
        print "!! WARNING !!";
    end if;
    setting1 := setting;
    LECdiscriminant(~setting1);
    return setting1`dLK;
end intrinsic;





/*******************************
 *  Correction Term
 *******************************/

// Berechnet den Korrekturterm $M_{L/K}$ der lokalen Epsilonkonstantenvermutung
// nach [Breuning PhD, 4.4].
intrinsic LECcorrectionTerm(setting::Rec) -> .
{ Compute the correction term for the local epsilon constant conjecture. }
    //return LECcorrectionTerm(Krel, setting`L, setting`G, setting`psiL, setting`P, setting`p);
    return LECcorrectionTerm(setting`Krel`QG, setting`psiL, setting`P);
end intrinsic;

intrinsic LECcorrectionTerm(~setting::Rec)
{ Compute the correction term for the local epsilon constant conjecture. }
    //setting`mLK := LECcorrectionTerm(Krel, setting`L, setting`G, setting`psiL, setting`P, setting`p);
    setting`mLK := LECcorrectionTerm(setting`Krel`QG, setting`psiL, setting`P);
end intrinsic;

intrinsic LECcorrectionTerm(QG::AlgGrp[FldRat], psi::Map, P::RngOrdIdl) -> AlgGrpElt
{ Compute the correction term for the local epsilon constant conjecture. }
    local G, L, OL, I, F, OF, H, res, psiH, PF, sigma, K, OK, PK, q,
          eI, eG, mLK;
    
    require Group(QG) eq Domain(psi) :
            "Group of Q[G] and Domain(psi) must be equal!";
    require Degree(Order(P)) eq #Group(QG) :
            "Group of Q[G] must be Galois group of L over its base field!";
    
    G := Domain(psi);
    L := Domain(psi(Id(G)));
    OL := RingOfIntegers(L);
    I := inertia_group(P, psi);
    vprintf epsConj, 3: "Inertia group #I=%o\n", #I;
    
    K := BaseField(L);
    OK := RingOfIntegers(K);
    PK := P meet OK;
    
    // Berechne den Frobenius von G/I
    F := fixed_field(psi, I, []);
    vprintf epsConj, 3: "Fixed field if I is %o\n", F;
    OF := RingOfIntegers(F);
    H, res := quo<G | I>;
    psiH := restrict_aut_to_quotient(psi, F, H, res);
    if Type(PK) eq RngInt then
        PF := Factorization((Generators(PK)[1])*OF)[1,1];
    else
        PF := Factorization(Parent(1*OF)!PK)[1,1];
    end if;
    
    sigma := frobenius_grp_elt(F, PF, H, psiH);
    // Representant in G
    sigma := (res^(-1))(sigma);
    
    if (AbsoluteDegree(K) eq 1) then
        q := Generators(PK)[1];
    else
        q := Index(OK, PK);
    end if;
    
    //Elemente in Q[G]
    //K0Rel := RelativeGroup(G, 3); 
    eI := 1/#I * &+([QG!g : g in I]);
    eG := 1/#G * &+([QG!g : g in G]);
    
    //vprint epsConj, 3: "Elemente", [* #quo<G|I> * eG, (1-1/q * QG!sigma) *eI , (1-QG!(sigma^(-1)))*eI *];
    //vprint epsConj, 3: "*Elemente (inv)", [* star(#quo<G|I> * eG) , star((1-1/q * QG!sigma) *eI) , star((1-QG!(sigma^(-1)))*eI)^(-1) *];
    
    mLK := star(#quo<G|I> * eG) * star((1-1/q * QG!sigma) *eI) * star((1-QG!(sigma^(-1)))*eI)^(-1);
    vprintf epsConj, 3: "mLK = %o\n", mLK;
    return mLK;
end intrinsic;



/* star(x::AlgGrpElt) -> AlgGrpElt
 * Apply the star-map as defined in [Breu PhD, 4.1.1] to x in Q[G], i.e.
 * read x=(x_i)\in\prod\C, replace x_i=0 by x_i=1 and return
 * resulting element in \C[G].
 * Uses functions from the RelAlgKTheory-package and therefore depends
 * on the attributes OC, H and X of Q[G] to be set, as it is done
 * by RelativeGroup(G,p).
 */
function star(x)
    assert(IstZentral(x));
    // x in zeta(Q[G])
    QG := Parent(x);
    // compute Tupel in prod K_i = zeta(Q[G])
    seq := PsiInv(x, QG`H, QG`OC);
    // write 1 instead of 0
    seq := [* s eq 0 select Parent(s)!1 else s : s in seq*];
    // compute element in Q[G]
    y := Psi(seq, QG`X , QG`H, QG) ;
    return y;
end function;





/*******************************
 *  Unramified Term
 *******************************/

//
intrinsic LECunramifiedTerm(~setting::Rec)
{ Compute the unramified term in N[G] for an extension L/Q, a prime p,
  psi: Gal(L/Q) -> \Aut(L), and N an unramified extension
  with [N:Q]=[L^ab/Q] (see [Bley/Breuning, 4.2.7]). }
    setting`uLK := LECunramifiedTerm(setting);
end intrinsic;
intrinsic LECunramifiedTerm(setting::Rec) -> AlgGrpElt
{ Compute the unramified term in N[G] for an extension L/Q, a prime p,
  psi: Gal(L/Q) -> \Aut(L), and N an unramified extension
  with [N:Q]=[L^ab/Q] (see [Bley/Breuning, 4.2.7]). }
    return LECunramifiedTerm(setting`psiL, setting`p, setting`N);
end intrinsic;

//
// Gegeben eine Erweiterung L/Q mit Galoisgruppe G, psi:G -> \Aut(L)
// und Erweiterung N/Q vom gleichen Grad wie L^{ab}/Q.
// Berechnet den unverzweigten Term U_{L/Qp}.
//
intrinsic LECunramifiedTerm(psi::Map, p::RngIntElt, N::FldNum) -> AlgGrpElt
{ Compute the unramified term in N[G] for an extension L/Q, a prime p,
  psi: Gal(L/Q) -> \Aut(L), and N an unramified extension
  with [N:Q]=[L^ab/Q] (see [Bley/Breuning, 4.2.7]). }

    local G, L, H, F, Fa, a, PP, psiF, res, phi, s, OF, PF, C, kappa, xi,
          N1, ON, PN, GN, psiN, NG, QGN;

    //assert(IsCyclic(N));

    G := Domain(psi);
    L := Domain(psi(G.1));
    // Maximal abelsche Teilerweiterung
    H := CommutatorSubgroup(G);
    F := FixedField(L, H);
    Fa := AbelianExtension(F);

    // Galoisgruppe GF von F/Q und psiF:GF -> Aut(F)
    GF, res := quo<G | H >;
    psiF := restrict_aut_to_quotient(psi, F, GF, res);

    //lokales Normrestsymbol (p, F/Q)
    //vprint epsConj, 2: "Berechne globales Ideal";
    //a := localNormResidueSymbolAsGlobalIdeal(Rationals()!p, Factorization(Conductor(Fa) meet Integers()), p*Integers());
    vprint epsConj, 3: "Compute Artin map";
    //phi := artinMap(a, F, GF, psiF);
    phi := localNormResidueSymbol(Rationals()!p, Fa, p*Integers(), psiF);

    // phi in GF, waehle Representant in G
    s := Order(phi);
    phi := (res^(-1))(phi);
    vprint epsConj, 5: "Using phi", phi;


    if (s eq 1) then
        vprint epsConj, 5: "phi has order 1";
        // phi = Id(G)
        // N1 = Q
        // xi = 1
        NG := GroupAlgebra(N, G);
        //u := &+([ N!N1!((psiN(f^i))(xi)) *NG!(phi^(-i))  : i in [0..s-1]]);
        u := N!1 *NG!(phi^(0));

    else
        //Unverzweigte Erweiterung initialisieren
        if s ne Degree(N) then
            vprintf epsConj, 5: "Subfield of degree %o\n", s;
            // Brauche ich hier zyklisch ???
            // teste mit IsCylic !!
            //N1 := [F : F in Subfields(N, s) |  IsCyclotomicPolynomial(DefiningPolynomial(F[1]))][1,1];
            N1 := [F : F in Subfields(N, s)][1,1];
            //N1 := Subfields(N,s)[1,1];
        else
            N1 := N;
        end if;

        vprintf epsConj, 3: "Using field of degree %o\n", Degree(N1);

        ON := RingOfIntegers(N1);
        PN := Factorization(p*ON)[1,1];
        GN, _, psiN := AutomorphismGroup(N1);
        psiN := map<Domain(psiN)->Codomain(psiN) | g:->psiN(g^-1)>;

        vprint epsConj, 3: "Compute Frobenius and integral normal basis element";
        // Frobenius von N/Q berechnen  
        f := frobenius_grp_elt(N1, PN, GN, psiN);
        vprintf epsConj, 5: "Using Frobenius %o\n", f;

        // Normalbasiselement
        //xi := NormalBasisElement(ON, psiN);
        // N ist zyklisch, da unverzweigt
        lst := Abelian_nb(DefiningPolynomial(N1));
        // --> [K, theta, Atheta, AssOrd, nb, h, QG]
        xi := QGAction(lst[7] ! ElementToSequence(lst[5]), lst[2], lst[6]);

        bool, m := IsIsomorphic(lst[1],N1);
        assert bool;
        xi := ON!m(xi);
        vprint epsConj, 5: "Using integral normal basis element %o\n", xi;


        // Tests
        vprint epsConj, 4: "Test results";
        // f ist Frobenius in Gal(N1/Q)
        R, projR := ResidueClassField(ON,PN);
        assert {projR(psiN(f)(b))-(projR(b))^p : b in Basis(ON)} eq {0};
        // phi|_F = (p, F/Qp)
        // d.h. psiF(phi) ist Frobenius, da (p, F/Qp)=Frob(F/Qp)^vp(p)=Frob(F/Qp)
        PF := Factorization(p*RingOfIntegers(F))[1,1];
        R, projR := ResidueClassField(RingOfIntegers(F),PF);
        m := psiF(res(phi));
        assert { projR(m(b))-projR(b)^p : b in Basis(RingOfIntegers(F))} eq {0};


        //u := [ <(psiN(f^i))(xi), phi^(-i)>  : i in [0..s-1]];
        NG := GroupAlgebra(N, G);
        //vprint epsConj, 3: "element", [ N!N1!((psiN(f^i))(xi)) *NG!(phi^(-i))  : i in [0..s-1]];
        u := &+([ N!N1!((psiN(f^i))(xi)) *NG!(phi^(-i))  : i in [0..s-1]]);
    end if;

    vprintf epsConj, 3: "uLK = %o\n", u;
    return u;
end intrinsic;


//
intrinsic pCubedLECunramifiedTerm(~setting::Rec)
{ Compute the unramified term in N[G] for an extension L/Q, a prime p,
  psi: Gal(L/Q) -> \Aut(L), and N an unramified extension
  with [N:Q]=[L^ab/Q] (see [Bley/Breuning, 4.2.7]). }
    setting`uLK, setting`N := pCubedLECunramifiedTerm(setting);
end intrinsic;

/*
intrinsic pCubedLECunramifiedTerm(setting::Rec) -> AlgGrpElt
{ Compute the unramified term in N[G] for an extension L/Q, a prime p,
  psi: Gal(L/Q) -> \Aut(L), and N an unramified extension
  with [N:Q]=[L^ab/Q] (see [Bley/Breuning, 4.2.7]). }
    return LECunramifiedTerm(setting`psiL, setting`p);
end intrinsic;
*/

//
// Gegeben eine Erweiterung L/Q mit Galoisgruppe G, psi:G -> \Aut(L)
// und Erweiterung N/Q vom gleichen Grad wie L^{ab}/Q.
// Berechnet den unverzweigten Term U_{L/Qp}.
//
intrinsic pCubedLECunramifiedTerm(setting :: Rec) -> AlgGrpElt, FldNum
{ Compute the unramified term in N[G] for an extension L/Q, a prime p,
  psi: Gal(L/Q) -> \Aut(L), and N an unramified extension
  with [N:Q]=[L^ab/Q] (see [Bley/Breuning, 4.2.7]). }
    
    local G, L, H, F, Fa, a, PP, psiF, res, phi, s, OF, PF, C, kappa, xi, 
          N1, ON, PN, GN, psiN, NG, QGN;
    
    //assert(IsCyclic(N));
    
    psi := setting`psiL; p := setting`p;  
    G := Domain(psi);
    L := Domain(psi(G.1));
    // Maximal abelsche Teilerweiterung
    H := CommutatorSubgroup(G);
    F := FixedField(L, H);
    Fa := AbelianExtension(F);
    
    // Galoisgruppe GF von F/Q und psiF:GF -> Aut(F)
    GF, res := quo<G | H >;
    psiF := restrict_aut_to_quotient(psi, F, GF, res);
    
    //lokales Normrestsymbol (p, F/Q)
    //vprint epsConj, 2: "Berechne globales Ideal";
    //a := localNormResidueSymbolAsGlobalIdeal(Rationals()!p, Factorization(Conductor(Fa) meet Integers()), p*Integers());
    vprint epsConj, 3: "Compute Artin map";
    //phi := artinMap(a, F, GF, psiF);
    phi := localNormResidueSymbol(Rationals()!p, Fa, p*Integers(), psiF);
    
    // phi in GF, waehle Representant in G
    phi := (res^(-1))(phi);
    vprint epsConj, 5: "Using phi", phi;
    s := Order(phi);
    if s eq 3 then
        N := FixedField(L, RamificationGroup(setting`P, 0, setting`psiL));   
    else
        N := Subfields(CyclotomicField(19), s)[1,1];
    end if;    

    ON := RingOfIntegers(N);
    PN := Factorization(p*ON)[1,1];
    GN, _, psiN := AutomorphismGroup(N);
    psiN := map<Domain(psiN)->Codomain(psiN) | g:->psiN(g^-1)>;
        
    vprint epsConj, 3: "Compute Frobenius and integral normal basis element";
    // Frobenius von N/Q berechnen  
    f := frobenius_grp_elt(N, PN, GN, psiN);
    vprintf epsConj, 5: "Using Frobenius %o\n", f;
        
    lst := Abelian_nb(DefiningPolynomial(N));
    // --> [K, theta, Atheta, AssOrd, nb, h, QG]
    xi := QGAction(lst[7] ! ElementToSequence(lst[5]), lst[2], lst[6]);
        
    bool, m := IsIsomorphic(lst[1],N);
    assert bool;
    xi := ON!m(xi);
    vprint epsConj, 5: "Using integral normal basis element %o\n", xi;
        
        
    // Tests
    vprint epsConj, 4: "Test results";
    // f ist Frobenius in Gal(N1/Q)
    R, projR := ResidueClassField(ON,PN);
    assert {projR(psiN(f)(b))-(projR(b))^p : b in Basis(ON)} eq {0};
    // phi|_F = (p, F/Qp)
    // d.h. psiF(phi) ist Frobenius, da (p, F/Qp)=Frob(F/Qp)^vp(p)=Frob(F/Qp)
    PF := Factorization(p*RingOfIntegers(F))[1,1];
    R, projR := ResidueClassField(RingOfIntegers(F),PF);
    m := psiF(res(phi));
    assert { projR(m(b))-projR(b)^p : b in Basis(RingOfIntegers(F))} eq {0};
       
        
    //u := [ <(psiN(f^i))(xi), phi^(-i)>  : i in [0..s-1]];
    NG := GroupAlgebra(N, G);
    //vprint epsConj, 3: "element", [ N!N1!((psiN(f^i))(xi)) *NG!(phi^(-i))  : i in [0..s-1]];
    u := &+([ N!((psiN(f^i))(xi)) *NG!(phi^(-i))  : i in [0..s-1]]);
   
    vprintf epsConj, 3: "uLK = %o\n", u;
    return u, N;
end intrinsic;

/* extension_for_unramified_term(L::FldNum, p::RngIntElt, G::Grp) -> FldNum
 * Computes a number field N of degree [N:Q]=[L^ab:Q] which is unramified at p;
 * N is computed as a subfield of a suitable cyclotomic field.
 * If the degree of N is greater than 4, only a heuristical method is implemented.
 */
function extension_for_unramified_term(L, P, G)
    local F, p, N;
    
    p := Generators(P meet Integers())[1];
    F  := FixedField(L, CommutatorSubgroup(G));
    //PF := Factorization(p*RingOfIntegers(F))[1,1];
    K := FixedField(L,G);
    // L/F hat Grad Degree(L)/Degree(F)
    // F/K hat Grad Degree(F)/Degree(K)
    
    // brauche unverzweigte Erweiterung N/Q vom gleichen Grad wie F/Q
    if InertiaDegree(P) mod Integers()!(Degree(F)/Degree(K)) eq 0 then
        // nehme Teilkoerper von InertiaField
        N := InertiaField(P);
        if Degree(N) eq Degree(F) then
            return N;
        else
            // Teilkoerper
            return Subfields(N, Degree(F))[1,1];
        end if;
        
    elif InertiaDegree(P) eq 2 and Degree(F) eq 4 then
        // versuche C2 in C4 einzubetten
        K := InertiaField(P);
        bool, g := embeddingC2C4(K : p := p);
        if bool then
            return NumberField(g);
        else
            // besser: C4 Erweiterung ohne Einbettung
            return unramified_cyclotomic_subfield(Degree(F), p);
        end if;
        
    else
        // erzeuge zyklotomischen Koerper
        vprintf epsConj, 3: "Using heuristic to find a suitable unramified extenssion!\n";
        return unramified_cyclotomic_subfield(Degree(F), p);
    end if;
end function;

/* Heuristically generates an extension of degree d which is unramified at p
 * as a subfield of a cyclotomic field
 */
function unramified_cyclotomic_subfield(d, p)
    local L, n;
    
    n := 0;
    repeat
        n := n+1;
        if n mod 500 eq 0 then
            if n eq 500 then
                print "!! WARNING !!";
                print "The heuristical seach for a cyclotomic field Q(zeta_n) which is unramified at p exceeded n = ", n;
            end if;
            print "... exceeded n = ", n;
        end if;
        if EulerPhi(n) mod d eq 0 then
            // can have subfield of degree d
            L := CyclotomicField(n);
            seq := DecompositionType(L,p);
            if seq[1,2] eq 1 then
                // unramified
                subL := [K[1] : K in Subfields(CyclotomicField(n),d) | 
                    #seq eq 1 and seq[1,2] eq 1 where seq is DecompositionType(K[1], p) ];
                if #subL gt 0 then
                    return subL[1];
                end if;
            end if;
        end if;
    until false;
end function;





/*******************************
 *  Cohomological Term
 *******************************/

/**
 * The main function for the computation of the cohomological term
 * is LECcohomologicalTerm.
 * It depends on several attributes in the LEC-record, as computed
 * by LECcreateRec and LECpreparations.
 * The algorithm first computes a cocycle for the local fundamental
 * class and then continues as in [Bley/Breuning, 4.2.4] by
 * computing the splitting module C(gamma), its projective resolution,
 * and finally the Q[G]-isomorphism between K+Q[G] and Q[G]^r from
 * [Bley/Burns, 3.3].
 */

// nach [Bley/Breuning, 4.2.4]
intrinsic LECcohomologicalTerm(~setting::Rec)
{ Compute the cohomological term for the local epsilon constant conjecture
  as desribed in [Bley/Breuning, 4.2.4]. }
    local c, gamma, gammaM, s, GG;
    
    setting`eLK := LECcohomologicalTerm(setting);
end intrinsic;

intrinsic LECcohomologicalTerm(setting::Rec) -> Rec
{ Compute the cohomological term for the local epsilon constant conjecture
  as desribed in [Bley/Breuning, 4.2.4]. }
    local s, X, c, gamma, gammaM, Cgamma, actionCgamma, actionCgammaMat, 
          zGBasisCgamma, K, psi_FCgamma, F,
          GG, g, G, mats, seq, i, h, M, actionFMat, actionF,
          CgammaQ, actionCgammaQMat, FQ, actionFQMat, zGBasisCgammaQ,
          actionCgammaQ, actionFQ,
          E_LK;
    
    vprintf epsConj, 2: "In LECcohommologicalTerm epsConj = %o\n", GetVerbose("epsConj");
    // return rec< K0RelElt | >;

    s := setting`prec;
    G := setting`G;
    GG := [g : g in setting`G];
    if (not assigned setting`gamma) then
        vprintf epsConj, 2: "Compute cocycle, precision %o\n", s;
        //c := cocycleLFC(setting, s+1);
        IndentPush();
        gamma := CocycleLFC(setting`LP, pAdicField(setting`LP), s : psi := setting`psiLP);
        IndentPop();
        //vprint epsConj, 1: "invertiert!";
        //gamma := map<Domain(c) -> setting`LP | x :-> valueModPn(c(x),setting`iota(setting`pi),s,setting`LP)^(-1) >;
        //vprint epsConj, 1: "nicht invertiert!";
        //gamma := map<Domain(c) -> setting`LP | x :-> valueModPn(c(x),setting`iota(setting`pi),s,setting`LP) >;
        
        //vprintf epsConj, 2: "Berechne alle Werte des Kozykels";
        //vprint epsConj, 3: "G-Aufzaehlung:", GG;
        gammaM := Matrix(Codomain(gamma), [[ gamma(g,h) : h in GG] : g in GG]);
        // ca 70 sec
        gamma := map< Domain(gamma) -> Codomain(gamma) | x :-> gammaM[Index(GG, x[1]), Index(GG, x[2])] >;
        
    else
        //vprint epsConj, 1: "Nutze gespeicherten Kozykel";
        gamma := setting`gamma;
    end if;
    
    vprint epsConj, 2: "Compute splitting module";
    Cgamma, actionCgammaMat := splitting_module_LPmul_modX(setting, gamma);
    vprint epsConj, 5: "C(gamma):",Cgamma;
    vprint epsConj, 5: "G-Wirkung auf C(gamma):", actionCgammaMat;
    
    actionCgamma := map< car<G,Cgamma> -> Cgamma | x:-> Cgamma!Vector(ElementToSequence(x[2]*actionCgammaMat[Index(GG,x[1])])) >;
    if GetVerbose("epsConj") ge 5 then
        // Teste G-Wirkung auf Cgamma
        //assert test_G_action(Cgamma, actionCgamma);
        assert test_G_action(actionCgamma);
    end if;
    
    vprint epsConj, 2: "Compute projective resolution";
    zGBasisCgamma, K, psi_FCgamma := proj_resolution_ZGr(Cgamma, actionCgamma);
    vprint epsConj, 5: "zGBasis C(gamma):", zGBasisCgamma;
    vprint epsConj, 5: "projektive Aufloesung, K:", K;
    
    // F = Z[G]^r
    F := Domain(psi_FCgamma);
    // Jetzt:  0 --> K --> F --> C(gamma) --> 0  exakt
    
    // G-Wirkung auf F
    G := setting`G;
    GG := [g : g in G];
    actionFMat := [];
    for g in GG do
        // Rechtswirkung:
        // seq := [Vector(Insert([0 : i in [2..#G]],s,1))   : s in [Index(GG, h*g) : h in GG]];
        // Linkswirkung:
        seq := [Vector(Insert([0 : i in [2..#G]],s,1))   : s in [Index(GG, g*h) : h in GG]];
        // Reihenfolge sehr wichtig !!!                              hier  ^^^
        M := VerticalJoin(seq);
        M := DiagonalJoin([M : i in [1..Degree(F)/#G]]);
        Append(~actionFMat, M);
    end for;
    //actionFMat := map< G -> GL(Degree(K), Rationals())  | g :-> mats[Index(GG, g)] >;
    actionF := map< car<G,F> -> F | x :-> x[2]*actionFMat[Index(GG,x[1])] >;
    
    if GetVerbose("epsConj") ge 5 then
        // Teste G-Vertraeglichkeit von  F --> Cgamma
        assert test_G_compatible(psi_FCgamma, actionF, actionCgamma, false, 0, 0);
    end if;
    
    // tensorieren mit Q
    CgammaQ, actionCgammaQMat := tensor_Zmod_with_Q(Cgamma, actionCgammaMat);
    FQ, actionFQMat := tensor_Zmod_with_Q(F, actionFMat);
    
    // Z[G]-Basis von C(gamma)_Q
    zGBasisCgammaQ := [ CgammaQ!Vector(ElementToSequence(b)[1..Degree(CgammaQ)]) : b in zGBasisCgamma];
    vprint epsConj, 5: "zGBasis C(gamma)_Q:", zGBasisCgammaQ;
    vprint epsConj, 5: "G-Wirkung auf C(gamma)_Q:", actionCgammaQMat;

    // maps from sequences
    actionCgammaQ := map< car<G,CgammaQ> -> CgammaQ | x:-> x[2]*actionCgammaQMat[Index(GG,x[1])] >;
    actionFQ      := map< car<G, FQ>     ->  FQ     | x:-> x[2]*actionFQMat[Index(GG,x[1])] >;
    
    if GetVerbose("epsConj") ge 5 then
        // Teste G-Vertraeglichkeit von  F_Q --> Cgamma_Q
        psi_Cgamma_CgammaQ := map< Cgamma -> CgammaQ | x :-> CgammaQ!ElementToSequence(x)[1..Degree(CgammaQ)] >;
        psi_FQCgammaQ := map< FQ -> CgammaQ | x :->  psi_Cgamma_CgammaQ( psi_FCgamma(F!x) ) >;
        assert test_G_compatible(psi_FQCgammaQ, actionFQ, actionCgammaQ, false, 0, 0);
    end if;
    
    // G-Wirkung auf K
    mats := [];
    for g in GG do
        M := VerticalJoin([Vector(Rationals(), Coordinates(K, K!actionF(g, b))) : b in Basis(K)]);
        Append(~mats, M);
    end for;
    actionKMat := map< G -> GL(Dimension(K), Rationals()) | g :-> mats[Index(GG,g)]>;
    //actionK := map< car<G, K> -> K | x :-> addProducts(Vector(Rationals(), Coordinates(K,x[2]))*actionKMat(x[1]), Basis(K)) >;
    actionK := map< car<G, K> -> K | x :-> &+([ a[i]*b[i] : i in [1..NumberOfColumns(a)] ]
        where a is Vector(Rationals(), Coordinates(K,x[2]))*actionKMat(x[1]) 
        where b is Basis(K)) >;
    
    // checks
    // K liegt im Kern  F --> C(gamma)
    assert &and([ psi_FCgamma(b) eq Zero(Cgamma) :  b in Basis(K)]);
    // K --> F ist G-aequivariant
    assert &and([ F!actionK(g,b) eq actionF(g, F!b) : g in G, b in Basis(K)]);
    // F --> C(gamma) ist G-aequivariant
    assert &and([  psi_FCgamma(actionF(g,b)) eq actionCgamma(g, psi_FCgamma(b)) : g in G, b in Basis(F) ]);
    // G-Wirkung auf K
    if GetVerbose("epsConj") ge 5 then
        //assert test_G_action(K, actionK);
        assert test_G_action(actionK);
    end if;
        
    // Element in K0Rel
    vprint epsConj, 2: "Compute element in relative K-group";
    // Rec-ZGModule erwartet Bilder in Spalten !
    actionKMatT := map< G -> GL(Dimension(K), Rationals()) | g :-> Transpose(actionKMat(g)) >;
    Krec := ZGModuleInit(IdentityMatrix(Rationals(), Dimension(K)), actionKMatT);
    QG := GroupAlgebra(Rationals(), G);
    // Z[G] als Z[G] Modul mit Linkswirkung:
    A := RegularModule(RegularRep(QG));
    // Z[G] als Z[G] Modul mit Rechtswirkung:
    // A := RegularModule(RegularRightRep(QG));
    KplusZG := ZGModuleDirectSum(Krec,A);
    
    // F als Rec-ZGModule
    Frec := A;
    for i in [2..Dimension(F)/#G] do
        Frec := ZGModuleDirectSum(Frec,A);
    end for;
    
    // Isomorphismus
    M := isomorphism_matrix_KQF(K, CgammaQ, actionCgammaQ, FQ, actionFQ, zGBasisCgammaQ);
    M := Transpose(M);
    
    if GetVerbose("epsConj") ge 5 then
        // Teste G-Vertraeglichkeit des Isomorphismus
        assert &and({M*KplusZG`phi(g)*KplusZG`hnf eq Frec`phi(g)*M*KplusZG`hnf : g in G});
    end if;
    
    // C_LK := 
    E_LK := rec< K0RelElt | >;
    E_LK`A := KplusZG;
    E_LK`B := Frec;
    E_LK`Theta := M;
    
    return E_LK;
end intrinsic;


intrinsic LEClattice(~setting : pCubed:=false)
{ Compute a few suitable lattices for the conjecture and an
  integer m such that the lattice includes P^m.
  Choose the (computationally) best one for further computations. }
    
    local theta, m, theta1, m1, i;
    
    if pCubed then
        print "In LEClattice: pCubed = true.";
        theta, m := LEClattice(setting`P, setting`pi, setting`psiL : pCubed:=true);
    else
        theta, m := LEClattice(setting`P, setting`pi, setting`psiL);
        for i in [2..6] do
            theta1, m1 := LEClattice(setting`P, setting`pi, setting`psiL);
            if (m1 lt m) then
                theta := theta1;
                m := m1;
            end if;
        end for;
     end if;
    
    setting`theta := theta;
    setting`prec := m;
end intrinsic;

intrinsic LEClattice(P::RngOrdIdl, pi::RngOrdElt, psi::Map : pCubed:=false) -> FldNumElt, RngIntElt
{ Given a prime ideal P of L with uniformizing element pi and automorphism
  map psi:G->Aut(L). Compute a generator theta of a suitable lattice and an
  integer m such that the lattice includes P^m.}
    
    if pCubed then
        return lec_lattice(P,pi,psi : pCubed:=true);
    else
        return lec_lattice(P,pi,psi);
    end if;
end intrinsic;

/* lec_lattice(P::RngOrdIdl, pi::RngOrdElt, psiL::Map : pCubed:=false) -> FldNumElt, RngIntElt
 * See LEClattice.
 */
function lec_lattice(P, pi, psiL : pCubed:=false)
    local L, OL, p, theta, v, v1, erz, x, M, ZpGtheta, k, m, M1, M1I, j, M2, T;
    
    OL := Order(P);
    L := NumberField(OL);
    p := Generator(P meet Integers());
    G := Domain(psiL);
    
    if pCubed then
        print "In lec_lattice: pCubed = true ------> compute theta and v.";
	theta := NormalBasisElement(OL, psiL);
	while Valuation(theta, P) lt p+2 do
	    theta := p*theta;
	end while;
	v := 0; 
	print "theta = ", theta;
        erz := [OL!(psiL)(g)(theta) : g in G];
        // !!! hier manchmal nicht in OL !?
        M := VerticalJoin( [Vector(ElementToSequence(x)) : x in erz ]);
        ZpGtheta := Lattice(M);
        assert Rank(ZpGtheta) eq Degree(L);
    else
        repeat
            theta, v := lattice_generator_theta(psiL, P, pi);
        
            // erzeuger des Gitters global
            erz := [OL!(psiL)(g)(theta) : g in G];
            // !!! hier manchmal nicht in OL !?
            M := VerticalJoin( [Vector(ElementToSequence(x)) : x in erz ]);
            ZpGtheta := Lattice(M);
        until Rank(ZpGtheta) eq Degree(L);
    end if;
    
    // finde m mit P^m in ZpGtheta
    // einfacher Ansatz:
    k := Index(StandardLattice(#G), ZpGtheta);
    m := Valuation(k*OL, P);
    //m := Valuation(k, setting`p)*RamificationIndex(setting`P);
    
    // kleinstes m
    // schreibe Basis von ZpGtheta in Matrix
    M1 := Matrix(Rationals(), [ElementToSequence(x) : x in erz]);
    M1I := M1^(-1);
    for j in [v+1..m] do
        m := j;
        // schreibe Basis von P^m in Matrix
        M2 := Matrix(Rationals(), [ElementToSequence(x) : x in Basis(P^j)]);
        // Basiswechsel durch T: M1*T=M2
        T := M1I*M2;
        // Elemente in T sollen nach Lokalisierung bei p ganz sein
        if not IsDivisibleBy(Denominator(T),p) then
            break;
        end if;
    end for;
    
    return theta, m;
end function;

/* Computes generator theta of the lattice for the local epsilon constant conjecture
 * by computing a normal basis element and multiplying it by a sufficient power pi^v.
 * Returns theta and v.
 */
function lattice_generator_theta(psi, P, pi)
    local OL, p, theta, v, v1;
    
    OL := Order(P);
    p := Generators(P meet Integers())[1];
    
    theta := NormalBasisElement(OL, psi);
    v := Valuation(theta, P);
    v1 := 1+Floor(RamificationIndex(P)/(p-1));
    v := Maximum(0,v1-v);
    theta := OL!(theta*(pi)^v);
    
    return theta, v;
end function;

/* splitting_module_LPmul_modX(setting::Rec, gamma::Map) -> ModTupRng, SeqEnum
 * Compute the splitting module C(gamma) for the given setting:
 * C=L_P^\times/exp(lattice) and gamma: G \times G -> L_P, then
 *     C(gamma) = B \times C = \langle b_sigma, sigma \neq Id \rangle \times \langle c_0=pi,c_1,\ldots c_n\rangle
 * with G-action
 *     G \times C(gamma) --> C(gamma)
 * The module C(gamma) a sequence of matrices representing the G-action are returned
 * (see [NSW, p.115]).
 */
function splitting_module_LPmul_modX(setting, gamma)
    local G, psiG, GG, piL,
          C, actCMat, phi_LP_C, M,
          g, bilder, sigma, tau, h, v, seq;
    
    G    := setting`G;
    psiG := setting`psiL;
    GG   := [g : g in setting`G];
    //OL   := MaximalOrder(setting`L);
    //pi   := setting`pi;
    piL  := setting`iota(setting`pi);
    
    // Berechne LP^\times/X
    C, actCMat, phi_LP_C := compute_LPmul_modX(setting);
    
    if GetVerbose("epsConj") ge 5 then
        assert &and({  Zero(C) eq C!(actCMat[j]*actCMat[i] - actCMat[Index(GG, GG[i]*GG[j])])[k]
            : i in [1..#G], j in [1..#G], k in [1..Dimension(C)]});
    end if;
    
    // C = <pi, s_1,\ldots s_n>
    // C(gamma) = C + B
    // B = <g, g in G>
    
    // G-Wirkung auf B
    GB := [];
    
    //printf "Berechne G-Wirkung im Splitting Modul";
    for g in GG do
        sigma := g;
        //printf ".";
        M := Vector([0]);
        // Wirkung auf b_tau in NSW :  sigma b_tau = b_{sigma*tau} - b_sigma + gamma(sigma,tau)
        //!!!!
        // Wirkung auf b_tau
        // Rechtswirkung:  sigma b_tau = b_{tau*sigma} - b_sigma + gamma(sigma,tau)
        // Linkswirkung:   sigma b_tau = b_{sigma*tau} - b_sigma + gamma(sigma,tau)
        for tau in [ h : h in GG | h ne Id(G)] do
            //Berechne gamma(g, tau)
            v := gamma(sigma,tau);
            // Rechtswirkung
            // if tau*sigma eq Id(G) then
            // Linkswirkung
            if sigma*tau eq Id(G) then
                // ACHTUNG: HIER ADDITIV ODER MULTIPLIKATIV ?????
                v := v * gamma(1,1);
            end if;
            if g eq Id(G) then
                v := v / gamma(1,1);
            end if;
            // Bild ist
            v := HorizontalJoin(
                   // Linkswirkung:
                   // eine 1 bei sigma*tau, eine -1 bei sigma und 0 bei allen anderen Elementen von G
                   Vector([ h eq sigma*tau select 1 else ( h eq sigma select -1 else 0  )  : h in GG | h ne Id(G)])
                   // Rechtswirkung:
                   // eine 1 bei tau*sigma, eine -1 bei sigma und 0 bei allen anderen Elementen von G
                   // Vector([ h eq tau*sigma select 1 else ( h eq sigma select -1 else 0  )  : h in GG | h ne Id(G)])
                   // und die Koeffizienten von v in C
                   , phi_LP_C(v)
                 );  // !! ZEILEN-VEKTOR !!
            if NumberOfColumns(M) eq 1 and NumberOfRows(M) eq 1 then
                // initialisieren
                M := Matrix(v);
            else
                M := VerticalJoin(M, v);
            end if;
        end for;
        
        Append(~GB, M);
    end for;
    //print "";
    
    // Jetzt ist GM eine Liste von Matrizen
    // jede dieser Matrizen M hat
    //    (#G-1) Zeilen
    //    (#G-1)+dim(C) Spalten
    // und stellt die G-Wirkung auf den Erzeugern von B dar
    
    // actCMat enthaelt die Matrizen fuer die G-Wirkung auf C
    Z := ZeroMatrix(Integers(), Dimension(C), (#G-1));
    GCgamma := [];
    for i in [1..#GG] do
        M := GB[i];
        // M ist (#G-1) x (#G-1)+dim(C) Matrix
        N := actCMat[i];
        // N ist dim(C) x dim(C)-Matrix
        // und wird "unten rechts" in die Matrix M eingefuegt
        // so dass M eine (#G-1)+dim(C) x (#G-1)+dim(C) Matrix wird
        
        // Z ist dim(C) x (#G-1) mit Nullen
        N := HorizontalJoin(Z,N);
        // jetzt hat N (#G-1)+dim(C) Spalten; links Nullen
        M := VerticalJoin(M,N);
        // jetzt ist M quadratisch
        Append(~GCgamma, M);
    end for;
    
    // RSpace Moduli
    seq := [0 : i in [1..#GG-1]] cat Moduli(C);
    Cgamma := RSpaceWithModuli(Integers(), seq);
    Cgamma`moduli := seq;
    
    if GetVerbose("epsConj") ge 5 then
        assert &and({  Zero(Cgamma) eq Cgamma!(GCgamma[j]*GCgamma[i] - GCgamma[Index(GG,GG[i]*GG[j])])[k]
            : i in [1..#G], j in [1..#G], k in [1..Dimension(Cgamma)]});
    end if;
    
    return Cgamma, GCgamma;
end function;


intrinsic LECcomputeLPmulModX(setting::Rec) -> ModTupRng, SeqEnum, Map
{Compute the module C=L_P^\times/X, X=exp(lattice) for the given setting
as well as a sequence of matrices representing the G-action and
a map L_P^\times -> C.}
    
    local V, M, proj;
    V, M, proj := compute_LPmul_modX(setting);
    return V, M, proj;
end intrinsic;

/* compute_LPmul_modX(setting::Rec) -> ModTupRng, SeqEnum, Map
 * Compute the module C=L_P^\times/X, X=exp(lattice), for the given setting
 * and a sequence of matrices representing the G-action and
 * a map L_P^\times -> C.
 */
function compute_LPmul_modX(setting)
    local G, psiG, psiL, GG, OL, pi, piL, iota, P,
          m, Q, pi_OL_Q, Qmal, i_Qmal_Q, phi_OL_Qmal, phi_Qmal_OL,
          expTheta, conjQ, S, pi_Qmal_S, actS, ZgenS, 
          M, k, g, bilder, seq, proj, V;
    
    G    := setting`G;
    psiG := setting`psiL;
    psiL := setting`psiLP;
    GG   := [g : g in setting`G];
    OL   := MaximalOrder(setting`L);
    pi   := setting`pi;
    piL  := setting`iota(setting`pi);
    iota := setting`iota;
    P    := setting`P;
    
    // X = exp(calL)
    // L_P^\times / X = ( L_P^\times / U^m ) / ( exp(L) / U^m )
    
    // ( L_P^\times / U^m ) = pi^\Z \times Q,
    // Q=(O_L / P^m)^\times
    m := setting`prec;
    // Erzeuge Q und Qmal
    Q, pi_OL_Q := quo<OL | P^m>;
    Qmal, i_Qmal_Q := UnitGroup(Q);
    phi_OL_Qmal := pi_OL_Q*(i_Qmal_Q^(-1));
    phi_Qmal_OL := i_Qmal_Q*(pi_OL_Q^(-1));
    
    // exp(calL) in Q wird erzeugt von exp(theta)
    // brauche exp(theta) nur bis zu einer gewissen Genauigkeit
    expTheta := (iota^(-1))(truncated_exp(iota(setting`theta),  truncated_exp_precision(setting) ));
    //expTheta := truncated_exp(setting`theta,  truncated_exp_precision(setting) );
    
    // expTheta und Konjugierte in Q lesen
    conjQ := [ phi_OL_Qmal( psiG(g)(expTheta) ) : g in GG];
    //H := sub<Qmal|conjQ>;
    S, pi_Qmal_S := quo<Qmal | sub<Qmal|conjQ> >;
    
    // Jetzt gilt: L_P^\times / X  =  pi^\Z \times S
    phi_OL_S := phi_OL_Qmal*pi_Qmal_S;
    phi_S_OL := phi_OL_S^(-1);
    
    // G-Wirkung auf S
    actS := map< car<G, S> -> S  |  x :-> phi_OL_S( psiG(x[1])( phi_S_OL(x[2]) ) ) >;
    if GetVerbose("epsConj") ge 5 then
        assert test_G_action(actS);
    end if;
    // Z-Erzeuger
    ZgenS := [S.i : i in [1..#Generators(S)] ];
    
    // G-Wirkung auf L_P^\times / X  als Matrizen
    M := [];
    for k in [1..#GG] do
        g := GG[k];
        bilder := [];
        // Wirkung auf pi lokal
        bild := psiL(g)(piL);
        seq := ElementToSequence( phi_OL_S((iota^(-1))(bild/piL)));
        // Wirkung auf pi global
        //bild := psiG(g)(pi);
        //num,den:=numden(bild/pi, P);
        //seq := ElementToSequence(phi_OL_S(num)-phi_OL_S(den));
        Append(~bilder, [1] cat seq);
        
        bilder := bilder cat [ [0] cat ElementToSequence(actS(g,s) ) : s in ZgenS];
        Append(~M ,  Matrix(Integers(), bilder) );
    end for;
    
    // L_P^\times / X  als RSpace
    V := RSpaceWithModuli(Integers(), [0] cat [Order(g) : g in ZgenS] );
    V`moduli := Moduli(V);
    // G-Wirkung testen
    // G-Wirkung von links, Matrizen operieren von rechts
    vprint epsConj, 5: "G-Action on LP/X: ",
        &and([   V!(M[j]*M[i] - M[Index(GG, GG[i]*GG[j])])[1] eq Zero(V) : i in [1..#G], j in [1..#G] ]);
    if GetVerbose("epsConj") ge 5 then
        assert &and([   V!(M[j]*M[i] - M[Index(GG, GG[i]*GG[j])])[1] eq Zero(V) : i in [1..#G], j in [1..#G] ]);
    end if;
    
    // Projektion (lokale Rechnung)
    proj := map< setting`LP -> V | x :-> 
        Vector([Valuation(x)] cat ElementToSequence(phi_OL_S((iota^(-1))(x/piL^Valuation(x)))))
    >;
    // Projektion (globale Rechnung)
    //print "global projection";
    //proj := map< setting`LP -> V | x :-> projLPV(x, setting, phi_OL_S) >;
    
    return V, M, proj;
end function;

/* compute_LPmul_modX(setting::Rec) -> ModTupRng, SeqEnum, Map
 * Compute the module C=L_P^\times/X, X=exp(lattice), for the given setting
 * and a sequence of matrices representing the G-action and
 * a map L_P^\times -> C.
 */
function compute_LPmul_modX_vals2(L, P, pi, psiG, iota, LP, psiL, theta, m)
    local G, GG, H, OL, piL,
          Q, pi_OL_Q, Qmal, i_Qmal_Q, phi_OL_Qmal, phi_Qmal_OL,
          expTheta, conjQ, S, pi_Qmal_S, actS, ZgenS, 
          M, k, g, bilder, seq, proj, V;
    
    G   := Domain(psiG);  GG := [g : g in G];
    H   := Domain(psiL);  HH := [g : g in H];
    OL  := MaximalOrder(L);
    piL := iota(pi);
    
    // X = exp(calL)
    // L_P^\times / X = ( L_P^\times / U^m ) / ( exp(L) / U^m )
    
    // ( L_P^\times / U^m ) = pi^\Z \times Q,
    // Q=(O_LP / P^m)^\times
    // Erzeuge Q und Qmal
    Q, pi_OL_Q := quo<OL | P^m>;
    Qmal, i_Qmal_Q := UnitGroup(Q);
    phi_OL_Qmal := pi_OL_Q*(i_Qmal_Q^(-1));
    phi_Qmal_OL := i_Qmal_Q*(pi_OL_Q^(-1));
    
    // exp(calL) in Q wird erzeugt von exp(theta)
    // brauche exp(theta) nur bis zu einer gewissen Genauigkeit
    expTheta := (iota^(-1))(truncated_exp(iota(theta),  truncated_exp_precision_vals(theta, m, P, pi) ));
    
    // expTheta und Konjugierte in Q lesen
    conjQ := [ phi_OL_Qmal( psiG(g)(expTheta) ) : g in HH];
    //H := sub<Qmal|conjQ>;
    S, pi_Qmal_S := quo<Qmal | sub<Qmal|conjQ> >;
    
    // Jetzt gilt: L_P^\times / X  =  pi^\Z \times S
    phi_OL_S := phi_OL_Qmal*pi_Qmal_S;
    phi_S_OL := phi_OL_S^(-1);
    
    // G-Wirkung auf S
    actS := map< car<G, S> -> S  |  x :-> phi_OL_S( psiG(x[1])( phi_S_OL(x[2]) ) ) >;
    // Z-Erzeuger
    ZgenS := [S.i : i in [1..#Generators(S)] ];
    
    // G-Wirkung auf L_P^\times / X  als Matrizen
    M := [];
    for k in [1..#HH] do
        g := HH[k];
        bilder := [];
        // Wirkung auf pi lokal
        bild := psiL(g)(piL);
        seq := ElementToSequence( phi_OL_S((iota^(-1))(bild/piL)));
        // Wirkung auf pi global
        //bild := psiG(g)(pi);
        //num,den:=numden(bild/pi, P);
        //seq := ElementToSequence(phi_OL_S(num)-phi_OL_S(den));
        Append(~bilder, [1] cat seq);
        
        bilder := bilder cat [ [0] cat ElementToSequence(actS(g,s) ) : s in ZgenS];
        Append(~M ,  Matrix(Integers(), bilder) );
    end for;
    
    // L_P^\times / X  als RSpace
    V := RSpaceWithModuli(Integers(), [0] cat [Order(g) : g in ZgenS] );
    V`moduli := Moduli(V);
    
    // Projektion (lokale Rechnung)
    proj := map< LP -> V |
        x :-> Vector([Valuation(x)] cat ElementToSequence(phi_OL_S((iota^(-1))(x/piL^Valuation(x))))),
        y :-> piL^yy[1]*iota(phi_S_OL(S!yy[2..#yy]))  where yy := ElementToSequence(y)
    >;
    // Projektion (globale Rechnung)
    //print "global projection";
    //proj := map< setting`LP -> V | x :-> projLPV(x, setting, phi_OL_S) >;
    
    return V, M, proj;
end function;

function compute_LPmul_modX_vals(L, P, pi, psiG, iota, LP, psiL, theta, m)
    local G, GG, H, OL, piL,
          Q, pi_OL_Q, Qmal, i_Qmal_Q, phi_OL_Qmal, phi_Qmal_OL,
          expTheta, conjQ, S, pi_Qmal_S, actS, ZgenS, 
          M, k, g, bilder, seq, proj, V;
    
    G   := Domain(psiG);  GG := [g : g in G];
    H   := Domain(psiL);  HH := [g : g in H];
    OL  := MaximalOrder(L);
    piL := iota(pi);
    
    // X = exp(calL)
    // L_P^\times / X = ( L_P^\times / U^m ) / ( exp(L) / U^m )
    
    // ( L_P^\times / U^m ) = pi^\Z \times Q,
    // Q=(O_LP / P^m)^\times
    // Erzeuge Q und Qmal
    Q, pi_OL_Q := quo<OL | P^m>;
    Qmal, i_Qmal_Q := UnitGroup(Q);
    phi_OL_Qmal := pi_OL_Q*(i_Qmal_Q^(-1));
    phi_Qmal_OL := i_Qmal_Q*(pi_OL_Q^(-1));
    
    // exp(calL) in Q wird erzeugt von exp(theta)
    // brauche exp(theta) nur bis zu einer gewissen Genauigkeit
    expTheta := (iota^(-1))(truncated_exp(iota(theta),  truncated_exp_precision_vals(theta, m, P, pi) ));
    
    // expTheta und Konjugierte in Q lesen
    conjQ := [ phi_OL_Qmal( psiG(g)(expTheta) ) : g in HH];
    //H := sub<Qmal|conjQ>;
    S, pi_Qmal_S := quo<Qmal | sub<Qmal|conjQ> >;
    
    // Jetzt gilt: L_P^\times / X  =  pi^\Z \times S
    phi_OL_S := phi_OL_Qmal*pi_Qmal_S;
    phi_S_OL := phi_OL_S^(-1);
    
    // G-Wirkung auf S
    actS := map< car<G, S> -> S  |  x :-> phi_OL_S( psiG(x[1])( phi_S_OL(x[2]) ) ) >;
    // Z-Erzeuger
    ZgenS := [S.i : i in [1..#Generators(S)] ];
    
    // G-Wirkung auf L_P^\times / X  als Matrizen
    M := [];
    for k in [1..#HH] do
        g := HH[k];
        bilder := [];
        // Wirkung auf pi lokal
        bild := psiL(g)(piL);
        seq := ElementToSequence( phi_OL_S((iota^(-1))(bild/piL)));
        // Wirkung auf pi global
        //bild := psiG(g)(pi);
        //num,den:=numden(bild/pi, P);
        //seq := ElementToSequence(phi_OL_S(num)-phi_OL_S(den));
        Append(~bilder, [1] cat seq);
        
        bilder := bilder cat [ [0] cat ElementToSequence(actS(g,s) ) : s in ZgenS];
        Append(~M ,  Matrix(Integers(), bilder) );
    end for;
    

    // L_P^\times / X  als RSpace
    V := RSpaceWithModuli(Integers(), [0] cat [Order(g) : g in ZgenS] );
    V`moduli := Moduli(V);
    
    // Projektion (lokale Rechnung)
    proj := map< LP -> V |
        x :-> Vector([Valuation(x)] cat ElementToSequence(phi_OL_S((iota^(-1))(x/piL^Valuation(x))))),
        y :-> piL^yy[1]*iota(phi_S_OL(S!yy[2..#yy]))  where yy := ElementToSequence(y)
    >;
    // Projektion (globale Rechnung)
    //print "global projection";
    //proj := map< setting`LP -> V | x :-> projLPV(x, setting, phi_OL_S) >;
    
    // L_P^\times / X  als abelian group 
    Y:=FreeAbelianGroup(#ZgenS+1);
    mmY := map< H -> Aut(Y) | g :-> hom< Y -> Y | 
        y :-> Y!Eltseq(Vector(Eltseq(y))*M[ Index(HH,g)]) > >;
    X, qX := quo<Y | [Order(ZgenS[i])* Y.(i+1) : i in [1..#ZgenS] ]>;
    mmX := map< H -> Aut(X) | g :-> hom< X -> X | x :-> qX( x@@qX @ mmY(g) ) > >; 
    
    // Projektion (lokale Rechnung)
    f := map< LP -> X |
      //x :->  qX(Y!([Valuation(x)] cat Eltseq(phi_OL_S((iota^(-1))(x/piL^Valuation(x)))))) ,
      x :->  qX(Y!([Valuation(x)] cat Eltseq( (x/piL^Valuation(x)) @@ iota @ phi_OL_S ))) ,
      y :->  piL^yy[1]*iota(phi_S_OL(S!yy[2..#yy])) where yy := Eltseq( y @@ qX )
    >;
    
    return X, mmX, f;
end function;


/* proj_resolution_ZGr(Y::ModTupRng, actY::Map) -> SeqEnum, ModTupRng, Map
 * Given a module Y of Type ModTupRng with G-action actY: G\times Y -> Y.
 * Compute a projective resolution of the form
 *      K --> Z[G]^r --> Y
 * Returns Z[G]-generators of Y, the module K and a map Z[G]^r --> Y
 */
function proj_resolution_ZGr(Y, actY)
    local G, GG, g, B, zGB, r, i, j, rows, M, mods, m, K, tors, d, gens;
    
    G := Component(Domain(actY),1);
    GG := [g : g in G];
    // Berechne ZG-Basis
    B := Basis(Y);
    // !!! auch einfache Summen ???
    // A := Basis(Cgamma) cat [a+b : a in Basis(Cgamma), b in Basis(Cgamma) | a ne b];
    zGB := zG_generators(Basis(Y), actY);
    //vprint epsConj, 3: "Z[G]-Basis:", zGB;
    // ZG-Rang
    r := #zGB;
    //vprint epsConj, 3: "Z[G]-Rank:", r;
    
    // Falls Z[G]^r=<z_1,\ldots,z_r>
    // Dann wird folgende Z-Basis von Z[G]^r zugrunde gelegt:
    // < g_1*z_1,\ldots, g_#G*z_1, \ldots\ldots g_1*z_r,\ldots, g_#G*z_r>
    //
    // Abbildung
    //   F=Z[G]^r  --> Y
    //   g_i*z_j   :-> actY(g_i,zGB_j)
    // Bestimme Matrix fuer Abbildung
    rows := [];
    for i in [1..r] do
        for j in [1..#GG] do
            g := GG[j];
            Append(~rows, actY(g, zGB[i]));
        end for;
    end for;
    M := VerticalJoin(rows);
    vprint epsConj, 5: "Abbildung F --> C(gamma):", M;
    
    // Abbildung
    F := RSpace(Integers(), #zGB*#G);
    psi_FY := map< F -> Y | x :-> Y!(x*M)>;
    // G-action
    //actionF := map< car<G,F> -> F | x:-> 
    //             
    //           >;
    
    // Moduli of Y
    tors := NumberOfRows(M);
    mods := Moduli(Y);
    for i in [1..#mods] do
        m := mods[i];
        if m ne 0 then
            M := VerticalJoin(M, Vector(Insert([0 : j in [2..#mods]], i, m)));
        end if;
    end for;
    // Torsion gives more rows
    tors := NumberOfRows(M)-tors;
    vprint epsConj, 3: "Torsion:", tors;
    
    
    // Compute Kernel
    K := Kernel(M);
    
    // we are not interested in the components that generated the torsion
    if tors gt 0 then
        // new degree of K
        d := Degree(K)-tors;
        vprint epsConj, 3: "Neuer Degree vom Kernel:", d;
        //print "Torsion: ", tors;
        //print "New Degree: ", d;
        gens := [];
        for g in Basis(K) do
            Append(~gens, Vector(ElementToSequence(g)[1..d]));
        end for;
        K := RSpaceWithBasis(gens);
    end if;
    
    return zGB, K, psi_FY;
end function;

/* zG_generators(Zgens::SeqEnum, psi::Map) -> SeqEnum
 * Given Z-generators Zgens of H and the G-action psi: (G, H) -> H.
 * Compute Z[G]-generators of H.
 */
function zG_generators(Zgens, psi)
  local G, H, subseqs, seq, S, num, count;
  
  G := Component(Domain(psi),1);
  H := Parent(Zgens[1]);
  for num in [Ceiling(#Zgens/#G)..#Zgens-1] do
    //subseqs := SubSequences(Zgens, num);
    subseqs := [Zgens[SetToSequence(s)] : s in Subsets({1..#Zgens}, num)];
    count := 0;
    for seq in subseqs do
      count := count+1;
      //if count gt 100 then
      //  break;
      //end if;
      S := sub<H | [ psi(g,x)   : x in seq, g in G]>;
      if {z in S : z in Zgens } eq {true} then
        return seq;
      end if;
    end for;
  end for;
  
  return Zgens;
end function;

/* tensor_Zmod_with_Q(R::ModTupRng, actRMat::SeqEnum) -> ModTupFld, SeqEnum
 * Given a Z-module R and a sequence of matrices representing the G-action.
 * Compute R_Q and its G-action.
 */
function tensor_Zmod_with_Q(R, actRMat)
    local moduli, free, V, M, actVMat, MatRing;
    
    moduli := [0 : i in [1..Dimension(R)]];
    if Type(R) eq ModTupRng then
        if assigned R`moduli then
            moduli := R`moduli;
        end if;
    end if;
    free := [[i,moduli[i]] : i in [1..#moduli] | moduli[i] eq 0];
    if free[#free,1] ne #free then
        error "Freier Anteil nicht vorne!";
    end if;
    
    // behalte den freien Anteil
    V := VectorSpace(Rationals(), #free);
    MatRing := MatrixRing(Rationals(), #free);
    actVMat := [ MatRing!Submatrix(M, 1,1, #free, #free) : M in actRMat];
    
    return V, actVMat;
end function;

/* isomorphism_matrix_KQF(K::ModTupRng, CgammaQ::ModTupRng, actionCgammaQ::Map, 
 *     FQ::ModTupRng, actionFQ::Map, zGBasisCgammaQ::SeqEnum) -> Mtrx
 * Given a projective resolution of the splitting module C(gamma)_Q as
 * K_Q --> F=Q[G]^r --> C(gamma)_Q, matrices representing the G-action on
 * C(gamma)_Q and F, and Z[G]-generators of C(gamma)_Q.
 * Compute the Q[G]-isomorphism
 *     K_Q + Q[G] --> K_Q + W_Q + Q --> K_Q + W_Q + C_Q --> K_Q + C(gamma)_Q --> F = Q[G]^r
 * See [Bley/Breuning, 4.2.4] and [Bley/Burns, 3.3].
 */
function isomorphism_matrix_KQF(K, CgammaQ, actionCgammaQ, FQ, actionFQ, zGBasisCgammaQ)
    local G, GG, g, QG, M1, i, 
          WQ, phi_CgammaQ_WQ, actionWQ, v, x, sigma, b, seq, j, M3,
          rho, M4;
    
    // GroupAlgebra
    G := Component(Domain(actionCgammaQ),1);
    GG := [g : g in G];
    //vprint epsConj, 3: "Isomorphismus, G-Aufzaehlung:", GG;

    QG := GroupAlgebra(Rationals(), G);
    
    // Schritt 1:
    // ============
    // Exakte Sequenz:   0 --> W --> Q[G] --aug--> Q --> 0
    // Isomorphismus:    Q[G] --> W_Q + Q
    // Matrix fuer den Isomorphismus
    M1 := HorizontalJoin(
        VerticalJoin(
            // erste Zeile -1/#G
            Vector(Rationals(), [-1/#G : i in [1..#G-1]]),
            // auf Subdiagonale 1-1/#G, sonst -1/#G
            IdentityMatrix(Rationals(), #G-1)+
            Matrix(Rationals(), #G-1,#G-1, [-1/#G : i in [1..(#G-1)*(#G-1)]])
        ),
        // letze Spalte jeweils 1
        Matrix(Rationals(), [[1] : i in [1..#G]])
    );
    
    
    // Schritt 2:
    // ============
    // Isomorphismus:    Q --v^(-1)--> C_Q
    // Matrix ist die Identitaet, da C_Q von pi erzeugt wird
    
    
    // Schritt 3:
    // ============
    // Exakte Sequenz:   C_Q --> C(gamma)_Q --> W_Q
    // Isomorphismus:    W_Q + C_Q --> C(gamma)_Q
    
    // W_Q = <g-1 : g in G>
    WQ := VectorSpace(Rationals(), #G-1);
    phi_CgammaQ_WQ := map< CgammaQ -> WQ | v :-> Submatrix(v, 1,1, 1, #G-1)>;
    // G-Rechtswirkung
    // actionWQ := map< car<G,WQ> -> WQ | x :-> 
    //    WQ!Vector(Remove(ElementToSequence(&+([x[2][i]*(QG!GG[i+1]-QG!Id(G)) : i in [1..#G-1]])*x[1]),1))
    // >;
    // G-Linkswirkung
    actionWQ := map< car<G,WQ> -> WQ | x :-> 
        WQ!Vector(Remove(ElementToSequence(x[1]*&+([x[2][i]*(QG!GG[i+1]-QG!Id(G)) : i in [1..#G-1]])),1))
    >;
    
    // Berechne Schnitt sigma: W_Q --> C(gamma)_Q
    // Q-linear
    // Q[G]-linear
    sigma := make_QG_linear(
        map< WQ -> CgammaQ | w :-> 
             CgammaQ!Vector(ElementToSequence(w) cat [0: i in [1..Dimension(CgammaQ)-Dimension(WQ)]])
           >,
        G, actionWQ, actionCgammaQ);
    
    // Matrix fuer Isomorphismus
    seq := [ sigma(b) :  b in Basis(WQ)];
    // dim(CgammaQ) = #G, dim(W)=#G-1
    // [0,...,0,1] anhaengen, fuer die Einbettung von C_Q in C(gamma)_Q
    Append(~seq, Vector([0 : j in [1..#G-1]] cat [1]));
    M3 := VerticalJoin(seq);
    
    
    // Schritt 4:
    // ============
    // Exakte Sequenz:   K --> F --> C(gamma)_Q
    // Isomorphismus:    K + C(gamma)_Q --> F
    
    // Schnitt rho: C(gamma)_Q --> F
    // Stelle ZG-Basis von C(gamma)_Q bezueglich Z-Basis dar
    vprint epsConj, 5: "Schnitt C(gamma)_Q --> F";
    vprintf epsConj, 5: "Berechne Loesung von M=\n%o\nb=\n%o\n",
        VerticalJoin([actionCgammaQ(g, b) : g in G, b in zGBasisCgammaQ]),
        Basis(CgammaQ);
    
    M := VerticalJoin(Solution(
            VerticalJoin([actionCgammaQ(g, b) : g in G, b in zGBasisCgammaQ]),
            Basis(CgammaQ)
         ));
    // Q[G]-linearer Schnitt
    rhoQ := map< CgammaQ -> FQ | x :-> FQ!(x*M) >;        
    if GetVerbose("epsConj") ge 5 then
        // Ist ein Schnitt?
        F := Matrix([actionCgammaQ(g, b) : g in G, b in zGBasisCgammaQ]);
        f := map< FQ -> CgammaQ | x :-> x*F >;
        assert &and({ f(rhoQ(b)) eq b  : b in Basis(CgammaQ) });
        // ist linear ?
        assert &and({ rhoQ(b1+b2) eq rhoQ(b1)+rhoQ(b2)  : b1 in Basis(CgammaQ), b2 in Basis(CgammaQ) });
        // ist G-vertraeglich ?
        //&and({ rhoQ(actionCgammaQ(g,b)) eq actionFQ(g, rhoQ(b)) : g in G, b in Basis(CgammaQ) });
        // ist f G-vertraeglich ?
        assert &and({ actionCgammaQ(g,f(b)) eq f(actionFQ(g, b)) : g in G, b in Basis(FQ) });
    end if;
    rho := make_QG_linear(rhoQ, G, actionCgammaQ, actionFQ);
    if GetVerbose("epsConj") ge 5 then
        // ist Schnitt ?
        assert &and({ f(rho(b)) eq b  : b in Basis(CgammaQ) });
        // ist linear ?
        assert &and({ rho(b1+b2) eq rho(b1)+rho(b2)  : b1 in Basis(CgammaQ), b2 in Basis(CgammaQ) });
        // ist G-vertraeglich ?
        assert &and({ rho(actionCgammaQ(g,b)) eq actionFQ(g, rho(b)) : g in G, b in Basis(CgammaQ) });
    end if;
    
    // Matrix darstellen
    seq := [];
    for b in Basis(K) do
        Append(~seq, FQ!b);
    end for;
    for b in Basis(CgammaQ) do
        Append(~seq, rho(b));
    end for;
    M4 := VerticalJoin(seq);
    
    // Matrizen multiplizieren
    // M1*M3 um Id auf K ergaenzen
    M := DiagonalJoin(IdentityMatrix(Rationals(), Dimension(K)), M1*M3);
    
    return M*M4;
end function;

/* make_QG_linear(m::Map, G::Grp, actionDom::Map, actionCodom::Map) -> Map
 * Given a Q-linear map m:A->B, a group G and a group action on the domain A
 * and the codomain B of m. Compute a Q[G]-linear map m': A->B.
 */
function make_QG_linear(m, G, actionDom, actionCodom)
    return map< Domain(m) -> Codomain(m) | 
                x :-> 1/#G* &+([  actionCodom(g, m(actionDom(g^(-1),x)))      : g in G])
              >; 
end function;

/* truncated_exp(alpha::., N::RngIntElt) -> .
 * Compute the exponential value of alpha truncated at the N's summand [Bley HS132, Remark 3.6].
 */
function truncated_exp(alpha, N)
    return &+( [  alpha^n/ Factorial(n) : n in [0..N] ] );
end function;

/* truncated_exp_precision(setting::Rec) -> RngIntElt
 * Compute the precision to which the exponential function values must be known to get
 * correct results in L_P/exp(lattice).
 */
function truncated_exp_precision(setting)
    m := setting`prec;
    n := Valuation(setting`theta, setting`P);
    N := Ceiling( m / (n - RamificationIndex(setting`P)/ (setting`p-1)  ) );
    // => N >=  m / (n - e(L/K) / (p-1) )
    return N;
end function;


function truncated_exp_precision_vals(theta, m, P, pi)
    p := Generator(P meet Integers());
    n := Valuation(theta, P);
    N := Ceiling( m / (n - RamificationIndex(P)/ (p-1)  ) );
    // => N >=  m / (n - e(L/K) / (p-1) )
    return N;
end function;



/*******************************
 *  Epsilon Constants
 *******************************/

// Berechne ( tau(L_P/Q_p, chi) )_chi


intrinsic pCubedLECepsilonConstant(~setting::Rec)
{ Compute epsilon constants as described in [Bley/Breuning, 2.5]. }

    local inclCE;

    if not &and({assigned setting`brauerInd, assigned setting`t, assigned setting`linChrs}) then
        pCubedLECprepareEps(~setting);
        //error "call LECprepareEps first";
    end if;
    if assigned setting`E then
        inclCE := complex_compatible_embedding(setting`Qmpt, setting`E, 100, true);
    else
        //inclCE := map< setting`Qmpt -> setting`Qmpt | x :-> x, y :-> y>;
        inclCE := IdentityFieldMorphism(setting`Qmpt);
    end if;

    setting`tLK := pCubed_lec_epsilon_constant(setting`psiL, setting`brauerInd, setting`linChrs, setting`p, setting`t, inclCE);
end intrinsic;



/* This is the same as Ruben's LECprepareEps but uses pCubedBrauerInductionDeg0. */
intrinsic pCubedLECprepareEps(~setting::Rec)
{ Compute Brauer inductions of all irreducible characters and the required
  precision t for the Galois Gauss sums, see [Bley/Breuning, Remark 2.7].
  Uses the special structure of characters of a group of order p^3.}

    local G, IrrG, linChrs, brauerInd, chi, t;

    //vprint epsConj, 1: "Berechne lineare Charaktere aller Untergruppen";
    G := setting`G;
    IrrG := CharacterTable(G);
    setting`IrrG := IrrG;
    //linChrs := linearCharactersWithFields(G, setting`L, setting`psiL);

    vprint epsConj, 2: "Compute Brauer induction in degree 0 for a group of order p^3";
    brauerInd := pCubedBrauerInductionDeg0(setting);

    vprintf epsConj, 2: "Compute abelian extensions";
    vtime   epsConj, 2: linChrs := abelianSituationForLinearCharacters(SequenceToList(&cat(brauerInd)), setting`L, setting`psiL);
    brauerInd := SequenceToList(brauerInd);

    vprintf epsConj, 2: "Precision for Epsilon Constants";
    t := comp_precision(brauerInd, linChrs, setting`p, setting`psiL);
    vprintf epsConj, 2: ": %o\n", t;

    setting`linChrs := linChrs;
    setting`brauerInd := brauerInd;
    setting`t := t;
    C := MinimalCyclotomicField({
                        RootOfUnity(Exponent(setting`G)),
                        RootOfUnity(setting`p^t) });
    if Degree(C) eq 1 then
        C := CyclotomicField(1);
    end if;
    setting`Qmpt := C;
end intrinsic;


intrinsic G63LECprepareEps(~setting::Rec)
{ Compute Brauer inductions of all irreducible characters and the required
  precision t for the Galois Gauss sums, see [Bley/Breuning, Remark 2.7].
  Uses the special structure of characters of a group of order 63.}

    local G, IrrG, linChrs, brauerInd, chi, t;

    //vprint epsConj, 1: "Berechne lineare Charaktere aller Untergruppen";
    G := setting`G;
    IrrG := CharacterTable(G);
    setting`IrrG := IrrG;
    //linChrs := linearCharactersWithFields(G, setting`L, setting`psiL);

    vprint epsConj, 2: "Compute Brauer induction in degree 0 for a group of order p^3";
    brauerInd := G63BrauerInductionDeg0(setting);

    vprintf epsConj, 2: "Compute abelian extensions";
    if #G eq 63 then
        vtime   epsConj, 2: linChrs := G63abelianSituationForLinearCharacters(SequenceToList(&cat(brauerInd)), setting`L, setting`psiL);
    else
        vtime   epsConj, 2: linChrs := abelianSituationForLinearCharacters(SequenceToList(&cat(brauerInd)), setting`L, setting`psiL);
    end if;
    brauerInd := SequenceToList(brauerInd);

    vprintf epsConj, 2: "Precision for Epsilon Constants";
    // t := comp_precision(brauerInd, linChrs, setting`p, setting`psiL);
    t := 1;
    vprintf epsConj, 2: ": %o\n", t;

    setting`linChrs := linChrs;
    setting`brauerInd := brauerInd;
    setting`t := t;
    C := MinimalCyclotomicField({
                        RootOfUnity(Exponent(setting`G)),
                        RootOfUnity(setting`p^t) });
    if Degree(C) eq 1 then
        C := CyclotomicField(1);
    end if;
    setting`Qmpt := C;
end intrinsic;

                             

/* comp_precision(brauerInd::List, linChrs::List, p::RngIntElt, psi::Map) -> RngIntElt
 * Given the Brauer inductions in degree zero, the abelian situations for all
 * linear characters, a prime p and psi:G->Aut(L).
 * Compute the precision for Galois Gauss Sums from the Brauer inductions
 * in degree zero.
 */
function comp_precision(brauerInd, linChrs, p, psi)
    local t, B, t1, Hphi, H, phi, N, M, PN, PM, H2, res, phi2, s, ram, bool;
    
    // berechne maximales t
    t := 0;
    for B in brauerInd do
        t1 := 0;
        for Hphi in B do
            H := Hphi[1];
            phi := Hphi[2];
            
            bool, N,M,_ := abelianSituationForLinearCharacter(linChrs, H, phi);
            assert bool;
            PM := Factorization(p*RingOfIntegers(M))[1,1];
            PN := Factorization(p*RingOfIntegers(N))[1,1];
            // H2 ist Galoisgruppe von N/M
            H2, res := quo<H | Kernel(phi)>;
            // phi2=phi ist Charater von H2
            phi2 := reduceCharacter(phi, H2, res);
            // s = Bewertung des Fuehrers
            
            ram := nontriv_ramGrps(PN, M, res, psi);
            s:=conductor(phi2, ram);
            t1 := Maximum(t1,Ceiling(s/ RamificationIndex(PM) ) );
        end for;
        t := Maximum(t, t1);
    end for;
    
    return t;
end function;

/* nontriv_ramGrps(PN::RngOrdIdl, M::Fld, res::Map, psi::Map) -> SeqEnum
 * Computations of non trivial ramification groups of PN in the following situation:
 * Let L/N/M be Galos extensions with H=\Gal(L/M), H_1 =\Gal(L/N), H_2 =\Gal(N/M),
 * res is the quotient map res: H -> H_2=H/H_1, and psi: H -> Aut(L/M). Then
 * the non-trivial ramification groups G_i<H_2 of PN are computed by using
 * the automorphisms from psi.
 */
function nontriv_ramGrps(PN, M, res, psi)
    local RamGroups, i, done, H, H2, ON, b;
    
    H := Domain(res);
    H2 := Codomain(res);
    ON := Order(PN);
    
    RamGroups := [];
    done := false;
    i := 0;
    while done eq false do
        Gi := ramGrp(PN, M, res, psi, i);
        /* Gi := sub< H2 | { res(sigma) : sigma in H |  
                            { Valuation(psi(sigma)(b)-b, PN) ge i+1 : b in Basis(ON)} eq {true}  
                        } >;
        */
        Append(~RamGroups, Gi);
        if #Gi eq 1 then
            done := true;
        end if;
        i := i+1;
    end while;
    
    return RamGroups;
end function;


function ramGrp(PN, M, res, psi, i)
    local G, I, B, g, is_in_Gi, j, a;

    H := Domain(res);
    H2 := Codomain(res);
    Igens := [One(H2)];
    I := sub<H2 | Igens >;
    B := Basis(Order(PN));
    cnt := 1;
    for g in H2 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(psi(g@(res^-1))(a) - a, PN) le i then
                    is_in_Gi := false;
                end if;
                j := j+1;
            end while;
            if is_in_Gi then
                Append(~Igens, g);
                I := sub<H2 | Igens >;
            end if;
            cnt := cnt+1;
        end if;
    end for;
    return I;
end function;


function lec_epsilon_constant(psi, brauerInd, linChrs, p,t, inclCE : pCubed:=false)
    
    local C, E, eps, B, tau, Hphi, H, phi;
        
    vprintf epsConj, 2: "Precision %o\n", t;
    vprintf epsConj, 2: "%o Characters\n", #brauerInd;
    
    C := Domain(inclCE);
    E := Codomain(inclCE);
    
    eps := [**];
    for i in [1..#brauerInd] do
        print "------------------------> lec_epsilon_constant: i = ", i;
        B := brauerInd[i];
        vprint epsConj, 2: "Compute Galois-Gaus sum for character"; // , IrrG[i];
        IndentPush();
        vprintf epsConj, 3: "%o linear characteres\n", #B;
        
        tau := E!1;
        cnt := 0;
        for Hphi in B do
            cnt := cnt+1; 
            // fuer jedes Paar (H, phi)
            H := Hphi[1];
            phi := Hphi[2];
            vprintf epsConj, 3: "character %o  ^  %o\n", phi, Hphi[3];
            IndentPush();
            //tau := tau*( ( galois_gauss_sum(setting, H, phi, inclCE) )^Hphi[3] );
            if pCubed and #Domain(psi) eq 27 then
                tau := tau*( ( pCubed_galois_gauss_sum(psi, linChrs, p, t, inclCE, H, phi) )^Hphi[3] );
            elif #Domain(psi) eq 63 then
                tau := tau*( ( G63_galois_gauss_sum(psi, linChrs, p, t, inclCE, H, phi) )^Hphi[3] );
            else   
                tau := tau*( ( galois_gauss_sum(psi, linChrs, p, t, inclCE, H, phi) )^Hphi[3] );
            end if;
            IndentPop();
        end for;
        Append(~eps, tau);
        IndentPop();
    end for;
    vprintf epsConj, 3: "tLK = %o\n", [* (inclCE^(-1))(x) : x in eps *];
    return eps;
end function;


/* This function is essentially the same as Ruben's lec_epsilon_constant with two  differences: 
   (a) We use pCubed_galois_gauss_sum instead of galois_gauss_sum.
   (b) At the end we compute \tau(\Qp, psi_2(chi) - a*chi)), where psi_2 denotes the second Adams operation.

   This function is used to compute Jacobi sums and the contribution of Gauss sums in the definition of \fra_{N/\Qp).
*/
function pCubed_lec_epsilon_constant(psi, brauerInd, linChrs, p,t, inclCE, IrrG, a)
    local C, E, eps, B, tau, Hphi, H, phi;

    vprintf epsConj, 2: "Precision %o\n", t;
    vprintf epsConj, 2: "%o Characters\n", #brauerInd;

    C := Domain(inclCE);
    E := Codomain(inclCE);

    gauss_sums := [**];
    for i in [1..#brauerInd] do
        B := brauerInd[i];
        vprint epsConj, 2: "Compute Galois-Gaus sum for character"; // , IrrG[i];
        IndentPush();
        vprintf epsConj, 3: "%o linear characteres\n", #B;

        tau := E!1;
        for Hphi in B do
            // fuer jedes Paar (H, phi)
            H := Hphi[1];
            phi := Hphi[2];
            vprintf epsConj, 3: "character %o  ^  %o\n", phi, Hphi[3];
            IndentPush();
            tau := tau*( ( pCubed_galois_gauss_sum(psi, linChrs, p, t, inclCE, H, phi) )^Hphi[3] );
            IndentPop();
        end for;
        Append(~gauss_sums, tau);
        IndentPop();
    end for;
    
    // print "gauss_sums = ", gauss_sums;
  
    /* Compute \tau(\Qp, psi_2(chi) - a*chi)) */
    eps := [**];
    for j:=1 to #IrrG do
        chi := IrrG[j];
        i:= [i : i in [1..#IrrG] | AdamsOperation(chi,2) eq IrrG[i]][1];
        eps[j] := gauss_sums[i]/gauss_sums[j]^a;
    end for;

    vprintf epsConj, 3: "tLK = %o\n", [* (inclCE^(-1))(x) : x in eps *];
    return eps;
end function;

/* This function is essentially the same as Ruben's lec_epsilon_constant with two  differences: 
   (a) We use G63_galois_gauss_sum instead of galois_gauss_sum.
   (b) At the end we compute \tau(\Qp, psi_2(chi) - a*chi)), where psi_2 denotes the second Adams operation.

   This function could also be used to compute Jacobi sums and the contribution of Gauss sums in the definition of \fra_{N/\Qp).
   However, the computation of the Gauss sums takes a long time. One should compute these Gauss sums and save them forever
   for any future computation and write a separate function which computes the part described in (b) above.
*/

function G63_lec_epsilon_constant(psi, brauerInd, linChrs, p,t, inclCE, IrrG, a)
    local C, E, eps, B, tau, Hphi, H, phi;

    vprintf epsConj, 2: "Precision %o\n", t;
    vprintf epsConj, 2: "%o Characters\n", #brauerInd;

    C := Domain(inclCE);
    E := Codomain(inclCE);

    gauss_sums := [**];
    for i in [1..#brauerInd] do
    // for i in [12] do
        print "i = ", i;
        B := brauerInd[i];
        vprint epsConj, 2: "Compute Galois-Gaus sum for character"; // , IrrG[i];
        IndentPush();
        vprintf epsConj, 3: "%o linear characteres\n", #B;

        tau := E!1;
        for Hphi in B do
            // fuer jedes Paar (H, phi)
            H := Hphi[1];
            phi := Hphi[2];
            vprintf epsConj, 3: "character %o  ^  %o\n", phi, Hphi[3];
            IndentPush();
            tau := tau*( ( G63_galois_gauss_sum(psi, linChrs, p, t, inclCE, H, phi) )^Hphi[3] );
            IndentPop();
        end for;
        Append(~gauss_sums, tau);
        IndentPop();
    end for;
    
    eps := [**];
    for j:=1 to #IrrG do
        chi := IrrG[j];
        i:= [i : i in [1..#IrrG] | AdamsOperation(chi,2) eq IrrG[i]][1];
        eps[j] := gauss_sums[i]/gauss_sums[j]^a;
    end for;

    vprintf epsConj, 3: "tLK = %o\n", [* (inclCE^(-1))(x) : x in eps *];
    return eps;
end function;



/* galois_gauss_sum(setting::Rec, H::Grp, phi::AlgChtrElt, E::FldNum) -> FldNumElt
 */
function galois_gauss_sum(psi, linChrs, p, t, inclCE, H, phi) //(setting, H, phi, inclCE)
    local C, E, bool, N, M, Na, PM, PN, H2, res, phi2, psi2,
          s, c, reps, StdChar, taux, x, a, psiLM;
    
    C := Domain(inclCE);
    E := Codomain(inclCE);
    
    // Koepersituation
    bool, N, M, Na := abelianSituationForLinearCharacter(linChrs, H, phi);
    assert bool;
    
    // Primideale
    PM  := Factorization(p*RingOfIntegers(M))[1,1];
    PN  := Factorization(p*RingOfIntegers(N))[1,1];
    // Charakter auf Faktorgruppe reduzieren
    kerPhi := Kernel(phi);
    H2, res := quo<H|kerPhi>;
    
    phi2 := reduceCharacter(phi, H2, res);
    // Abbildung G--> Aut einschraenken
    psi2 := restrict_aut_to_quotient(psi, N, H2, res);
    
    // Fuehrer berechnen
    vprintf epsConj, 4: "Compute conductor\n";
    IndentPush();
    s := conductor(phi2, nontriv_ramGrps(PN, M, res, psi));
    if M eq Rationals() then
        if s eq 0 then
            c := MaximalOrder(M)!1;
        else
            c := Generator(PM^s);
        end if;
    else
        c := PrimitiveElement(PM^Valuation(Different(MaximalOrder(M))*PM^s,PM));
    end if;
    vprintf epsConj, 4: "conductor s=%o, generator c=%o\n", s, c;
    IndentPop();
    
    // Representanten von (OM/PM^s)^\times
    reps := representants(PM, s);
    // Additiver Charakter (Werte in E)
    StdChar := standard_add_char(PM, t, C, E);
    
    // fuer jeden Representanten Artin-Symbol und dann
    // Summand in Galois-Gauss-Summe berechnen
    vprintf epsConj, 4: "%o representant(s) ", #reps;

    
    taux := E!0;
    /* The local symbol for c is computed only once to save time.*/
    b := localNormResidueSymbol(M!(c), Na, PM, psi2 : NoArchPlaces:=IsOdd(#H2));
    for x in reps do
        vprintf epsConj, 4:  ".";
        a := localNormResidueSymbol(M!(x), Na, PM, psi2 : NoArchPlaces:=IsOdd(#H2)) / b;
        taux := taux + inclCE(phi2(a))*StdChar(M!(x/c));
    end for;
    vprintf epsConj, 4: "\n";
    
    return taux;
end function;

/* This is essentially Ruben's galois_gauss_sum.
   In the final loop we compute the Artin symbol of c only once. This saves time.
*/
function pCubed_galois_gauss_sum(psi, linChrs, p, t, inclCE, H, phi) //(setting, H, phi, inclCE)
    local C, E, bool, N, M, Na, PM, PN, H2, res, phi2, psi2,
          s, c, reps, StdChar, taux, x, a, psiLM;

    C := Domain(inclCE);
    E := Codomain(inclCE);

    // Koepersituation
    bool, N, M, Na := abelianSituationForLinearCharacter(linChrs, H, phi);
    assert bool;

    OM := MaximalOrder(M);

    // Primideale
    PM  := Factorization(p*RingOfIntegers(M))[1,1];
    PN  := Factorization(p*RingOfIntegers(N))[1,1];
    // Charakter auf Faktorgruppe reduzieren
    kerPhi := Kernel(phi);
    H2, res := quo<H|kerPhi>;

    phi2 := reduceCharacter(phi, H2, res);
    // Abbildung G--> Aut einschraenken
    psi2 := restrict_aut_to_quotient(psi, N, H2, res);

    // Fuehrer berechnen
    vprintf epsConj, 4: "Compute conductor\n";
    IndentPush();
    s := conductor(phi2, nontriv_ramGrps(PN, M, res, psi));
    if M eq Rationals() then
        if s eq 0 then
            c := MaximalOrder(M)!1;
        else
            c := Generator(PM^s);
        end if;
    else
       c := PrimitiveElement(PM^Valuation(Different(MaximalOrder(M))*PM^s,PM));
    end if;
    vprintf epsConj, 4: "conductor s=%o, generator c=%o\n", s, c;
    IndentPop();

    reps := representants(PM, s);
    
    // Additiver Charakter (Werte in E)
    StdChar := standard_add_char(PM, t, C, E);

    // fuer jeden Representanten Artin-Symbol und dann
    // Summand in Galois-Gauss-Summe berechnen
    vprintf epsConj, 4: "%o representant(s) ", #reps;

    taux := E!0;
    b := localNormResidueSymbol(M!(c), Na, PM, psi2 : NoArchPlaces:=IsOdd(#H2));
    for x in reps do
        vprintf epsConj, 4:  ".";
        // a := localNormResidueSymbol(M!(x/c), Na, PM, psi2);
        a := localNormResidueSymbol(M!(x), Na, PM, psi2 : NoArchPlaces:=IsOdd(#H2)) / b;
        taux := taux + inclCE(phi2(a))*StdChar(M!(x/c));
    end for;
    vprintf epsConj, 4: "\n";

    return taux;
end function;


/* This is essentially Ruben's galois_gauss_sum or pCubed_galois_gauss_sum.

   Certain improvements of performance are implemented so to speed up the computation.
   Essentially we use a function FastLocalNormResidueSymbol which, in addition to compute the local norm
   residue symbol, also returns computed data such as the values g(b) for b in Basis(ON) and g in Gal(N/M) and also a list of already
   computed Frobenius elements. This data is then used in later calls of FastLocalNormResidueSymbol.

   A clean and probably much faster way to compute the local norm residue symbol is to implement one of the follwoing methods
       (a) The Dokchitser's approach to compute (conjugacy classes of) Frobenius elements.
       (b) Use MAGMA's built-in reciprocity map. Probably one has to normalize this map so that it is compatible with everything else.
*/

function G63_galois_gauss_sum(psi, linChrs, p, t, inclCE, H, phi) //(setting, H, phi, inclCE)
    local C, E, bool, N, M, Na, PM, PN, H2, res, phi2, psi2,
          s, c, reps, StdChar, taux, x, a, psiLM;

    C := Domain(inclCE);
    E := Codomain(inclCE);

    // Koepersituation
    bool, N, M, Na := abelianSituationForLinearCharacter(linChrs, H, phi);
    assert bool;

    OM := MaximalOrder(M);

    // Primideale
    PM  := Factorization(p*RingOfIntegers(M))[1,1];
    PN  := Factorization(p*RingOfIntegers(N))[1,1];
    // Charakter auf Faktorgruppe reduzieren
    kerPhi := Kernel(phi);
    H2, res := quo<H|kerPhi>;

    phi2 := reduceCharacter(phi, H2, res);
    // Abbildung G--> Aut einschraenken
    psi2 := restrict_aut_to_quotient(psi, N, H2, res);

    // Fuehrer berechnen
    vprintf epsConj, 4: "Compute conductor\n";
    IndentPush();
    s := conductor(phi2, nontriv_ramGrps(PN, M, res, psi));
    if M eq Rationals() then
        if s eq 0 then
            c := MaximalOrder(M)!1;
        else
            c := Generator(PM^s);
        end if;
    else
       c := PrimitiveElement(PM^Valuation(Different(MaximalOrder(M))*PM^s,PM));
    end if;
    vprintf epsConj, 4: "conductor s=%o, generator c=%o\n", s, c;
    IndentPop();

    reps := representants(PM, s);

    // Additiver Charakter (Werte in E)
    StdChar := standard_add_char(PM, t, C, E);

    // fuer jeden Representanten Artin-Symbol und dann
    // Summand in Galois-Gauss-Summe berechnen
    vprintf epsConj, 4: "%o representant(s) ", #reps;
    
    H := Domain(psi2);
    N := Domain(Codomain(psi2)[1]);
    /* We use a multiple of the conductor: by construction we know that at most the primes above 2,7,11,13 ramify and only 7 ramifies wildly (and weakly).*/
    // F := Factorization(Conductor(Na));
    F := G63Factorization(M);    
    B := Basis(MaximalOrder(N));
    index_set := [<h, b> : h in H, b in B];
    precomp_data := [N!0 : x in index_set];
    data := [* B, index_set, precomp_data, [], [**] *];

    taux := E!0;
    
    b, data := FastLocalNormResidueSymbol(M!(c), F, PM, psi2, data);
    for x in reps do
        // print "x = ", x;
        vprintf epsConj, 4:  ".";
        a1, data := FastLocalNormResidueSymbol(M!(x), F, PM, psi2, data ); a := a1 / b;
        taux := taux + inclCE(phi2(a))*StdChar(M!(x/c));
    end for;
    vprintf epsConj, 4: "\n";

    return taux;
end function;

intrinsic G63Factorization(M :: FldRat) -> SeqEnum
{}
    OM := MaximalOrder(M);
    P := Factorization(7*OM)[1,1];
    fac := Factorization(2*11*13*OM);
    return [<P, 2>] cat [<f[1], 1> : f in fac];
end intrinsic;

intrinsic G63Factorization(M :: FldNum) -> SeqEnum
{}
    OM := MaximalOrder(M);
    P := Factorization(7*OM)[1,1];
    fac := Factorization(2*11*13*OM);
    return [<P, 2>] cat [<f[1], 1> : f in fac];
end intrinsic;



/* representants(PM::., s::RngIntElt) -> SeqEnum
 * Computes representatives of the multiplicative group (OM/PM^s)^\times.
 */
function representants(PM, s)
    if Type(PM) eq RngInt then
        OM := Integers();
    else
        OM := Order(PM);
    end if;

    Q, m1     := quo<OM | PM^s>;
    Qmult, m2 := MultiplicativeGroup(Q);
    if #Qmult eq 1 then
        return [1];
    else
        return [ (m1^(-1))(m2(q)) : q in Qmult];
    end if;
end function;



/* standard_add_char(PM::RngOrdIdl/RngInt, t::RngIntElt, C::FldCyc, E::FldNum) -> Map
 * Compute the standard additive Character for \zeta_(p^t) in E, where C is
 * a cyclotomic extension s.t. \zeta_(p^t)\in C\subset E (see [Bley/Breuning, 2.5]).
 */
function standard_add_char(PM, t, C, E)
    local M, p, MC, iota, Zp, Qp, zeta, addChar;
    
    if Type(PM) eq RngInt then
        M := Rationals();
        p := Generator(PM);
    else
        M := NumberField(Order(PM));
        p := Generator(PM meet Integers());
    end if;
    
    MC, iota := Completion(M, PM : Precision:=60);
    Zp := pAdicRing(p);
    Qp := pAdicField(p);
    zeta := root_of_unity(p^t, C, E);
    
    if Type(PM) eq RngInt then
        addChar := map<M -> E | z :->  zeta^(
            [a : a in [0..p^t-1] | (z - a/p^t) in Zp ][1]) >;
    else
        addChar := map<M -> E | z :->  zeta^(
            [a : a in [0..p^t-1] | Trace(iota(z), Qp) - a/p^t in Zp ][1] ) >;
    end if;
    
    return addChar;
end function;


/* root_of_unity(n::RngIntElt, C::FldCyc, E::FldNum) -> FldNumElt
 * Given an integer n, a cyclotomic field C and a number field E,
 * s.t. zeta_n\subset C\subset E.
 * Compute the n-th root of unity zeta_n in E, whose first complex
 * embedding is exp(2 pi i/n).
 */
function root_of_unity(n, C, E)
    CC := ComplexField();
    zeta := E!C!RootOfUnity(n,C);
    expWert := Exp(2*Pi(CC)*CC.1/n);
    //expWert := ChangePrecision(C!expWert, Precision(expWert)-1);
    //print expWert;
    
    for k in [1..n] do
        //print [k, C!Conjugates(zeta^k)[1] ];
        if Abs(CC!Conjugates(zeta^k)[1] - expWert) lt 10^(-Precision(CC)+1) then
            return zeta^k;
        end if;
    end for;
    
    error "Error while computing root of unity.";
    return E!0;
end function;

/* Computes the twisted unramified characteristic frc_{E/Qp} as an element of the centre. 
 * This should work for an arbitrary extension. 
 */
function pCubed_lec_unram_char(setting)
    
    QG := setting`Krel`QG;
    G := setting`G;
    p := setting`p;
    OL := setting`OL;
    h := setting`psiL;

    P := Factorization(p*OL)[1,1];
    I := RamificationGroup(P, 0, h);
    eI := SubgroupIdempotent(QG, I);
    eG := SubgroupIdempotent(QG, G);
    Frob := QG ! Frobenius(OL, h, p) ;
    eins := One(QG);
    // eta := (eins - eI) + Frob*eI; /* This is the "old" twisted unramified characteristic. */
    eta := (eins - eI) + Frob^-1 * eI;

    incl_id := Coercion(setting`Qmpt, setting`Qmpt);
    nr := NewtonReducedNormEG(eta, incl_id);
    return nr;
end function;






/**************************************************
 *         Relative algebraic K-Groups
 **************************************************/

/* K0RelGrp wie im RelAlgKTheory Paket.
 * Allerdings wird statt QG allgemeiner EG verwendet.
 */
K0RelGrp := recformat<DT : GrpAb,
                      PrimeIdeals : List,
                      PrimeElts : List,
                      OCmodG : GrpAb,
                      m : List,
                      Gp : List,
                      f : Map,
                      p : RngIntElt,
                      EG : AlgGrp
                    >;

/* Input: K0RelGrp-Record K0(Z_p[G], Q_p) und eine bei p voll zerlegte
 * Erweiterung E/Q.
 *
 * Berechnet ein K0RelGrp-Record, welches das Rechnen in K0(Z_p[G], E_Q)
 * erlaubt, sowie einen Isomorphismus zwischen den beiden relativen K-Gruppen.
 *
 * Ziel ist es, ein Element (alpha_chi)\in \prod EK_i als Element in 
 * K0(Z_p[G], Q_p) auffassen zu koennen.
 */
intrinsic extendGlobalField(Krel::Rec, P::RngOrdIdl : Precision := 100) -> Rec
{}
    local OE, G, EG, H;
    
    require Krel`p eq Generators(P meet Integers())[1] : 
            "Prime Ideal does not divide prime";
    
    // Ein Primideal P in E ueber p waehlen
    OE := Order(P);
    E  := NumberField(OE);
    
    // Charakterkoerper E(chi) bilden und Primideale ueber P
    // und Primelemente bestimmen
    charFields  := [* *];
    PrimeIdeals := [* *];
    PrimeElts   := [* *];
    Gp := [* *];
    H := [* *];
    for i in [1..#Krel`QG`H] do
        Qchi := Krel`QG`H[i,1];
        // Konstuiere Echi
        f := DefiningPolynomial(Qchi);
        if Type(f) eq SeqEnum then
            assert #f eq 1;
            f := f[1];
        end if;
        
        Echi := ext< E| f>;
        Echi := AbsoluteField(Echi);
        
        OEchi := RingOfIntegers(Echi);
        
        // Einbettung Qchi->Echi
        iota := complex_compatible_embedding(Qchi, Echi, Precision, false);
        incl := complex_compatible_embedding(E, Echi, Precision, false);
        
        // Primideal(e) ueber P finden
        PEchi := [ v[1] : v in Factorization(ideal< OEchi | { incl(x) : x in Generators(P) } >)];
        
        // Berechne alle Primideale ueber p
        ideals := [ v[1] : v in Factorization(Krel`p * OEchi) ];
        // nur die, die auch ueber P liegen
        ideals := [ v : v in ideals | v in PEchi ];
        
        // Berechne Ideal ueber gp
        gp := ideal<OEchi | { iota(x) : x in Generators(Krel`Gp[i]) } >;
        // aber nur den Anteil ueber P
        factors := Factorization(gp);
        if #factors ne 0 then
            gp := &*( [ v[1]^v[2] : v in factors  | v[1] in PEchi  ]  );
            //gp := factors[1,1]^factors[1,2];
        else
            // Qchi = Q
            // => gp = OEchi
        end if;
        
        // Primelemente
        primes := Uniformizer(ideals, gp);
        
        Append(~charFields,  Echi);
        Append(~PrimeIdeals, ideals);
        Append(~PrimeElts,   primes);
        Append(~Gp, gp);
        Append(~H, [* Echi, Krel`QG`H[i,2], iota, OEchi *]);
    end for;
    
    
    EG := GroupAlgebra(E, Group(Krel`QG));
    //K1AmodF := Oumf(QG, ZG, F, Ideals,p);
    OCmodG, m, pIdeals := residueClassGroup(charFields, Gp, Krel`p);
    //OCmodG, m, pIdeals := ResidueClassGroup(QG, Ideals,p);
    //DT:= TorsionSubgroup(QG, OCmodG, m, K1AmodF);
    
    KrelE := rec<K0RelGrp |>;
    // KrelE`DT
    KrelE`PrimeIdeals := PrimeIdeals;
    KrelE`PrimeElts := PrimeElts;
    KrelE`OCmodG := Krel`OCmodG;
    KrelE`m := m;
    KrelE`Gp := Gp;
    KrelE`f := Krel`f;   // DT`f
    KrelE`p := Krel`p;
    KrelE`EG := EG;
    KrelE`EG`H := H;
    
    return KrelE;
end intrinsic;

intrinsic residueClassGroup(charFields::List, Ideals::List, p::RngIntElt) -> GrpAb, List, List
{ Computes the residue class group as ResidueClassGroup from RelAlgKTheory for
  character fields different from Q(chi). }
    local i, ordA, m, QR, A, f, pIdeals;

    ordA := [];
    m := [* *];
    pIdeals := [* *];

    for i:=1 to #Ideals do
        Id := pPrimaryPart(Ideals[i], p);
        QR := quo< MaximalOrder(charFields[i]) | Id >;
        A, f := MultiplicativeGroup(QR);
        Append(~m, f^-1);
        Append(~pIdeals, Id);
        if NumberOfGenerators(A) eq 0 then
            Append(~ordA, [ 1 ]);
        else
            Append(~ordA, [Order(A.j) : j in [1..NumberOfGenerators(A)]]);
        end if;
    end for;
    //print &cat(ordA);
    OCmodG := AbelianGroup( &cat(ordA) );

    return OCmodG, m, pIdeals;
end intrinsic;

intrinsic reducedNorm(K::Rec, TOmega::Rec) -> List
{ Compute the reduced norm of an element of the relative K-group.
  See first part of K0RelLog from the RelAlgKTheory package. }
    local QG, BasisA, BasisB, S, nr;

    QG := K`QG;

    BasisA := LocalBasis(TOmega`A, K`p);
    BasisB := LocalBasis(TOmega`B, K`p);

    S :=  QGMatrix(QG, TOmega`Theta, TOmega`A, BasisA, TOmega`B, BasisB);
    S :=  Transpose(S);
    /* nr:= ReducedNorm(S); */
    nr:= NewtonReducedNorm(S);
    return nr;
end intrinsic;

intrinsic mapToRelKGroup(Krel::Rec, nr::List) -> .
{ Given a tupel representing an element in a relative K-group,
  e.g. from reducedNorm. Map the tupel of norms to the relative K-group.
  See second part of K0RelLog from the RelAlgKTheory package. }
    local IdealPart, QG, i, j, e, row, v; 
    
    if "EG" in Names(Krel) then
        EG := Krel`EG;
    else
        EG := Krel`QG;
    end if;
    
    IdealPart := [* 1*h[4] : h in EG`H *];
    for i:=1 to #nr do
        for j:=1 to #Krel`PrimeIdeals[i] do
            e := Valuation(nr[i], Krel`PrimeIdeals[i][j]);
            IdealPart[i] := IdealPart[i]*Krel`PrimeIdeals[i][j]^e;
            nr[i] := nr[i]*Krel`PrimeElts[i][j]^(-e);
        end for;
    end for;
    
    row := [];
    for i:=1 to #nr do
        //print "numden";
        beta, gamma := numden(nr[i], Krel`PrimeIdeals[i,1]);
        //print "eltseq";
        v := ElementToSequence( (Krel`m[i])(beta) - ((Krel`m[i])(gamma)) );
        if #v eq 0 then
            v := [0];
        end if;
        row cat:= v;
    end for;
    
    return [* IdealPart, Krel`f(Krel`OCmodG!row) *];
end intrinsic;

/* numden(xi::FldNumElt, P::RngOrdIdl) -> RngOrdElt, RngOrdElt
 * Copy of numden from the RelAlgKTheory package
 * allowing xi to be an element of the field, not only the order.
 */
function numden(xi, P)
    local OK, F, IdA, IdB, f, b, beta, gamma, p, j;

    assert Valuation(xi, P) eq 0;
    OK := Order(P);

    /* Compute the ideal numerator and denominator of xi*OK. Could be replaced by IdealNumDen. */
    F := Factorization(xi*OK);
    IdA := 1*OK;
    IdB := 1*OK;
    for f in F do
        if f[2] lt 0 then
            IdB := IdB * f[1]^(-f[2]);
        end if;
        if f[2] gt 0 then
            IdA := IdA * f[1]^f[2];
        end if;
    end for;

    b, beta := IsPrincipal(IdA);
    if b then
        /* Ideal numerator and denominator are already principal. */
        gamma := beta / xi;
        return beta, gamma;
    end if;

    /* Find a prime Q such that IdA*Q is principal. This could be improved a lot. This was
       probably never tested, since in most applications in our context the rings of integers
       have class number 1. 
    */
    p := 2;
    found := false;
    while not found do
        F := Factorization(p*OK);
        j := 1;
        while not found and j le #F do
            if F[j][1] ne P then
                b, beta := IsPrincipal(F[j][1]*IdA);
                if b then
                    found := true;
                    gamma := beta / xi;
                end if;
             end if;
             j := j+1;
        end while;
        p :=NextPrime(p);
     end while;
     return beta, gamma;
end function;

/* Gegeben lambda in L[G]\subset\C[G].
 * Berechne die reduzierte Norm
 *      nr(lambda) \in \prod_{\chi\in Irr(G)} \C
 * [Bley/Wilson, 3.3]
 * Hierbei sei C ein zyklischer Koerper, so dass chi(g)\in C fuer alle chi und g,
 * und E sei eine Erweiterung vom Kompositum LC, so dass nr(lambda)_i \in E
 */
intrinsic NewtonReducedNormEG(lambda :: AlgGrpElt, incl::Map) -> List
{}
    local G, IrrG, N, h, chi, n, sigma, s, k;
    
    G := Group(Parent(lambda));
    IrrG := CharacterTable(G);
    //print IrrG[2];
    
    N := [* *];
    for chi in IrrG do
        n := Integers() ! Degree(chi);
        // Fortsetzung von chi auf E[G] auf lambda anwenden
        sigma := [lin_ext_char(chi, lambda, incl)];
        s := [ sigma[1] ];
        for k:=2 to n do
            Append(~s, lin_ext_char(chi, lambda^k, incl));
            Append(~sigma, (-1)^(k+1) * (s[k] + &+[ (-1)^l * s[k-l]*sigma[l] : l in [1..k-1] ]) / k);
        end for;
        Append(~N, sigma[n] );
   end for;
   return N;
end intrinsic;

/* lin_ext_char(chi::AlgChtrElt, lambda::AlgGrpElt, incl::Map) -> FldCycElt
 * Applies the character chi on an element lambda of the group ring R[G]
 * by lifting the character-map to R via the inclusion incl.
 * See LinExtChar from the RelAlgKTheory package.
 */
function lin_ext_char(chi, lambda, incl)
    local G, g;
    G := Group(Parent(lambda));
    return &+[Coefficient(lambda, g) * incl(Domain(incl)!chi(g)) : g in G];
end function;

intrinsic '+'(X::List, Y::List) -> List
{ Compute the sum of to elements in K0(Zp[G], Qp). }
    local Z;
    
    require #X eq 2 and #Y eq 2 and
            Type(X[2]) eq GrpAbElt and Type(Y[2]) eq GrpAbElt and
            Type(X[1]) eq List and Type(Y[1]) eq List and
            &and({Type(X[1,i]) in {RngOrdIdl, RngOrdFracIdl} : i in [1..#X[1]]}) and
            &and({Type(Y[1,i]) in {RngOrdIdl,RngOrdFracIdl} : i in [1..#Y[1]]}) :
            "Erwartet: Listen X und Y, die Elemente in K0(Z_p[G], Q_p) repraesentieren!";
            
    Z := [* [* X[1,i]*Y[1,i] : i in [1..#X[1]]  *], X[2]+Y[2]   *];
    return Z;
end intrinsic;

intrinsic '-'(X::List, Y::List) -> List
{ Compute the sum of to elements in K0(Zp[G], Qp). }
    local Z;

    require #X eq 2 and #Y eq 2 and
            Type(X[2]) eq GrpAbElt and Type(Y[2]) eq GrpAbElt and
            Type(X[1]) eq List and Type(Y[1]) eq List and
            &and({Type(X[1,i]) in {RngOrdIdl, RngOrdFracIdl} : i in [1..#X[1]]}) and
            &and({Type(Y[1,i]) in {RngOrdIdl,RngOrdFracIdl} : i in [1..#Y[1]]}) :
            "Erwartet: Listen X und Y, die Elemente in K0(Z_p[G], Q_p) repraesentieren!";

    Z := [* [* X[1,i]/Y[1,i] : i in [1..#X[1]]  *], X[2]-Y[2]   *];
    return Z;
end intrinsic;


intrinsic RelativeGroupC(G::Grp, p::RngIntElt) -> Rec
{ Compute the relative K-group and make its embeddings compatible
  with the complex embedding. }
    local Krel, H;
    
    Krel := RelativeGroup(G,p);
    
    for i in [1..#Krel`QG`H] do
        H := Krel`QG`H[i];
        iota := complex_compatible_embedding(Codomain(H[3]), Domain(H[3]), 100, false);
        
        Krel`QG`H[i,3] := iota^(-1);
    end for;
    
    return Krel;
end intrinsic;





/**************************************************
 *                Number Fields
 **************************************************/


/* Let L and E be number fields such that L can be embedded into E.
 * Compute the embedding of L into E which is compatible to the embeddings
 * of L and E into C.
 */
function complex_compatible_embedding(L, E,  precision, coercion)
    local CC, B, rows, b, rts, r, M, Id, T, i;
    
    CC := ComplexField(precision);
    assert #Generators(L) eq 1;
    b := L.1;
    if coercion then
        rts := Roots(MinimalPolynomial(b), L);
        rts := [E!r[1] : r in rts];
    else
        rts := Roots(MinimalPolynomial(b), E);
        rts := [r[1] : r in rts];
    end if;
    rts := [r : r in rts |  Abs(CC!Conjugates(r)[1]-CC!Conjugates(b)[1]) lt 10^(-precision+1) ];
    r := rts[1];
    
    return embedding(L,E,r);
end function;

/*
 * Given two absolute number fields L \subset E and a root r of the
 * defining polynomial of L in E.
 * Create the embedding map L --> E, L.1 :-> r and its inverse.
 */
function embedding(L, E, r)
//intrinsic constructInvertibleEmbedding(L::FldNum, E::FldNum, r::FldNumElt) -> Map
    if not (BaseField(L) eq Rationals() and BaseField(E) eq Rationals()) then
        error "Absolute fields required!";
    end if;
    
    rows := [];
    for i in [0..Degree(L)-1] do
        Append(~rows, ElementToSequence(r^i) );
    end for;
    M := Matrix(rows);
    Id := IdentityMatrix(Rationals(), NumberOfRows(M) );
    T := Transpose(Matrix(Solution(Transpose(M), [Id[i]: i in  [1..NumberOfRows(Id)]])));
    assert M*T eq Id;
    
    return hom< L -> E | x :-> E!ElementToSequence(Vector(ElementToSequence(x))*M), y :-> L!inverse_embedding(M, T, y)   >;
end function;

/*
 * Let M be the matrix with rows corresponding to the embeddings of the basis of L
 * into E. And let T be such that M T=Id, i.e. T represents the inverse of the
 * embedding as far as possible.
 * If y lies in the subspace generated by the rows of M, y T is computed.
 * Otherwise an error is thrown.
 */
function inverse_embedding(M, T, y)
    row  := Vector(ElementToSequence(y));
    if Rank(M) eq Rank(VerticalJoin(M, Matrix(row))) then
        return ElementToSequence(row*T);
    else
        error "Coercion not possible";
    end if;
end function;

/* root_with_complex_embedding(f::RngUPolElt, K::FldNum, x::FldNumElt, precision::RngIntElt) -> FldNumElt
 * Compute the zero of f in K, whose first complex embedding is closest to the
 * first complex embedding of x.
 */
function root_with_complex_embedding(f, K, x, precision)    
    local CC, y, rts, z, m, i;
    
    y := Conjugates(x)[1];
    CC := Parent(y);
    rts := Roots(f, K);
    rts := [r[1] : r in rts |  Abs(CC!Conjugates(r[1])[1]-y) lt 10^(-precision+1) ];
    if #rts eq 1 then
        z := rts[1];
    else
        m,i := Minimum([ Abs(CC!Conjugates(r)[1]-y) : r in rts]);
        z := rts[i];
    end if;
    
    return z;
end function;

/* ramification_group(P::RngOrdIdl, psi::Map, s::RngIntElt) -> Grp
 * For an ideal P of L, psi: Gal(L)-> Aut(L), compute a higher ramification group of P.
 * If the base field of L is Q, this uses Magma's internal function RamificationGroup.
 */
function ramification_group(P, psi, s)
    local OL, R;
    
    if Type(BaseRing(Order(P))) eq RngInt then
        return RamificationGroup(P,s);
    end if;
    
    // v(sigma(a)-a) >= s+1
    // bzw. sigma(a)-a = 0 in OK/P^(s+1)
    
    OL := Order(P);
    R := quo<OL | P^(s+1) >;
    G := Domain(psi);
    
    return sub<G | [ h : h in G | &and([R!(psi(h)(x)-x) eq Zero(R) : x in Basis(OL)])  ] >;
end function;

/* inertia_group(P::RngOrdIdl, psi::Map) -> Grp
 * For an ideal P of L, psi: Gal(L)-> Aut(L), compute the inertia group of P.
 * If the base field of L is Q, this uses Magma's internal function RamificationGroup.
 */
function inertia_group(P, psi)
    return ramification_group(P, psi, 0);
end function;

/* fixed_field(psi::Map, H::Grp, subfields::SeqEnum) -> FldNum
 * Returns the subfield of Domain(psi) which ist fixed by H.
 * If not given, all subfields are computed first.
 */
function fixed_field(psi, H, subfields)
    local G;
    
    G := Domain(psi);
    L := Domain(psi(Id(G)));
    
    if #H eq #G then
        return BaseField(L);
    else
        if #subfields eq 0 then
            subfields := Subfields(L);
        end if;
        return [ K[1] : K in subfields | 
                    Degree(K[1]) eq Degree(L)/#H and  
                    &and([  psi(sigma)(K[2](x)) eq K[2](x) : x in  Generators(K[1]), sigma in H ])
               ][1];
    end if;
end function;

/* frobenius_grp_elt(L::FldNum/FldRat, PL::RngOrdIdl/RngInt, G::Grp, psi::Map) -> GrpElt
 * Given a number field extension L/K with Galois group G, the
 * automorphism map psi:G->Aut(L/K) and a prime ideal PL.
 * Compute an element in G representing the Frobenius 
 * automorphism of PL.
 */
function frobenius_grp_elt(L, PL, G, psi)
    local OL, OK, PK, B, b, g, f, q, Q, s1;
    
    if Type(L) eq FldRat then
        assert #G eq 1;
        return Id(G);
    end if;
    
    OL := Order(PL);
    B := Basis(OL);
    
    OK := RingOfIntegers(BaseField(L));
    PK := PL meet OK;
    
    Q, s1 := quo<OL | PL>;
    
    if OK eq Integers() then
        q := #ResidueClassField(Integers(), PK);
    else
        q := Index(OK,PK);
    end if;
    // Suche g mit
    // Valuation(g(b)-b^q) > 0 
    for g in G do
        for b in B do
            if s1(psi(g)(b)-b^q) ne Zero(Q) then
                continue g;
            end if;
        end for;
        //durchgekommen
        f := g;
    end for;
    
    return f;    
end function;


/* restrict_aut_to_quotient(psi::Map, K::FldNum/FldRat, H::Grp, res::Map) -> .
 * Given the automorphism map psi:G->Aut(L), K\subset L and a factor H
 * of U<G with canonical projection res:U-->H.
 * Compute psi|_H:H->Aut(K).
 */
function restrict_aut_to_quotient(psi, K, H, res)
    local Hs, Auts;
    Hs := [h : h in H];
    if Type(K) eq FldRat then
        assert(#H eq 1);
        Auts := [ map< K -> K | x :-> x> ] ;
    else
        Auts := [ restricted_aut(psi((res^(-1))(h)),K) : h in Hs];
    end if;
    return map<H -> Auts | x:-> Auts[Index(Hs,x)] >;
end function;

/* Given an automorphism sigma of L and K\subset L.
 * Compute sigma|_K.
 */
function restricted_aut(sigma, K)
    //if not IsSubfield(K,Codomain(sigma)) then
    //    error "K must be a subfield of codomain of sigma";
    //end if;
    return hom<K-> K | x:-> sigma(x), y:-> (sigma^(-1))(y)>; 
end function;










/**************************************************
 *                Local Fields
 **************************************************/

/* completion_with_prec(L::FldNum, P::RngOrdIdl, psi::Map : Precision := 30) -> FldPad, Map, Map
 * Given a number field L, a prime ideal P and the automorphism map psi:G->Aut(L).
 * Computes the completion L_P,iota:=Completion(L,P) with given precision and
 * the local automorphism map psi_{L_P}:G->Aut(L_P) such that 
 *       psi(g)(b)= iota^(-1) ( psi_{L_P}(g)(iota(b)) )
 * holds for any generator b of L, every g in G.
 * Returns L_P, iota and psi_{L_P}.
 */
function completion_with_prec(L, P, psi, precision)
    local prec, min, err, compatible;
    
    if Generator(P meet Integers()) eq 2 then
        prec := 3*Maximum(precision,30);
    else
        prec := 2*Maximum(precision,30);
    end if;
    
    repeat
        err := false;
        compatible := false;
        try
            //print "completion precision ", prec;
            vprint epsConj, 5: "compute completion", prec;
            LP, iota := Completion(L, P : Precision:=prec);
            autLP := Automorphisms(LP, pAdicField(LP));
            _, psiLP := localized_automorphism_group(psi, P, iota, autLP);
            
            vprint epsConj, 5: "test compatibility";
            compatible := test_G_compatible(iota, psi, psiLP, true, 0, precision);
            if compatible then
                compatible := test_G_compatible((iota)^(-1), psiLP, psi, false, P, precision);
            end if;
            //time min := Minimum(test_G_compatible(iota, psi, psiLP, true, 0)
            //    join test_G_compatible((iota)^(-1), psiLP, psi, false, P));
            //compatible := (min ge precision);
            
        catch e
            //print e`Object;
            err := true;
        end try;
        
        prec := 2*prec;
        if err then
            continue;
        end if;
    until (not err) and compatible;
    
    return LP, iota, psiLP;
end function;

/* localized_automorphism_group(m::Map, P::RngOrdIdl, iota::Map, AutLoc::SeqEnum) -> GrpPerm, Map
 * Given an automorphism map m:G->Aut(L), a prime ideal P with inclusion
 * iota: L -> L_P, and a list Aut(L_P) of automorphisms of L_P.
 * Compute the local Galois group H and the automorphism map m':H->Aut(L_P).
 */
function localized_automorphism_group(m, P, iota, AutLoc)
    local G, H, f, L, OL, Rts, RtsLok, i,j,prec,index, z, y, S;
    
    G := Domain(m);
    // Untergruppe von G
    //H := DecompositionGroup(P);
    H := [g : g in G | &and([  m(g)(x) in P   : x in Generators(P)]) ];
    H := sub< G | H >;
    
    // Wenn G und H gleich sind, kann es sein, dass Magma die Gruppen
    // unterschiedlich aufzaehlt. D.h.
    // G eq H liefert true und 
    // [g : g in G] eq [g : g in H] liefert false
    
    // dieses Verhalten ist nicht ganz nachvollziehbar und wird
    // hiermit umgangen
    if G eq H then
        H := G;
    end if;
    
    L := Domain(m(H.1));
    OL := Domain(AutLoc[1]);
    Rts := Roots(DefiningPolynomial(L), L);
    //RtsLok := Roots(ChangePrecision(Polynomial(OL, ElementToSequence(f)),OL`DefaultPrecision));
    RtsLok := Roots(DefiningPolynomial(L), OL);
    
    assert #Rts eq #RtsLok;
    
    // Zuordnung:    globale Nst x <-> lokale Nst y
    z := [];
    for i in [1..#Rts] do
        S := [ Valuation(RtsLok[j,1] - OL!iota(Rts[i,1])) : j in [1..#RtsLok] ];
        prec, index := Maximum(S);
        z := z cat [index];
    end for;
    //print z;
    
    // Zuordnung:    g in AutLoc <-> index von g(RtsLok[1]) in RtsLok
    y := [];
    for i in [1..#AutLoc] do
        S := [ Valuation(AutLoc[i](RtsLok[1,1]) - RtsLok[j,1] ) : j in [1..#RtsLok] ];
        //print S;
        prec, index := Maximum(S);
        y := y cat [index];
    end for;
    //print y;
    
    // Zuordnung:    Index globale Nst x  <->  Index von g in AutLoc, so dass g(RtsLok[1])=y
    z := [ Index(y, z[i]) : i in [1..#z] ];
    
    return H, map< H -> AutLoc | x :-> local_map(x, m, AutLoc, Rts, z) >;
end function;

/* Computes the image of g for the map from G to the automorphisms
 * of the local field.
 */
function local_map(g, m, HomL, Rts, z)
//localMap(g::., m::Map, HomL::SeqEnum, Rts::SeqEnum, z::SeqEnum) -> Map
    // Nehme die globale Nst x0, die auf die erste lokale Nst abbildet
    first := Index(z,1);
    x := m(g)(Rts[first,1]);
    // Finde Index der globalen Nst x, so dass g(x0)=x
    S := [ x- Rts[i,1] : i in [1..#Rts] ];
    j := Index(S, 0);
    // Der Index der lokalen Abb, die das gleiche tut, steht in z
    return HomL[z[j]];
end function;

/* Computes the coefficients of the element x in L, according to the basis
 * of L/K.
 */
function elem_to_seq(x, K)
    if RingOfIntegers(Parent(x)) cmpeq RingOfIntegers(K) then
        return [x];
    end if;
    y := ElementToSequence(x);
    while not RingOfIntegers(Parent(y[1])) eq RingOfIntegers(K) do
        y := &cat([Coefficients(y[j]) : j in [1..#y]]);
    end while;
    return y;
end function;





/**************************************************
 *                   Tests
 **************************************************/


/* test_G_compatible(phi::Map, actDom::Map, actCodom::Map : modP := false, prime := 0) -> BoolElt/Set
 * Given a map phi and group actions actDom:G \times D -> D and actCodom: H \times C -> C
 * on its domain D and codomain C.
 * Test whether phi is compatible with group actions.
 * The groups G and H must either be equal or one of them is a subgroup of the other one.
 * If modP is false, just true/false is returned.
 * Otherwise, the algorithm works locally and returns valuations at the given prime.
 */
function test_G_compatible(phi, actDom, actCodom, modP, prime, min)
    local D, B, gens, actD, actB, seq, U;

   // print "In test_G_compatible with modP = ", modP;
    
    if Type(prime) eq RngOrdIdl then
        modP := true;
    end if;
    
    D := Domain(phi);
    B := Codomain(phi);
    
    if Type(Domain(actDom)) ne SetCart then
        G := Domain(actDom);
        actD := map< car<G, D> -> D | x :-> actDom(x[1])(x[2]) >;
    else
        G := Component(Domain(actDom),1);
        actD := actDom;
    end if;
    
    if Type(Domain(actCodom)) ne SetCart then
        H := Domain(actCodom);
        //assert G eq Domain(actCodom);
        actB := map< car<H, B> -> B | x :-> actCodom(x[1])(x[2]) >;
    else
        H := Component(Domain(actCodom),1);
        //assert G eq Component(Domain(actCodom),1);
        actB := actCodom;
    end if;
    
    if G eq H then
        // groups equal
        U := G;
    else
        // take the smaller group
        if #H lt #G then
            U := H;
        else
            U := G;
        end if;
        // and make sure the elements can be read in the other group
        assert &and([x in G and x in H :x in U]);
    end if;
    
    
    if Type(D) in {RngOrd, FldNum, ModTupRng, ModTupFld} then
        gens := Basis(D);
    elif Type(D) in {FldPad, RngPad} then
        gens := basis_pad(D, pAdicField(D));
    else
        print "not yet implemented: Generators/Basis for ", Type(D);
        try
            gens := Basis(D);
        catch e
            gens := Generators(D);
        end try;
    end if;
    
    if Type(B) in {FldPad,RngPad} then
        val := func< x | Valuation(x) >;
    elif Type(B) in {FldNum,RngOrd} then
        if Type(prime) ne RngOrdIdl then
            error "Prime Ideal for Valuation needed!";
        end if;
        val := func< x | Valuation(x, prime) >;
    end if;
    
    if modP then
        vals := {};
        for x in gens, sig in U do
            s := phi(actD(sig, x)) - actB(sig, phi(x));
            v := val(s);
            vals := vals join {v};
            if min eq 0 and v lt min then
                break;
            end if;
        end for;
        if min eq 0 then
            return vals;
        else
            return &and({v ge min : v in vals });
        end if;
        /*seq := [ phi(actD(sig, x)) - actB(sig, phi(x)) : x in gens, sig in Generators(U)];
        if Type(B) in {FldPad,RngPad} then
            return {Valuation(x) : x in seq};
        elif Type(B) in {FldNum,RngOrd} then
            if Type(prime) ne RngOrdIdl then
                error "Prime Ideal for Valuation needed!";
            end if;
            return {Valuation(x, prime) : x in seq};
        else
            error "not yet implemented: Valuation";
        end if;*/
    else
        //seq := [ [* sig, x, B!(phi(actD(sig, x)) - actB(sig, phi(x))) *] : x in gens, sig in G ];
        //print seq;
        return &and({ phi(actD(sig, x)) eq actB(sig, phi(x)) : x in gens, sig in U});
    end if;
end function;

/*
 * Computes a basis of L/K from the products of powers of the primitive
 * elements of the subextensions.
 * If L has basis (1, a, a^2) over E and E has basis (1, b) over K,
 * the output is (1,b,  a,a b,  a^2,a^2 b).
 */
function basis_pad(L, K)
    if (BaseRing(L) eq K) then
        n := Degree(L, K);
        return [L.1^(i-1) : i in [1..n]];
    else
        n := Degree(L); // ueber BaseRing
        B := basis_pad(BaseRing(L), K);
        return &cat([ [ L.1^(i-1)*b : b in B] : i in [1..n]   ]);
    end if;
end function;

/* test_G_action(m::Map) -> BoolElt
 * Test whether the given map m:G\times A -> A is a G-action.
 */
function test_G_action(m)
    local G, B, C, b, sig, tau;
    
    G := Component(Domain(m),1);
    C := Component(Domain(m),2);
    
    if Type(C) in {GrpAb} then
        B := Generators(C);
    elif Type(C) in {FldNum, RngOrd, ModTupRng} then
        B := Basis(C);
    else
        try
            B := Basis(C);
        catch e
            B := Generators(C);
        end try;
    end if;
    return &and({ m(sig*tau,b) eq m(sig,m(tau,b))  : b in B, sig in G, tau in G});
end function;



/* test_cocycle(gamma::Map, actCodom::Map)
 * Given a map gamma and a G-action actCodom: G\times C -> C on the
 * codomain C of gamma.
 * Test whether gamma is a (multiplicative) cocyle and return either
 * true/false or a set of valuations.
 */
function test_cocycle(gamma, actCodom)
    local G, actB, seq, B;
    
    G := Component(Domain(gamma),1);
    B := Codomain(gamma);
    if Type(Domain(actCodom)) ne SetCart then
        assert G eq Domain(actCodom);
        actB := map< car<G, B> -> B | x :-> actCodom(x[1])(x[2]) >;
    else
        actB := actCodom;
    end if;
    
    if Type(B) in {RngPad, FldPad} then
        return { Valuation(gamma(sig*tau, rho)*gamma(sig,tau)-
                 actB(sig,gamma(tau,rho))*gamma(sig,tau*rho))
                : sig in G, tau in G, rho in G};
        
    elif Type(B) in {ModTupRng} then
        seq := [ gamma(sig*tau, rho)+gamma(sig,tau) -
               ( actB(sig,gamma(tau,rho))+gamma(sig,tau*rho) )
               : sig in G, tau in G, rho in G];
        
        return &and([ B!ElementToSequence(x) eq Zero(B) : x in seq ]);
    else
        seq := [ gamma(sig*tau, rho)*gamma(sig,tau) -
               ( actB(sig,gamma(tau,rho))*gamma(sig,tau*rho) )
               : sig in G, tau in G, rho in G];
        
        return &and([ x eq Zero(B) : x in seq ]);
    end if;
end function;

intrinsic QmptIdempotents(setting :: Rec) -> SeqEnum
{Computes the Q-rational idempotents}

    local G, K, i, chi, psi, Kchi, i1, i2, v, e, E;

    idem := [];
    G := setting`G;
    QmptG := GroupAlgebra(setting`Qmpt, G);

    Qmpt := setting`Qmpt;
    for chi in setting`IrrG do
        v := [Qmpt ! chi(g^-1) : g in G];
        v := [Degree(chi) / #G * a : a in v];
        e := QmptG ! v;
        Append(~idem, e);
    end for;
    return idem;
end intrinsic;



intrinsic ProjectionToClassgroup(omega :: List, Krel :: Rec, cl :: Rec) -> GrpAbElt
{}

    local QG, rho, A, index, ps, f, TOmega, p, K, t, off, gamma, i, g, Oki, CoD, D, w, eta, y, IdB;

    // print "From ProjectionToClassgroup";

    K := Krel;
    QG := K`QG;
    p := K`p;
    log := cl`cl ! 0;
    t := ElementToSequence( ((K`f)^-1)(omega[2]) );
    off := 0;
    IdB := [* *];
    for i:=1 to #QG`H do
        // print "i = ", i;
        /* Compute gamma_i as in Bley/Wilson, Section 5. */
        g := K`m[i];
        OKi := Order(cl`Ideals[i]);
        CoD := Codomain(g);
        D := Domain(g);
        w := [ t[off + i] : i in [1..Ngens(CoD)] ];
        off := off + Ngens(CoD);

        eta := OKi ! (g^-1)(CoD!w);
        if eta eq (OKi ! 0) then
            eta := (OKi ! 1);
            assert #CoD eq 1;
        end if;
        y := &*[ K`PrimeElts[i][j]^Valuation(omega[1][i], K`PrimeIdeals[i][j]) : j in [1..#K`PrimeIdeals[i]] ];
        gamma := y * eta;

        /* Set up data to do the weak approximation step. */
        fac := Factorization(cl`Ideals[i]);
        alpha := [ OKi ! 1 : f in fac ];
        for j:=1 to #fac do
            if pPrimaryPart(fac[j][1], p) eq fac[j][1] then
                alpha[j] := gamma;
            end if;
        end for;
        beta := Approximation(alpha, fac, cl`Ideals[i],  cl`InfinitePlaces[i]);

        /* Finally compute the ideal whose class in the ray classgroup represents [P_p,id, Z_pG^rk(P)]. */
        Append(~IdB, pPrimaryPart(gamma*OKi, p) * beta);
    end for;

    /* Now read IdB (which comes componentwise) as an element of the locally free class group. */
    cl_elt := [];
    for i:=1 to #IdB do
        NumB, DenB := IdealNumDen(IdB[i]);
        v := ElementToSequence( (cl`m[i])(NumB) - (cl`m[i])(DenB) );
        if #v eq 0 then
            v := [0];
        end if;
        cl_elt cat:= v;
    end for;

    log := cl`f (cl`rcgp ! cl_elt);
    return log;
end intrinsic;



intrinsic Reduce(alpha :: RngOrdElt, N :: RngIntElt) -> RngOrdElt
{}
    coeff := [Integers() ! c : c in Eltseq(alpha)];
    coeff := [c mod N : c in coeff];
    beta := Parent(alpha) ! coeff;
    return beta;
end intrinsic;



intrinsic RelativeEmbeddingInverse(t::FldNumElt, M::FldNum, E::FldNum, Qmpt::FldCyc, zeta::FldNumElt) -> FldCycElt
{}
    local d, c, A, b, row, i, k, j, y, alpha;

    d := [Eltseq(b) : b in Eltseq(t)];
    c := [[Eltseq(a) : a in Eltseq(zeta^j)] : j in [0..Degree(Qmpt)-1]];

    A := ZeroMatrix(Rationals(), Degree(M)*Degree(E), Degree(Qmpt));
    b := ZeroMatrix(Rationals(), Degree(M)*Degree(E), 1);

    row := 0;
    for i:=0 to Degree(E)-1 do
        for k:=0 to Degree(M)-1 do
	    row := row+1;
	    for j:=0 to Degree(Qmpt)-1 do
	        A[row, j+1] := c[j+1,i+1,k+1];
	    end for;
	    b[row, 1] := d[i+1, k+1];
	end for;
    end for;
    y := Solution(Transpose(A), Transpose(b));
    alpha := &+[y[1,i+1]*Qmpt.1^i : i in [0..Degree(Qmpt)-1]];
    return alpha;
end intrinsic;
