/*
 * artin.m
 * (c) 2008-2009 Ruben Debeerst
 *
 * Methods to compute the Artin map
 *
 * Depends On:
 *   RelAlgKTheory.m (Approximation, MakeTotallyPositive, RealSign, Uniformizer)
 *
 *
 * Methods provided:
 * 
 * localNormResidueSymbol
 *     Compute the norm residue symbol (x, N/M) for a local field
 *     extension N/M. 
 *
 * localNormResidueSymbolAsGlobalIdeal
 *     Computes a global ideal a such that its artin symbol (a, L/K)
 *     represents the local norm residue symbol (x, N/M).
 *
 * globalArtinSymbol
 *     Compute the norm residue symbol (x, L/K) for a global field
 *     extension L/K.
 *     
 */


declare verbose artin, 3;
SetVerbose("artin", 0);

// local functions
forward loc_normres_symb, 
        loc_normres_symb_glob_idl, loc_normres_symb_glob_idl_Z,
        glob_artin_symb, glob_artin_symb_Q;
forward approximation_Z, fast_decomposition_group, decomposition_group, primitive_crt,
        primitive_crt_seq, sort_factorization, sort_factorization_Z,
        test_congruence;


intrinsic localNormResidueSymbol(x::., N::FldNum, M::FldNum, PM::.) -> GrpElt, FldAb
{ Let N/M be a global abelian extension and x in M^\times and PM an ideal of M,
  such that there is just one prime ideal PN in N above PM.
  Compute the local norm residue symbol (x, N_PN/M_PM) in Gal(N/M).
  The extension N/M can also be given as abelian field. If the map
  psiH: Gal(N/M) -> Aut(N/M) is not given, the automorphism group of N/M is
  computed first.}
    
    local N2, H, psiH, Na, a;
    require (Type(x) eq FldNumElt and Type(PM) eq RngOrdIdl) or
            (Type(x) eq FldRatElt and Type(PM) eq RngInt) :
            "Expected: x in M and PM ideal in M.";
    
    if BaseField(N) eq M then
        N2 := N;
    else
        vprint artin, 1: "Compute relative field";
        N2 := ext<M|Factorization(PolynomialRing(M)!DefiningPolynomial(N))[1,1]>;
    end if;
    vprint artin, 1: "Compute automorphism group";
    H,_, psiH := AutomorphismGroup(N2,M);
    Na := AbelianExtension(N2);
    
    a := loc_normres_symb(M!x, Na, PM, psiH);
    return a, Na;
end intrinsic;

intrinsic localNormResidueSymbol(x::., Na::FldAb, PM::.) -> GrpElt
{ Let N/M be a global abelian extension and x in M^\times and PM an ideal of M,
  such that there is just one prime ideal PN in N above PM.
  Compute the local norm residue symbol (x, N_PN/M_PM) in Gal(N/M).
  The extension N/M can also be given as abelian field. If the map
  psiH: Gal(N/M) -> Aut(N/M) is not given, the automorphism group of N/M is
  computed first.}
    
    local N,M, H, psiH, a;
    require (Type(x) eq FldNumElt and Type(PM) eq RngOrdIdl) or
            (Type(x) eq FldRatElt and Type(PM) eq RngInt) :
            "Expected: x in M and PM ideal in M.";
    
    N := NumberField(Na);
    M := BaseField(Na);
    H,_, psiH := AutomorphismGroup(N,M);
    return loc_normres_symb(x, Na, PM, psiH);
end intrinsic;

intrinsic localNormResidueSymbol(x::., Na::FldAb, PM::., psiH::Map : NoArchPlaces:=false) -> GrpElt
{ Let N/M be a global abelian extension and x in M^\times and PM an ideal of M,
  such that there is just one prime ideal PN in N above PM.
  Compute the local norm residue symbol (x, N_PN/M_PM) in Gal(N/M).
  The extension N/M can also be given as abelian field. If the map
  psiH: Gal(N/M) -> Aut(N/M) is not given, the automorphism group of N/M is
  computed first.}
    
    local N, H, a, F;


    // print "In localNormResidueSymbol with x = ", x;
    require (Type(x) eq FldNumElt and Type(PM) eq RngOrdIdl) or
            (Type(x) eq FldRatElt and Type(PM) eq RngInt) :
            "Expected: x in M and PM ideal in M.";
    return loc_normres_symb(x, Na, PM, psiH : NoArchPlaces:=NoArchPlaces);
