/*
 * globalrep.m
 * (c) 2008-2011 Ruben Debeerst
 *
 * Methods to compute global representations of local Galois extensions
 * heuristically.
 *
 *
 * Methods provided:
 *
 * GlobalRepresentations
 * GlobalRepresentationsJR
 *     Find global representations for local extensions in the database
 *     of Kluener/Malle, either starting with specific Galois groups and
 *     using the command AllExtensions or starting with a list of local polynomials
 *     from the database by Jones/Roberts.
 *
 * kluenersMallePols
 *     Send an internet-request to retrieve polynomials with given
 *     Galois group from the database by Klueners/Malle.
 *
 * jonesRobertsPols
 *     Send an internet-request to retrieve polynomials
 *     with given Galois group.
 *
 * genericC4Pol
 * genericD4Pol
 * randomD4Pol
 *     Return a polynomial of Galois Group C4 or D4.
 *
 * embeddingC2C4
 *     Embedd a field of Galois group C2 into a field
 *     of Galois group C4, if possible.
 *
 *
 * helper functions:
 *
 * allExtensionsForGroup
 *     Compute all local extension with given Galois group.
 *
 * sumOfSquares
 *     Write an integer as a sum of squares, if possible.
 *
 * isUndecomposed
 *     Check whether the number field generated by f is undecomposed at p.
 *
 */



// local functions
forward automorphism_group;
forward find_global_representations;


declare verbose GlobalRepresentations, 1;

/**************************************************
 *           Global representations
 **************************************************/

intrinsic GlobalRepresentations(G::Grp, p::RngIntElt : JR := false, candlist := []) -> .
{ Given a Galois group G and a prime p or a list of tuples <G,p>.
  For each tuple compute all local extensions of degree #G of Qp,
  and search for global representations using the database
  by Klueners/Malle.
  Also shows corresponding polynomials from the database by 
  Jones/Roberts if JR is set to true.
}
    return GlobalRepresentations([<G,p>] : JR:=JR, candlist:=candlist)[1];
end intrinsic;

