One shortcut you could use for computing the level 17 form you link to would be the following. There are exactly 8 Eisenstein series of weight 1 for $\Gamma_1(17)$ and they are all given by completely explicit $q$-expansion formulae (involving sums of divisors, etc). The pairwise products of these series span the space of forms of weight 2. So if $F$ is the form from your link, you can find some finite collection of weight 1 Eisenstein series $f_i, g_i$ and constants $\lambda_i$ such that
$$ F = \sum_i \lambda_i f_i g_i .$$
So if you can compute the $f_i$ and $g_i$ up to the $q^N$ term, which just amounts to computing the divisors of the integers up to $N$ (about 30 seconds in Sage for $N = 10^6$ on my machine), and then evaluate the pairwise products to the same precision (again, this should take a couple of minutes at most in that sort of range), then this will allow you to compute $F$ up to precision $N$.
EDIT. I implemented this algorithm (in Sage), mostly out of curiosity to see how it would perform. See here for the code. It turns out that it doesn't perform too badly: it took 38 minutes to compute the first million coefficients. As an illustration, if $p = 1,000,003$ is the next prime after 1 million, then
$$ a_p(F) = 277 \zeta_8^3 + 277 \zeta_8^2 - 109 \zeta_8 + 109. $$