#############################################################################
####
##
#W  SplittField.gi            RADIROOT package                Andreas Distler
##
##  Installs the functions to compute the splitting field of a polynomial
##
#H  $Id:$
##
#Y  2004
##


#############################################################################
##
#M  SplittingField( <f> )
##
##  Returns the smallest field, that contains the roots of the irreducible,
##  rational polynomial <f> 
##
InstallMethod( SplittingField, "rational polynomials",
[ IsUnivariatePolynomial ],
function( f )

    if ForAny( CoefficientsOfUnivariatePolynomial(f),
               c -> not c in Rationals ) then TryNextMethod( ); fi;

    return RR_Zerfaellungskoerper( f, rec( roots := [ ],
                                           degs := [ ],
                                           coeffs := [ ],
                                           K := FieldByMatrices([ [[ 1 ]] ]),
                                           H := Rationals ) ).H;

end );

#############################################################################
##
#F  RR_BegleitMatrix( <f>, <A> )
##
##  Computes the companionmatrix of the polynomial <f> with respect to
##  the field generated by the matrix <A>
##
InstallGlobalFunction( RR_BegleitMatrix, function( f, A )
    local matrix, coeff, blockmat, deg, i, k, l;

    deg := Degree(f);
    coeff := CoefficientsOfUnivariatePolynomial(f);
    matrix := NullMat( deg*Size(A), deg*Size(A), Rationals);

    # Letzte Zeile erstellen
    for i in [ 1..deg ] do
        # Matrix, die dem i-ten Koeffizienten entspricht
        blockmat := -RR_RootInK( A, coeff[i] );
        for k in [1..Size(A)] do
            for l in [1..Size(A)] do
                matrix[(deg-1)*Size(A)+k][(i-1)*Size(A)+l] := blockmat[k][l];
            od;
        od;
    od;

    # Nebendiagonale mit Einsen fuellen
    for i in [1..(deg-1)] do
        for k in [1..Size(A)] do 
            matrix[(i-1)*Size(A)+k][i*Size(A)+k] := 1;
        od;
    od;
    	
    return matrix;
end );


#############################################################################
##
#F  RR_BlowUpMat, function( <mat>, <n> )
##
##  Computes a matrix that is <n>-times bigger than <mat> and has
##  <mat> on the <n> blocks with size of <mat> at the diagonal
##
InstallGlobalFunction( RR_BlowUpMat, function( mat, n )
    local i, j, k, Mat;

    Mat := NullMat( n * Size(mat), n * Size(mat), Rationals );
    # auf den n Diagonalbloecken werden die Nullen durch mat ersetzt
    for i in [ 1..n ] do
        for j in [ 1..Size( mat ) ] do
            for k in [ 1..Size( mat ) ] do
                Mat[ (i-1)*Size(mat)+j ][ (i-1)*Size(mat)+k ] := mat[j][k];
            od;
        od;
    od;

    return Mat;
end );


#############################################################################
##
#F  RR_MatrixField( <f>, <mat> )
##
##  Returns the matrixfield that arises from adjoining a root of the
##  polynomial <f> to the matrixfield generated by <mat>
##
InstallGlobalFunction( RR_MatrixField, function( f, mat )
    local A, B;

    # mat als Matrix in Oberkoerper
    # mat wird deg(f)-mal auf die Diagonale gesetzt 
    A := RR_BlowUpMat( mat, Degree(f) );

    # Begleitmatrix von f bzgl. des durch mat erzeugten Koerpers
    B := RR_BegleitMatrix( f, mat );

    return FieldByMatricesNC( [A, B] );
end );


#############################################################################
##
#F  RR_RootInH( <erw>, <a> )
##
##  The record <erw> contains two isomorphic fields. One generated
##  with AlgebraicExtension and the other as matrixfield. Both are
##  defined by a primitive element. This function transfers the
##  matrix <a> to it's isomorphic symbolic represenation
##
InstallGlobalFunction( RR_RootInH, function( erw, a )
    local coeff, bas, n;

    n := Size(PrimitiveElement(erw.K));

    # Basis: 1, primEl, ... , primEl^(n-1) in Matrix umwandeln
    bas := Basis( erw.K, List( [1..n], i -> PrimitiveElement(erw.K)^(i-1) ) );

    # Koeffizienten bzgl. der standardisierten Basis
    coeff := Coefficients( bas, a );

    return Sum( [ 1..n ], i -> coeff[i] * PrimitiveElement(erw.H)^(i-1) );
end ); 
    

#############################################################################
##
#F  RR_RootInK( <primEl>, <coeff> )
##
##  Does the inverse of RR_RootInH; the fieldelement given symbolic
##  by it's external representation <coeff> is transfered in a matrix
##  of the field generated by <primEl>
##
InstallGlobalFunction( RR_RootInK, function( primEl, coeff )
    local i, mat;

    mat := NullMat( Size(primEl), Size(primEl), Rationals );
    for i in [1..Size(primEl)] do
        mat := mat +  ExtRepOfObj(coeff)[i] * primEl^(i-1);
    od;

    return mat;
end );


