intrinsic MyFibo(n :: RngIntElt) -> RngIntElt
{Meine Fibonaccizahlen.}
    if n eq 0 then return 0; end if;
    if n eq 1 then return 1; end if;
    a := 0; b := 1;
    for i:=2 to n do
        c := b;
        b := a+b;
        a := c;
    end for;
    return b;
end intrinsic;


intrinsic MyPrimes(N :: RngIntElt) -> SeqEnum
{Sequenz der Primzahlen kleiner oder gleich N.}
    p := 2;
    P := [];
    while p le N do
        Append(~P, p);
        p := NextPrime(p);
    end while;
   return P;
end intrinsic;

intrinsic MyPrime(N :: RngIntElt, a :: RngIntElt, m :: RngIntElt) -> RngIntElt
{Die kleinste Primzahl p >= N mit p = a mod m}
    p := NextPrime(N-1);
    found := p mod m eq a;
    while not found do
        p := NextPrime(p);
        found := p mod m eq a;
    end while;
    return p;
end intrinsic;

intrinsic MyDn(n :: RngIntElt) -> GrpPerm
{Diedergruppe der Ordnung 2n.}
    S := Sym(n);
    sigma := S ! ([2..n] cat [1]);
    tau := One(S);
    if IsEven(n) then
        for i:=1 to Integers()! (n/2) do
            tau := tau * S!(i, n-i+1);
        end for;
    else
        for i:=2 to Integers()! ((n+1)/2) do
            tau := tau * S!(i, n-i+2);
        end for;
    end if;
    return sub<S | [sigma, tau]>;
end intrinsic;


intrinsic MyQ8() -> GrpPerm
{Quaternionengruppe der Ordnung 8.}

    /* Zuerst als Matrixgruppe. */
    ZX<x> := PolynomialRing( Integers());
    F<i> := NumberField(x^2+1);
    M := GL(2,F);
    I := M![i,0,0,-i];
    J := M![0,-1,1,0];
    Q := MatrixGroup<2,F|[I,J]>;

    /* und jetzt als Permutationsgruppe */
    Qelts := [g : g in Q];
    v := [Index(Qelts, I*g) : g in Qelts];
    S8 := Sym(8);
    sigma := S8!v;
    w := [Index(Qelts, J*g) : g in Qelts];
    tau := S8!w;
    H := sub<S8 | [sigma, tau]>;
    return H;
end intrinsic;

intrinsic MyIsInvertible(a :: RngIntElt, m :: RngIntElt) -> BoolElt, RngIntElt
{Ueberprueft, ob a modulo m invertierbar ist und berechnet gegebenenfalls das Inverse.}
    d,p,q := XGCD(a,m);
    if d ne 1 then
        return false, 0;
    else
        return true, p mod m;
    end if;
end intrinsic;


intrinsic MyPrimRoot(F :: FldFin) -> FldFinElt
{Berechnet einen Erzeuger von F^x.}
    U, j := UnitGroup(F);
    return j(U.1);
end intrinsic;


intrinsic ComputeGDB(q :: RngIntElt) -> GrpMat, GrpMat, GrpMat
{Berechnet G = Gl_2(F_q) und zwei spezielle Untergruppen von G.}
    F := GF(q);
    G := GL(2, F);
    w := MyPrimRoot(F);
    D := sub<G | [ G![w,0,0,1], G![1,0,0,w] ]>;
    B := sub<G | [ G![1,1,0,1] ]>;
    return G, D, B;
end intrinsic;



intrinsic Nebenklasse(g :: GrpPermElt, U :: GrpPerm) -> SetEnum
{Berechnet die Linksnebenklasse gU.}
    return {g*u : u in U};
end intrinsic;

intrinsic Nebenklasse(U :: GrpPerm, g :: GrpPermElt) -> SetEnum
{Berechnet die Linksnebenklasse gU.}
    return {u*g : u in U};
end intrinsic;


intrinsic Linksnebenklassen(G :: GrpPerm, U :: GrpPerm) -> SeqEnum, SeqEnum
{Berechnet G/U.}
    Gset := Set(G);
    R := [Id(G)];
    Q := Nebenklasse(Id(G), U);
    Nk := [ Q ];
    while #Q lt #G do
        g := [g : g in Gset diff Q][1];
        gU := Nebenklasse(g, U);
        Q := Q join gU;
        Append(~R, g);
        Append(~Nk, gU);
    end while;
    return Nk, R;
end intrinsic;

intrinsic Rechtsnebenklassen(G :: GrpPerm, U :: GrpPerm) -> SeqEnum, SeqEnum
{Berechnet U\G.}
    Gset := Set(G);
    R := [Id(G)];
    Q := Nebenklasse(U, Id(G));
    Nk := [ Q ];
    while #Q lt #G do
        g := [g : g in Gset diff Q][1];
        gU := Nebenklasse(U, g);
        Q := Q join gU;
        Append(~R, g);
        Append(~Nk, gU);
    end while;
    return Nk, R;
end intrinsic;

intrinsic IstNormal(G :: GrpPerm, U :: GrpPerm) -> BoolElt
{Entscheidet, ob U normal in G ist.}
    for g in Generators(G) do
        if Nebenklasse(g,U) ne Nebenklasse(U,g) then
	    return false;
	end if;
    end for;
    return true;
end intrinsic;

