The derivatives of neg_binomial_lpmf are incorrect for large alpha

Description

The neg_binomial_lpmf(n, alpha, beta) is approximated by Poisson for alpha > 1e10. While this approximation is OK for the density, the derivatives get messed up.

Example

Analytically, the derivative wrt. alpha for n = 0 is log(beta / (1 + beta)). So, for beta = 0.8 the derivative should be roughly -0.81093 irrespective of alpha. However, when we set alpha = 1e10 + 1, the current implementation evaluates the gradient as -1.25.

Similarly (and perhaps more worrying) the derivative wrt. beta for n = 0 is (alpha * (1 + beta) *(-(beta/(1 + beta)^2) + 1/(1 + beta)))/beta, so for beta = 0.8 and alpha = 1e10 + 1 the analytic derivative is 6.94e9 while the current implementation in Stan returns 1.56e10.

(see bottom of https://www.wolframcloud.com/obj/martin.modrak/Published/NegBinomial_test.nb for Mathematica code I used to obtain those).

Expected Output

On a branch (https://github.com/martinmodrak/math/blob/bugfix/1578-neg_binomial_derivatives_at_cutoff/test/unit/math/rev/prob/neg_binomial_test.cpp) I wrote a test that compares the derivatives to those computed by complex step. The test is derived from the tests implemented in #1497 for neg_binomial_2_lpmf. A good question is: "Is the complex step derivative accurate for large alpha?" Unfortunately Mathematica (at least the free cloud version) refuses to compute exact derivatives for n > 0 and alpha = 1e10, but the complex step results match exact solutions for n = 0 and the test passes for values just below the alpha cutoff, so I believe, this is a true failure. I would expect the test to pass.

Current Version:

v3.0.0

Start with the neg_binomial_lpmf implementation and the test at test/unit/math/rev/prob/neg_binomial_test.cpp on the referenced branch. Run that test around the alpha cutoff and compare the alpha and beta derivatives with the complex-step results; done means the large-alpha derivative checks pass.