Counting points on the Klein quartic In Moreno's book "Algebraic Curves over Finite Fields", he mentions the following in passing with no further comments ($K$ denotes the Klein quartic defined by $X^3 Y + Y^3 Z + Z^3 X = 0$):

The Jacobian of $K$ is a product of three elliptic curves all isogenous to the elliptic curve with complex multiplication over the field $Q(\sqrt{-7})$. This last fact implies the existence of a formula for counting the number of rational points on $K$ over the field $\mathbb{F}_p$ which depends on how the prime $p$ splits in the ring $\mathbb{Z}[\sqrt{-7}]$.

Does anyone know what he is referring to? The Klein quartic is the same as the modular curve $X(7)$ which has genus 3, so I'd be very interested in such a formula!
EDIT: After some computation, it appears that the number of points on $X(7)$ over $\mathbf{F}_q$ is precisely $q+1$ as long as $q \not\equiv 1 \bmod 7$ (*). Otherwise, the behaviour of the error term $a_q = q+1-\#X(7)(\mathbf{F}_q)$ is kind of complicated, but it appears to be constant if we restrict attention to primes $q$ lying in certain quadratic progressions, for example when $q$ is of the form $28n^2 - 28n + 43$ my data suggests that we always have $a_q = -12$. Does anyone have any idea what's going on here? Note that when $q \equiv 1 \bmod 7$ then $a_q/3 = b_q$ where $b_n$ are the Fourier coefficients of the unique cusp form of weight 2 for the congruence subgroup $\Gamma_0(49)$.
Furthermore, what would the moduli interpretation be of statement (*) above? We know $X(7)$ has 24 cusps, but the only time there are elliptic curves defined over $\mathbf{F}_q$ with complete rational 7-torsion is when $q \equiv 1 \bmod 7$ and the trace of Frobenius is 2 mod 7... so why should there be precisely $q-23$ noncuspidal points on $X(7)$ when $q \not\equiv 1 \bmod 7$?

$\def\FF{\mathbb{F}}$I'm switching to another answer to discuss relations with $X_0(49)$ and complex multiplication. Again, Elkies notes are a good reference and $\zeta$ is a primitive $7$-th root of unity. 
Let $\Delta \subset PSL_2(\FF_7)$ be the group of diagonal matrices. Let $Y = \Delta \backslash K$. Then $Y$ is the modular curve corresponding to the arithmetic subgroup $\left( \begin{smallmatrix} a & 7b \\ 7c & d \end{smallmatrix} \right)$. Conjugating by $\left( \begin{smallmatrix} 7 & 0 \\ 0 & 1 \end{smallmatrix} \right)$ turns this group into $\left( \begin{smallmatrix} a & b \\ 49c & d \end{smallmatrix} \right)$ which is $\Gamma_0(49)$, so $Y \cong X_0(49)$. On the level of moduli, when our base field contains a $7$-th root of unity, $K$ parametrizes elliptic curves with a basis for $E[7]$ (up to a condition on the determinant) and $Y$ parametrizes curves with a chosen splitting of $E[7]$ into two cyclic subgroups. See the penultimate paragraph of Elkies Section 4.2 for why $Y \cong X_0(49)$ from a moduli perspective.
We can also use this perspective to explain why $Y$ has CM by $(1+\sqrt{-7})/2$. 
The curve $Y$ has genus $1$; make it into an elliptic curve by choosing the cusp at $\infty$ as the origin. Let $y \in Y$ and let $z \in K$ be a preimage of $y$. Write $\pi$ for the map $K \to Y$. Define 
$$\alpha(y) = \pi \left( \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 2 \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 4 \\ 0 & 1 \end{pmatrix} z \right)$$
and
$$\beta(y) = \pi \left( \begin{pmatrix} 1 & 3 \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 5 \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 6 \\ 0 & 1 \end{pmatrix} z \right)$$
where the sums are in the group law on $Y$.
Note that replacing $z$ by a different preimage $\left( \begin{smallmatrix} a & 0 \\ 0 & a^{-1} \end{smallmatrix}  \right) z$ would just turn $\alpha(y)$ into
$$\pi \left( \begin{pmatrix} a & a^{-1} \\ 0 & a^{-1} \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} a & 2 a^{-1} \\ 0 & a^{-1} \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} a & 4 a^{-1} \\ 0 & a^{-1} \end{pmatrix} z \right)=$$
$$\pi \left( \begin{pmatrix} 1 & a^{-2} \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 2 a^{-2} \\ 0 & 1 \end{pmatrix} z \right)+ \pi \left( \begin{pmatrix} 1 & 4 a^{-2} \\ 0 & 1 \end{pmatrix} z \right) = \alpha(y).$$
To get to the second line, we used that $\pi$ is the quotient by $\Delta$, so $\pi\left( \left( \begin{smallmatrix} a & 0 \\ 0 & a^{-1} \end{smallmatrix} \right) w \right) = \pi(w)$. So $\alpha$ and $\beta$ are well defined endomorphisms of $Y$.
We claim that $\alpha+\beta+1=0$ and $\alpha \beta = \beta \alpha = 2$. So $\alpha$ and $\beta$ generate a quadratic ring, with $\alpha$ and $\beta$ corresponding to $\frac{-1 \pm \sqrt{-7}}{2}$.
We verify that $\alpha+\beta+1=0$: Writing this out, we must show that 
$$\sum_{b=0}^6 \pi\left( \begin{pmatrix} 1 & b \\ 0 & 1 \end{pmatrix} z \right) =0.$$
Write $\psi$ for the map $Y \to X_0(7)$ which forgets one of the two summands into which $Y$ decomposes $E[7]$. Then the above sum is $\sum_{\psi(y')= \psi(y) } y'$. But $X_0(7)$ has genus $0$, so the class of $\sum_{\psi(y') = w} y'$ is independent of the choice of point $w \in X_0(7)$, and taking $y$ to be the cusp shows that this constant value is $0$.
We sketch the verification of $\alpha \beta = \beta \alpha = 2$. Multiplying out both sides, we must show
$$3 y + \sum_{b=1}^6 \pi\left( \begin{pmatrix} 1 & b \\ 0 & 1 \end{pmatrix} z \right) = 2y.$$
Cancelling $2y$ from both sides, this turns into our previous verification.
Of course, all of this is just using that $\mathbb{Z}[\zeta]$ acts on the Jacobian of $K$ and that $\zeta+\zeta^2+\zeta^4 = \frac{-1+\sqrt{-7}}{2}$, but I think it is fun to see it in terms of matrices.
