To add to Aşağı Güzdək's nice answer, there is something you can do also when the genus is greater than zero which is also in some sense an exact formula. Namely, if $\Gamma$ is one of $\Gamma_0(N), \Gamma_1(N)$ or $\Gamma(N)$ [or probably some other congruence subgroups as well, but I am not familiar with this story], and $Y_\Gamma$ (resp. $X_\Gamma$) is the corresponding open (resp. compact) modular curve, then one has $$ |Y_\Gamma(\mathbf{F}_p)| = p+1-\mathrm{Tr}(\mathrm{Frob}_p|E_2(\Gamma))-\mathrm{Tr}(\mathrm{Frob}_p|S_2(\Gamma)),$$ where $E_k(\Gamma)$ (resp. $S_k(\Gamma)$) is the Galois representation attached to Eisenstein series (resp. cusp forms) of weight $k$ for $\Gamma$. The Galois representation $E_2(\Gamma)$ is given by the $\ell$-adic cohomology of $X^\infty_\Gamma = X_\Gamma \setminus Y_\Gamma$, hence the appearance of the number of cusps in the previous answer. In general one can for a given $\Gamma$ and p figure out which of the cusps are defined over what fields by looking explicitly at the modular interpretation of $X_\Gamma^\infty$ described in the paper by Deligne & Rapoport.

