We describe the practical implementation of an average polynomial-time algorithm for counting points on superelliptic curves defined over $mathbb Q$ that is substantially faster than previous approaches. Our algorithm takes as input a superelliptic curves $y^m=f(x)$ with $mge 2$ and $fin mathbb Z[x]$ any squarefree polynomial of degree $dge 3$, along with a positive integer $N$. It can compute $#X(mathbb F_p)$ for all $ple N$ not dividing $mmathrm{lc}(f)mathrm{disc}(f)$ in time $O(md^3 Nlog^3 Nloglog N)$. It achieves this by computing the trace of the Cartier--Manin matrix of reductions of $X$. We can also compute the Cartier--Manin matrix itself, which determines the $p$-rank of the Jacobian of $X$ and the numerator of its zeta function modulo~$p$.