// Given a group G, elements Hgens generating a subgroup H and
// a list gList of elements, check that
// gList contains exactly one representative of each left coset
// gH.
function CheckRepresentativesForLeftQuotient(G, Hgens, gList)
    Gperm, iso := PermutationGroup(G);
    H := sub<Gperm | [ iso(h) : h in Hgens ]>;
    gListPerm := [ iso(g) : g in gList ];
    Gsize := Order(Gperm);

    return
        (Gsize eq #gList * Order(H)) and
        (Gsize eq #{ g * h : g in gListPerm, h in H });
end function;

function CheckConjugatorsAndDehnfillingsForDifferentCusps(
                                        BI, P, conjugatorsAndDehnfillings)
    for i in [1..#P] do
        conjugators := [ c[1] : c in conjugatorsAndDehnfillings[i] ];
        if not CheckRepresentativesForLeftQuotient(BI, P[i], conjugators) then
            return false;
        end if;
    end for;
    return true;
end function;

function ComputePeripheralCurves(P, nklTriples, conjugatorsAndDehnfillings)
    return [ c[1] *
          ((P[i][1]^nklTriples[i][1])^c[2][1] *
           (P[i][1]^nklTriples[i][2]*P[i][2]^nklTriples[i][3])^c[2][2]) *
          c[1]^-1
          : c in conjugatorsAndDehnfillings[i],
            i in [1..#P] ];
end function;

function ComputeBIandNI(Bianchi, P, nklTriples)
    PI := [ p : p in [ P[i][1]^nklTriples[i][1],
                       P[i][1]^nklTriples[i][2] * P[i][2]^nklTriples[i][3] ],
                i in [1..#P] ];
    BI := quo<Bianchi | PI>;
    NI := NormalClosure(Bianchi, sub<Bianchi | PI>);
    NI := Rewrite(Bianchi, NI : Simplify := false);

    return BI, NI;
end function;

function VerifyLink(
    Bianchi, P, nklTriples, size, conjugatorsAndDehnfillings)

    BI, NI := ComputeBIandNI(Bianchi, P, nklTriples);

    if not (Order(BI) eq size) then
        return "Unexcepted size of B(I)";
    end if;

    if not CheckConjugatorsAndDehnfillingsForDifferentCusps(
                                BI, P, conjugatorsAndDehnfillings) then
        return "Not exactly one curve for each cusp.";
    end if;

    Q := quo<NI |
        ComputePeripheralCurves(P, nklTriples, conjugatorsAndDehnfillings) >;
    Q := ReduceGenerators(Q);

    if not (Order(Q) eq 1) then
        return "Peripheral curves do not kill fundamental group";
    end if;

    numComp := #[ c : c in cs, cs in conjugatorsAndDehnfillings ];

    return IntegerToString(numComp) cat "-Link";
end function;

function Symmetrize(
            g, n, conjugatorsAndDehnfillings, conjugatorsAndDehnfillingsSym)
    return [ conjugatorsAndDehnfillings[j] cat
                 [ <g^i * c[1], c[2]>
                   : i in [0..n-1], c in conjugatorsAndDehnfillingsSym[j] ]
             : j in [1..#conjugatorsAndDehnfillings] ];
end function;

function Expand(conjugators, dehnFillings)
    return [ [ <conjugators[i], d[i]> : i in [1..#conjugators] ]
               : d in dehnFillings ];
end function;

function VerifyLink2(Bianchi, P, exponents, size)
    Pcurve := [ P[i][1]^exponents[i][1] * P[i][2]^exponents[i][2]
                : i in [1..#P] ];
    G := quo<Bianchi | Pcurve>;
    if not (Order(G) eq size) then
        return "Unexcepted size of quotient.";
    end if;

    return "Link";
end function;