The trace of Frobenius on the space of cusp forms is the really interesting one. In many cases this can be computed for a given $\Gamma$ and $p$ by looking in tables like the Modular Forms Database of William Stein, or asking SAGE to spit them out for you (see e.g. http://www.sagemath.org/doc/reference/sage/modular/modform/cuspidal_submodule.html) -- for a prime p, the trace of $\mathrm{Frob_p}$ coincides with the trace of the Hecke operator $T_p$, which can be read off from the pth Fourier coefficients of a basis of normalized eigenforms for the space of cusp forms you are looking at.

It seems the only built-in functionality in SAGE is for the groups $\Gamma_0(N)$ and $\Gamma_1(N)$. But perhaps one can still do something for small N, e.g. using that conjugation by the matrix with entries $(0,-1;N,0)$ takes $\Gamma(N)$ to $\Gamma_0(N^2) \cap \Gamma_1(N)$ -- this tells you how to get cusp forms for $\Gamma(N)$ both from $\Gamma_0(N^2)$ and $\Gamma_1(N)$, and by the genus formulae for $\Gamma(N)$ you can tell if this has produced all of them.

To add to Aşağı Güzdək's nice answer, there is something you can do also when the genus is greater than zero which is also in some sense an exact formula. Namely, if $\Gamma$ is one of $\Gamma_0(N), \Gamma_1(N)$ or $\Gamma(N)$ [or probably some other congruence subgroups as well, but I am not familiar with this story], and $Y_\Gamma$ (resp. $X_\Gamma$) is the corresponding open (resp. compact) modular curve, then one has $$ |Y_\Gamma(\mathbf{F}_p)| = p+1-\mathrm{Tr}(\mathrm{Frob}_p|E_2(\Gamma))-\mathrm{Tr}(\mathrm{Frob}_p|S_2(\Gamma)),$$ where $E_k(\Gamma)$ (resp. $S_k(\Gamma)$) is the Galois representation attached to Eisenstein series (resp. cusp forms) of weight $k$ for $\Gamma$. The Galois representation $E_2(\Gamma)$ is given by the $\ell$-adic cohomology of $X^\infty_\Gamma = X_\Gamma \setminus Y_\Gamma$, hence the appearance of the number of cusps in the previous answer. In general one can for a given $\Gamma$ and p figure out which of the cusps are defined over what fields by looking explicitly at the modular interpretation of $X_\Gamma^\infty$ described in the paper by Deligne & Rapoport.

The trace of Frobenius on the space of cusp forms is the really interesting one. In many cases this can be computed for a given $\Gamma$ and $p$ by looking in tables like the Modular Forms Database of William Stein, or asking SAGE to spit them out for you (see e.g. http://www.sagemath.org/doc/reference/sage/modular/modform/cuspidal_submodule.html) -- for a prime p, the trace of $\mathrm{Frob_p}$ coincides with the trace of the Hecke operator $T_p$, which can be read off from the pth Fourier coefficients of a basis of normalized eigenforms for the space of cusp forms you are looking at.

It seems the only built-in functionality in SAGE is for the groups $\Gamma_0(N)$ and $\Gamma_1(N)$. But perhaps one can still do something for small N, e.g. using that conjugation by the matrix with entries $(0,-1;N,0)$ takes $\Gamma(N)$ to $\Gamma_0(N^2) \cap \Gamma_1(N)$ -- this tells you how to get cusp forms for $\Gamma(N)$ both from $\Gamma_0(N^2)$ and $\Gamma_1(N)$, and by the genus formula for $\Gamma(N)$ you can tell if this has produced all of them.

To add to Aşağı Güzdək's nice answer, there is something you can do also when the genus is greater than zero which is also in some sense an exact formula. Namely, if $\Gamma$ is one of $\Gamma_0(N), \Gamma_1(N)$ or $\Gamma(N)$ [or probably some other congruence subgroups as well, but I am not familiar with this story], and $Y_\Gamma$ (resp. $X_\Gamma$) is the corresponding open (resp. compact) modular curve, then one has $$ |Y_\Gamma(\mathbf{F}_p)| = p+1-\mathrm{Tr}(\mathrm{Frob}_p|E_2(\Gamma))-\mathrm{Tr}(\mathrm{Frob}_p|S_2(\Gamma)),$$ where $E_k(\Gamma)$ (resp. $S_k(\Gamma)$ is the Galois representation attached to Eisenstein series (resp. cusp forms) of weight $k$ for $\Gamma$. The Galois representation $E_2(\Gamma)$ is given by the $\ell$-adic cohomology of $X^\infty_\Gamma = X_\Gamma \setminus Y_\Gamma$, hence the appearance of the number of cusps in the previous answer. In any case one can for any $p$ and $\Gamma$ give a bound for this trace by the number of cusps.

The trace of Frobenius on the space of cusp forms is the really interesting one. In many cases this can be computed for a given $\Gamma$ and $p$ by looking in tables like the Modular Forms Database of William Stein, or asking SAGE to spit them out for you (see e.g. http://www.sagemath.org/doc/reference/sage/modular/modform/cuspidal_submodule.html) -- for a prime p, the trace of $\mathrm{Frob_p}$ coincides with the trace of the Hecke operator $T_p$, which can be read off from the pth Fourier coefficients of a basis of normalized eigenforms for the space of cusp forms you are looking at.

To add to Aşağı Güzdək's nice answer, there is something you can do also when the genus is greater than zero which is also in some sense an exact formula. Namely, if $\Gamma$ is one of $\Gamma_0(N), \Gamma_1(N)$ or $\Gamma(N)$ [or probably some other congruence subgroups as well, but I am not familiar with this story], and $Y_\Gamma$ (resp. $X_\Gamma$) is the corresponding open (resp. compact) modular curve, then one has $$ |Y_\Gamma(\mathbf{F}_p)| = p+1-\mathrm{Tr}(\mathrm{Frob}_p|E_2(\Gamma))-\mathrm{Tr}(\mathrm{Frob}_p|S_2(\Gamma)),$$ where $E_k(\Gamma)$ (resp. $S_k(\Gamma)$) is the Galois representation attached to Eisenstein series (resp. cusp forms) of weight $k$ for $\Gamma$. The Galois representation $E_2(\Gamma)$ is given by the $\ell$-adic cohomology of $X^\infty_\Gamma = X_\Gamma \setminus Y_\Gamma$, hence the appearance of the number of cusps in the previous answer. In general one can for a given $\Gamma$ and p figure out which of the cusps are defined over what fields by looking explicitly at the modular interpretation of $X_\Gamma^\infty$ described in the paper by Deligne & Rapoport.

The trace of Frobenius on the space of cusp forms is the really interesting one. In many cases this can be computed for a given $\Gamma$ and $p$ by looking in tables like the Modular Forms Database of William Stein, or asking SAGE to spit them out for you (see e.g. http://www.sagemath.org/doc/reference/sage/modular/modform/cuspidal_submodule.html) -- for a prime p, the trace of $\mathrm{Frob_p}$ coincides with the trace of the Hecke operator $T_p$, which can be read off from the pth Fourier coefficients of a basis of normalized eigenforms for the space of cusp forms you are looking at.

It seems the only built-in functionality in SAGE is for the groups $\Gamma_0(N)$ and $\Gamma_1(N)$. But perhaps one can still do something for small N, e.g. using that conjugation by the matrix with entries $(0,-1;N,0)$ takes $\Gamma(N)$ to $\Gamma_0(N^2) \cap \Gamma_1(N)$ -- this tells you how to get cusp forms for $\Gamma(N)$ both from $\Gamma_0(N^2)$ and $\Gamma_1(N)$, and by the genus formulae for $\Gamma(N)$ you can tell if this has produced all of them.