end intrinsic;

// [Exact Algo. (alt), sec 2.3, S. 5]
intrinsic localNormResidueSymbolAsGlobalIdeal(alpha::FldNumElt, F::SeqEnum, PK::RngOrdIdl : NoArchPlaces:=false) -> RngOrdIdl
{ Given the factorization F of the Artin conductor of an abelian extension L/K,
  an element alpha in K and an ideal PK of K sucht that there is just one prime
  ideal PL of L above PK.
  Compute an ideal a of K such that the global Artin symbol (a, L/K) is equal to
  the local norm residue symbol (alpha, L_PL/K_PK). }
    
    require alpha in FieldOfFractions(Order(PK)) :
            "The element alpha and the ideal PK must have the same parent!";
    return loc_normres_symb_glob_idl(alpha, F, PK : NoArchPlaces:=NoArchPlaces);
end intrinsic;

intrinsic localNormResidueSymbolAsGlobalIdeal(alpha::FldRatElt, F::SeqEnum, PK::RngInt : NoArchPlaces:=false) -> RngOrdIdl
{ Given the factorization F of the Artin conductor of an abelian extension L/K,
  an element alpha in K and an ideal PK of K sucht that there is just one prime
  ideal PL of L above PK.
  Compute an ideal a of K such that the global Artin symbol (a, L/K) is equal to
  the local norm residue symbol (alpha, L_PL/K_PK). }

    require alpha in Rationals() :
            "The element alpha must be a rational number!";
    return loc_normres_symb_glob_idl_Z(alpha, F, PK);
end intrinsic;

// [Exact Algo. (alt), sec 2.3, S. 6]
intrinsic globalArtinSymbol(a::RngOrdFracIdl, psi::Map) -> GrpElt
{ For an abelian extension L/K, an ideal a in K and psi: Gal(L/K) --> Aut(L).
  Compute the Artin symbel (a, L/K) in Gal(L/K). }
    return glob_artin_symb(a, psi);
end intrinsic;

// Wie globalArtinSymbol fuer K=\Q
intrinsic globalArtinSymbol(a::RngInt, psi::Map) -> GrpElt
{ For an abelian extension L/K, an ideal a in K and psi: Gal(L/K) --> Aut(L).
  Compute the Artin symbel (a, L/K) in Gal(L/K). }
    return glob_artin_symb_Q(a, psi);
end intrinsic;