#############################################################################
##
#F  RR_Zerfaellungskoerper( <poly>, <erw> )
##
##  Computes the splitting field of the polynomial <poly>. In the
##  record <erw> the field is stored as matrix field as well as in a
##  symbolic represenation generated by
##  AlgebraicExtension. The roots of <poly> are also stored.  
##
InstallGlobalFunction( RR_Zerfaellungskoerper, function( poly, erw )
    local matA,matB,faktoren,i,f,minpol,roots,primEl;

    if Degree( poly ) = 1 then 
        erw.roots := [ [ ], [ ] ];
        return erw;
    fi;
    roots := [ ];
 
    # Schleife bis poly in Linearfaktoren zerfaellt
    while Length(erw.roots) + Length(roots) + 1 < Degree(poly) do
	
        # poly ueber neuem H faktorisieren	
	faktoren := FactorsPolynomialKant( poly, erw.H );;
	Info( InfoRadiroot, 4, "    Factorization of polynomial:\n",
              faktoren );
    	f := faktoren[ Length( faktoren ) ];
        if Degree( f ) = 1 then break; fi;

        roots := RR_Roots( [ erw.roots, roots, 
                             List( Filtered( faktoren, f -> Degree( f ) = 1 ), 
                                   f -> 
                                   -CoefficientsOfUnivariatePolynomial(f)[1])],
                           erw );

	erw.K := RR_MatrixField( f, PrimitiveElement( erw.K ) );
	Add( erw.degs, Degree(f) );
        SetDegreeOverPrimeField( erw.K, Product( erw.degs ));
        Info( InfoRadiroot, 2,"    Degree of the extension: ", Degree(f) );

	matA := GeneratorsOfField( erw.K )[ 1 ];;
	matB := GeneratorsOfField( erw.K )[ 2 ];;

	# Wurzelliste aktualisieren
	for i in [ 1..Length(erw.roots) ] do
	    erw.roots[i] := RR_BlowUpMat( erw.roots[i], Degree( f ) );
	od;
	for i in [ 1..Length(roots) ] do
	    roots[i] := RR_BlowUpMat( roots[i], Degree( f ) );
	od;
        if IsBound( erw.unity ) then 
            erw.unity := RR_BlowUpMat( erw.unity, Degree( f ) );
        fi;

        Info( InfoRadiroot, 4, "Adjoined root:\n", matB );
	Add( erw.roots, matB );

	# Primitives Element suchen
        Info( InfoRadiroot, 2, "    Searching for a primitive element" );
        primEl := Sum([1..Length(erw.roots)], i -> i * erw.roots[i]);
        if IsBound( erw.unity ) then
            primEl := primEl+erw.unity;
        fi;
        minpol := MinimalPolynomial( Rationals, primEl );
        if Degree( minpol ) = Product( erw.degs ) then
            SetPrimitiveElement( erw.K, primEl );
            SetDefiningPolynomial( erw.K, minpol );
            Add( erw.coeffs, Length( erw.roots ) );
        else
            for i in [ Minimum( 2 * [ Length( erw.degs )-1, 1 ] )..99 ] do
                minpol := MinimalPolynomial( Rationals, i * matA + matB );
                if Degree( minpol ) = Product( erw.degs ) then
                    SetPrimitiveElement( erw.K, i * matA + matB );
                    SetDefiningPolynomial( erw.K, minpol );
                    erw.coeffs := Flat( [ i * erw.coeffs, 1 ] ); 
		    break;
	        fi;
            od;
        fi;

        erw.H := AlgebraicExtension( Rationals, minpol );
        Info( InfoRadiroot, 3, "    ", minpol, " is defining polynomial.");
    od;
    erw.roots := [ Concatenation( erw.roots, roots ),
                   List( faktoren, f -> -Value( f, 0 ) ) ];
    if IsBound( erw.unity ) then
        erw.degs := erw.degs{[ 2..Length(erw.degs) ]};
    fi;

    Info( InfoRadiroot, 3, "        Composition of the primitive element: ",
                           erw.coeffs );

    return erw;
end );


#############################################################################
##
#F  RR_SplittField( <poly>, <m> )
##
##  Calls the function RR_Zerfaellungskoerper for the polynomial <poly> with
##  special initial values. The splitting field is constructed over a
##  cyclotomic field.
##
InstallGlobalFunction( RR_SplittField, function( poly, m )
    local erw, cyclopol;

    cyclopol := CyclotomicPolynomial(Rationals, m);

    erw := rec( roots := [ ], degs := [ Degree(cyclopol) ], coeffs := [ ],
                K := RR_MatrixField( cyclopol, [[ 1 ]]),
                H := AlgebraicExtension( Rationals, cyclopol ));
    erw.unity := PrimitiveElement( erw.K );

    erw := RR_Zerfaellungskoerper( poly, erw );;

    return erw;
end );


#############################################################################
##
#E












