Numerical issues in neg_binomial_lccdf if beta is very large

Description

When the second argument to neg_binomial_lccdf is very large, the function can return -Inf even if the answer is finite.

Example

We take advantage of the fact the neg_binomial_lccdf(0|a, b) == log(1 - exp(neg_binomial_lpmf(0| a, b))
The formula then simplifies to neg_binomial_lccdf(0|a, b) == log1m( (b / (b + 1)) ^ a). Now for large b b/(b +1) is numerically exactly 1, but we can further rewrite as neg_binomial_lccdf(0|a, b) == log1m_exp(-a * log1p(1/b)).

So consider the following Stan program:

transformed data {
  real alpha = 1;
  real beta = 1e18;
  print(neg_binomial_lccdf(0 | alpha, beta));
  print(log1m_exp(-alpha * log1p(1/beta)));
}

Which outputs

-inf
-41.4465

Expected Output

-41.4465
-41.4465

Note that the two versions start to diverge only at around beta > 1e14, so this is unlikely to be super important in practice.

Current Version:

v3.3.0

Start at the neg_binomial_lccdf implementation and reproduce the Stan program from the issue with beta set to 1e18. The work is done when the function returns the finite value shown in the expected output, with coverage for this large-beta case.