/*
 * brauer.m
 * (c) 2008-2011 Ruben Debeerst
 *
 * This file includes some algorithms written by W. Bley
 * (see bottom of file). They were included to get a
 * Magma-package which works independently.
 *
 * Methods to compute with local and global Brauer groups.
 *
 * Methods provided:
 *
 * LocalBrauerGroup
 *     Compute a module cohomologically isomorphic to L_v^\times
 *     and the corresponding cohomology groups. Also compute
 *     the local fundamental class if required.
 *
 * FrobeniusAutomorphism
 * FrobeniusEquation
 *     Implementation of the constructive proof of
 *     Neukirch, V, Lemma 2.1.
 *
 * LocalFundamentalClassDirect
 * LocalFundamentalClassSerre
 *     Methods to compute the local fundamental class as described
 *     in ...
 *
 *
 * GlobalCocycleInvariants
 * GlobalCocycle
 *     Methods to compute in global Brauer groups as described
 *     in Chapter ...
 * 
 * RightToLeftAction
 * LeftToRightAction
 * RightToLeftCocycle
 * LeftToRightCocycle
 *     Turn a right G-action into a left action and vice versa.
 *     Similarly, for cocycles.
 *
 */

// verbose comments
declare verbose Brauer, 2;
declare verbose CocycleLFC, 2;
//SetVerbose("CocycleLFC", 0);
declare verbose FrobEq, 3;
//SetVerbose("FrobEq", 1);




// forward local functions
forward local_brauer_group,
        completion_with_prec, lattice, compute_LPmul_modX;
forward lattice_generator_theta, lattice_absolute, lattice_relative;


// for Frobenius Automorphism
forward frobenius_automorphism, frobenius_automorphism_grpelt,
        continue_frobenius, continued_automorphism_image;
// for Frobenius Equation
forward frobenius_equation, frobenius_equation_OK, frobenius_equation_cont;
// for Local Fundamental Class
forward cocycle_lfc, cocycle_lfc_G, restrictAutomorphism,
        generates_maximal_unram, galois_act_L_tensor_Knr,
        apply_map, is_zero, is_equal, find_power, precompute_map, 
        tupelQuotient, tupelProduct, tupelInverse,
        continuations_with_unram_restriction;
// to compute cocycle invariant
forward cocycle_invariant;
// for local fields
forward automorphism_group, generators_pad, field_tower,
        elem_to_seq, primitive_element, maximal_unram_subext_simple;
forward local_map, localized_automorphism_group;
forward test_G_compatible, basis_pad;

forward local_fundamental_class_direct;

forward global_cocycle_invariants, global_cocycle;



// local Brauer group
locBrGrp := recformat<
    L    : FldNum,        // the corresponding global field
    P    : RngOrdIdl,     // with prime ideal
    p    : RngIntElt,
    //m    : RngIntElt      // precision needed
    
    M    : GrpAb,        // the module
    actM : Map,          // group action on M
    qM   : Map,          // projection map L^\times ->> M
    theta: RngOrdElt,
    
    C    : ModCoho,      // local cohomology group structure
                         // computed for M and the subgroup H
    f1   : Map,          // corresponding map of M into the local structure
    lfc  : ModTupRngElt  // local fundamental class
>;

declare attributes FldNum: localBrauerGroups;


intrinsic LocalBrauerGroup(L::FldNum, p::RngIntElt : autMap := 0, lfc := false) -> Rec
{ Compute the local cohomology group at p: H^2(G_P,L_P^\times).
Returns a record containing:
  the cohomology structure C,
  the module M=L_P^\times/X s.t. H^2(G_P,L_P^\times)=H^2(G_P,M),
  the group action on M,
  a map qM:L->M,
  a map f1 from M onto the internal module,
  and the representant of the local fundamental class (if required).
}
    P := Factorization(p*RingOfIntegers(L))[1,1];
    return LocalBrauerGroup(L,P : autMap := autMap, lfc := lfc);
end intrinsic;

