When modeling real-world data, we face the challenge of determining which probability distribution best represents the data. To address this intricate problem, we rely on goodness-of-fit tests. However, when the data come from a bivariate negative binomial distribution, the literature reveals no existing goodness-of-fit test for this distribution. For this reason, in this article, we propose and study a computationally convenient goodness-of-fit test for the bivariate negative binomial distribution. This test is based on a bootstrap approximation and a parallelization strategy. To this end, we use a reparameterization technique based on the probability generating function and a Cram & eacute;r-von Mises-type statistic. From the simulation studies, we conclude that the results converge to the established nominal levels as the sample size increases, and in all cases considered, the parametric bootstrap method provides an accurate approximation of the null distribution of the statistic we propose. Additionally, we verify the power of the proposed test, as well as its application to five real datasets. To accelerate the massive computational work, we employ the parallelization strategy that, according to Novoa-Mu & ntilde;oz (2024), was the most efficient among the techniques he analyzed.