gamma_p and gamma_q gradient problems

Description

gamma_p doesn't work great with vars and gamma_q doesn't work great with fvar<double>

Here is test code:

#include <test/unit/math/test_ad.hpp>
#include <limits>

TEST(mathMixScalFun, gamma_q) {
  using stan::math::var;
  using stan::math::fvar;
  using stan::math::gamma_q;
  using stan::math::gamma_p;
  using stan::math::gamma_lccdf;
  using stan::math::log;
  using stan::math::log1m;

  var yv = 12;
  var alphav = 3;
  var betav = 3;

  fvar<double> yf(12.0, 1.0);
  fvar<double> alphaf(3.0, 1.0);
  fvar<double> betaf(3.0, 1.0);

  var lp1 = log1m(gamma_p(alphav, betav * yv));
  fvar<double> lp2 = log(gamma_q(alphaf, betaf * yf));

  lp1.grad();

  std::cout << "Computed: " << std::endl;
  std::cout << "lp1.val(): " << lp1.val() << std::endl;
  std::cout << "lp2.val(): " << lp2.val() << std::endl;
  std::cout << "y.adj(): " << yv.adj() << std::endl;
  std::cout << "alpha.adj(): " << alphav.adj() << std::endl;
  std::cout << "beta.adj(): " << betav.adj() << std::endl;
  std::cout << "lp2.d_: " << lp2.d_ << std::endl;
}

TEST(mathMixScalFun, gamma_q_ref) {
  using stan::math::var;
  using stan::math::fvar;
  using stan::math::gamma_q;
  using stan::math::gamma_p;
  using stan::math::gamma_lccdf;
  using stan::math::log;
  using stan::math::log1m;

  var yv = 12;
  var alphav = 3;
  var betav = 3;

  fvar<double> yf(12.0, 1.0);
  fvar<double> alphaf(3.0, 1.0);
  fvar<double> betaf(3.0, 1.0);

  var lp1 = gamma_lccdf(yv, alphav, betav);
  fvar<double> lp2 = gamma_lccdf(yf, alphaf, betaf);

  lp1.grad();

  std::cout << "Reference: " << std::endl;
  std::cout << "lp1.val(): " << lp1.val() << std::endl;
  std::cout << "lp2.val(): " << lp2.val() << std::endl;
  std::cout << "y.adj(): " << yv.adj() << std::endl;
  std::cout << "alpha.adj(): " << alphav.adj() << std::endl;
  std::cout << "beta.adj(): " << betav.adj() << std::endl;
  std::cout << "lp2.d_: " << lp2.d_ << std::endl;
}

The output for gamma_p and gamma_q are (computing the same function, lp1 is with gamma_p and lp2 is with gamma_q):

Computed: 
lp1.val(): -29.4707
lp2.val(): -29.4706
y.adj(): 0
alpha.adj(): 0
beta.adj(): 0
lp2.d_: -6.70499e+12

And the reference values are here:

Reference: 
lp1.val(): -29.4706
lp2.val(): -29.4706
y.adj(): -2.83796
alpha.adj(): 2.68924
beta.adj(): -11.3518
lp2.d_: -11.5005

(edit: updated output)

Current Version:

v3.3.0

Start with the supplied mathMixScalFun gamma_q and gamma_q_ref reproducers, focusing on the gamma_p, gamma_q, and gamma_lccdf entry points. Run the examples and compare reverse-mode adjoints and forward-mode derivatives for y, alpha, and beta. Done means gamma_p and gamma_q agree with gamma_lccdf for both values and derivatives.