/* This is essentially the same as Ruben's function glob_artin_symb. However, we speed up the computation by using
   precomputed data: (i) We save all computed values g(b) for g in Gal(N/M) and b in Basis(ON).
                     (ii) We save all computed Frobenius elements.

   One should try other, probably faster methods as described in the comments of G63_galois_gauss_sums.

   We also use a new function fast_decomposition_group instead of decomposition_group.
*/
intrinsic FastLocalNormResidueSymbol(a :: FldNumElt , F :: SeqEnum, PM ::RngOrdIdl , psiH :: Map, data :: List) -> GrpPermElt, List
{}

    B := data[1]; index_set := data[2]; precomp_data := data[3]; list_of_p := data[4]; list_of_frob := data[5];
    a := localNormResidueSymbolAsGlobalIdeal(a, F, PM : NoArchPlaces:=true);
    vprint artin, 1: "Compute global artin symbol for ideal";

    H := Domain(psiH);
    N := Domain(Codomain(psiH)[1]);

    ON := MaximalOrder(Parent(B[1]));
    OK := Order(a);

    h := Id(H);
    vprintf artin, 2: "Factorization of ideal";
    factors := Factorization(a);
    //print factors;
    vprintf artin, 2: ": %o factors\n", #factors;
    for fac in factors do
        vprintf artin, 3: "factor: %o", fac;
        p := fac[1];
        factors := Factorization(ideal<ON | {ON!b : b in Generators(p) } >);
        P := factors[1,1];
        Q, s1 := quo<ON | P>;

        // time GD1 := decomposition_group(P,psiH);
        ind := Index(list_of_p, p);
        if ind ne 0 then        
            frob := list_of_frob[ind];
        else
            GD := fast_decomposition_group(P,psiH); 
            q := Index(OK,p);

            candidates := {g : g in GD};
            j := 2;
            while #candidates ne 1 do
                b := B[j];
                cand := candidates;
                for g in candidates do
                    ind := Index(index_set, <g,b>);
                    if precomp_data[ind] ne 0 then
                        y := precomp_data[ind];
                    else
		        y := psiH(g)(b);
		        precomp_data[ind] := y;
		    end if;
                    if s1(y)-s1(b)^q ne Zero(Q) then
                        cand:= cand diff {g};
                    end if;
                end for;
                candidates := cand;
                j := j+1;
            end while;
            frob := SetToSequence(candidates)[1];
	    Append(~list_of_p, p); Append(~list_of_frob, frob);
	end if;
        h := h*(frob^fac[2]); 
    end for;

    vprintf artin, 2: "\n";
    return h, [* B, index_set, precomp_data, list_of_p, list_of_frob *];
end intrinsic;


/* Same as above for M being the rationals.
*/
intrinsic FastLocalNormResidueSymbol(a :: FldRatElt , F :: SeqEnum, PM :: RngInt , psiH :: Map, data :: List) -> GrpPermElt, List
{}

    B := data[1]; index_set := data[2]; precomp_data := data[3]; list_of_p := data[4]; list_of_frob := data[5];
    a := localNormResidueSymbolAsGlobalIdeal(a, F, PM : NoArchPlaces:=true);
    vprint artin, 1: "Compute global artin symbol for ideal";

    H := Domain(psiH);
    N := Domain(Codomain(psiH)[1]);

    ON := MaximalOrder(Parent(B[1]));
    OK := Order(a);

    h := Id(H);
    vprintf artin, 2: "Factorization of ideal";
    factors := Factorization(a);
    //print factors;
    vprintf artin, 2: ": %o factors\n", #factors;
    for fac in factors do
        vprintf artin, 3: "factor: %o", fac;
        p := fac[1];
        factors := Factorization(ideal<ON | {ON!b : b in Generators(p) } >);
        P := factors[1,1];
        Q, s1 := quo<ON | P>;

        ind := Index(list_of_p, p);
        if ind ne 0 then
	    frob := list_of_frob[ind];
        else
            GD := fast_decomposition_group(P,psiH);
            q := Abs(Generator(p));
            candidates := {g : g in GD};
            j := 2;
            while #candidates ne 1 do
                b := B[j];
                cand := candidates;
                for g in candidates do
                    ind := Index(index_set, <g,b>);
                    if precomp_data[ind] ne 0 then
                        y := precomp_data[ind];
                    else
                        y := psiH(g)(b);
                        precomp_data[ind] := y;
                    end if;
                    if s1(y)-s1(b)^q ne Zero(Q) then
                        cand:= cand diff {g};
                    end if;
                end for;
                candidates := cand;
                j := j+1;
            end while;
            frob := SetToSequence(candidates)[1];
            Append(~list_of_p, p); Append(~list_of_frob, frob);
        end if;
        h := h*(frob^fac[2]);
    end for;
    vprintf artin, 2: "\n";
    return h, [* B, index_set, precomp_data, list_of_p, list_of_frob *];
end intrinsic;


/***************************************************
 *               Local functions
 ***************************************************/
 
/* loc_normres_symb(x::., Na::FldAb, PM::., psiH::Map) -> GrpElt
 * See localNormResidueSymbol.
 */
