This is a question about computing coefficients of modular forms, in particular about $E_4^4$, or $E_4^n$. Just a disclaimer, I am quite a beginner in this subject, and I am still slowly learning about modular forms of level $> 1$ and with non-trivial characters, so hopefully the answers will try to stick to level $1$ modular forms :)
To be exact, I want to compute $[q^k]E_4^n$, where the factorisation of $k$ is known and quite smooth. Therefore, I want to decompose the expressions into sum of those that have multiplicative coefficients, such as $\Delta$ (where $\tau(n)$ is multiplicative) or $E_k$ (where $\sigma_{k - 1}(n)$ is multiplicative), from which I can compute them in $\tilde{O}\left(\min_{p \mid n} p\right)$ (or $\log$ of that according to literature). Here is my attempt.
For the first expression, I know that $E_4^2 = E_8$, so the coefficients will equal the convoluting sum $[q^n](E_4^4) = [q^n]E_8^2 = \sum_{k = 1}^{n - 1} \sigma_7(k) \sigma_7(n - k)$. However, $E_8^2 \in \mathcal{M}_{16}$ which is two dimensional, so we can't just write $E_8^2 = cE_{16}$. Instead, we have $E_8^2 = \frac{1}{1911}\left(5258E_{16} - 3617E_4E_{12}\right)$, or other equivalent expressions in any basis of $\mathcal{M}_{16}$. Of course, we run into the same problem, as the coefficients of $E_4E_{12}$ is not multiplicative.
From what I know, eigenforms should be what I am looking for, since they are eigenvectors of the Hecke operators, and such coefficients are multiplicative. Therefore, I can first find an eigenbasis of $\mathcal{M}_{16}$, which is $\{\Delta E_4, E_{16}\}$, then write $E_4^4$ as a linear combination of them, getting $E_4^4 = \frac{3456000}{3617}\Delta E_4 + E_{16}$. Since $E_4 \Delta$ is an eigenform, by theory of Hecke operators, (1) coefficients are multiplicative and (2) we have $T_{p^{e + 2}} = T_pT_{p^{e + 1}} - p^{16 - 1}T_{p^e}$. This is verified by the following Sage script:
sage: M = ModularFormsRing(1)
sage: E4 = M(ModularForms(1, 4).0)
sage: Delta = M(CuspForms(1, 12).0)
sage: coefs = list((E4 * Delta).qexp(101))
sage: for i in range(1, 101):
....: for j in range(1, 101):
....: if i * j >= 100: break
....: if gcd(i, j) > 1: break
....: assert coefs[i] * coefs[j] == coefs[i * j], f"{i} {j}"
....: for p in prime_range(2, 100):
....: for e in range(2, floor(log(100, p)) + 1):
....: assert coefs[p^e] == coefs[p] * coefs[p^(e - 1)] - p^15 * coefs[p^(e - 2)]
However, this already took quite some work, and often times the eigenbasis is not as simple. For example, (if I am not mistaken) the eigenbasis of $\mathcal{M}_{24}$ is given by $\{E_{24}, \Delta E_{12} \pm \left(540 \pm 12\sqrt{144169}\right)\Delta^2\}$, which is slightly annoying to deal with. Moreover, even after I obtain such an eigenbasis of $\mathcal{M}_{4n}$ (which is a degree $8$ number field when $n = 25$) I still have to write $E_4^n$ as a linear combination of them, then also compute the first coefficients of them (over the number field), and finally combine them.
Another direction I thought about would be some kind of recursive representation. However, I am not aware of any such ideas, e.g. can I represent the eigenbasis of $\mathcal{M}_{4k + 12}$ in terms of that of $\mathcal{M}_{4k}$? It seems quite difficult, at least numerically I don't see any relation between the following polynomials (the modulus of the number field):
\begin{align*} P_{12}(x) &= x\\ P_{24}(x) &= x^2 - 1080x - 20468736 \\ P_{36}(x) &= x^3 - 139656x^2 - 59208339456x - 1467625047588864 \\ P_{48}(x) &= x^4 - 5785560x^3 - 467142374034432x^2 + 1426830562183253852160x + 3297913828840214320807673856 \end{align*}
(These are the number fields that the coefficients of the new forms live in, I am not sure if that is what I want here...)
Hence, my question is, is there a simpler way that I am missing?
Thanks a lot in advance!
Thanks to @Aphelli for the hint / answer! I will elaborate a bit here (for my own reference):
I was actually pretty close to the solution at the end, the problem was how to polish the idea of "writing as linear combination of eigenforms". For any cuspidal space $S_{2k}$, by Maeda's conjecture there exists a cusp form $g(q) \in S_{2k}$ such that $S_{2k}$ is spanned by $\{g(q), g^{\sigma}(q), \cdots, g^{\sigma^{d - 1}}(q)\}$, where $d = \dim S_{2k}$ and $\sigma$ is Galois conjugation. In Sage, this is obtained as:
The key is then to note that if say $E_4^{12} - E_{48} = c_1 g + c_2 g^{\sigma} + c_3 g^{\sigma^2} + c_4 g^{\sigma^4} \in \mathbb{Q}[q]$, then by replacing $c_2, c_3, c_4$ with their appropriate conjugation we get
$$ E_4^{12} - E_{48} = c_1g + (c_2g)^{\sigma} + (c_3g)^{\sigma^2} + (c_4g)^{\sigma^3} $$
Since the left coefficients are all rational, if we apply $\sigma^k$ to both sides we get
\begin{align*} \underbrace{\sigma(E_4^{12} - E_{48})}_{\text{fixed}} &= (c_1g)^{\sigma} + (c_2g)^{\sigma^2} + (c_3g)^{\sigma^3} + (c_4g) \\ \sigma^2(E_4^{12} - E_{48}) &= (c_1g)^{\sigma^2} + (c_2g)^{\sigma^3} + (c_3g) + (c_4g)^{\sigma} \\ \sigma^3(E_4^{12} - E_{48}) &= (c_1g)^{\sigma^3} + (c_2g) + (c_3g)^{\sigma} + (c_4g)^{\sigma^2} \\ \end{align*}
Adding up the columns just gives $4(E_4^{12} - E_{48}) = \mathrm{Tr}(c_1g + c_2g + c_3g + c_4g)$, i.e. $E_4^{12} - E_{48} = \mathrm{Tr}(cg)$
This is elementary Galois theory but it went over my head for a second :)
Hopefully this code is useful for someone in the future!