Optimizing GAP Code - decomposing Regular Characters into Irreducible Characters & Projecting $\mathbb C^n$ onto an Invariant Subspace

Question

I'm trying to write some GAP code for a couple equations I'm trying to solve. So far I have working code, but I don't have much (or really any) experience with GAP, so I'm hoping my naive translation might be improved. As is I won't be able to solve the problems I need it for, as the input will be way too big.

First, I'm trying to take the character $\chi$ of the regular representation $\rho$ of a group $G$, and find its decomposition into a sum of irreducible characters, i.e. find $m_1,m_2,...,m_k$ such that $\chi = m_1\chi_1 + m_2\chi_2 + ... + m_3\chi_3$. For this I have the formula:

$$m_i = \frac{1}{|G|}\sum_{g \in G} \chi_i(g)\chi(g^{-1})$$

This corresponds to a decomposition of $\rho$ as $\rho = \rho_1^{(m_1)} \oplus \rho_2^{(m_2)} \oplus ... \oplus \rho_k^{(m_k)}$. Once I compute all of these, I want to find the matrices

$$P_i = \frac{m_i \chi_i(1)}{|G|}\sum_{g \in G} \chi_i(g^{-1})\rho(g)$$

Here $P_i$ will be the orthogonal projection of $\mathbb C^n$ onto $V_i$, the invariant subspace for $\rho_i^{(m_i)}$. My ultimate goal is to find bases for each of these invariant subspaces, which I can get easily from each $P_i$.

Now, my current GAP code to compute all of this is:

cclasses := ConjugacyClasses(G);;
irrs := Irr(G);;
reg_chars := RegularCharacters(G,dim);

m := List(irrs, i -> 0);;
for i in [1..Size(irrs)] do
    for g in G do
        irr_index := Position(cclasses, ConjugacyClass(G,g));
        reg_index := Position(cclasses, ConjugacyClass(G,Inverse(g)));
        m[i] := m[i] + irrs[i][irr_index]*reg_chars[reg_index];
    od;
    m[i] := m[i] / Order(G);
od;

P := List(irrs, i -> NullMat(dim,dim));;
for i in [1..Size(irrs)] do
    if m[i] <> 0 then
        for g in G do
            irr_index := Position(cclasses, ConjugacyClass(G,Inverse(g)));
            reg_rep := PermutationMat(g,dim);
            P[i] := P[i] + irrs[i][irr_index]*reg_rep;
        od;
        P[i] := m[i]*irrs[i][1]/Order(G) * P[i];
    fi;
od;

where RegularCharacters is a function

RegularCharacters := function(G,dim)
    local cclasses, ch_list;
    cclasses := ConjugacyClasses(G);
    ch_list := List([1..NrConjugacyClasses(G)],
                    i -> TraceMat(PermutationMat(Representative(cclasses[i]),dim)));
    return ch_list;
end;

Are there any ways to optimize this so I can run it for large groups (i.e. hopefully up to around the size of $S_{20}$ or so) within a reasonable time? (Note: I'm aware there are faster methods for $S_{n}$, I'm interested in other groups that may not have a simpler way to do this).

I've taken this from T. Rehn, Exploring Core Points for Fun and Profit, section 5.5.1. Afterwards, there is another algorithm that should be more suitable for large groups, but despite implementing it (as far as I can tell) correctly, I can't manage to get correct solutions from it. I and another student have both attempted to solve it by hand, and both arrived at the same (incorrect) answer. If you know of any other references to a better algorithm for this problem, please let me know!

Thank you!