intrinsic GlobalRepresentations(list::SeqEnum[Tup] : JR := false, candlist := []) -> .
{ Given a Galois group G and a prime p or a list of tuples <G,p>.
  For each tuple compute all local extensions of degree #G of Qp,
  and search for global representations using the database
  by Klueners/Malle.
  Also shows corresponding polynomials from the database by 
  Jones/Roberts if JR is set to true.
}
    allpols := [];
    allExt := [];
    lastext := <0,0>;
    for x in list do
        G := x[1];
        p := x[2];
        vprintf GlobalRepresentations, 1: "G = %o, p = %o\n", IdentifyGroup(G), p;
        vprint  GlobalRepresentations, 1: "====================";
        if lastext ne <#G,p> then
            // alle lokalen Erweiterungen eines bestimmten Grades
            vprintf GlobalRepresentations, 1: "Compute extensions of degree %o of Q_%o ... ", #G,p;
            vtime   GlobalRepresentations, 1: allExt := AllExtensions(pAdicField(p,100), #G : Galois);
            vprintf GlobalRepresentations, 1: "%o extensions\n", #allExt;
            vprintf GlobalRepresentations, 1: "Find Galois extensions ... ", #G,p;
            vtime   GlobalRepresentations, 1: allExt := [E : E in allExt | #Automorphisms(E,pAdicRing(E)) eq Degree(E,pAdicRing(E)) ];
            vprintf GlobalRepresentations, 1: "%o extensions\n", #allExt;
            vprint  GlobalRepresentations, 1: "--------------------";
            lastext := <#G, p>;
        end if;
        if IsCyclic(G) then
            vprint  GlobalRepresentations, 1: "cyclic group";
        elif IsAbelian(G) then
            vprint  GlobalRepresentations, 1: "abelian group";
        end if;
        
        vprintf GlobalRepresentations, 1: "Select extensions for group ... ";
        vtime   GlobalRepresentations, 1: ext := allExtensionsForGroup(G, p : ext := allExt, otherAutGrp);
        vprintf GlobalRepresentations, 1: "%o extensions\n", #ext;
        if #ext gt 0 then
            pols := GlobalRepresentations(ext, G : JR:=JR, candlist:=candlist);
            vprint  GlobalRepresentations, 1: [ Type(f) eq SeqEnum select f[1] else f : f in pols];
        else
            pols := [];
        end if;
        Append(~allpols, pols);
        vprint  GlobalRepresentations, 1: "";
    end for;
    return allpols;
end intrinsic;

intrinsic GlobalRepresentations(ext::SeqEnum, G::Grp : JR := false, candlist := [] ) -> SeqEnum
{ Given a list of local extensions which have Galois group G.
  Search for global representations using the database by Klueners/Malle.
  Also shows corresponding polynomials from the database by 
  Jones/Roberts if JR is set to true.
}
    Zx<x> := PolynomialRing(Integers());
    if #ext gt 0 then
        p := Prime(ext[1]);
        polsJR := [];
        if #G le 11 and JR then
            vprint GlobalRepresentations, 1: "Corresponding Jones/Roberts polynomials";
            polsJR := jonesRobertsPols(#G,p);
            polsJR := [ [ f : f in polsJR | #Roots(f,E) eq #G]   : E in ext ];
            if &and([ #f eq 1: f in polsJR]) then
                polsJR :=  [f[1] : f in polsJR];
            end if;
            vprint  GlobalRepresentations, 1: polsJR;
            
        elif #G gt 11 then
            vprint  GlobalRepresentations, 1: "No corresponding Jones/Roberts polynomials for degree > 11";
        end if;
        // Polynome fuer Gruppe
        if #candlist eq 0 then
            vprintf GlobalRepresentations, 1: "Search for Klueners/Malle polynomials ... ";
            vtime   GlobalRepresentations, 1: polsG := kluenersMallePolsG(G);
            polsG := [ f : f  in polsG | isUndecomposed(f,p)  ];
        else
            vprint GlobalRepresentations, 1: "Using candidate polynomials ";
            polsG := candlist;
        end if;
        pols := find_global_representations(ext, polsG);
        if #polsJR eq #pols then
            return [ [ Zx!pols[i], Zx!polsJR[i] ] : i in [1..#pols]];
        elif JR then
            return [ [ Zx!pols[i], Zx!0 ] : i in [1..#pols]];
        else
            return pols;
        end if;
    end if;
    return [];
end intrinsic;



intrinsic GlobalRepresentationsJR(pols::List, p::RngIntElt) -> .
{ Given a list of polynomials in format of the database by Jones/Roberts,
  representing extensions of degree n of Qp.
  For each Galois group of this degree, select corresponding polynomials
  from the list and search for global representations using the database
  by Klueners/Malle.
}
    return GlobalRepresentationsJR(pols, Degree(pols[1,6]), p);
end intrinsic;

intrinsic GlobalRepresentationsJR(pols::List, n::RngIntElt, p::RngIntElt) ->.
{ Given a list of polynomials in format of the database by Jones/Roberts,
  representing extensions of degree n of Qp.
  For each Galois group of this degree, select corresponding polynomials
  from the list and search for global representations using the database
  by Klueners/Malle.
}
    allPols := [];
    Zx<x> := PolynomialRing(Integers());
    for i in [1..NumberOfSmallGroups(n)] do
        vprintf GlobalRepresentations, 1: "JR group id <%o,%o>\n", n,i;
        vprint  GlobalRepresentations, 1: "====================";
        vprintf GlobalRepresentations, 1: "Compute extensions for group ... ";
        polsGlocal := [s : s in pols | <u[1],u[2]> eq <n,i>  where u := s[7,3] ];
        vprint  GlobalRepresentations, 1: [ f[6] :f in polsGlocal];
        vtime   GlobalRepresentations, 1: ext := [
            ChangePrecision(Completion(L, Factorization(p*RingOfIntegers(L))[1,1]),300)
            where L := NumberField(f) where f := s[6] : s in polsGlocal];
        vprintf GlobalRepresentations, 1: "%o extensions\n", #ext;
        if #ext gt 0 then
            vprintf GlobalRepresentations, 1: "Identify Galois group ... ";
            G := AutomorphismGroup(ext[1], pAdicField(ext[1]));
            vprint  GlobalRepresentations, 1: IdentifyGroup(G);
            vprintf GlobalRepresentations, 1: "Search for Klueners/Malle polynomials ... ";
            vtime   GlobalRepresentations, 1: polsGglobal := kluenersMallePolsG(G);
            vprintf GlobalRepresentations, 1: "Select undecomposed extensions ... ";
            vtime   GlobalRepresentations, 1: polsGglobal := [ f : f  in polsGglobal | isUndecomposed(f,p)  ];
            vprintf GlobalRepresentations, 1: "Match with extensions ... ";
            vtime   GlobalRepresentations, 1: extpols := find_global_representations(ext, polsGglobal);
            vprint  GlobalRepresentations, 1: extpols;
            
            Append(~allPols, [ [ Zx!extpols[i], Zx!polsGlocal[i,6] ] : i in [1..#polsGlocal]] );
        else
            Append(~allPols, []);
        end if;
    end for;
    return allPols;
end intrinsic;

/* find_global_representations(localExt::SeqEnum, globalPols::SeqEnum) -> SeqEnum
 * Given a list of local extensions and a list of polynomials, returns a list
 * of polynomials (subset of globalPols) which represent the given extensions
 * (as far as possible).
 */
function find_global_representations(localExt, globalPols)
    local globPols, R, E, g;
    
    R := PolynomialRing(Rationals());
    globPols := [];
    for E in localExt do
        for g in globalPols do
            if #Roots(R!g, E) eq Degree(g) then
                Append(~globPols, g);
                continue E;
            end if;
        end for;
        Append(~globPols, R!0);
    end for;
    
    return globPols;
end function;


/**************************************************
 *           All local extensions
 **************************************************/

intrinsic allExtensionsForGroup(G::., p::RngIntElt : precision := 100 , ext := [], otherAutGrp := false) -> SeqEnum
{ Compute all extensions of Qp usgin AllExtensions and select those which
  have the given Galois group. If a list ext of extensions is given,
  this list is being searched for suitable extensions.
  
  If otherAutGrp is true, a replacement for AutomorphismGroup is
  used to compute the automorphism group. }
  
    local Qp, grp, GId, n, j, i;
    
    if Type(G) ne Tup then
        GId := IdentifyGroup(G);
        n := #G;
    else
        GId := G;
        n := GId[1];
    end if;
    
    Qp := pAdicRing(p,precision);
    if #ext eq 0 then
        // all extension of degree #G
        ext := AllExtensions(Qp, n : Galois := true);
    end if;
    // only Galois extensions
    ext := [E : E in ext | #Automorphisms(E,pAdicRing(E)) eq Degree(E,pAdicRing(E)) ];
    // group identifiers
    if otherAutGrp then
        // use this for local extensions with infinite precision
        grp := [ IdentifyGroup(automorphism_group(E,pAdicRing(E))) : E in ext];
    else
        grp := [ IdentifyGroup(AutomorphismGroup(E)) : E in ext];
    end if;
    // select matching extensions
    ext := [ ext[j] :  j in [ i :  i in [1..#grp] | grp[i] eq GId ] ];
    // remove isomorphic extensions
    // removeIsomorphic(ext, Qp);
    return ext;
end intrinsic;



/**************************************************
 *           Database requests
 **************************************************/

intrinsic kluenersMallePols(d::RngIntElt, t::RngIntElt) -> SeqEnum
{ Get all polynomials of degree d with Galois group identifier <d,t>
  from the database by Klueners/Malle. Note that the identifier of
  Magma does not always agree with the identifier of Klueners/Malle.
  Depends on an internet connection and the Unix wget command. }
  
    local x, pols, i;
    
    x := Pipe(Sprintf("wget -q -O - \"http://reh.math.uni-duesseldorf.de/cgi-klueners/groups3.pl?deg=%o&t=%o&zer=&bet=&teil=&noteil=&maxteil=&numprime=\" | grep \"td>x\" | sed 's/.*td>x/x/g; s/<.*//g; s/^x/,1*x/; s/-x/-1*x/g; s/+x/+1*x/g';", d, t), "r");
    
    x := Pipe(Sprintf("wget -q -O - \"http://reh.math.uni-duesseldorf.de/cgi-klueners/groups3.pl?deg=%o&t=%o\""
       * " | grep \"td>x\" "        // Zeilen mit Polynom isolieren
       * " | sed 's/.*td>x/x/g; "   // unwichtige Zeichen loeschen
               * "s/<.*//g;"        // dto
               * "s/^x/,1*x/;"      // vorne Komma und 1-Faktor einfuegen
               * "s/-x/-1*x/g;"     // -x durch -1*x ersetzen
               * "s/+x/+1*x/g';"    // +x durch +1*x ersetzen
        , d, t), "r");
    
    x := Split(x, ",");
    pols := [];
    for y in x[2..#x] do
        // nullter Koeffizient
        bool, _, match := Regexp("([+-][0-9]+)\n$", y);
        if bool then
            coeff := [StringToInteger(match[1])];
        else
            bool, _, match := Regexp("^([0-9]+)\n$", y);
            if bool then
                coeff := StringToInteger(match[1]);
            else 
                coeff := [0];
            end if;
        end if;
        // erster Koeffizient
        bool, _, match := Regexp("([+-][0-9]+)\\*x[^\\^]", y);
        if bool then
            Append(~coeff, StringToInteger(match[1]));
        else
            bool, _, match := Regexp("^([0-9]+)\\*x[^\\^]", y);
            if bool then
                Append(~coeff, StringToInteger(match[1]));
            else 
                Append(~coeff, 0);
            end if;
        end if;
        // weitere Koeffizienten
        for i in [2..d] do
            bool, _, match := Regexp(Sprintf("([+-][0-9]+)\\*x\\^%o[^0-9]", i), y);
            if bool then
                Append(~coeff, StringToInteger(match[1]));
            else
                bool, _, match := Regexp(Sprintf("^([0-9]+)\\*x\\^%o[^0-9]", i), y);
                if bool then
                    Append(~coeff, StringToInteger(match[1]));
                else 
                    Append(~coeff, 0);
                end if;
            end if;
        end for;
        
        Append(~pols, coeff);
    end for;
    
    //PrintFile("kluenersDBPols.tmp", x);
    //load "kluenersDBPols.tmp";
    //Pipe("rm " cat tmpfile, "r");
    
    return [ Polynomial(y) : y in pols ];
end intrinsic;

intrinsic kluenersMallePolsG(G::Grp) -> SeqEnum
{ Get all polynomials with Galois group G from the database by Klueners/Malle.
  Depends on an internet connection and the Unix wget command. }
    polsG := kluenersMallePols(#G,IdentifyGroup(G)[2]);
    i := 0;
    while not IdentifyGroup(AutomorphismGroup(NumberField(polsG[1]))) eq IdentifyGroup(G) do
        i := i+1;
        vprintf GlobalRepresentations, 1: "Test Group Identifier <%o,%o>\n", #G, i;
        polsG := kluenersMallePols(#G,i);
    end while;
    
    return polsG;
end intrinsic;


intrinsic jonesRobertsPols(n::RngIntElt, p::RngIntElt) -> SeqEnum
{ Get polynomials generating all extensions of degree n of Qp
  from the database by Jones/Roberts.
  Depends on an internet connection and the Unix wget command. }

    x := Pipe(Sprintf("wget -q -O - \"http://hobbes.la.asu.edu/LocalFields/basic-table.cgi?degree=%o&prime=%o\""
       * " | grep \"<\\!-- a\" "    // Zeilen mit Polynom isolieren
       * " | sed 's/.*|.*|//g; "    // unwichtige Zeichen loeschen
               * "s/-->//g;"        // dto
               * "s/ //g;"          // dto
               * "s/^x/,1*x/;"      // vorne Komma und 1-Faktor einfuegen
               * "s/-x/-1*x/g;"     // -x durch -1*x ersetzen
               * "s/+x/+1*x/g"    // +x durch +1*x ersetzen
        * "'", n, p), "r");
    
    x := Split(x, ",");
    pols := [];
    for y in x[2..#x] do
        // nullter Koeffizient
        bool, _, match := Regexp("([+-][0-9]+)\n$", y);
        if bool then
            coeff := [StringToInteger(match[1])];
        else
            bool, _, match := Regexp("^([0-9]+)\n$", y);
            if bool then
                coeff := StringToInteger(match[1]);
            else 
                coeff := [0];
            end if;
        end if;
        // erster Koeffizient
        bool, _, match := Regexp("([+-][0-9]+)\\*x[^\\^]", y);
        if bool then
            Append(~coeff, StringToInteger(match[1]));
        else
            bool, _, match := Regexp("^([0-9]+)\\*x[^\\^]", y);
            if bool then
                Append(~coeff, StringToInteger(match[1]));
            else 
                Append(~coeff, 0);
            end if;
        end if;
        // weitere Koeffizienten
        for i in [2..n] do
            bool, _, match := Regexp(Sprintf("([+-][0-9]+)\\*x\\^%o[^0-9]", i), y);
            if bool then
                Append(~coeff, StringToInteger(match[1]));
            else
                bool, _, match := Regexp(Sprintf("^([0-9]+)\\*x\\^%o[^0-9]", i), y);
                if bool then
                    Append(~coeff, StringToInteger(match[1]));
                else 
                    Append(~coeff, 0);
                end if;
            end if;
        end for;
        
        Append(~pols, coeff);
    end for;

    return [ Polynomial(y) : y in pols ];
    
end intrinsic;






/**************************************************
 *           Generic Polynomials
 **************************************************/

intrinsic genericC4Pol(s::FldRatElt, t::FldRatElt) -> RngUPolElt
{ Returns the generic C4-Polynomial for s and t from [Jensen et al.: Generic Polynomials, Cor. 2.2.6].
  The given polynomial generates a C4-extension if s!=0 and 1+t^2 is not a square. }
    return Polynomial(Rationals(), [s^2*t^2*(1+t^2), 0, -2*s*(1+t^2), 0, 1]);
end intrinsic;

intrinsic genericC4Pol(s::FldNumElt, t::FldNumElt) -> RngUPolElt
{ Returns the generic C4-Polynomial for s and t from [Jensen et al.: Generic Polynomials, Cor. 2.2.6].
  The given polynomial generates a C4-extension if s!=0 and 1+t^2 is not a square. }
    return Polynomial(Parent(s), [s^2*t^2*(1+t^2), 0, -2*s*(1+t^2), 0, 1]);
end intrinsic;

intrinsic genericD4Polynomial(a::., b::. ) -> RngUPolElt
{ If b and b(a^2-4b) are both not square, the polynomial
  f=b X^4+a X^2+1 which generates a D4 extension is returned.
  Otherwise an error occurs.
  See [Jensen et al.: Generic Polynomials, Cor. 2.2.4]. }
    
    local f, K;
    
    require Type(a) in [FldNumElt, FldRatElt] and
            Type(b) in [FldNumElt, FldRatElt] and
            Type(a) eq Type(b) and
            b in Parent(a):
            "Bad argument types";
    
    if IsSquare(b) then
        error "b is a square!";
    end if;
    if IsSquare(b*(a^2-4*b)) then
        error "b*(a^2-4*b) is a square!";
    end if;
    
    K := Parent(a);
    f := Polynomial(K, [b, 0, a, 0, 1]);
    if not IsIrreducible(f) then
        error "The polynomial is not irreducible!";
    end if;
    
    return f;
end intrinsic;

intrinsic randomD4Polynomial(K::., bound::RngIntElt : maxTries := 5) -> RngUPolElt
{ Computes a random polynomial generating a D_4 extension over K. }
    
    local d, f;
    
    require Type(K) in [FldNum, FldRat] :
            "Bad argument types\nNumber Field or Rational Field required";
    
    d := Degree(K);
    
    for i in [1..maxTries] do
        a := K![Random(2*bound)-bound : i in [1..d]];
        b := K![Random(2*bound)-bound : i in [1..d]];
        try
            f := genericD4Polynomial(a,b);
            return f;
        catch e
            f := 0;
        end try;
    end for;
    error "Couldn't find D4-polynomial in", maxTries, "tries!";
end intrinsic;




/**************************************************
 *           Embedding Problems
 **************************************************/

intrinsic embeddingC2C4(K::FldNum : p:= 0 ) -> BoolElt, RngUPolElt
{ Computes a generating polynomial for a C4-Extension L/Q which includes K/Q, [K:Q]=2.
  If p is specified, L will be unramified and undecomposed at p.
  L can either be created as absolute field over Q or relative over K.
  [Jensen et al.: Generic Polynomials, Thm. 2.2.5] }
  
    local f, c, a, M1, OM1, S, z, G, psi, sigma, w, g;
    
    require AbsoluteDegree(K) eq 2 :
            "K must be of degree 2";
    
    if p ne 0 then
        // ensure that K is undecomposed unramified at p
        require #DecompositionType(K,p) eq 1 and DecompositionType(K,p)[1,2] eq 1 :
                "K must be undecomposed and unramified at p";
    end if;
    
    f := DefiningPolynomial(K);
    // Schreibe M=K(sqrt(a))
    c := Coefficient(f, 1);
    if c ne 0 then
        // x -> x-c/2
        f := Evaluate(f, Polynomial(Rationals(), [-c/2,1]));
    end if;
    // Nenner im Polynom loswerden
    c := Denominator(Coefficient(f,0));
    if c ne 1 then
        bool, src := IsSquare(c);
        if bool then
            c := src;
        end if;
        // x -> x/c
        f:=c^2*Evaluate(f,Polynomial(Rationals(), [0, 1/c]) );
    end if;
    // Jetzt gilt: f=X^2-a
    a := Integers()!-Coefficient(f,0);
    
    bool, x, y := sumOfSquares(a);
    if not bool then
        return false, _;
    end if;
    
    // a = x^2 + y^2
    // a/x^2 = 1+(y/x)^2 = 1 + t^2
    t := y/x;
    i := 0;
    
    repeat
        // Zufallszahl in -5..5, -50..50, ...
        i := i+1;
        s := Random( (i div 10 + 1) * 10 );
        s := s-(i div 10 + 1) * 5;
        if s eq 0 then
            s := 1;
        end if;
        g := genericC4Pol(Rationals()!s,Rationals()!t);
        bool := IsIrreducible(g);
        if p ne 0 and bool then
            // test if undecomposed and unramified
            L := NumberField(g);
            bool := #DecompositionType(L,p) eq 1 and DecompositionType(L,p)[1,2] eq 1;
        end if;
    until bool;
    
    // make monic and integral
    c := LCM([Denominator(x) : x in Coefficients(g)]);
    if c ne 1 then
        bool, src := IsSquare(c);
        if bool then
            c := src;
            bool, src := IsSquare(c);
            if bool then
                c := src;
            end if;
        end if;
        g := c^4 * Evaluate(g, Polynomial(Rationals(), [0, 1/c]));
    end if;
    
    return true, g;
end intrinsic;

intrinsic sumOfSquares(a::RngIntElt) -> BoolElt, RngIntElt, RngIntElt
{ Returns true if a can be written as a=x^2+y^2 for x,y\in\Z, and also x,y if so. }
    local bool, x, y, i;
    
    if a lt 0 then
        return false, _, _;
    elif a eq 0 then
        return true, Integers()!0, Integers()!0;
    end if;
    
    for i in [1..Ceiling(SquareRoot(a))] do
        bool, y := IsSquare(a-i^2);
        if bool then
            x := i;
            return bool,x,y;
            break;
        end if;
    end for;
    
    return false, _, _;
end intrinsic;



/**************************************************
 *                   Others
 **************************************************/

intrinsic isUndecomposed(f::RngUPolElt, p::RngIntElt) -> BoolElt
{ Check whether the number field generated by f is undecomposed at p. }
    return (#DecompositionType(NumberField(f), p) eq 1);
end intrinsic;



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;

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;


