We develop an efficient numerical scheme to solve accurately the set of nonlinear integral equations derived previously in [A. Saichev and D. Sornette, J. Geophys. Res. 112, B04313 (2007)], which describes the distribution of interevent times in the framework of a general model of earthquake clustering with long memory. Detailed comparisons between the linear and nonlinear versions of the theory and direct synthetic catalogs show that the nonlinear theory provides an excellent fit to the synthetic catalogs, while there are significant biases resulting from the use of the linear approximation. We then address the suggestions proposed by some authors to use the empirical distribution of interevent times to obtain a better determination of the so-called clustering parameter. Our theory and tests against synthetic and empirical catalogs find a rather dramatic lack of power for the distribution of interevent times to distinguish between quite different sets of parameters, casting doubt on the usefulness of this statistic for the specific purpose of identifying the clustering parameter.