beta_binomial_lpmf is unstable for large shape arguments

Description

As the shape arguments for the beta part get large (roughly > 1e14) the beta_binomial_lpmf function starts to return problematic values.

Example

beta_binomial_lpmf(500, 1000,  0.5 * 1e10,  0.5 * 1e10)
// -3.67992
beta_binomial_lpmf(500, 1000,  0.5 * 1e13,  0.5 * 1e13)
// -3.6802
beta_binomial_lpmf(500, 1000,  0.5 * 1e19,  0.5 * 1e19)
// 689.467

Expected Output

The values should converge towards binomial_lpmf(500, 1000, 0.5) = -3.67992.

Current Version:

v4.1.0

Proposed solution

It is possible that a similar trick to the one used in lbeta would be needed - compute separately a stirling approximation to the gamma functions involved, simplify analytically and then add corrections via lgamma_stirling_diff.

I'll try to start a PR with some tests demonstrating the problematic behaviour

Start by reproducing the large-shape beta_binomial_lpmf examples and compare them with binomial_lpmf. Read the lbeta implementation in stan/math/prim/fun/lbeta.hpp and the use of lgamma_stirling_diff, then add regression tests for the reported cases and verify convergence toward the expected binomial value.