// *********************** Package rjf.m ****************************
//
// Contains functions for constructing the Real Jordan Form of a 
// matrix.

intrinsic Re(v::ModTupFldElt) -> ModTupFldElt
    {Return the real part of a vector}

    return VectorSpace(RealField(), 
                       Degree(v))![Re(v[i]) : i in [1..Degree(v)]];
end intrinsic;

intrinsic Im(v::ModTupFldElt) -> ModTupFldElt
    {Return the imaginary part of a vector}

    return VectorSpace(RealField(), 
                       Degree(v))![Im(v[i]) : i in [1..Degree(v)]];
end intrinsic;

intrinsic LambdaMatrix(lambda::.) -> AlgMatElt
    {Given a real or complex number lambda return a matrix used
     in the construction of the RealJordanBlock(lambda)}

    return IsReal(lambda) 
               select MatrixAlgebra(RealField(), 1)![lambda]
               else   MatrixAlgebra(RealField(), 2)![a, b, -b, a]
                          where a := Re(lambda)
                          where b := Im(lambda);
end intrinsic;

intrinsic RealJordanBlock(lambda::., k::RngIntElt) -> AlgMatElt
    {Return the RealJordanBlock of lambda of dimension k i.e.
     J(lambda, k) if lambda is real and K(lambda, k) otherwise}

    M := MatrixAlgebra(RealField(), k);
    N := M!0;
    for i in [1 .. k-1] do
        N[i+1, i] := 1;
    end for;
    return TensorProduct(Ik, LambdaMatrix(lambda)) +
           TensorProduct(N, I)
               where Ik := M!1
               where I := MatrixAlgebra(RealField(), 
                                        IsReal(lambda)
                                            select 1
                                            else   2
                                        )!1;
end intrinsic;

intrinsic JoinJordanBlocks(blocks::List) -> ModMatFldElt
    {Return the Jordan matrix with Jordan Blocks blocks}

    J := blocks[1];
    Remove(~blocks, 1);
    for B in blocks do
        J := VerticalJoin(HorizontalJoin(J, OmJnB),
                          HorizontalJoin(OmBnJ, B))
                 where OmJnB := RMatrixSpace(RealField(), mJ, nB)!0
                 where OmBnJ := RMatrixSpace(RealField(), mB, nJ)!0
                 where mJ := Nrows(J)
                 where nJ := Ncols(J)
                 where mB := Nrows(B)
                 where nB := Ncols(B);
    end for;
    return J;
end intrinsic;

intrinsic RealJordanForm(X::., p::RngIntElt) 
                         -> ModMatFldElt, ModMatFldElt
    {Given a sequence or matrix X return the Real Jordan Form J of X
     and a matrix U such that J = U^-1 X U, where calculations in the
     Real Field are done to precision p}

    if p eq 0 then
        R := RealField();
        C := ComplexField();
    else
        R := RealField(p);
        C := ComplexField(p);
    end if;
    if Type(X) eq SeqEnum then
        isSquare, n := IsSquare(#X);
        if not isSquare then
            error "Argument 1 is not a square matrix";
        else
	    delete isSquare;
        end if;
    else
        n := Nrows(X);
    end if;
    // Magma considers matrices to act on the right ...
    // hence we need to transpose X and consider the vectors
    // we get as multiplying on the left.
    M := MatrixAlgebra(R, n);
    X := Transpose(M!X);
    I := M!1;
    eigenvalues := [e : e in Eigenvalues(X)
                      | IsReal(e[1]) or Im(e[1]) gt 0];
    eigenvalues := Sort([e : e in eigenvalues | IsReal(e[1])],
                        func<a, b | b[1] - a[1]>) cat
                   Sort([e : e in eigenvalues | not IsReal(e[1])],
                        func<a, b | Norm(b[1]) - Norm(a[1])>);
    eigenbasis := [];
    jordanblocks := [* *];
    jordandescription := "";
    AssertAttribute(FldPr, "OutputPrecision", 5);
    for lambdai in eigenvalues do
        V := VectorSpace(IsReal(lambdai[1])
                             select R
                             else   C,
                         n);
        XmlambdaI := X - lambdai[1] * I;
        B := I;
        geneigenspaces := [ sub<V|> ];  // 0-dimensional subspace of V,
                                        // used as a sentinel
        while lambdai[2] gt Dimension(geneigenspaces[#geneigenspaces]) do
            B *:= XmlambdaI;
            geneigenspaces cat:= [ NullSpace(B) ];
        end while;
        eigendepth := #geneigenspaces;
        eigenchains := [ [u] : u in Basis(geneigenspaces[eigendepth])
                             | u notin geneigenspaces[eigendepth - 1] ];
        while eigendepth gt 2 do
            for i in [1..#eigenchains] do
                Append(~eigenchains[i], 
                       eigenchains[i, #eigenchains[i]] * XmlambdaI);
            end for;
            eigendepth -:= 1;
            eigenchains cat:= 
                [ [u] : u in Basis(geneigenspaces[eigendepth])
                      | u notin geneigenspaces[eigendepth - 1] and
                        u notin sub<V | [chain[#chain] : chain in
                                                         eigenchains]> ];
        end while;
        for chain in eigenchains do
            Append(~jordanblocks, RealJordanBlock(lambdai[1], #chain)); 
            jordandescription 
                cat:= IsReal(lambdai[1])
                          select Sprintf(" (+) J_%o(%o)",
                                         #chain, lambdai[1])
                          else   Sprintf(" (+) K_%o(%o, %o)",
                                         #chain, Re(lambdai[1]),
                                         Im(lambdai[1]));
        end for;
        eigenbasis 
            cat:= &cat (IsReal(lambdai[1])
                           select eigenchains
                           else   [ &cat [ [ Re(chaini), Im(chaini)]
                                             : chaini in chain]
                                           : chain in eigenchains ]);
    end for;
    jordandescription := Substring(jordandescription, 6,
                                   #jordandescription - 5);
    return JoinJordanBlocks(jordanblocks),     // Real Jordan Form J
           Transpose(VerticalJoin(eigenbasis)),// Transformation Matrix U
                                               // such that J = U^-1 X U
           jordandescription;                  // J(lambda, k) (+) ...
end intrinsic;

intrinsic RealJordanForm(X::.) -> ModMatFldElt, ModMatFldElt
    {Given a sequence or matrix X return the Real Jordan Form J of X
     and a matrix U such that J = U^-1 X U}

    return J, U, S where J, U, S := RealJordanForm(X, 0);
end intrinsic;
// *********************** End of rjf.m *****************************