function loc_normres_symb(x, Na, PM, psiH : NoArchPlaces:=false)
    local H, N, F, a;
    
    H := Domain(psiH);
    if Type(Codomain(psiH)) eq SeqEnum then
        N := Domain(Codomain(psiH)[1]);
    else
        assert Type(Codomain(psiH)) eq PowMapAut;
        N := Domain(Codomain(psiH));
    end if;
    vprint artin, 1: "Factorization of conductor";
    F := Factorization(Conductor(Na));
    vprint artin, 1: "Compute ideal representing local norm residue symbol";
    a := localNormResidueSymbolAsGlobalIdeal(x, F, PM : NoArchPlaces:=NoArchPlaces);
    vprint artin, 1: "Compute global artin symbol for ideal";
    a := globalArtinSymbol(a, psiH);
    return a;
end function;

/* loc_normres_symb_glob_idl(alpha::FldNumElt, F::SeqEnum, PK::RngOrdIdl) -> RngOrdIdl
 * See localNormResidueSymbolAsGlobalIdeal.
 */
function loc_normres_symb_glob_idl(alpha, F, PK : NoArchPlaces:=false)
    local PP, P, pi, e, xi, a, OK, K, I;
    // PK ist erster Faktor in dieser Liste (moeglicherweise mit Exponent 0).
    PP := sort_factorization(F, PK);
    
    OK := Order(PK);
    K := FieldOfFractions(OK);
    //K := Parent(alpha);
    
    e := Valuation(alpha, PK);
    if NoArchPlaces then
        ArchPlaces := [];
    else
        ArchPlaces := RealPlaces(K);
    end if;
    
    if #PP gt 1 then
        pi := primitive_crt_seq([P[1] : P in PP]);
        //pi := Uniformizer([P[1] : P in PP], &*([P[1] : P in PP]))[1];
        if Signature(K) ne 0 and not NoArchPlaces then
            pi := MakeTotallyPositive(K!pi, PP[1,1]^2* &*([P[1] : P in PP]) );
        end if;
        //print [Valuation(pi, P[1]) : P in PP];
        // Val_PK(pi)=1, sonst Val_P(pi)=0
        // und pi is totally positive
        xi := Approximation([K!(alpha/pi^e)] cat [K!(1/pi^e)  : i in [1..#PP-1] ],
                             [<P[1], Max(1,P[2]) > : P in PP], &*([P[1]^Max(1,P[2]) : P in PP]),
                             ArchPlaces);
        // jetzt ist xi/pi^e totally positive
        
    else
        pi := PrimitiveElement(PP[1,1]);
        if Signature(K) ne 0 and not NoArchPlaces then
            pi := MakeTotallyPositive(K!pi, PP[1,1]^2);
        end if;
        xi := Approximation([K!(alpha/pi^e)], [<PP[1,1],Max(1,PP[1,2])>], PP[1,1]^Max(1,PP[1,2]), ArchPlaces);
    end if;

    assert xi in OK;
    
    // Bedingungen testen
    // print "signs2: ", RealSign(xi/pi^e);
    if not NoArchPlaces and (Signature(K) ne 0) and (not {r gt 0 : r in RealEmbeddings(xi/pi^e)} eq {true}) then
        error "Real Embeddings sind nicht alle positiv!";
    end if;
    assert test_congruence(OK!xi, pi^e/alpha, PK, Max(1,PP[1,2]));
    for P in Remove(PP,1) do
        assert test_congruence(OK!xi, pi^e, P[1], P[2]);
    end for;
    // xi prim zu PK
    assert Valuation(xi, PK) eq 0;
    
    Factors := Factorization(pi*OK);
    a := (OK!xi) * &*([1*OK] cat [P[1]^(-e*P[2]) : P in Factors | not P[1] eq PK]);
    
    // a prim zu PK
    assert OK!1 in a + PK;
    
    return a;
end function;

/* loc_normres_symb_glob_idl_Z(alpha::FldRatElt, F::SeqEnum, PK::RngInt) -> RngOrdIdl
 * See localNormResidueSymbolAsGlobalIdeal; for ideals of Z.
 */
function loc_normres_symb_glob_idl_Z(alpha, F, PK)
    local PP, a, primes, pi, num, den, Q, s1, Qmal, s2;
    
    PP := sort_factorization_Z(F, PK);
    primes := [Generators(P[1])[1] : P in PP];
    
    e := Valuation(alpha, PK);
    pi := primes[1];
    if pi lt 0 then
        pi := -pi;
    end if;
    // print "Before approximation_Z.";
    xi := approximation_Z([Rationals()!(alpha/pi^e)] cat [Rationals()!(1/pi^e)  : i in [1..#PP-1] ],
                          [<P[1], Max(1,P[2]) > : P in PP], &*([P[1]^Max(1,P[2]) : P in PP]));
    assert test_congruence(Integers()!xi, pi^e/alpha, PK, Max(1,PP[1,2]));
    for P in Remove(PP,1) do
        assert test_congruence(Integers()!xi, pi^e, P[1], P[2]);
    end for;
    
    return (Integers()!xi)*Integers();
end function;

/* glob_artin_symb(a::RngOrdFracIdl, psi::Map) -> GrpElt
 * See globalArtinSymbol.
 */
function glob_artin_symb(a, psi)
    local G, L, OK, OL, factors, fac, g, B, p, P, GD, q, h, H, frob, Q, s1;
    
    G := Domain(psi);
    //L := Domain(Codomain(psi));
    if Type(Codomain(psi)) eq SeqEnum then
        L := Domain(Codomain(psi)[1]);
    else
        assert Type(Codomain(psi)) eq PowMapAut;
        L := Domain(Codomain(psi));
    end if;

    OL := RingOfIntegers(L);
    OK := Order(a);
    B := Basis(OL);
    
    h := Id(G);
    vprintf artin, 2: "Factorization of ideal";
    factors := Factorization(a);
    //print factors;
    vprintf artin, 2: ": %o factors\n", #factors;
    for fac in factors do
        vprintf artin, 3: "factor: %o", fac;
        //print "====\nFactor", fac;
        p := fac[1];
        factors := Factorization(ideal<OL | {OL!b : b in Generators(p) } >);
        P := factors[1,1];
        Q, s1 := quo<OL | P>;
        GD := decomposition_group(P,psi); /*<--*/
        q := Index(OK,p);
        //print q;
candidates := {g : g in GD};
j := 2;
while #candidates ne 1 do
    b := B[j];
    cand := candidates;
    for g in candidates do
        if s1(psi(g)(b))-s1(b)^q ne Zero(Q) then
            cand:= cand diff {g};
        end if;
    end for;
    candidates := cand;
    j := j+1;
end while;
frob := SetToSequence(candidates)[1];

/*
        for g in GD do
            for b in B do
                //if Valuation(psi(g)(b)-b^q, P) eq 0 then
                if s1(psi(g)(b))-s1(b)^q ne Zero(Q) then
                   continue g;
                end if;
            end for;
            //durchgekommen
            frob := g;
            break;
        end for;
        Cputime() - start;
        //print frob;
        assert &and([s1(psi(frob)(b))-s1(b)^q eq Zero(Q) : b in B]);
        //H := [g : g in GD | {Valuation(psi(g)(b)-b^q, P) ge 1 : b in B} eq {true}];
        //print [*#H, [g^fac[2] : g in H ]   *];
*/
        h := h*(frob^fac[2]);
    end for;
    vprintf artin, 2: "\n";
    //print times;
    return h;
end function;

/* glob_artin_symb_Q(a::RngInt, psi::Map) -> GrpElt
 * See globalArtinSymbol; for ideals of Z.
 */
function glob_artin_symb_Q(a, psi)
    local G, L, OL, factors, fac, g, B, p, P, GD, q, h, H, Q, s1;
    G := Domain(psi);
    //L := Domain(Codomain(psi));
    if Type(Codomain(psi)) eq SeqEnum then
        L := Domain(Codomain(psi)[1]);
    else
        assert Type(Codomain(psi)) eq PowMapAut;
        L := Domain(Codomain(psi));
    end if;
    OL := RingOfIntegers(L);
    B := Basis(OL);
    
    h := Id(G);
    vprintf artin, 2: "Factorization of ideal";
    factors := Factorization(a);
    vprintf artin, 2: ": %o factor(s)\n", #factors;
    for fac in factors do
        vprintf artin, 3: "factor: %o\n", fac;
        p := Generators(fac[1])[1];
        P := Factorization(p*OL)[1,1];
        //GD := DecompositionGroup(P);
        vprint artin, 3: "Compute decomposition group";
        GD := decomposition_group(P,psi);
        q := #ResidueClassField(Integers(), fac[1]);
        Q, s1 := quo<OL | P>;
        vprint artin, 3: "Find Frobenius";
        H := [g : g in GD | {s1(psi(g)(b)) eq s1(b)^q : b in B} eq {true}];
        assert #H eq 1;
        h := h*(H[1]^fac[2]);
    end for;
    vprint artin, 2: "\n";
    return h;
end function;





/***************************************************
 *              Helper functions
 ***************************************************/

/* Test the congruence xi=rest mod P^e during the computation of the local norm residue symbol.
 */
function test_congruence(xi, rest, P, e)
    if Type(P) eq RngInt then
        OK := Integers();
    else
        OK := Order(P);
    end if;
    K := FieldOfFractions(OK);
    Q, s1 := quo< OK | P^e >;
    Qmal, s2 := MultiplicativeGroup(Q);
    if Type(P) eq RngInt then
        num := Numerator(  K!(rest^(-1)));
        den := Denominator(K!(rest^(-1)));
    else
        num, den := numden(K!(rest^(-1)), P);
    end if;
    //print (s2^(-1))(s1(xi)) + (s2^(-1))(s1(num)) - (s2^(-1))(s1(den)) eq Zero(Qmal);
    
    return (s2^(-1))(s1(xi)) + (s2^(-1))(s1(num)) - (s2^(-1))(s1(den)) eq (s2^(-1))(s1(OK!1));
end function;

/* sort_factorization(Factors::SeqEnum, PK::RngOrdIdl) -> SeqEnum
 * Sort the factors of a factorization, s.t. the first prime ideal appearing
 * in the list is PK (possibly with exponent 0).
 */
function sort_factorization(Factors,PK)
    PP := [];
    for P in Factors do
        if P[1] eq PK then
            Insert(~PP, 1, P);
        else
            Append(~PP, P);
        end if;
    end for;
    if not PP[1,1] eq PK then
        Insert(~PP, 1, <PK, 0>);
    end if;
    
    return PP;
end function;

/* sort_factorization_Z(Factors::SeqEnum, PK::RngInt) -> SeqEnum
 * As above for prime ideals PK of Z.
 */
function sort_factorization_Z(Factors, PK)
    PP := [];
    for P in Factors do
        Q := P[1] meet Integers();
        if Q eq PK then
            Insert(~PP, 1, <Q, P[2]>);
        else
            Append(~PP, <Q, P[2]>);
        end if;
    end for;
    if not PP[1,1] eq PK then
        Insert(~PP, 1, <PK, 0>);
    end if;
    
    return PP;
end function;

/* approximation_Z(alpha::SeqEnum, fac::SeqEnum, IdG::RngInt) -> RngCycElt
 * Algorithm as in the RelAlgKTheory package for ideals of Z.
 */
function approximation_Z(alpha, fac, IdG)
    OK := Integers();
    uniformizer := [Generators(f[1])[1] : f in fac];
    pi := &*[ uniformizer[j]^( -Valuation(alpha[j], fac[j][1]) ) : j in [1..#fac] ] / One(OK);
    c := [];
    for j:=1 to #fac do
       Q, s1 := quo< OK | fac[j][1]^fac[j][2] >;
       Qmal, s2 := MultiplicativeGroup(Q);
       num := Numerator(alpha[j]*pi);
       den := Denominator(alpha[j]*pi);
       Append(~c, OK ! s2( (s2^(-1))(s1(den)) - (s2^(-1))(s1(num)) ) ); 
    end for;

    y := ChineseRemainderTheorem(c, [uniformizer[i]^fac[i,2] : i in [1..#fac]]);
    beta := y * pi;
    if beta lt 0 then
        beta := -beta;
    end if;
    return beta;
end function;

/* primitive_crt(P1::RngOrdIdl, P2::RngOrdIdl) -> RngOrdElt
 * Given two prime ideals P1 and P2. Compute a uniformizing element of P1
 * which ist not in P2.
 */
function primitive_crt(P1, P2)
    OL := Order(P1);
    x := PrimitiveElement(P1);
    // Valuation_P1(x)=1
    if not Valuation(x,P2) eq 0 then
        // Bestimme y in P1^2 mit y notin P2
        y := CRT(P1^2,P2, OL!0, OL!1);
        // Valuation_P1(y)=2
        // Valuation_P2(y)=0
        x := x-y;
    end if;
    
    return x;
end function;

/* primitive_crt_seq(P::SeqEnum) -> RngOrdElt
 * Given prime ideals [P1,P2,...,Pn]. Compute a uniformizing element of P1
 * which ist not in P2,...,Pn.
 */
function primitive_crt_seq(P)
    OL := Order(P[1]);
    x := PrimitiveElement(P[1]);
    // Valuation_P1(x)=1
    // x jetzt so um y abaendern, dass Bewertungen an anderen Stellen
    // auch stimmen
    
    // y hat bei P1 Bewertung >1, damit 
    // die Bewertung von x-y weiterhin 1 ist
    X := [P[1]^2];
    R := [OL!0];
    
    for i in [2..#P] do
        if Valuation(x, P[i]) eq 0 then
            // das soll so bleiben
            // y in P[i], also mit Bewertung >0 waehlen
            // dann hat x-y immer noch Bewertung 0
            Append(~X, P[i]);
            Append(~R, OL!0);
            
        else
            // y mit Bewertung 0 waehlen,
            // dann hat x-y auch Bewertung 0
            Append(~X, P[i]);
            Append(~R, OL!1);
        end if;
    end for;
    
    y := CRT(R, X);
    return x-y;
end function;

/* 
 * Let K/k be an extension of global fields with Galois group G
 * and psi: G --> Aut(K/k).
 * Compute the decomposition group of the prime ideal PK
 * as subgroup of G.
 * 
 */
function decomposition_group(PK, psi)
//intrinsic decompGroup(PN::RngOrdIdl, H::Grp, psi::Map) -> Grp
    G := Domain(psi);
    return sub< G | {sigma : sigma in G |
                        &and([ psi(sigma)(b) in PK : b in Generators(PK)])
                    } >;
end function;


/* 
 * Let K/k be an extension of global fields with Galois group G
 * and psi: G --> Aut(K/k).
 * Compute the decomposition group of the prime ideal PK
 * as subgroup of G.
 * 
 * Could still be improved by computing the lattice of subgroups first and 
 * then work subgroup by subgroup going from large groups to small groups.
 * This somehow generalizes the approach taken for cyclic groups taken below.
 * 
 */
function fast_decomposition_group(PK, psi)
    G := Domain(psi);
    if IsPrime(#G) then
        sigma := G.1;
	if &and[ psi(sigma)(b) in PK : b in Generators(PK)] then
	    return G;
	else
	    return sub<G | [One(G)] >;
	end if;
    else
        D := sub<G | [Id(G)] >;
        Elts := {g : g in G | g ne One(G)};
        while #Elts ne 0 do
            sigma := [g : g in Elts][1];
            if &and[ psi(sigma)(b) in PK : b in Generators(PK)] then
                D := sub<G | Generators(D) join {sigma}>;
                Elts := Elts diff {g : g in D};
            else
                Elts := Elts diff {sigma};
            end if;
        end while;
        return D;
    end if;
end function;