intrinsic LocalBrauerGroup(L::FldNum, P::RngOrdIdl : autMap := 0, lfc := false) -> Rec
{ Compute the local cohomology group at p: H^2(G_P,L_P^\times).
Returns a record containing:
  the cohomology structure C,
  the module M=L_P^\times/X s.t. H^2(G_P,L_P^\times)=H^2(G_P,M),
  the group action on M,
  a map qM:L->M,
  a map f1 from M onto the internal module,
  and the representant of the local fundamental class (if required).
}

    require IsAbsoluteField(L) and IsPrime(P) :
            "Absolute field L and prime ideal P required.";
    require IsNormal(L):
            "The field must be normal.";
    
    if not assigned L`localBrauerGroups then
        L`localBrauerGroups := [**];
    end if;
    p := Generator(P meet Integers());
    i := 0;
    for locBr in L`localBrauerGroups do
        i := i+1;
        if locBr`p cmpeq p then
            if not lfc or assigned locBr`lfc then
                return locBr;
            else
                // local fundamental class has to be computed
                // remove old record
                Remove(~L`localBrauerGroups, i);
            end if;
        end if;
    end for;
    
    if autMap cmpeq 0 then
        _,_,autMap := AutomorphismGroup(L);
        autMap := map< Domain(autMap) -> Codomain(autMap) | x :-> autMap(x^(-1)) >;
    end if;
    locBr := local_brauer_group(L,P,autMap, lfc);
    
    Append(~L`localBrauerGroups, locBr);
    return locBr;
end intrinsic;


intrinsic LocalFundamentalClassDirect(L::FldPad, n::RngIntElt) -> Map
{Compute the cocyle representing the local fundamental class of L/Qp
up to the given precision.}
    return local_fundamental_class_direct(L, n);
end intrinsic;


intrinsic LocalFundamentalClassSerre(L::FldPad, K::FldPad, steps::RngIntElt : psi:= 0, NoNorm := false) -> Map
{Compute the cocyle representing the local fundamental class of L/K
up to the given precision. Optionally, one can pass the map psi:G->Aut(L/K)
computed by AutomorphismGroup(L,K). If the parameter NoNorm=true, 
the algorithm without norm equations is used.}
    return cocycle_lfc_G(RingOfIntegers(L),RingOfIntegers(K),psi,steps,NoNorm);
end intrinsic;

intrinsic LocalFundamentalClassSerre(L::RngPad, K::RngPad, steps::RngIntElt : psi:= 0, NoNorm := false) -> Map
{Compute the cocyle representing the local fundamental class of L/K
up to the given precision. Optionally, one can pass the map psi:G->Aut(L/K)
computed by AutomorphismGroup(L,K). If the parameter NoNorm=true, 
the algorithm without norm equations is used.}
    return cocycle_lfc_G(L,K,psi,steps,NoNorm);
end intrinsic;

intrinsic CocycleLFC(L::FldPad, K::FldPad, steps::RngIntElt : psi:= 0, NoNorm := false) -> Map
{Compute the cocyle representing the local fundamental class of L/K
up to the given precision. Optionally, one can pass the map psi:G->Aut(L/K)
computed by AutomorphismGroup(L,K). If the parameter NoNorm=true, 
the algorithm without norm equations is used.}
    return cocycle_lfc_G(RingOfIntegers(L),RingOfIntegers(K),psi,steps,NoNorm);
end intrinsic;

intrinsic CocycleLFC(L::RngPad, K::RngPad, steps::RngIntElt : psi:= 0, NoNorm := false) -> Map
{Compute the cocyle representing the local fundamental class of L/K
up to the given precision. Optionally, one can pass the map psi:G->Aut(L/K)
computed by AutomorphismGroup(L,K). If the parameter NoNorm=true, 
the algorithm without norm equations is used.}
    return cocycle_lfc_G(L,K,psi,steps,NoNorm);
end intrinsic;


intrinsic GlobalCocycleInvariants(L::FldNum, gamma::Map) -> SeqEnum
{ Compute the invariants of the cocycle gamma. }
    return global_cocycle_invariants(L, gamma);
end intrinsic;

intrinsic GlobalCocycle(L::FldNum, locCond::SeqEnum) -> Map
{ Computes a global cocycle in H^2(G,L^\times) respecting the given
  local conditions. These must be given in the form <p,inv_p>
  with inv_p in 1/[L_P:Q_p] Z and sum(inv_p)=0. }
    
    require &+([x[2] : x in locCond ]) in Integers():
            "sum of invariants must be an integer!";
    require &and([IsPrime(x[1]) : x in locCond]) :
            "only specify invariants for prime numbers!";
    // denominators must divide local degree
    require &and([ LocalDegree(L,x[1])/Denominator(x[2]) in Integers() : x in locCond ] ) :
            "denominators of invariants must divide local degrees!";
    
    return global_cocycle(L, locCond);
end intrinsic;


intrinsic FrobeniusAutomorphism(L::FldPad, K::FldPad) -> Map
{ Computes the Frobenius automorphism of K, restricted to L/K.}
    return frobenius_automorphism(L,K);
end intrinsic;
  
intrinsic FrobeniusAutomorphism(L::RngPad, K::RngPad) -> Map
{ Computes the Frobenius automorphism of K, restricted to L/K.}
    return frobenius_automorphism(L,K);
end intrinsic;

intrinsic FrobeniusAutomorphism(psi::Map) -> GrpElt
{ Given the automorphism map psi: G-> Aut(L) of an extension
  L/K. Compute the group element that represents the Frobenius
  automorphism. }
    return frobenius_automorphism_grpelt(psi);
end intrinsic;

intrinsic FrobeniusEquation(c::RngPadElt, precision::RngIntElt) -> RngPadElt, Map
{Solves the equation x^(phi-1)=c, c in OE^\times, up to the given precision,
where phi is the Frobenius automorphism of OK.
The solution x and the automorphism phi are returned.
If a sequence C of elements is given, a sequence of solutions is returned.
If OK is not given, OK=OE. Otherwise, OE must be an extension of OK.
Note, that whenever the norm of c over OK is not 1, this can generate
huge extensions of OE.}
    
    local x, frob;
    x, frob := frobenius_equation(c, precision);
    return x, frob;
end intrinsic;

intrinsic FrobeniusEquation(C::SeqEnum[RngPadElt], precision::RngIntElt) -> SeqEnum[RngPadElt], Map
{Solves the equation x^(phi-1)=c, c in OE^\times, up to the given precision,
where phi is the Frobenius automorphism of OK.
The solution x and the automorphism phi are returned.
If a sequence C of elements is given, a sequence of solutions is returned.
If OK is not given, OK=OE. Otherwise, OE must be an extension of OK.
Note, that whenever the norm of c over OK is not 1, this can generate
huge extensions of OE.}
    
    local x, frob, c, X;
    
    c := C[1];
    q := #ResidueClassField(Parent(c));
    x, frob := frobenius_equation(c, precision);
    X := [x];
    
    for i in [2..#C] do
        O := Parent(x);
        c := C[i];
        x, frob := frobenius_equation_OK(O!c, precision, frob, q);
        X := X cat [x];
    end for;
    
    return X, frob;
end intrinsic;

intrinsic FrobeniusEquation(c::RngPadElt, precision::RngIntElt, OK::RngPad) -> RngPadElt, Map
{Solves the equation x^(phi-1)=c, c in OE^\times, up to the given precision,
where phi is the Frobenius automorphism of OK.
The solution x and the automorphism phi are returned.
If a sequence C of elements is given, a sequence of solutions is returned.
If OK is not given, OK=OE. Otherwise, OE must be an extension of OK.
Note, that whenever the norm of c over OK is not 1, this can generate
huge extensions of OE.}
    
    local q, frob, x;
    
    q := #ResidueClassField(OK);
    frob :=  map<OK -> OK | x :-> x>;
    frob := continue_frobenius(frob, Parent(c), q);
    x, frob := frobenius_equation_OK(c, precision, frob, q);
    return x, frob;
end intrinsic;

intrinsic FrobeniusEquation(C::SeqEnum[RngPadElt], precision::RngIntElt, OK::RngPad) -> SeqEnum[RngPadElt], Map
{Solves the equation x^(phi-1)=c, c in OE^\times, up to the given precision,
where phi is the Frobenius automorphism of OK.
The solution x and the automorphism phi are returned.
If a sequence C of elements is given, a sequence of solutions is returned.
If OK is not given, OK=OE. Otherwise, OE must be an extension of OK.
Note, that whenever the norm of c over OK is not 1, this can generate
huge extensions of OE.}
    
    local q, x, frob, c, X;
    
    q := #ResidueClassField(OK);
    frob :=  map<OK -> OK | x :-> x>;
    frob := continue_frobenius(frob, Parent(C[1]), q);
    
    c := C[1];
    x, frob := frobenius_equation_OK(c, precision, frob, q);
    X := [x];
    
    for i in [2..#C] do
        x, frob := frobenius_equation_OK(Parent(x)!C[i], precision, frob, q);
        X := X cat [x];
    end for;
    
    return X, frob;
end intrinsic;

intrinsic RightToLeftCocycle(m::Map) -> Map
{ Create a left-cocycle from a right-cocycle and vice versa. }
    require Type(Domain(m)) eq SetCart :
            "Map on Cartesian Product required.";
    require NumberOfComponents(Domain(m)) eq 2 :
            "TwoCocycle required.";
    return map< Domain(m) -> Codomain(m) | x :-> m( < x[2]^(-1), x[1]^(-1) >) >;
end intrinsic;

intrinsic LeftToRightCocycle(m::Map) -> Map
{ Create a left-cocycle from a right-cocycle and vice versa. }
    return RightToLeftCocycle(m);
end intrinsic;

intrinsic RightToLeftAction(m::Map) -> Map
{ From a map representing the action of a group acting from right,
  create a corresponding map for the group action from left
  and vice versa. }
    require Type(Domain(m)) eq GrpPerm:
            "Permutation group required.";
    return map< Domain(m) -> Codomain(m) | x :-> m( x^(-1) ) >;
end intrinsic;

intrinsic LeftToRightAction(m::Map) -> Map
{ From a map representing the action of a group acting from right,
  create a corresponding map for the group action from left
  and vice versa. }
    return RightToLeftAction(m);
end intrinsic;










/**************************************************
 *          Local Brauer group
 **************************************************/

function local_brauer_group(L, P, psi, computeLFC)
    local G, pi, theta, m, LP, iota, H, psiL, lfc,
          V, M, proj, mm, HH, g, v, Y, mmY, y, X, qX, mmX, x,
          C, f1, f2, lfc2, H2;
    
    // compute prime ideal
    OL := RingOfIntegers(L);
    p := Generator(P meet Integers());
    vprint Brauer, 1: "Computing cohomology at", p;
    IndentPush();
    
    // Global Galois group
    G := Domain(psi);
    
    // Compute lattice
    pi := UniformizingElement(P);
    theta, m := lattice(P, pi, psi);
    vprint Brauer, 1: "lattice precision:" , m;
    vprint Brauer, 2: "theta: ", theta;
    
    // Localization with enough precision
    vprint Brauer, 1: "compute completion, prec=", 2*m+2;
    vtime Brauer, 1: LP, iota, psiL := completion_with_prec(L, P, psi, 2*m+2);
    // and Galois group
    //H, psiL := localized_automorphism_group(psi, P, iota, Automorphisms(LP, pAdicField(LP)));
    H := Domain(psiL);
    
    // compute V=LP^\times/X
    vprint Brauer, 1: "compute module";
    X, mmX, qX := compute_LPmul_modX(L, P, pi, psi, iota, LP, psiL, theta, m);
    // make this a right action
    mmX := map< H -> Aut(X) | g :-> mmX(g^(-1)) >; 
    // map of L onto X
    qM := iota*qX;
    vprint Brauer, 2: "Dimension will be", #Invariants(X);
    
    // compute cohomology
    vprint Brauer, 1: "compute cohomology group";
    C := CohomologyModule(H, X, mmX);
    H2 := CohomologyGroup(C,2);
    // compute mappings to and from internal module V
    // f1: X -> V
    f1 := map< X -> RSpace(Integers(),Dimension(C)) |
        x:-> Eltseq(x),
        y:-> X!Eltseq(y)
    >;
    
    // compute local fundamental class
    // with values in L^\times
    if computeLFC then
        vprint Brauer, 1: "compute cocycle, prec =",m;
        lfc := cocycle_lfc_G(LP, pAdicField(LP), psiL, m, false);
        // read lfc in cohomology group
        vprint Brauer, 2: "identify cocycle";
        lfc2 := func< x | lfc(x[2]^(-1),x[1]^(-1)) @ qX @ f1 >;
        g := IdentifyTwoCocycle(C, lfc2);
        
        IndentPop();
        return rec< locBrGrp |
            L := L, P := P, p := p,
            M := X, actM := mmX, qM := qM,
            C:=C,f1:=f1, lfc := g,
            theta:=theta
        >;
    else
        IndentPop();
        return rec< locBrGrp |
            L := L, P := P, p := p,
            M := X, actM := mmX, qM := qM,
            C:=C,f1:=f1,
            theta:=theta
        >;
    end if;
end function;

/* Given a bunch of local and global data, this function computes
 * the module M=LP^\times/X, X=exp(Z[G]theta), as an abelian groups
 * as well as the G-automorphisms on M and a projection mappings
 * from LP^\times onto M.
 */
function compute_LPmul_modX(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));
    
    // Compute the exponential value of alpha truncated at the N's summand [Bley HS132, Remark 3.6].
    truncated_exp := func< alpha, N | &+( [  alpha^n/ Factorial(n) : n in [0..N] ] ) >;
    // Compute the precision to which the exponential function values must be known
    // to get correct results in L_P/exp(lattice).
    p := Generator(P meet Integers());
    N := Ceiling( m / (Valuation(theta, P) - RamificationIndex(P)/ (p-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),  N ));
    
    // 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 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;

function compute_LPmul_modX_global(L, P, pi, psiG, theta, m)
    G   := Domain(psiG);  GG := [g : g in G];
    OL  := MaximalOrder(L);
    
    // 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));
    
    
    // Compute the exponential value of alpha truncated at the N's summand [Bley HS132, Remark 3.6].
    truncated_exp := func< alpha, N | &+( [  alpha^n/ Factorial(n) : n in [0..N] ] ) >;
    // Compute the precision to which the exponential function values must be known
    // to get correct results in L_P/exp(lattice).
    p := Generator(P meet Integers());
    N := Ceiling( m / (Valuation(theta, P) - RamificationIndex(P)/ (p-1)  ) );
    
    // exp(calL) in Q wird erzeugt von exp(theta)
    // brauche exp(theta) nur bis zu einer gewissen Genauigkeit
    expTheta := truncated_exp(theta,  N );
    n := Numerator(expTheta);
    d := Denominator(expTheta);
    if Valuation(n, P) + Valuation(OL!d, P) gt 0 then
        print "numden!!!";
        n,d := numden(expTheta, P);
    end if;
    
    // expTheta und Konjugierte in Q lesen
    conjQ := [ phi_OL_Qmal( psiG(g)(n) ) - phi_OL_Qmal( psiG(g)(d) ) : g in GG];
    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..#GG] do
        g := GG[k];
        bilder := [];
        // 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 abelian group 
    Y:=FreeAbelianGroup(#ZgenS+1);
    mmY := map< G -> Aut(Y) | g :-> hom< Y -> Y | 
        y :-> Y!Eltseq(Vector(Eltseq(y))*M[ Index(GG,g)]) > >;
    X, qX := quo<Y | [Order(ZgenS[i])* Y.(i+1) : i in [1..#ZgenS] ]>;
    mmX := map< G -> Aut(X) | g :-> hom< X -> X | x :-> qX( x@@qX @ mmY(g) ) > >; 
    
    // Projektion (globale Rechnung)
    f := map< L -> X |
      x :->  qX(Y!([v] cat Eltseq( phi_OL_S(n) - phi_OL_S(d) )))
            where n,d := numden(x/pi^v, P)  where v := Valuation(x, P) ,
      y :->  pi^yy[1]*phi_S_OL(S!yy[2..#yy]) where yy := Eltseq( y @@ qX )
    >;
    
    return X, mmX, f;
end function;



/**************************************************
 *               Lattice
 **************************************************/

/* lattice(P::RngOrdIdl, pi::RngOrdElt, psi::Map) -> RngOrdElt, RngIntElt
 * Given an ideal P=<pi> of OL and an automorphism map psi:G->Aut(OL),
 * this function computes an element theta,
 * such that L=Z[G]*theta is a full projective lattice upon which
 * the exponential function is injective,
 * and an integer m such that P^m subset L.
 */
function lattice(P, pi, psi)
    local OL;
    OL := Order(P);
    if Degree(OL) eq AbsoluteDegree(OL) then
        return lattice_absolute(P,pi,psi);
    else
        return lattice_relative(P,pi,psi);
    end if;
end function;

/* lattice_absolute(P::RngOrdIdl, pi::RngOrdElt, psiL::Map) -> RngOrdElt, RngIntElt
 * See lattice.
 */
function lattice_absolute(P, pi, psiL)
    local OL, p, theta, v, v1, erz, x, M, ZpGtheta, k, m, M1, M1I, j, M2, T;
    
    OL := Order(P);
    p := Generator(P meet Integers());
    G := Domain(psiL);
    rand := false;
    
    repeat
        theta, v := lattice_generator_theta(psiL, P, pi, rand);
        rand := true;
        // erzeuger des Gitters global
        erz := [OL!(psiL)(g)(theta) : g in G];
        M := VerticalJoin( [Vector(ElementToSequence(x)) : x in erz ]);
        ZpGtheta := Lattice(M);
    until Rank(ZpGtheta) eq Degree(OL);
    
    // finde m mit P^m in ZpGtheta
    // einfacher Ansatz:
    k := Index(StandardLattice(Degree(OL)), ZpGtheta);
    m := Valuation(k*OL, P);
    //m := Valuation(k, setting`p)*RamificationIndex(setting`P);
    
    if m eq 0 then
        return theta, v+1;
    end if;
    
    // 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;

/* lattice_relative(P::RngOrdIdl, pi::RngOrdElt, psiL::Map) -> RngOrdElt, RngIntElt
 * See lattice.
 */
function lattice_relative(P, pi, psiL)
    local OL,OK,p,G,theta,b,erz,g,M,x,y,ZpGtheta,k,m,M1,M1I,M2,j,T;
    
    OL := Order(P);
    OK := BaseRing(OL);
    assert(BaseRing(OK) cmpeq Integers());
    p := Generator(P meet Integers());
    G := Domain(psiL);
    rand := false;
    
    repeat
        theta, v := lattice_generator_theta(psiL, P, pi, rand);
        print theta;
        rand := true;
        // erzeuger des Gitters global
        erz := [OL!(psiL)(g)(theta) : g in G];
        // multiplizieren mit Basis von OK
        erz := [x*(OL!y) : x in erz, y in Basis(OK)];
        M := VerticalJoin( [Vector(&cat([ Eltseq(y) : y in Eltseq(x)])) : x in erz ]);
        ZpGtheta := Lattice(M);
    until Rank(ZpGtheta) eq AbsoluteDegree(OL) and &and([x in Integers() : x in Eltseq(M) ]);
    
    // finde m mit P^m in ZpGtheta
    // einfacher Ansatz:
    print M;
    k := Index(StandardLattice(AbsoluteDegree(OL)), ZpGtheta);
    m := Valuation(k*OL, P);
    
    if m eq 0 then
        return theta, v+1;
    end if;
    
    // kleinstes m
    // schreibe Basis von ZpGtheta in Matrix
    M1 := Matrix(Rationals(), M);
    M1I := M1^(-1);
    for j in [v+1..m] do
        m := j;
        // schreibe Basis von P^m in Matrix
        M2 := Matrix(Rationals(),  &cat([
            [ &cat([Eltseq(y) : y in Eltseq(b*x)]) : x in Basis(OK)]
            : b 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;

// only for L absolute / Q
function find_included_P_power_absolute(P, psiL, theta)
    p := Generator(P meet Integers());
    G := Domain(psiL);
    OL := Order(P);
    v := Valuation(theta, P);
    
    erz := [OL!(psiL)(g)(theta) : g in G];
    M := VerticalJoin( [Vector(ElementToSequence(x)) : x in erz ]);
    ZpGtheta := Lattice(M);
    k := Index(StandardLattice(AbsoluteDegree(OL)), ZpGtheta);
    m := Valuation(k*OL, P);
    
    if m eq 0 then
        return v+1;
    end if;
    
    // 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 m;
end function;

/* Computes generator theta of the lattice 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, rand)
    local OL, p, theta, v, v1;
    
    OL := Order(P);
    p := Generators(P meet Integers())[1];
    
    theta := NormalBasisElement(OL, psi : rand := rand);
    v := Valuation(theta, P);
    v1 := 1+Floor(AbsoluteRamificationIndex(P)/(p-1));
    v := Maximum(0,v1-v);
    theta := OL!(theta*(pi)^v);
    
    return theta, v;
end function;


intrinsic AbsoluteRamificationIndex(P::RngOrdIdl) -> RngIntElt
{ Compute the absolute ramification index of P. }
    local O, Q, e;
    O := Order(P);
    Q := P;
    e := RamificationIndex(P);
    while O cmpne Integers() do
        O := BaseRing(O);
        Q := Q meet O;
        e *:= RamificationIndex(Q);
    end while;
    return e;
end intrinsic;















/**************************************************
 *          Local Fundamental Class
 **************************************************/


function local_fundamental_class_direct(L, n)
    LN := ext<L | RamificationIndex(L,pAdicField(L))>;
    printf "Prec(Qp)=%o, deg(LN)=%o, Prec(LN)=%o\n",
        Precision(pAdicField(L)), AbsoluteDegree(LN), Precision(LN);
    
    // dont change precision, or save "old" field and 
    // read values there at the end
    //ChangePrecision(~LN, 100); // ???
    L := BaseField(LN);
    Qp := pAdicField(LN);
    
    try
        Gamma, psi := AutomorphismGroup(LN,Qp);
    catch e
        print "An error occured within the AutomorphismGroup command:", e`Object;
        print "This seems to happen sometimes.";
        print "Restart Magma and try again.";
        error "Could not compute automorphism group.";
        //time Gamma, psi,_ := automorphismGroup(RingOfIntegers(LN),RingOfIntegers(Qp));
    end try;
    
    print "Computing lfc in unramified extension and its inflation";
    OLN := RingOfIntegers(LN);
    RLN := ResidueClassField(OLN);
    
    // find Frobenius map
    // there will be multiple continuations to LN; take one of them
    frobs := [sig : sig in Gamma | &and([ RLN!OLN!(LN!OLN!x @ psi(sig)) eq x^Prime(L) : x in Basis(RLN)])];
    F := frobs[1];
    // compute orbits and construct the inflation of the local fundamental class
    H := sub< Gamma| F>;
    Q, q := quo< Gamma | H>;
    bahnen := [ [ (sig @@ q)*F^i : sig in Q] : i in [0..5] ];
    lfc := map< car<Gamma, Gamma> -> LN | 
        x :-> ( Index([ x[1] in b : b in bahnen], true) + Index([ x[2] in b : b in bahnen], true) - 2 lt 6
                select LN!1 else LN! Prime(L)  )
    >;
    
    t := Cputime();
    print "Computing cohomology of compositum, deg =", AbsoluteDegree(LN);
    // Compute M=LN^\times/ P^n
    time LNmul, mLNmul_LN := MultiplicativeGroup(LN : Prec := n); // prec = ???
    M := LNmul; mM_LN := mLNmul_LN;
    /*basis := Basis(LN);
    basis := [ b @ psi(g) : b in basis, g in Gamma ];
    pi := UniformizingElement(LN);
    basis2 := [(1+pi^n*x) : x in basis];
    time M, tmp := quo< LNmul | [x @@ mLNmul_LN : x in basis2] >;
    mM_LN := Inverse(tmp)*mLNmul_LN;
    // Compute Galois-action on M
    time psiM := map<Gamma -> Aut(M) | x :-> hom<M -> M | [M.i @ mM_LN @ (x@psi) @@ mM_LN : i in [1..Ngens(M)]]>  >;*/
    psiM := map< Gamma -> Aut(LNmul) | sig :-> hom< LNmul -> LNmul | [ LNmul.i @ mLNmul_LN @ psi(sig) @@ mLNmul_LN : i in [1..Ngens(LNmul)] ] > >;
    
    
    GammaAll := [g : g in Gamma];
    print "Compute Galois-action on module";
    time psiMall := [psiM(g) : g in GammaAll];
    psiMfast := map< Gamma -> Aut(M) | g :-> psiMall[Index(GammaAll,g)] >;
    // Compute cohomology
    CM := CohomologyModule(Gamma, M, psiMfast );
    print "Time for cohomology of LN:", Cputime(t);
    
    //print "Compute H^2(LN)";
    //time _:=CohomologyGroup(CM, 2);
    print "Compute H^1(LN)";
    time _:=CohomologyGroup(CM, 1);
    // Construct map from LN^\times to the cohomology module (via M)
    f1 := map< M -> RSpace(Integers(),Dimension(CM)) |
            x:-> Eltseq(x),
            y:-> M!Eltseq(y)
        >;
    f1 := Inverse(mM_LN)*f1;
    
    
    print "Computing cohomology of L";
    // Compute M_L= L^\times/P^n
    N := 3*n;
    G, psiG := AutomorphismGroup(L,Qp);
    //psiG := map< Domain(psiG) -> Codomain(psiG) | x :-> psiG(x^(-1)) >;
    Lmul, mLmul_L := MultiplicativeGroup(L : Prec := n);   // prec = ???
    ML := Lmul; mML_L := mLmul_L;
    /*
    basis := Basis(L);
    pi := UniformizingElement(L);
    basis2 := [(1+pi^N*x) : x in basis];
    ML, tmp := quo< Lmul | [x @@ mLmul_L : x in basis2] >;
    mML_L := Inverse(tmp)*mLmul_L;
    */
    // Compute Galois-action on M
    psiML := map< G -> Aut(ML) | x :-> hom<ML -> ML | [ML.i @ mML_L @ (x@psiG) @@ mML_L : i in [1..Ngens(ML)]]>  >;
    Gall := [g : g in G];
    psiMLall := [psiML(g) : g in Gall];
    psiMLfast := map< G -> Aut(ML) | g :-> psiMLall[Index(Gall,g)] >;
    // Compute cohomology
    CML := CohomologyModule(G, ML, psiMLfast );
    _:=CohomologyGroup(CML, 2);
    // Construct map from L^\times to the cohomology module (via M_L)
    f2 := map< ML -> RSpace(Integers(),Dimension(CML)) |
            x:-> Eltseq(x),
            y:-> ML!Eltseq(y)
        >;
    f2 := Inverse(mML_L)*f2;
    
    print "Comparing cocycles in compositum";
    // construct all cocyles in H^2(G, M_L)
    H2L := CohomologyGroup(CML, 2);
    cocyclesL := [ map< car<G,G> -> L | x:-> x @ TwoCocycle(CML, c) @@ f2> : c in H2L];
    bahnen := [[h : h in Gamma | psi(h)(LN!L.1) eq psiG(g)(L.1)] : g in Gall];
    // compute inflation to H^2(Gamma, M)
    cocyclesLN := [ func< x | 
         f1(c(< Gall[Index([ x[1] in b : b in bahnen], true)],
                    Gall[Index([ x[2] in b : b in bahnen], true)]>))
        >  : c in cocyclesL ];
    // identify these cocyles and compare with local fundamental class
    lfcLN := func< x | x @ lfc @ f1>;
    //lfcID := IdentifyTwoCocycle(CM, lfcLN);
    //ind := [i : i in [1..#cocyclesLN] | IdentifyTwoCocycle(CM, cocyclesLN[i]) eq lfcID ];
    dif := [ func< x | f(x)- lfcLN(x)> : f in cocyclesLN ];
    for i in [1..#dif] do
        if IsTwoCoboundary(CM, dif[i]) then
            return cocyclesL[i];
        end if;
    end for;
    
    //if #ind ne #H2L/#G then
        error "Should be more precise in computation of cocyles of L!";
    //end if;
    //return cocyclesL[ind[1]];
end function;



/**************************************************
 *      Computation of Frobenius Automorphism 
 **************************************************/

/*
 * Compute the Frobenius automorphism of K, restricted to L/K.
 */
function frobenius_automorphism(L, K)
    phi := map<K->K | x :-> x>;
    if L eq K then
        return phi;
    else
        q := #ResidueClassField(RingOfIntegers(K));
        return continue_frobenius(phi, L, q);
    end if;
end function;

/* 
 * Given the automorphism map psi: G-> Aut(L) of a p-adic extension
 * L/K. Compute a group element that represents the Frobenius
 * automorphism. 
 */
function frobenius_automorphism_grpelt(psi)
    local G, L, K, B;
    G := Domain(psi);
    L := Domain(psi(G!1));
    B := [*L.1*];
    // find base field K of the automorphism group and 
    // generate list of generating elements of L/K.
    K := BaseRing(L);
    if Type(L) eq RngPad then
        Qp := pAdicRing(L);
    else
        Qp := pAdicField(L);
    end if;
    while Degree(L,K) lt #G do
        Append(~B, K.1);
        K := BaseRing(K);
        if K cmpeq Qp and Degree(L,K) lt #G then
            error "Couldn't find fixed field of automorphism group in field tower of L/Qp!";
        end if;
    end while;
    
    OK := RingOfIntegers(K);
    q := #ResidueClassField(OK);
    
    OL := RingOfIntegers(L);
    R := ResidueClassField(OL);
    
    for g in G do
        sigma := psi(g);
        if &and([ R!sigma(b) - (R!b)^q eq R!0 : b in B  ]) then
            return g;
        end if;
    end for;
    
    error "No Frobenius element found!";
end function;

/* Given the Frobenius automorphism of K, restricted to E,
 * and the number q of elements in the residue class group of K.
 * Compute the Frobenius automorphism of K, restricted to L/K.
 */
function continue_frobenius(frob, L, q)

    local K, B, b, f, rts, cont, x, r, m, i, cand;
    //require (Type(L) eq FldPad) or 
    //        (Type(L) eq RngPad): 
    //        "Bad argument types\np-Adic Rings or Fields required";
    
    
    K := Domain(frob);
    
    B := Reverse(generators_pad(L,K));
    for b in B do
        // continuations of frob with values in L
        f := MinimalPolynomial(b);
        rts := {r[1] : r in Roots( Polynomial([ frob(x) :  x in Coefficients(f)]), L) };
        assert #rts eq Degree(f);
        cont := [ map< Parent(b) -> L | 
            x :-> &+([ frob(elem[i])*r^(i-1) : i in [1..#elem]] 
            where elem is elem_to_seq(x,BaseRing(Parent(b)))) >  
            : r in rts 
        ];
        
        for m in cont do
            cand := true;
            for i in [0..Degree(Parent(b))-1] do
                if Valuation(m(b^i)-(b^i)^q) eq 0 then
                    cand := false;
                    break;
                end if;
            end for;
            if cand then
                // found a lift of the frobenius
                frob := m;
                break;
            end if;
        end for;
        
        if not cand then
            error "No Frobenius found!";
        end if;
    end for;
    
    return frob;
end function;

/*
 * Given x=\sum a_i p^i, compute \sum m(a_i) b^i.
 */
function continued_automorphism_image(x, m, b)
    if Type(x) in {RngPadElt, FldPadElt} then
        if x eq ChangePrecision(Zero(Parent(x)),Precision(x)) then
            return ChangePrecision(Zero(Parent(x)),Precision(x));
        end if;
    end if;
    y := ElementToSequence(x);
    return &+([ m(y[i])*b^(i-1)  :  i in [1..#y] ]);
end function;




/**************************************************
 *          Solving Frobenius Equations
 **************************************************/

/* 
 * Given c in OK. Let phi be the Frobenius of \tilde K/K.
 * Solve x^{phi-1}=c up to the given precision.
 */
function frobenius_equation(c, precision)
    OK := Parent(c);
    frob := map<OK -> OK | x :-> x>;
    q := #ResidueClassField(OK);
    x, frob := frobenius_equation_OK(c, precision, frob, q);
    return x, frob;
end function;


/* 
 * Given c in OE, E/K unramified, the Frobenius phi of \tilde K/K and 
 * the characteristic q of the residue class field of K.
 * Solve x^{phi-1}=c up to the given precision.
 */
function frobenius_equation_OK(c, precision, frob, q)
    
    O := Domain(frob);
    pi := UniformizingElement(O);
    F := ResidueClassField(O);
    R<x> := PolynomialRing(F);
    pp := Factorization(x^q-F!c*x);
    f := pp[2,1];
    
    if Degree(f) gt 1 then
        vprintf FrobEq, 1: "-%o-", Degree(f);
        // Erzeuge Erweiterung in der f eine Nullstelle hat
        O  := ext<O | Degree(f)>;
        pi := UniformizingElement(O);
        F  := ResidueClassField(O);
    else
        //print "Keine Koerpererweiterung notwendig";
    end if;
    
    R<x> := PolynomialRing(O);
    rts := Roots(R!f);
    x := rts[1,1];
    if Degree(f) gt 1 then
        frob := continue_frobenius(frob, O, q);
    end if;
    
    // Weitere Schritte
    x, frob := frobenius_equation_cont(c, precision, x, 1, frob, q);
    
    return x, frob;
end function;

/*
 * Given c in OE, E/K unramified, the Frobenius phi of \tilde K/K on OE, 
 * the characteristic q of the residue class field of K.
 * Solve x^{phi-1}=c up to the given precision, staring with the solution x up to precision r.
 */
function frobenius_equation_cont(c, precision, x, r, frob, q)
    
    xi := [x];
    O := Parent(x);
    pi := UniformizingElement(O);
    F := ResidueClassField(O);
    
    // initialisiere Werte a,b fuer den naechsten Schritt
    a := O!(c*x/frob(x));
    bool, b := IsExactlyDivisible(a-1, pi); // ^ r  !!!
    if not bool then
        error "nicht teilbar 1!";
    end if;

    // weitere Schritte
    for j in [r+1..precision] do
        R<X> := PolynomialRing(F);
        pp := Factorization(X^q-X-F!O!b);
        f := pp[1,1];
        
        // Falls noetig, erzeuge Erweiterung
        if Degree(f) eq 1 then
            //print "Keine Erweiterung notwendig in Schritt", i;
        else
            //print "Erweiterung vom Grad", Degree(f);
            vprintf FrobEq, 1: "-%o-", Degree(f);
            O := ext<O | Degree(f)>;
        end if;
        pi := UniformizingElement(O);
        F  := ResidueClassField(O);
        rts := Roots(f, O);
        y := rts[1,1];
        
        lastfrob := frob;
        if Degree(f) gt 1 then
            frob := continue_frobenius(frob, O, q);
        end if;
        
        z := 1+y*pi^(j-1);
        // if not (Valuation(frob(z)/z-a[i]) ge #xi[i]+1) then
        //     print "\nFehler: Die Loesung in Schritt", i, "erfuellt die Valuation-Bedingung nicht!\nTeilergebnis wird ausgegeben.\n";
        //     return [ &*(xi[i])  : i in [1..#C]], lastfrob;
        // end if;
        xi := (xi cat [z]);
        
        if j lt precision then
            // initialisiere Werte a,b fuer naechsten Schritt
            a := (O!a) * z/frob(z);
            bool,bb:=IsExactlyDivisible(O!(a-1),pi^j);
            if bool then
                b := O!bb;
            else
                error "nicht teilbar 2!";
            end if;
        end if;
    end for;        
    
    return &*(xi), frob;
end function;





/**************************************************
 *     Compute Local Fundamental Class
 **************************************************/

/*
 * Compute the cocyle representing the local fundamental class of L/K
 * up to the given precision. If the parameter NoNorm = true, 
 * the algorithm without norm equations is used (not recommended).
 */
function cocycle_lfc_G(L, K, psi, precision, NoNorm)
    local steps, G, OL, Zp, Qp, pi, pi_sigma, pisigmapi, g, i, u_sigma, phi,
          e, OL2ub, OL2, OL1, L1, K1, incl, GG, AutL, AutL1, sigma, psi1,
          OLmal, m, u, bool, gamma1;
    
    // Initialization
    if psi cmpeq 0 then
        _,psi,_ := automorphismGroup(L,K);
        psi := map<Domain(psi) -> Codomain(psi) | g :-> psi(g^(-1))>;
    end if;
    
    steps := precision+2;
    G := Domain(psi);
    OL := RingOfIntegers(L);
    Zp := pAdicRing(RingOfIntegers(OL));
    Qp := FieldOfFractions(Zp);
    d := InertiaDegree(OL,Zp);
    
    // unverzweigt
    if d eq Degree(L,K) then
        pi := UniformizingElement(K);
        phi := FrobeniusAutomorphism(L, K);
        phi := [g : g in G | &and([ is_equal(psi(g)(b), phi(b)) : b in generators_pad(L,K) ]) ];
        if #phi eq 0 then
            // ungenauer
            phi := FrobeniusAutomorphism(L, K);
            phi := [g : g in G | &and([ Valuation(psi(g)(b) - phi(b)) ge L`DefaultPrecision*95/100 : b in generators_pad(L,K) ]) ];
            //print d;
            //print FrobeniusAutomorphism(L, K)(L.1);
            //print [psi(g)(L.1) : g in G ];
            //print G;
        end if;
        phi := phi[1];
        GG := [phi^i : i in [0..#G-1]];
        return map< car<G,G> -> L | x :-> (Index(GG, x[1])+Index(GG,x[2])-2 lt #G select L!1 else pi) >;
    end if;
    
    
    // Maximale Rechengenauigkeit haengt ab von
    // - Genauigkeit von L
    // - Genauigkeit der Automorphismen von L
    // - im Fall NoNorm=true: Genauigkeit von L mit bounded precision.
    if steps gt L`DefaultPrecision then
        error "Precision of p-adic field L not high enough to compute the cocycle!";
    end if;
    if steps gt Minimum([ Precision(r[1]) : r in Roots(DefiningPolynomial(L),L)]) then
        error "Precision of automorphisms of L not high enough to compute the cocyle!";
    end if;
    
    if NoNorm then
        pi := UniformizingElement(L);
        // Frobenius-Gleichungen loesen
        pi_sigma  := [ OL!psi(g)(pi) : g in G];
        pisigmapi := [ OL!(pi_sigma[i]/pi) : i in [1..#pi_sigma]];
        vprint CocycleLFC, 1: "Solve Frobenius equations without norms";
        vtime  CocycleLFC, 1: u_sigma, phi := FrobeniusEquation(pisigmapi, steps);
        
        // Erweiterungskoerper
        OL2 := Parent(u_sigma[1]);
        // Galois-Action on \prod OL2
        frobAction, GAction, frobIndex := galois_act_L_tensor_Knr(OL, OL2, psi, phi);
        
    else
        
        if Precision(L) eq Infinity() then
            e := RamificationIndex(L,K);
            OL2ub := UnramifiedExtension(OL,e);
            vprint CocycleLFC, 1: "Switch to bounded precision, prec = ", OL2ub`DefaultPrecision;
            OL2 := ChangePrecision(OL2ub, OL2ub`DefaultPrecision);
            OL1 := BaseRing(OL2);
            L1 := FieldOfFractions(OL1);
            //L1 := ChangePrecision(L, L`DefaultPrecision);
            //K1 := BaseField(L1);
            //K1 := ChangePrecision(K, K`DefaultPrecision);
            K1 := pAdicField(L1);
            
            incl := Coercion(L,L1);
            G := Domain(psi);
            GG := [g : g in G];
            AutL := [psi(g) : g in GG];
            AutL1 := [ (incl^(-1))*sigma*incl  : sigma in AutL];
            psi1 := map< G -> AutL1 | g :-> AutL1[Index(GG,g)] >;
            
            if steps gt Minimum([ Precision(r[1]) : r in Roots(DefiningPolynomial(L1),L1)]) then
            //if steps gt Minimum([ Precision(psi1(g)(L1.1)) : g in G]) then
                error "Precision of automorphisms of L (bounded) not high enough to compute the cocyle!";
            end if;
            
        else
            L1 := L;
            K1 := K;
            e := RamificationIndex(L1,K1);
            OL1 := RingOfIntegers(L1);
            OL2 := ext<OL1 |e>;
        end if;
        
        vprintf CocycleLFC, 1: "Compute UnitGroup... ";
        vtime   CocycleLFC, 1: OLmal, m := UnitGroup(OL1);
        
        u := UniformizingElement(K1)/UniformizingElement(L1)^e;
        vprintf CocycleLFC, 1: "Solve Norm equation... ";
        vtime   CocycleLFC, 1: bool, gamma1 := NormEquation(OL2,m,OL1!u);
        
        assert bool;
        pi:= gamma1*UniformizingElement(L1);
        
        if Precision(L) eq Infinity() then
            vprint CocycleLFC, 1: "Switch back to unbounded precision";
            pi := OL2ub!pi;
            OL2 := OL2ub;
        else
            pi := OL2!pi;
        end if;
        // Habe OL2 schon und muss die richtigen Forsetzungen haben
        // Galois-Action on \prod OL2
        phi := FrobeniusAutomorphism(OL2, OL);
        frobAction, GAction, frobIndex := galois_act_L_tensor_Knr(OL, OL2, psi, phi);
        pi_sigma := [GAction(g,<pi : i in [1..d]>)[1] : g in G];
        pisigmapi := [ OL2!(pi_sigma[i]/pi) : i in [1..#pi_sigma]];
        vprintf CocycleLFC, 1: "Solve Frobenius equations... ";
        vtime  CocycleLFC, 1: u_sigma, phi := FrobeniusEquation(pisigmapi, steps, OL);
    end if;
    
    if GetVerbose("CocycleLFC") ge 2 then
        vprint CocycleLFC, 2: "Test FrobeniusEquation result";
        assert &and({Valuation(phi(u_sigma[i])/u_sigma[i] - pisigmapi[i]) ge steps : i in [1..#u_sigma]});
    end if;
    
    // Kozykel
    d := InertiaDegree(OL,Zp);
    L2 := FieldOfFractions(OL2);
    prodL2 := Domain(frobAction);
    prodOL2 := CartesianProduct([OL2 : y in [1..d] ]);
    
    tup_sigma := [];
    GG := [g : g in G];
    for g in GG do
        ind := Index(GG,g);
        frobIdx := frobIndex[ind];
        if frobIdx eq 0 then
            frobIdx := d;
        end if;
        Append(~tup_sigma, prodOL2! < u_sigma[ind] * ( i le frobIdx select 1 else  pi_sigma[ind] ) : i in [1..d] > );
    end for;
    
    c := map<car<G,G> -> prodL2   | x  :->   // x = <x[1], x[2]>
                tupelQuotient(
                    tupelProduct(
                        GAction(x[1], tup_sigma[Index(GG, x[2])]),
                        tup_sigma[Index(GG, x[1])]
                    ),
                    tup_sigma[Index(GG, x[1]*x[2])]
                )
            >;
    
    // Testen
    c := precompute_map(c);
    // Komponenten gleich modulo pi^(precision+1)
    assert Minimum([ Minimum([Valuation(y[1]-y[i]) : i in [1..#y]])
        where y is c(x) :  x in Domain(c)]) ge (precision+1);
    if Degree(Codomain(c)[1]) ge 2 then
        // erste Komponente in L modulo pi^(precision+1)
        assert Minimum([ Minimum([ Valuation(z) : z in ElementToSequence(y[1])[2..Degree(Parent(y[1]))]])
            where y is c(x)  :  x in Domain(c)]) ge (precision+1);
    end if;
    
    // Wegen der Theorie muessen die Werte in einem Tupel c(x) alle gleich sein,
    // d.h. c(x) in L. Anschließend noch invertieren um die lokale Fundamentalklasse zu
    // erhalten.
    
    // Man verliert insgesamt maximal 2 Stellen Genauigkeit, weil im inversen Bild von c ein
    // Produkt sigma(b_tau)*b_sigma und die Komponenten von b_sigma jeweils
    // maximal Bewertung 1 haben.
    gamma := map< Domain(c) -> FieldOfFractions(L) | x :->  ( elem_to_seq(c(x)[1], L)[1] )^(-1) >;
    
    return gamma;
end function;

/* Given an automorphism sigma of L and K\subset L. 
 * Return the sigma|_K.
 */
function restrictAutomorphism(sigma, K)
    incl := Coercion(K, Domain(sigma));
    return hom<K-> K | x:-> sigma(x), y:-> (sigma^(-1))(y)>; 
end function;

/* galois_act_L_tensor_Knr(OL::RngPad, OL2::RngPad, psi::Map, phi::Map) -> Map, Map, SeqEnum
 * Let L/Qp be a p-adic field with Galois group G and inertia degree d,
 * such that L is an extensions of its maximal unramfied subextension.
 * Furthermore, let L2 be unramified over L, psi:G->Aut(L/Qp) a map
 * as computed by AutomorphismGroup(L,Qp) and phi the Frobenius
 * automorphism of L2/L.
 * If K2 denotes the maximal unramified extension in L2/Qp, this function
 * computes the actions of phi and G on L\tensor_Qp K2 = \prod_d L2
 */
function galois_act_L_tensor_Knr(OL, OL2, psi, phi)
    local OK, Zp, G, GG, g, d, sigmaHut, frobIndex, L2, prodL2, 
        frobeniusMap, Gaction, x;
    
    
    // Test whether algorithm is applicable
    OK := maximal_unram_subext_simple(OL);
    Zp := BaseRing(OK);
    
    if not generates_maximal_unram(OL2, OL, OK) then
        error "The maximal unramified extension in OL2/Zp cannot be deduced from OL2/OL!";
    end if;
    
    // Initialization
    G := Domain(psi);
    GG := [g : g in G];
    d := InertiaDegree(OL,Zp);
    
    // Compute i such that sigma_OK = phi^i for all sigma in G
    // and extensions sigmaHut of sigma such that sigmaHut^(-1)=phi^(-i) on K2
    sigmaHut, frobIndex := continuations_with_unram_restriction(OL, G, psi, OL2);
    
    // Frobenius automorphism on \prod L2
    L2 := FieldOfFractions(OL2);
    prodL2 := CartesianProduct([L2 : y in [1..d] ]);
    frobeniusMap := map< prodL2 -> prodL2 | x :->  < i eq d select phi(x[1]) else x[i+1]  : i in [1..d]> >;
    
    // action of G on \prod L2
    Gaction := map< car<G, prodL2> -> prodL2 | x :-> 
        apply_map( (frobeniusMap^((d-frobIndex[Index(GG,x[1])]) mod d))( x[2] ), sigmaHut[Index(GG, x[1])]) >;
    
    return frobeniusMap, Gaction, frobIndex;
end function;

/* apply_map(y::., m::Map) -> .
 * Apply m to all elements of y, which may be of Type Tup or SeqEnum.
 */
function apply_map(y, m)
    local n, i;
    
    case Type(y):
        when Tup:
            n := NumberOfComponents(Parent(y));
        when SeqEnum:
            n := #y;
    else:
      error "Typ nicht unterstuetzt.";
    end case;
    
    for i in [1..n] do
        y[i]:= m(y[i]);
    end for;
    return y;
end function;

/* tupelQuotient(t1::Tup, t2::Tup) -> Tup
 * Compute the quotient t_1/t_2 for t_1,t_2 of type Tup.
 */
function tupelQuotient(t1, t2)
    return < t1[i]/t2[i]  : i in [1..#t1]>;
end function;

/* tupelProduct(t1::Tup, t2::Tup) -> Tup
 * Compute the product t_1 t_2 for t_1,t_2 of type Tup.
 */
function tupelProduct(t1, t2)
    return < t1[i]*t2[i]  : i in [1..#t1]>;
end function;

/* tupelInverse(t::Tup) -> Tup
 * Compute the inverse 1/t for t of type Tup.
 */
function tupelInverse(t)
    return < t[i]^(-1) : i in [1..#t] >;
end function;

/* continuations_with_unram_restriction(OL::RngPad, G::Grp, psi_OL_Zp::Map, OL2::RngPad) -> SeqEnum, SeqEnum, SeqEnum
 * Given an unramified extension OL2/OL and the Galois group G of OL/Zp
 * as well as a map psi_OL_Zp: G -> Aut(OL, Zp).
 * Compute a continuation of each automorphism to OL2 and its restriction
 * to the unramified subextension of OL2 as a power of the frobenius
 * automorphism.
 * These continuations are selected in such a way, that the inverse of its
 * restriction to the unramified subextension are low powers of the
 * frobenius automorphism.
 */
function continuations_with_unram_restriction(OL, G, psi_OL_Zp, OL2)
    local Hom_OL_Zp, OK, Zp, sigma, sig, f,
          sigmaHut, frobIndex, sigmaKnrExponentInv;
    
    //require Type(Domain(psi_OL_Zp(G.1))) in {RngPad, FldPad} :
    //        "Bad argument types\nAutomorphisms of p-adic rings/fields needed.";
    
    Hom_OL_Zp := Codomain(psi_OL_Zp);
    OK := maximal_unram_subext_simple(OL);
    Zp := BaseRing(OK);
    
    if not generates_maximal_unram(OL2, OL, OK) then
        print "ACHTUNG: Fuer den angegebenen Koerperturm ist der schnelle Algorithmus nicht anwendbar!";
        print "Benutze den langsamen Algorithmus via primitiven Elementen.";
        
        error "Could not compute 'unramified behaviour' of extension";
    end if;
    
    // Jetzt ist folgendes bekannt:
    // Seien E_i die Zwischenkoerper von OL2/OL mit E_0=OL und E_n=OL2.
    // Dann wird die maximal unverzweigte Teilerweiterung Knr von OL2/OK
    // erzeugt von {E_1.1, E_2.1,...,E_n.1} und Knr(E_0.1)=OL2.
    
    // die Automorphismen von OL/Zp
    sigma := [psi_OL_Zp(g) : g in G];
    if Type(Domain(psi_OL_Zp(G.1))) eq FldPad then
        // Automorphismen auf Ganzheitsring einschraenken
        sigma := [map<OL -> OL | x :-> sig(x) > : sig in sigma];
    end if;
    
    d := Degree(OK);
    d2 := InertiaDegree(OL2, Zp);
    
    if d eq 1 then
        // Ausgangssituation voll verzweigt
        
        frobIndex := [Zero(Integers()) : sig in sigma];
        sigmaKnrExponentInv := frobIndex;
        
        // setze sigma fort mit Identitaet auf unverzweigten Teil
        sigmaHut := [];
        for i in [1..#sigma] do
            sig := sigma[i];
            B := generators_pad(OL2,OL);
            for r in Reverse(B) do
                sig := map<Parent(r) -> Parent(r) | x:-> continued_automorphism_image(x, sig, r) >;
            end for;
            Append(~sigmaHut, sig);
        end for;
        
    elif OL eq OK then
        // Ausgangssitugation unverzweigt
        
        f := FrobeniusAutomorphism(OK, Zp);
        frobIndex := [find_power(sig, f, OK.1 , InertiaDegree(OK, pAdicRing(OK))) : sig in sigma];
        fInv := inverseAutomorphism(f);
        
        sigmaHut := sigma;
        sigmaKnrExponentInv := [(d-i) mod d : i in frobIndex];
        
    else
        // Verhalten von sigma\in\Aut(OL,Zp) auf OK bestimmen
        f := FrobeniusAutomorphism(OK, Zp);
        frobIndex := [find_power(sig, f, OK.1 , d) : sig in sigma];
        fInv := inverseAutomorphism(f);
        
        // Frobenius und Inverse auf OL fortsetzen
        f := continue_frobenius(f, OL, Prime(Zp));
        fInv := inverseAutomorphism(f, fInv);
        
        // Frobenius und Inverses auf OL2 fortsetzen,
        // dabei jeweils die Bilder f^(-i)(b) merken, wobei i=0..d-1 und b
        // Erzeuger der Erweiterung OL2/OL
        //vprint cohomTerm, 3: "Berechne Frobenius und sein Inverses";
        B := Reverse(generators_pad(OL2,OL));
        fB := [**];
        for b in B do
            f := continue_frobenius(f, Parent(b), Prime(Zp));
            fInv := inverseAutomorphism(f, fInv);
            seq := [b];
            for i in [1..d-1] do
                Append(~seq, fInv(seq[#seq]));
            end for;
            Append(~fB, seq);
        end for;
        
        //vprint cohomTerm, 3: "Erzeuge Fortsetzungen mit gewuenschtem Verhalten";
        // setze nun sigma forst
        // Falls sigma=f^j auf OK, j=0..d-1, und [OL2:OL]=m
        // Dann gilt fuer die Fortsetzungen sigmaHut|_Knr = F^(k*d+j)
        // Wir wollen die Fortsetzungen sigmaHut, so dass
        // F^(d-j)*sigmaHut|_Knr = id|_Knr
        // Wir waehlen also die Bilder zu Finv^(d-j)
        sigmaHut := [];
        for i in [1..#sigma] do
            sig := sigma[i];
            idx := frobIndex[i];
            // Liste f^(-idx)(b) fuer die Erzeuger von OL2/OL
            B := [* seq[((d-idx) mod d) + 1]   : seq in fB *];
            // sigma schrittweise fortsetzen
            for r in B do
                sig := map<Parent(r) -> Parent(r) | x:-> continued_automorphism_image(x, sig, r) >;
            end for;
            Append(~sigmaHut, sig);
        end for;
        
        // so konstruiert, dass das Inverse auf Knr einfach zu berechnen ist
        sigmaKnrExponentInv := frobIndex;
    end if;
    
    return sigmaHut, frobIndex, sigmaKnrExponentInv;
end function;

// simple function to check whether an element is zero according to its precision
is_zero  := func< x   | x eq ChangePrecision(Zero(Parent(x)),Precision(x)) >;
is_equal := func< x,y | ChangePrecision(x,m) eq ChangePrecision(y,m)
                        where m := Minimum(Precision(x), Precision(y)) >;

/* find_power(g::Map, f::Map, X::SeqEnum : max := 0) -> RngIntElt
 * Returns integer k such that g=f^k for the elements X.
 * If one does not specify a maximum for k this can run forever. 
 */
function find_power(g, f, x, max)
    local p, m, fX, gX, X;
    
    if Type(x) eq SeqEnum then
        X := x;
    else
        X := [x];
    end if;
    
    
    m := (max eq 0 select Infinity() else max );
    p := 0;
    fX := X;
    gX := [g(x) : x in X];
    
    while p le m do
        if &and([ is_zero(gX[i]-fX[i]) : i in [1..#fX]]) then
            return p;
        end if;
        fX := [f(x) : x in fX];
        p := p+1;
    end while;
    
    error "Maximum power exceeded!!!";
end function;

/* precompute_map(m::Map) -> Map
 * Compute all values of m and return a new map producing
 * the same values. If the definition of the map m was very
 * complicated, this can speed things up significantly.
 */
function precompute_map(m)
    local domSeq, x, img;
    
    domSeq := [x : x in Domain(m)];
    img := [m(x) : x in domSeq];
    
    return map< Domain(m) -> Codomain(m) | x :->  img[Index(domSeq,x)] >;
end function;








/**************************************************
 *            Global Brauer groups
 **************************************************/

/* function global_cocycle_invariants(L::FldNum, gamma::Map) -> SeqEnum
 */
function global_cocycle_invariants(L, gamma)
    //print "Compute prime ideals";
    OL := RingOfIntegers(L);
    values := Setseq({gamma(g) : g in Domain(gamma)});
    // prime numbers for which the valuation is non-trivial
    // including ramified primes
    primes := [ p :
        p in Setseq(Seqset(
            //&cat([[ f[1] : f  in Factorization(Generator(v*OL meet Integers()))] : v in values])
            &cat([[ Generator(f[1] meet Integers()) : f in Factorization(v*OL)] : v in values])
            cat
            [f[1] : f in Factorization(Discriminant(OL))]
        ))
        | LocalDegree(L,p) ne 1
    ];
    // select a prime ideal above each prime
    // primes := [ Factorization(p*OL)[1,1] : p in primes ];
    
    // compute invariants
    invar := [];
    for p in primes do
        locBrGrp := LocalBrauerGroup(L,p : lfc);
        c := IdentifyTwoCocycle(locBrGrp`C, func< x|  gamma(x[2]^(-1),x[1]^(-1)) @ locBrGrp`qM @ locBrGrp`f1 >);
        inv := [i : i in [1..#Group(locBrGrp`C)] |  i*locBrGrp`lfc eq c ][1]/#Group(locBrGrp`C);
        if not inv in Integers() then
            Append(~invar, <p,inv>);
        end if;
    end for;
    
    return invar;
end function;

/* function global_cocycle(L::FldNum, locCond::SeqEnum) 
 */
function global_cocycle(L, locCond)
    local OL, G, psi, g, x, y, 
          S, Sall, US, mU, s, actUS,
          maxX, maxY, P, locBr, indX, alpha, cond, b, newvars,
          i, sig, tau, rho, p, inv, r, H, HH, c, e, beta, indY, cc, seq,
          numnewvars, M, N, sol, N1, t;
    
    t := Cputime();
    OL := RingOfIntegers(L);
    G, _, psi := AutomorphismGroup(L);
    psi := map< Domain(psi) -> Codomain(psi) | g :-> psi(g^(-1)) >;
    
    // conditions for ramified primes
    locCond cat:= [<x[1],0> : x in Factorization(Discriminant(OL)) | not x[1] in {y[1] : y in locCond} ];
    
    S := [x[1] : x in locCond];
    print "compute S-units S=", S;
    //Sall := Setseq( &join({{ x[1] : x in Factorization(Generator(P meet Integers())*OL)} : P in S}));
    Sall := Setseq( &join({{ x[1] : x in Factorization(p*OL)} : p in S}));
    US, mU := SUnitGroup(Sall);
    GG := [sig : sig in G];
    s := #Generators(US);
    print "compute S-unit action";
    // G-actions on US
    actUS := [ SUnitAction(mU, psi(g), Sall ) : g in GG];
    
    print "Time for set S:", Cputime(t);
    t := Cputime();
    
    print "compute local Brauer groups";
    
    // initialize local Brauer groups
    maxY := 0;
    for P in S do
        locBr :=LocalBrauerGroup(L, Integers()!P : lfc);
        maxY +:= #Invariants(locBr`M)*#DecompositionGroup(locBr`P);
    end for;
    
    // produces index for x_(sig,tau,i) in vector
    // (x_(1,1,1)...x_(1,1,s),x_(1,G.1,1)...,x_(1,G.1,s),......)
    indX := func< sig,tau,i | (Index(GG, sig)-1)*#GG*s + (Index(GG,tau)-1)*s+i>;
    maxX := (#G^2)*s;
    // the elements alpha_(sig,i,j) in Z
    // representing the G-action on US
    alpha := func< sig, i,j |  Eltseq(actUS[Index(GG,sig)](US.i))[j] >;
    
    print "Time for local Brauer groups:", Cputime(t);
    
    
    
    t := Cputime();
    // produce conditions
    cond    := []; // sequences with coefficients of x_(sig,tau,i)
    b       := []; // right side vector
    newvars := []; // orders for new variables
    
    print "create conditions for cocycle";
    // cocycle-conditions
    for i in [1..s], sig in G, tau in G, rho in G do
        // x_(sig*tau,rho,i)+x_(sig,tau,i)
        // - sum_j alpha_(sig,j,i)*x_(tau,rho,j) - x_(sig,tau*rho,i) 
        // - ord(US.i) * some_new_variable  == 0
        seq := [0 : i in [1..maxX+maxY]];
        seq[indX(sig*tau, rho, i)] +:=1;
        seq[indX(sig, tau, i)] +:=1;
        for j in [1..s] do
            seq[indX(tau, rho, j)] -:= alpha(sig,j,i);
        end for;
        seq[indX(sig, tau*rho, i)] -:=1;
        if Order(US.i) eq 0 then
            cond    cat:= [seq];
            b       cat:= [0];
            newvars cat:= [0];
        else
            // conditions with new variables at top
            cond    := [seq] cat cond;
            b       := [0] cat b;
            newvars := [Order(US.i)] cat newvars;
        end if;
    end for;
    
    
    // conditions for prime ideals
    // - compute local fundamental class c
    loc := locCond[1];
    posY := maxX;
    for loc in locCond do
        p := Integers()!loc[1];
        inv := loc[2];
        
        locBr:=LocalBrauerGroup(L, p : lfc);
        P := locBr`P;
        r := #Invariants(locBr`M);
        H := DecompositionGroup(P);
        HH := [h : h in H];
        // multiple of local fundamental class as cocycle
        c := inv*LocalDegree(locBr`P) * locBr`lfc;
        c := TwoCocycle(locBr`C, c);
        // cocycle left action
        c := map< car<H,H> -> locBr`M | x :-> c(<x[2]^(-1), x[1]^(-1)>) @@ locBr`f1 >;
        
        // read generators of US in M
        // function for their coefficients
        e := func< k,i | Eltseq(US.k @ mU @ locBr`qM)[i] >;
        // precompute values of the function e
        // reduced computation time by 60% !!
        eDom := [[k,i] : k in [1..s], i in [1..r] ];
        eImg := [e(x[1],x[2]) : x in eDom];
        e := func< k,i | eImg[Index(eDom,[k,i])] >;
        
        // M = prod <w_i>
        // rewrite sig(w_i) in terms of w_i
        // function for their coefficients
        beta := func< sig, i, j | Eltseq(locBr`M.i @ locBr`actM(sig))[j] >;
        
        // can differ by a coboundary y
        // write y(sig) = prod w_i^y_(sig,i) with indeterminants y_(sig,i)
        // function for the index of these indeterminants
        indY := func< sig, i | (Index(HH, sig)-1)*r + i + posY >;
        
        print "create conditions for p =", p, ", r =", r;
        
        // conditions
        for sig in H, tau in H do
            // compute the elements gamma_j in Z
            // representing c(sig,tau)
            cc := Eltseq(c(sig,tau));
            
            for i in [1..r] do
                // sum_k e_(k,i)*x(sig,tau,k) 
                // - sum_j beta_(sig,j,i)*y(sig,j)  - y_(sig,i) + y_(sig*tau,i)
                // - ord(w_i)*new_var
                // == c_i
                seq := [0 : i in [1..maxX+maxY]];
                for k in [1..s] do
                    seq[indX(sig, tau, k)] +:= e(k,i);
                end for;
                for j in [1..r] do
                    seq[indY(tau, j)] -:= beta(sig,j,i);
                end for;
                seq[indY(sig, i)] -:= 1;
                seq[indY(sig*tau, i)] +:= 1;
                
                if Order(locBr`M.i) eq 0 then
                    cond    cat:= [seq];
                    b       cat:= [ cc[i] ];
                    newvars cat:= [0];
                else
                    // conditions with new variables at top
                    cond    := [seq] cat cond;
                    b       := [ cc[i] ] cat b;
                    newvars := [ Order(locBr`M.i) ] cat newvars;
                end if;
            end for;
        end for;
        
        // add to posY for y-variables
        posY +:= #H*r;
    end for;
    
    print "create conditions for normalized cocycle";
    
    // normalized cocycle c(sig,1)=c(1,sig)=1
    for sig in G, i in [1..s] do
        seq := [0 : i in [1..maxX+maxY]];
        seq[indX(sig,1, i)] +:= 1;
        cond    cat:= [seq];
        b       cat:= [0];
        newvars cat:= [0];
        
        seq := [0 : i in [1..maxX+maxY]];
        seq[indX(1,sig, i)] +:= 1;
        cond    cat:= [seq];
        b       cat:= [0];
        newvars cat:= [0];
    end for;
    
    print "create matrix";
    
    // create matrix
    numnewvars := #newvars -#[0 : x in newvars | x eq 0 ];
    // first conditions with extra variable
    M := Matrix(cond[1..numnewvars]);
    // coefficient of extra variable
    M := HorizontalJoin(M, DiagonalMatrix(newvars[1..numnewvars]));
    // latter conditions without extra variables
    M := VerticalJoin(M, Matrix(
        [ seq cat [0 :i in [1..numnewvars] ] : seq in cond[numnewvars+1..#cond] ]
    ));
    print "Time for conditions:", Cputime(t);
    
    // compute solution
    print "solve equation with", NumberOfRows(M),"equations in", NumberOfColumns(M),"variables";
    time sol,N := Solution(Transpose(M),Vector(b));
    
    // a short solution with LLL
    //N1,S := LLL(VerticalJoin(sol, Matrix([ Eltseq(n) : n in Generators(N)])));
    //if Index(Eltseq(Transpose(S)[1]),1) ne 0 then
    //    sol := N1[Index(Eltseq(Transpose(S)[1]),1)];
    //end if;
    print "find short solution";
    time N1,S := LLL(VerticalJoin(
        Vector(Eltseq(sol)[1..maxX]),
        Matrix([ Eltseq(n)[1..maxX] : n in Generators(N)])
    ));
    if Index(Eltseq(Transpose(S)[1]),1) ne 0 then
        sol := N1[Index(Eltseq(Transpose(S)[1]),1)];
    else
        seq := Eltseq(Transpose(S)[1]);
        //print #seq;
        sol := Eltseq(Solution(Transpose(Matrix([seq])), Matrix([[1]])));
        //print #sol;
        //print N1[91];
        sol := &+([Eltseq(sol)[i]*N1[i] : i in [1..#Eltseq(sol)]]);
    end if;
    
    // write as cocycle
    c := map< car<G,G> -> L | x :-> mU(&*([ sol[indX(x[1],x[2],i)]*US.i  : i in [1..s]] )) >;
    
    return c;
end function;











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

// InertiaDegree(L,K) eq Degree(L,K)
intrinsic isUnramified(L::FldPad, K::FldPad) -> BoolElt
{ Test whether L/K is unramified. }
    return InertiaDegree(L,K) eq Degree(L,K);
end intrinsic;
intrinsic isUnramified(L::RngPad, K::RngPad) -> BoolElt
{ Test whether L/K is unramified. }
    return InertiaDegree(L,K) eq Degree(L,K);
end intrinsic;

// RamificationIndex(L,K) eq Degree(L,K)
intrinsic isTotallyRamified(L::FldPad, K::FldPad) -> BoolElt
{ Test whether L/K is totally ramified. }
    return RamificationIndex(L,K) eq Degree(L,K);
end intrinsic;
intrinsic isTotallyRamified(L::RngPad, K::RngPad) -> BoolElt
{ Test whether L/K is totally ramified. }
    return RamificationIndex(L,K) eq Degree(L,K);
end intrinsic;

intrinsic LocalDegree(P::RngOrdIdl) -> RngIntElt
{ return InertiaDegree*RamificationDegree}
    return InertiaDegree(P)*RamificationDegree(P);
end intrinsic;

intrinsic LocalDegree(L::FldNum, p::RngIntElt) -> RngIntElt
{ return InertiaDegree*RamificationDegree}
    return Degree(L)/#Decomposition(L,p);
end intrinsic;


/* 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 := Maximum(precision,100);
    else
        prec := Maximum(precision,30);
    end if;
    
    repeat
        err := false;
        compatible := false;
        try
            //print "compute completion", prec;
            LP, iota := Completion(L, P : Precision:=prec);
            autLP := Automorphisms(LP, pAdicField(LP));
            _, psiLP := localized_automorphism_group(psi, P, iota, autLP);
            /*H := [g : g in Domain(psi) | &and([  psi(g)(x) in P   : x in Generators(P)]) ];
            H := sub< Domain(psi) | H>;
            HH := [h : h in H];
            maps := [map< LP -> LP | x:-> x @@ iota @ psi(h) @ iota> :  h in HH];
            psiLP := map< H -> maps | h :-> maps[Index(HH,h)] >;
            */
            //print "test compatibility";
            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;




intrinsic automorphismGroup(OL::RngPad, OK::RngPad) -> GrpPerm, Map, SeqEnum
{ Compute the automorphism group of a local field extension;
  similar to the Magma function AutomorphismGroup, but also
  works, if the fields are of infinite precision. }
    return automorphism_group(OL, OK);
end intrinsic;

/* See automorphismGroup.
 */
function automorphism_group(OL, OK)
    local p, f, rts, r, Aut, permut, i, g, G, H, x, psi;
    
    p, f := primitive_element(OL, OK);
    
    rts := Roots(f, OL);
    rts := [r[1] : r in rts];
    //assert #rts eq Degree(OL, OK);
    
    Aut := Automorphisms(OL, OK);
    
    //G := PermutationGroup< {1..Degree(OL,OK)} | [[Index(rts,r) : r in  [g(r) : r in rts] ] : g in Aut] >;
    permut := [   [   [i : i in [1..#rts] | Valuation(rts[i]-r) eq Minimum(Precision(r),Precision(rts[i])) ][1]
             : r in  [g(r) : r in rts] ]
         : g in Aut ];
    G := PermutationGroup< {1..#rts} | permut >;
    H := [G!x : x in permut];
    
    psi := map< G -> Aut | x :-> Aut[Index(H,x)] >;
    
    return G, psi, Aut;
end function;

intrinsic inverseAutomorphism(aut::Map) -> Map
{ Compute the inverse of the given automorphism.
  If inv is given, the inverse is computed as a continuation of inv if possible. }
    local K;
    
    K := pAdicRing(Domain(aut));
    return inverseAutomorphism(aut, hom<K-> K| x:-> x>);
end intrinsic;

intrinsic inverseAutomorphism(aut::Map, inv::Map) -> Map
{ Compute the inverse of the given automorphism.
  If inv is given, the inverse is computed as a continuation of inv if possible. }
    local L, K, B, invAut, b, cont, i, c;
    
    L := Domain(aut);
    K := Domain(inv);
    
    invAut := inv;
    B := Reverse(generators_pad(L,K));
    for b in B do
        if (aut(b) eq b) then
            invAut := map<Parent(b) -> Parent(b) | x :-> continued_automorphism_image(x, invAut, b) >;
        else
            f := MinimalPolynomial(b);
            rts := {r[1] : r in Roots( Polynomial([ invAut(x) :  x in Coefficients(f)]), L) };
            assert #rts eq Degree(f);
            cont := [ map< Parent(b) -> L | 
                x :-> &+([ invAut(elem[i])*r^(i-1) : i in [1..#elem]] 
                where elem is elem_to_seq(x,BaseRing(Parent(b)))) >  
                : r in rts 
            ];
            _, i := Maximum([Valuation(aut(c(b))-b) :  c in cont]);
            invAut := cont[i];
        end if;
    end for;
    
    return invAut;
end intrinsic;

/*
 * Creates a list of generators of L/K.
 * If L=L_n/L_{n-1}/.../L_0=K, this returns the list [L_n.1,L_{n-1}.1,...,L_0.1].
 */
function generators_pad(L, K)
    G := [*  *];
    E := L;
    while not E eq K do
        if E.1 eq Zero(E) then
            Append(~G, One(E));
        else
            Append(~G, E.1);
        end if;
        E := BaseRing(E);
        //Insert(~G,1, E.1);
    end while;
    return G;
end function;

/* generates_maximal_unram(OL2::RngPad, OL::RngPad, OK::RngPad) -> BoolElt
 * Let OK be the maximal unramified extension of Zp in OL and OL2 an
 * extension of OL.
 * Test whether the maximal unramified extension OK2 in OL2/Zp can
 * be deduced from the extension OL2/OL, i.e. whether the coefficients of
 * the defining polynomial of OL2/OL are in OK.
 */
function generates_maximal_unram(OL2, OL, OK)
    local fields, E, d;
    
    if Degree(OL2, OL) eq 1 then
        return true;
    end if;
    
    //assert BaseRing(OL) eq OK;
    
    // erzeuge Koerperturm
    fields := Reverse(field_tower(OL2, OL));
    
    // OL selbst wollen wir nicht
    Remove(~fields,1);
    
    d := Degree(OL, OK);
    for E in fields do
        f := DefiningPolynomial(E);
        coeff := Coefficients(f);
        // Koeffizienten in Basisdarstellung
        coeff := [elem_to_seq(c, OK) : c in coeff];
        // jetzt darf jeweils nur jeder d-te Eintrag ungleich Null sein
        // entferne die Eintraege i*(d-1)+1
        for x in coeff do
            c := x;
            for i in [0..#c/d-1] do
                //print "entferne",  i*(d-1)+1;
                Remove(~c, i*(d-1)+1);
            end for;
            if SequenceToSet(c) ne {ChangePrecision(Zero(OK), Minimum([Precision(y) : y in c]))} then
                return false;
            end if;
        end for;
    end for;
    
    return true;
end function;

/* 
 * Compute the list of extensions from K to L, i.e. 
 * if L=E_1/E_2/.../E_n=K, return [*E_1,...,E_n*].
 */
function field_tower(L, K)
    lst := [* *];
    E := L;
    while not E eq K do
        Append(~lst, E);
        if E eq BaseRing(E) then
            error "K is not a subfield of L";
        end if;
        E := BaseRing(E);
    end while;
    Append(~lst, K);
    return lst;
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;

/*
 * Given L/Qp, where L is either totally ramified over Qp, unramified over Qp or
 * a totally ramified over an unramified extension of Qp.
 * Computes the maximal unramified extension of Qp in L.
 */
function maximal_unram_subext_simple(OL)
    local Zp, OK;
    
    if Type(OL) eq RngPad then
        Zp := pAdicRing(OL);
    else
        Zp := pAdicField(OL);
    end if;
    if isUnramified(OL,Zp) then
        OK := OL;
    elif isTotallyRamified(OL,Zp) then
        OK := Zp;
    else
        OK := OL;
        while AbsoluteDegree(OK) gt AbsoluteInertiaDegree(OK) do
            OK := BaseRing(OK);
        end while;
        assert AbsoluteDegree(OK) eq AbsoluteInertiaDegree(OK);
    end if;
    // habe nun Erweiterungen OL/OK/Zp 
    // mit OL/OK voll verzweigt und OK/Zp unverzweigt
    assert isTotallyRamified(OL, OK);
    assert isUnramified(OK, Zp);
    // OK soll "einfache" Erweiterung sein
    assert Degree(OK) eq Degree(OK, Zp);
    
    return OK;
end function;

/*
 * Compute a primitive element of the extension L/K heuristically; may fail!
 */
function primitive_element(L, K)
	//print "primitive_element";
    E := BaseRing(L);
    if E eq K then
        return L.1, MinimalPolynomial(L.1);
    else
        alpha := L.1;
        //print("Berechne primitives Element rekursiv");
        beta := primitive_element(E,K);
        //print beta;
        n := Degree(L,K);
        for j in [1..10] do
            gamma := Random(1,10)*alpha + beta;
            //print "berechne MiPo";
            f := MinimalPolynomial(gamma,K);
            if Type(K) eq RngPad then
                f := PolynomialRing(FieldOfFractions(K))!f;
            end if;
            try
                //print "faktorisiere quadratfrei";
                //f := SquareFreeFactorization(f)[1,1];
                //print "faktorisiere";
                if Degree(Factorization(f)[1,1]) eq n then
                    return gamma, f;
                end if;
            catch e
                gamma := 0;
            end try;
        end for;
        error "Did not find a primitive element!";
    end if;
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 fales, 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)
    local D, B, gens, actD, actB, seq, U;
    
    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 modP then
        seq := [ phi(actD(sig, x)) - actB(sig, phi(x)) : x in gens, sig in 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;



























/**************************************************
 *         Algorithms by Bley
 **************************************************/

intrinsic NormalBasisElement(OL :: RngOrd, h :: Map : rand := false) -> RngOrdElt
{ Given an order OL and a automorphism map h:G->Aut(OL),
  this function returns a normal basis element for OL. }
    local G, b, D, found;
    
    G := Domain(h);

    if not rand then
        for b in Basis(OL) do
            D := Matrix([ ElementToSequence( h(g)(b) ) : g in G ]);
            if Determinant(D) ne 0 then
                return OL!b;
            end if;
        end for;
    end if;
    
    found := false;
    while  not found  do
        b := OL ! [ Random(3) : i in [1..#G] ];
	D := Matrix([ ElementToSequence( h(g)(b) ) : g in G ]);
        if Determinant(D) ne 0 then
	   return OL!b;
	end if;
    end while;

    if not found then
        error "ERROR: No normal basis element found!!!";
    end if;
end intrinsic;

intrinsic numden(xi :: FldOrdElt, P :: RngOrdIdl) -> RngOrdElt, RngOrdElt
{}

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