intrinsic IstNormalFast(G :: GrpPerm, U :: GrpPerm) -> BoolElt
{Entscheidet, ob U normal in G ist.}
    for g in Generators(G), u in Generators(U) do
        if not (g^(-1)*u*g in U) then
            return false;
        end if;
    end for;
    return true;
end intrinsic;

            
intrinsic IstZentral(g :: GrpElt) -> BoolElt
{Entscheidet, ob g mit allen Gruppenelementen verstauscht.}
    for h in Generators(Parent(g)) do
        if h*g ne g*h then
	    return false;
	end if;
    end for;
    return true;
end intrinsic;


intrinsic Zentrum(G :: Grp) -> Grp
{Bestimmt das Zentrum von G.}
    return sub<G | [g : g in G | IstZentral(g)]>;
end intrinsic;


intrinsic RandomPoly(FpX :: RngUPol, n :: RngIntElt) -> RngUPolElt
{Zufaelliges normiertes Polynom vom Grad n ueber dem endlichen Koerper F.}
    x := FpX.1;
    Fp := CoefficientRing(FpX);
    q := x^n;
    for i:=0 to n-1 do
        q := q + Random(Fp)*x^i;
    end for;
    return q;
end intrinsic;


intrinsic IrrPoly(FpX :: RngUPol, n) -> RngUPolElt
{Sucht ein irreduzibles Polynom vom Grad n.}
    q := RandomPoly(FpX, n);
    cnt := 1;
    while not IsIrreducible(q) do
        q := RandomPoly(FpX, n);
        cnt := cnt + 1;
    end while;
    return q, cnt;
end intrinsic;


intrinsic PrimWurz(K :: RngUPolRes) -> RngUPolResElt, RngIntElt
{Sucht heuristisch eine Primitivwurzel in K.}
    N := #K - 1;
    P := PrimeDivisors(N);
    found := false;
    cnt := 0;
    while not found do
        w := Random(K);
	cnt := cnt + 1;
	if w ne Zero(K) then
	    found := true;
	    j := 1;
	    while j lt #P and found do
	        M := Integers() ! (N/P[j]);
		found := found and ( w^M ne One(K) );
		j := j+1;
	    end while;
	end if;
    end while;
    return w, cnt;
end intrinsic;


intrinsic MyQuadMipo(alpha :: FldNumElt) -> RngUPolElt
{}
    QX<x> := PolynomialRing(Rationals());
    K := Parent(alpha);
    d := Integers() ! (K.1^2);
    coeff := Eltseq(alpha);
    a := coeff[1];
    b := coeff[2];
    if b eq 0 then
        return x - a;
    else
        return x^2 - 2*a*x + (a^2 - b^2*d);
    end if;    
end intrinsic;

intrinsic MyBiquadMipo(alpha :: FldNumElt, om1 :: FldNumElt, om2 :: FldNumElt) -> RngUPolElt
{}
    QX<x> := PolynomialRing(Rationals());
    L := Parent(alpha);
    LX<z> := PolynomialRing(L);

    coeff := Eltseq(alpha);
    r := coeff[1];
    s := coeff[2];
    coeff := Eltseq(r);
    a1 := coeff[1];
    b1 := coeff[2];
    coeff := Eltseq(s);
    a2 := coeff[2];
    b2 := coeff[2];
    a2 := coeff[1];

    d1 := Integers() ! ( om1^2 );
    d2 := Integers() ! ( om2^2 );
    d3 := d1*d2;

    
    if a2 eq 0 and b1 eq 0 and b2 eq 0 then
            return x - a1;
    end if;

    if a2 eq 0 and b2 eq 0 then
        return x^2 - 2*a1*x + (a1^2 - b1^2*d1);
    end if;

    if b1 eq 0 and b2 eq 0 then
        return x^2 - 2*a1*x + (a1^2 - a2^2*d2);
    end if;

    if b1 eq 0 and a2 eq 0 then
        return x^2 - 2*a1*x + (a1^2 - b2^2*d3);
    end if;

    
    alpha1 := (a1-b1*om1) + (a2-b2*om1)*om2;
    alpha2 := (a1+b1*om1) - (a2+b2*om1)*om2;
    alpha3 := (a1-b1*om1) - (a2-b2*om1)*om2;

    h := QX ! ( (z - alpha) * (z - alpha1) * (z - alpha2) * (z - alpha3) );
    return h;
end intrinsic;

intrinsic WhichAut(G :: SeqEnum, i :: RngIntElt, j :: RngIntElt, p :: RngIntElt, m :: RngIntElt) -> Map
{}
    F := Domain(G[1]);
    FX<t> := PolynomialRing(F);
    f := t^p - m;
    nullst := Roots(f); 
    omega := nullst[1,1];
    g := &+[t^i : i in [0..p-1]];
    nullst := Roots(g); 
    zeta := nullst[1,1];

    for g in G do
        if g(zeta) eq zeta^i and g(omega) eq zeta^j*omega then
            return g;
        end if;
    end for;
end intrinsic;


intrinsic FinFldSubfield(K :: FldFin, f :: RngIntElt) -> Map
{}
    fac := Factorization(#K);
    p := fac[1,1]; d := fac[1,2];

    Kmal, j := UnitGroup(K);
    alpha := j(Kmal.1);
    n := Integers()!((p^d-1)/(p^f-1));                                                                                                                    
    beta := alpha^n;
    mipo := MinimalPolynomial(beta);

    return beta, mipo;
end intrinsic;

