  In this paper, we present a probabilistic algorithm to compute the number of
$\mathbb{F}_p$-points of modular curve $X_1(n)$. Under the Generalized Riemann
Hypothesis(GRH), the algorithm takes
$\textrm{O}(n^{56+\delta+\epsilon}\log^{9+\epsilon} p)$ bit operations, where
$\delta$ is an absolute constant and $\epsilon$ is any positive real number. As
an application, we can compute $#X_1(17)(\mathbb{F}_p)\textrm{mod} 17$ for huge
primes $p$. For example, we have
$#X_1(17)(\mathbb{F}_{10^{1000}+1357})\textrm{mod} 17=3$.
