
intrinsic FindAllRels(ff::[RngElt], NecessaryMonomials, PossibleMonomials) -> SeqEnum
    {}
    _, JacRels := JacobianRelations(ff);
    
    M := Universe(JacRels);
    F := BaseRing(M);

    // We need to apply vector space routines:
    require IsField(F): "Base ring must be a field";
    
    MonList := [M|];
    for i in [1..16] do
        MonList cat:= [ M.i*M.j : j in [i..16]];
    end for;
    
    V := VectorSpace(F, #MonList);
    
    RelsMat := [];
    
    for i in [1..#JacRels] do
        Mons := Monomials(JacRels[i]);
        Coeffs := Coefficients(JacRels[i]);
        RelsMat[i] := &+[Coeffs[j]*V.(Index(MonList,Mons[j])) : j in [1..#Mons]];
    end for;
    
    RelsMat := Matrix(RelsMat);
    
    // We first list all the monomials that can appear in our relation:    
    WantedMonsVectors := [ V.Index(MonList, mon) : mon in PossibleMonomials];
        
    // We then create the corresponding space...
    WantedSpace := sub < V | WantedMonsVectors >;
    // ... and intersect it with the relations:
    WantedSpace :=  WantedSpace meet Image(RelsMat);
        
    // We now list all the monomials (one of which _must_ appear);
    NeededMons := [ V.Index(MonList, mon) : mon in NecessaryMonomials];
        
    B := Basis(WantedSpace);
        
    PossibleRelations := {};
        
    for i in [1..#B] do
        for j in NeededMons do
            if InnerProduct( j , B[i]) ne 0 then
                sol := Solution(RelsMat, B[i]);
                sol := Eltseq(sol);
                PossibleRelations join:= {&+[sol[i]*JacRels[i] : i in [1..#JacRels]]};
            end if;
        end for;
    end for;
    
    // Have we found a relation?
    if not IsEmpty(PossibleRelations) then
        return SetToSequence(PossibleRelations);
    else 
        return false;
    end if;
end intrinsic;


// JACOBIAN SOLVING RELATIONS

intrinsic JacSolveRels(ff::[RngElt], AlreadyKnown, Wanted) -> SeqEnum
    {}
    _, JacRels := JacobianRelations(ff);
    
    M := Universe(JacRels);
    F := BaseRing(M);

    require IsField(F): "Base ring must be a field";
    
    MonList := [M|];
    for i in [1..16] do
        MonList cat:= [ M.i*M.j : j in [i..16]];
    end for;
    
    V := VectorSpace(F, #MonList);
    
    RelsMat := [];
    
    for i in [1..#JacRels] do
        Mons := Monomials(JacRels[i]);
        Coeffs := Coefficients(JacRels[i]);
        RelsMat[i] := &+[Coeffs[j]*V.(Index(MonList,Mons[j])) : j in [1..#Mons]];
    end for;
    
    RelsMat := Matrix(RelsMat);
    
    // We first list all the monomials that can appear in our relation:
    WantedMons := SetToSequence({ M.(i+1)*M.(j+1) : i in AlreadyKnown, j in AlreadyKnown});
    WantedMons cat:= [ M.(i+1)*M.(Wanted+1) : i in AlreadyKnown];
    WantedMons cat:= [ M.(i+1)*M.(Wanted + 1) : i in AlreadyKnown];
        
    WantedMonsVectors := [ V.Index(MonList, mon) : mon in WantedMons];
        
    // We then create the corresponding space...
    WantedSpace := sub < V | WantedMonsVectors >;
    // ... and intersect it with the relations:
    WantedSpace :=  WantedSpace meet Image(RelsMat);
        
    // We now need to find an element in the WantedSpace that
    // contains a term of the form S.CurrentIndex*S.j for some j > CurrentIndex.
        
    NeededMons := [ M.(Wanted+1)*M.(i+1) : i in AlreadyKnown];
        
    B := Basis(WantedSpace);
        
    PossibleRelations := [];
        
    for i in [1..#B] do
        for j in NeededMons do
            if InnerProduct( V.Index(MonList, j), B[i]) ne 0 then
                sol := Solution(RelsMat, B[i]);
                sol := Eltseq(sol);
                PossibleRelations cat:= [&+[sol[i]*JacRels[i] : i in [1..#JacRels]]];
            end if;
        end for;
    end for;
    
    if not IsEmpty(PossibleRelations) then
        // We've found a relation:
        // "false" because we don't need to evaluate a square root to solve:
        return PossibleRelations, false;
    end if;
    
    // We are here if, regrettably(!), we need to evaluate a square root:
        
    // We add the square to both wanted and needed:
    WantedMonsVectors cat:= [ V.Index(MonList, M.(Wanted+1)*M.(Wanted+1))];
    NeededMons cat:= [ M.(Wanted+1)*M.(Wanted+1)];
        
    WantedSpace := Image(RelsMat) meet sub < V | WantedMonsVectors >;
        
    B := Basis(WantedSpace);
        
    PossibleRelations := [];
        
    for i in [1..#B] do
        for j in NeededMons do
            if InnerProduct( V.Index(MonList, j), B[i]) ne 0 then
                sol := Solution(RelsMat, B[i]);
                sol := Eltseq(sol);
                PossibleRelations cat:= [&+[sol[i]*JacRels[i] : i in [1..#JacRels]]];
            end if;
        end for;
    end for;

    if not IsEmpty(PossibleRelations) then
        return PossibleRelations, true;
    else
        printf "Failed to find a relation for s%o\n", Wanted;
        return false;
    end if;

end intrinsic;

// JACOBIAN SOLVING SEQUENCE

intrinsic JacobianSolvingSequence(ff::[RngElt], AlreadyKnown, Wanted) -> SeqEnum
    {}

    ReturnSequence := [];
    // We return a sequence in which each entry has the form:
    // < index, relation, true/false>;
    // Where index is the s_i that it is possible to solve for, relation is the 
    // relation from which it is possible to derive s_i and true/false
    // indicates whether or not it is necessary to evaluate a square root:

    AK := AlreadyKnown;
    for x in [1..#Wanted] do
        CurrentIndex := Wanted[x];
        Rels, NeedSquareRoot := JacSolveRels(ff, AK, CurrentIndex);
        Sizes := [ 1/#Coefficients(Rels[i]) : i in [1..#Rels]];
        _, RelIndex := Maximum(Sizes);
        ReturnSequence[x] := <CurrentIndex, Rels[RelIndex], NeedSquareRoot>;
        AK cat:= [Wanted[x]];
    end for;
        
    return ReturnSequence;
    
end intrinsic;

intrinsic x5Rel(ff::[RngElt]) -> SeqEnum
    {}
    
    _, Jac := JacobianRelations(ff);
        
    M := Universe(Jac);
    N<Y00,Y01,Y02,Y03,Y04,Y05,Y06,Y07,Y08,Y09,Y10,Y11,Y12,Y13,Y14,Y15>
        := FieldOfFractions(M);
    
    Expand := [N.i : i in [1..16]];
    
    rels := JacSolveRels(ff, [15,14,13,12,11,10,5], 0);
    pos0 := [];
    for i in [1..#rels] do
        top, bottom := Explode(Coefficients(rels[i],M.1));
        pos0[i] := -(top/bottom);
    end for;
    
    rels := JacSolveRels(ff, [15,14,13,12,11,10,5], 3);
    pos3 := [];
    for i in [1..#rels] do
        top, bottom := Explode(Coefficients(rels[i],M.4));
        pos3[i] := -(top/bottom);
    end for;
    
    rels := JacSolveRels(ff, [15,14,13,12,11,10,5], 4);
    pos4 := [];
    for i in [1..#rels] do
        top, bottom := Explode(Coefficients(rels[i], M.5));
        pos4[i] := -(top/bottom);
    end for;
    
    Rels := {};
    
    rels5 := JacSolveRels(ff, [15, 14, 13, 12, 11, 10, 3, 4], 5);
    
    for rel0 in pos0, rel3 in pos3, rel4 in pos4, rel in rels5 do
        Expand[1] := rel0;
        Expand[4] := rel3;
        Expand[5] := rel4;
        Rels join:= {Numerator(Evaluate(rel, Expand))};
    end for;
    
    return Rels;
    
end intrinsic;

