Computing millions of coefficients of non self-dual modular forms To test some conjectures made by some colleagues, I need to compute millions of coefficients of non self-dual modular forms, preferably in low weight (say 2 or 3).  A form such as this.
For elliptic curves over $\mathbb{Q}$ I can compute millions of coefficients and for newforms given as eta products I can, as well.  One example of a non self-dual form we can do is to compute a self-dual form and twist by a character, but we are looking for forms that don't arise in the way.
In Sage and MAGMA I can only compute hundreds of thousands of coefficients of the form linked to above.  
 A: You might want to try Brandt matrices and lattice methods instead. If you can write your function as a theta series, then the computation of the terms can be done very quickly through multiplication of power series. See http://projecteuclid.org/download/pdf_1/euclid.jmsj/1240433920 for an introduction to the theory, and note that $\sum q^{ax^2+by^2+cz^2}=\sum q^{ax^2} \cdot \sum q^{by^2} \cdot \sum q^{cy^2}$, and the product can be evaluated very quickly.
There is a complication given by the fact you want a character, but you can twist theta series to get them to transform with a particular character.
A: 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. $$
