diff --git a/stan/math/prim/prob/chi_square_cdf.hpp b/stan/math/prim/prob/chi_square_cdf.hpp
index 7255340570a..c00d39cd1ab 100644
--- a/stan/math/prim/prob/chi_square_cdf.hpp
+++ b/stan/math/prim/prob/chi_square_cdf.hpp
@@ -53,6 +53,7 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
   // Explicit return for extreme values
@@ -63,16 +64,17 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
     }
   }
 
+  VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
+      tgamma_vec(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
-
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return alpha_dbl = value_of(nu_vec[i]) * 0.5;
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      digamma_vec(size_nu);
+  if (!is_constant_all<T_y, T_dof>::value) {
+    for (size_t i = 0; i < size_nu; i++) {
+      const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+      tgamma_vec[i] = tgamma(half_nu_dbl);
+      if (!is_constant_all<T_dof>::value) {
+        digamma_vec[i] = digamma(half_nu_dbl);
+      }
     }
   }
 
@@ -83,23 +85,21 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
       continue;
     }
 
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return alpha_dbl = value_of(nu_vec[n]) * 0.5;
-    const T_partials_return beta_dbl = 0.5;
-
-    const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl);
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[n]);
+    const T_partials_return half_y_dbl = 0.5 * value_of(y_vec[n]);
+    const T_partials_return Pn = gamma_p(half_nu_dbl, half_y_dbl);
 
     cdf *= Pn;
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] += beta_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge1_.partials_[n] += 0.5 * exp(-half_y_dbl)
+                                          * pow(half_y_dbl, half_nu_dbl - 1)
+                                          / (tgamma_vec[n] * Pn);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
           -= 0.5
-             * grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
+             * grad_reg_inc_gamma(half_nu_dbl, half_y_dbl, tgamma_vec[n],
                                   digamma_vec[n])
              / Pn;
     }
@@ -111,7 +111,7 @@ return_type_t<T_y, T_dof> chi_square_cdf(const T_y& y, const T_dof& nu) {
     }
   }
   if (!is_constant_all<T_dof>::value) {
-    for (size_t n = 0; n < stan::math::size(nu); ++n) {
+    for (size_t n = 0; n < size_nu; ++n) {
       ops_partials.edge2_.partials_[n] *= cdf;
     }
   }
diff --git a/stan/math/prim/prob/chi_square_lccdf.hpp b/stan/math/prim/prob/chi_square_lccdf.hpp
index b6f1e19c08d..cccb6b57e59 100644
--- a/stan/math/prim/prob/chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/chi_square_lccdf.hpp
@@ -55,53 +55,52 @@ return_type_t<T_y, T_dof> chi_square_lccdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
   for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (value_of(y_vec[i]) == 0) {
+    const T_partials_return y_dbl = value_of(y_vec[i]);
+    if (y_dbl == 0) {
       return ops_partials.build(0.0);
     }
+    if (y_dbl == INFTY) {
+      return ops_partials.build(NEGATIVE_INFTY);
+    }
   }
 
+  VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
+  VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
+      tgamma_vec(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
-
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return alpha_dbl = value_of(nu_vec[i]) * 0.5;
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    half_nu[i] = half_nu_dbl;
+    if (!is_constant_all<T_y, T_dof>::value) {
+      tgamma_vec[i] = tgamma(half_nu[i]);
+    }
+    if (!is_constant_all<T_dof>::value) {
+      digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
 
   for (size_t n = 0; n < N; n++) {
-    // Explicit results for extreme values
-    // The gradients are technically ill-defined, but treated as zero
-    if (value_of(y_vec[n]) == INFTY) {
-      return ops_partials.build(negative_infinity());
-    }
-
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return alpha_dbl = value_of(nu_vec[n]) * 0.5;
-    const T_partials_return beta_dbl = 0.5;
-
-    const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_dbl);
+    const T_partials_return half_y_dbl = 0.5 * value_of(y_vec[n]);
+    const T_partials_return Pn = gamma_q(half_nu[n], half_y_dbl);
 
     ccdf_log += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] -= beta_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge1_.partials_[n] -= 0.5 * exp(-half_y_dbl)
+                                          * pow(half_y_dbl, half_nu[n] - 1)
+                                          / (tgamma_vec[n] * Pn);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
           += 0.5
-             * grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
+             * grad_reg_inc_gamma(half_nu[n], half_y_dbl, tgamma_vec[n],
                                   digamma_vec[n])
              / Pn;
     }
diff --git a/stan/math/prim/prob/chi_square_lcdf.hpp b/stan/math/prim/prob/chi_square_lcdf.hpp
index 3361c126dd1..da4df4877c9 100644
--- a/stan/math/prim/prob/chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/chi_square_lcdf.hpp
@@ -55,53 +55,52 @@ return_type_t<T_y, T_dof> chi_square_lcdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
   for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (value_of(y_vec[i]) == 0) {
-      return ops_partials.build(negative_infinity());
+    const T_partials_return y_dbl = value_of(y_vec[i]);
+    if (y_dbl == 0) {
+      return ops_partials.build(NEGATIVE_INFTY);
+    }
+    if (y_dbl == INFTY) {
+      return ops_partials.build(0.0);
     }
   }
 
+  VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
+  VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
+      tgamma_vec(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
-
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return alpha_dbl = value_of(nu_vec[i]) * 0.5;
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    half_nu[i] = half_nu_dbl;
+    if (!is_constant_all<T_y, T_dof>::value) {
+      tgamma_vec[i] = tgamma(half_nu[i]);
+    }
+    if (!is_constant_all<T_dof>::value) {
+      digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
 
   for (size_t n = 0; n < N; n++) {
-    // Explicit results for extreme values
-    // The gradients are technically ill-defined, but treated as zero
-    if (value_of(y_vec[n]) == INFTY) {
-      return ops_partials.build(0.0);
-    }
-
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return alpha_dbl = value_of(nu_vec[n]) * 0.5;
-    const T_partials_return beta_dbl = 0.5;
-
-    const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl);
+    const T_partials_return half_y_dbl = 0.5 * value_of(y_vec[n]);
+    const T_partials_return Pn = gamma_p(half_nu[n], half_y_dbl);
 
     cdf_log += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] += beta_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge1_.partials_[n] += 0.5 * exp(-half_y_dbl)
+                                          * pow(half_y_dbl, half_nu[n] - 1)
+                                          / (tgamma_vec[n] * Pn);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
           -= 0.5
-             * grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
+             * grad_reg_inc_gamma(half_nu[n], half_y_dbl, tgamma_vec[n],
                                   digamma_vec[n])
              / Pn;
     }
diff --git a/stan/math/prim/prob/chi_square_log.hpp b/stan/math/prim/prob/chi_square_log.hpp
index 71dd2bc8af0..38e31e48ed5 100644
--- a/stan/math/prim/prob/chi_square_log.hpp
+++ b/stan/math/prim/prob/chi_square_log.hpp
@@ -8,25 +8,7 @@ namespace stan {
 namespace math {
 
 /** \ingroup prob_dists
- * The log of a chi-squared density for y with the specified
- * degrees of freedom parameter.
- * The degrees of freedom parameter must be greater than 0.
- * y must be greater than or equal to 0.
- *
- \f{eqnarray*}{
- y &\sim& \chi^2_\nu \\
- \log (p (y \, |\, \nu)) &=& \log \left( \frac{2^{-\nu / 2}}{\Gamma (\nu / 2)}
- y^{\nu / 2 - 1} \exp^{- y / 2} \right) \\
- &=& - \frac{\nu}{2} \log(2) - \log (\Gamma (\nu / 2)) + (\frac{\nu}{2} - 1)
- \log(y) - \frac{y}{2} \\ & & \mathrm{ where } \; y \ge 0 \f}
- *
  * @deprecated use <code>chi_square_lpdf</code>
- * @param y A scalar variable.
- * @param nu Degrees of freedom.
- * @throw std::domain_error if nu is not greater than or equal to 0
- * @throw std::domain_error if y is not greater than or equal to 0.
- * @tparam T_y Type of scalar.
- * @tparam T_dof Type of degrees of freedom.
  */
 template <bool propto, typename T_y, typename T_dof>
 return_type_t<T_y, T_dof> chi_square_log(const T_y& y, const T_dof& nu) {
diff --git a/stan/math/prim/prob/chi_square_lpdf.hpp b/stan/math/prim/prob/chi_square_lpdf.hpp
index 3362f8bb306..7a8a55259f7 100644
--- a/stan/math/prim/prob/chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/chi_square_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
@@ -60,64 +61,59 @@ return_type_t<T_y, T_dof> chi_square_lpdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
-  for (size_t n = 0; n < stan::math::size(y); n++) {
+  for (size_t n = 0; n < size_y; n++) {
     if (value_of(y_vec[n]) < 0) {
       return LOG_ZERO;
     }
   }
 
-  VectorBuilder<include_summand<propto, T_y, T_dof>::value, T_partials_return,
-                T_y>
-      log_y(size(y));
-  for (size_t i = 0; i < stan::math::size(y); i++) {
-    log_y[i] = log(value_of(y_vec[i]));
-  }
-
+  VectorBuilder<true, T_partials_return, T_y> log_y(size_y);
   VectorBuilder<include_summand<propto, T_y>::value, T_partials_return, T_y>
-      inv_y(size(y));
-  for (size_t i = 0; i < stan::math::size(y); i++) {
+      inv_y(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    const T_partials_return y_dbl = value_of(y_vec[i]);
+    log_y[i] = log(y_dbl);
     if (include_summand<propto, T_y>::value) {
-      inv_y[i] = 1.0 / value_of(y_vec[i]);
+      inv_y[i] = inv(y_dbl);
     }
   }
 
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      lgamma_half_nu(size(nu));
+      lgamma_half_nu(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_half_nu_over_two(size(nu));
-
-  for (size_t i = 0; i < stan::math::size(nu); i++) {
-    T_partials_return half_nu = 0.5 * value_of(nu_vec[i]);
-    if (include_summand<propto, T_dof>::value) {
+      digamma_half_nu(size_nu);
+  if (include_summand<propto, T_dof>::value) {
+    for (size_t i = 0; i < size_nu; i++) {
+      T_partials_return half_nu = 0.5 * value_of(nu_vec[i]);
       lgamma_half_nu[i] = lgamma(half_nu);
-    }
-    if (!is_constant_all<T_dof>::value) {
-      digamma_half_nu_over_two[i] = digamma(half_nu) * 0.5;
+      if (!is_constant_all<T_dof>::value) {
+        digamma_half_nu[i] = digamma(half_nu);
+      }
     }
   }
 
   for (size_t n = 0; n < N; n++) {
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return half_y = 0.5 * y_dbl;
     const T_partials_return nu_dbl = value_of(nu_vec[n]);
-    const T_partials_return half_nu = 0.5 * nu_dbl;
+    const T_partials_return half_nu_m1 = 0.5 * nu_dbl - 1;
+
     if (include_summand<propto, T_dof>::value) {
       logp -= nu_dbl * HALF_LOG_TWO + lgamma_half_nu[n];
     }
-    logp += (half_nu - 1.0) * log_y[n];
-
     if (include_summand<propto, T_y>::value) {
-      logp -= half_y;
+      logp -= 0.5 * value_of(y_vec[n]);
     }
+    logp += half_nu_m1 * log_y[n];
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] += (half_nu - 1.0) * inv_y[n] - 0.5;
+      ops_partials.edge1_.partials_[n] += half_nu_m1 * inv_y[n] - 0.5;
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
-          -= HALF_LOG_TWO + digamma_half_nu_over_two[n] - log_y[n] * 0.5;
+          -= 0.5 * (LOG_TWO + digamma_half_nu[n] - log_y[n]);
     }
   }
   return ops_partials.build(logp);
diff --git a/stan/math/prim/prob/gamma_cdf.hpp b/stan/math/prim/prob/gamma_cdf.hpp
index 1d8c412091e..fbae3d3086d 100644
--- a/stan/math/prim/prob/gamma_cdf.hpp
+++ b/stan/math/prim/prob/gamma_cdf.hpp
@@ -32,9 +32,8 @@ namespace math {
  * @param y A scalar variable.
  * @param alpha Shape parameter.
  * @param beta Inverse scale parameter.
- * @throw std::domain_error if alpha is not greater than 0.
- * @throw std::domain_error if beta is not greater than 0.
  * @throw std::domain_error if y is not greater than or equal to 0.
+ * @throw std::domain_error if alpha or beta is not greater than 0.
  */
 template <typename T_y, typename T_shape, typename T_inv_scale>
 return_type_t<T_y, T_shape, T_inv_scale> gamma_cdf(const T_y& y,
@@ -42,6 +41,7 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_cdf(const T_y& y,
                                                    const T_inv_scale& beta) {
   using T_partials_return = partials_return_t<T_y, T_shape, T_inv_scale>;
   using std::exp;
+  using std::pow;
   static const char* function = "gamma_cdf";
   check_positive_finite(function, "Shape parameter", alpha);
   check_positive_finite(function, "Inverse scale parameter", beta);
@@ -60,29 +60,30 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_cdf(const T_y& y,
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_inv_scale> beta_vec(beta);
+  size_t size_y = stan::math::size(y);
+  size_t size_alpha = stan::math::size(alpha);
   size_t N = max_size(y, alpha, beta);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
+  for (size_t i = 0; i < size_y; i++) {
     if (value_of(y_vec[i]) == 0) {
       return ops_partials.build(0.0);
     }
   }
 
-  using std::exp;
-  using std::pow;
-
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
+  VectorBuilder<!is_constant_all<T_y, T_shape, T_inv_scale>::value,
+                T_partials_return, T_shape>
+      tgamma_vec(size_alpha);
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
-
-  if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < stan::math::size(alpha); i++) {
+      digamma_vec(size_alpha);
+  if (!is_constant_all<T_y, T_shape, T_inv_scale>::value) {
+    for (size_t i = 0; i < size_alpha; i++) {
       const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      tgamma_vec[i] = tgamma(alpha_dbl);
+      if (!is_constant_all<T_shape>::value) {
+        digamma_vec[i] = digamma(alpha_dbl);
+      }
     }
   }
 
@@ -96,36 +97,36 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_cdf(const T_y& y,
     const T_partials_return y_dbl = value_of(y_vec[n]);
     const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
     const T_partials_return beta_dbl = value_of(beta_vec[n]);
-
-    const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl);
-
+    const T_partials_return beta_y = beta_dbl * y_dbl;
+    const T_partials_return Pn = gamma_p(alpha_dbl, beta_y);
+    const T_partials_return rep_deriv = is_constant_all<T_y, T_inv_scale>::value
+                                            ? 0
+                                            : exp(-beta_y)
+                                                  * pow(beta_y, alpha_dbl - 1)
+                                                  / (tgamma_vec[n] * Pn);
     P *= Pn;
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] += beta_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge1_.partials_[n] += beta_dbl * rep_deriv;
     }
     if (!is_constant_all<T_shape>::value) {
       ops_partials.edge2_.partials_[n]
-          -= grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
+          -= grad_reg_inc_gamma(alpha_dbl, beta_y, tgamma_vec[n],
                                 digamma_vec[n])
              / Pn;
     }
     if (!is_constant_all<T_inv_scale>::value) {
-      ops_partials.edge3_.partials_[n] += y_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge3_.partials_[n] += y_dbl * rep_deriv;
     }
   }
 
   if (!is_constant_all<T_y>::value) {
-    for (size_t n = 0; n < stan::math::size(y); ++n) {
+    for (size_t n = 0; n < size_y; ++n) {
       ops_partials.edge1_.partials_[n] *= P;
     }
   }
   if (!is_constant_all<T_shape>::value) {
-    for (size_t n = 0; n < stan::math::size(alpha); ++n) {
+    for (size_t n = 0; n < size_alpha; ++n) {
       ops_partials.edge2_.partials_[n] *= P;
     }
   }
diff --git a/stan/math/prim/prob/gamma_lccdf.hpp b/stan/math/prim/prob/gamma_lccdf.hpp
index bd3a78854ce..4d25bff1a44 100644
--- a/stan/math/prim/prob/gamma_lccdf.hpp
+++ b/stan/math/prim/prob/gamma_lccdf.hpp
@@ -46,6 +46,7 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_inv_scale> beta_vec(beta);
+  size_t size_alpha = stan::math::size(alpha);
   size_t N = max_size(y, alpha, beta);
 
   // Explicit return for extreme values
@@ -56,16 +57,18 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
     }
   }
 
+  VectorBuilder<!is_constant_all<T_y, T_shape, T_inv_scale>::value,
+                T_partials_return, T_shape>
+      tgamma_vec(size_alpha);
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
-
-  if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < stan::math::size(alpha); i++) {
+      digamma_vec(size_alpha);
+  if (!is_constant_all<T_y, T_shape, T_inv_scale>::value) {
+    for (size_t i = 0; i < size_alpha; i++) {
       const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      tgamma_vec[i] = tgamma(alpha_dbl);
+      if (!is_constant_all<T_shape>::value) {
+        digamma_vec[i] = digamma(alpha_dbl);
+      }
     }
   }
 
@@ -73,32 +76,32 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lccdf(const T_y& y,
     // Explicit results for extreme values
     // The gradients are technically ill-defined, but treated as zero
     if (value_of(y_vec[n]) == INFTY) {
-      return ops_partials.build(negative_infinity());
+      return ops_partials.build(NEGATIVE_INFTY);
     }
 
     const T_partials_return y_dbl = value_of(y_vec[n]);
     const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
     const T_partials_return beta_dbl = value_of(beta_vec[n]);
-
-    const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_dbl);
-
+    const T_partials_return beta_y = beta_dbl * y_dbl;
+    const T_partials_return Pn = gamma_q(alpha_dbl, beta_y);
+    const T_partials_return rep_deriv = is_constant_all<T_y, T_inv_scale>::value
+                                            ? 0
+                                            : exp(-beta_y)
+                                                  * pow(beta_y, alpha_dbl - 1)
+                                                  / (tgamma_vec[n] * Pn);
     P += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] -= beta_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge1_.partials_[n] -= beta_dbl * rep_deriv;
     }
     if (!is_constant_all<T_shape>::value) {
       ops_partials.edge2_.partials_[n]
-          += grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
+          += grad_reg_inc_gamma(alpha_dbl, beta_y, tgamma_vec[n],
                                 digamma_vec[n])
              / Pn;
     }
     if (!is_constant_all<T_inv_scale>::value) {
-      ops_partials.edge3_.partials_[n] -= y_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge3_.partials_[n] -= y_dbl * rep_deriv;
     }
   }
   return ops_partials.build(P);
diff --git a/stan/math/prim/prob/gamma_lcdf.hpp b/stan/math/prim/prob/gamma_lcdf.hpp
index cbce5a4e045..0cced485288 100644
--- a/stan/math/prim/prob/gamma_lcdf.hpp
+++ b/stan/math/prim/prob/gamma_lcdf.hpp
@@ -46,26 +46,29 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lcdf(const T_y& y,
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_inv_scale> beta_vec(beta);
+  size_t size_alpha = stan::math::size(alpha);
   size_t N = max_size(y, alpha, beta);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
   for (size_t i = 0; i < stan::math::size(y); i++) {
     if (value_of(y_vec[i]) == 0) {
-      return ops_partials.build(negative_infinity());
+      return ops_partials.build(NEGATIVE_INFTY);
     }
   }
 
+  VectorBuilder<!is_constant_all<T_y, T_shape, T_inv_scale>::value,
+                T_partials_return, T_shape>
+      tgamma_vec(size_alpha);
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
-
-  if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < stan::math::size(alpha); i++) {
+      digamma_vec(size_alpha);
+  if (!is_constant_all<T_y, T_shape, T_inv_scale>::value) {
+    for (size_t i = 0; i < size_alpha; i++) {
       const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      tgamma_vec[i] = tgamma(alpha_dbl);
+      if (!is_constant_all<T_shape>::value) {
+        digamma_vec[i] = digamma(alpha_dbl);
+      }
     }
   }
 
@@ -79,26 +82,26 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lcdf(const T_y& y,
     const T_partials_return y_dbl = value_of(y_vec[n]);
     const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
     const T_partials_return beta_dbl = value_of(beta_vec[n]);
-
-    const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_dbl);
-
+    const T_partials_return beta_y = beta_dbl * y_dbl;
+    const T_partials_return Pn = gamma_p(alpha_dbl, beta_y);
+    const T_partials_return rep_deriv = is_constant_all<T_y, T_inv_scale>::value
+                                            ? 0
+                                            : exp(-beta_y)
+                                                  * pow(beta_y, alpha_dbl - 1)
+                                                  / (tgamma_vec[n] * Pn);
     P += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n] += beta_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge1_.partials_[n] += beta_dbl * rep_deriv;
     }
     if (!is_constant_all<T_shape>::value) {
       ops_partials.edge2_.partials_[n]
-          -= grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_dbl, gamma_vec[n],
+          -= grad_reg_inc_gamma(alpha_dbl, beta_y, tgamma_vec[n],
                                 digamma_vec[n])
              / Pn;
     }
     if (!is_constant_all<T_inv_scale>::value) {
-      ops_partials.edge3_.partials_[n] += y_dbl * exp(-beta_dbl * y_dbl)
-                                          * pow(beta_dbl * y_dbl, alpha_dbl - 1)
-                                          / tgamma(alpha_dbl) / Pn;
+      ops_partials.edge3_.partials_[n] += y_dbl * rep_deriv;
     }
   }
   return ops_partials.build(P);
diff --git a/stan/math/prim/prob/gamma_lpdf.hpp b/stan/math/prim/prob/gamma_lpdf.hpp
index 0babd07d6fd..6871538b65f 100644
--- a/stan/math/prim/prob/gamma_lpdf.hpp
+++ b/stan/math/prim/prob/gamma_lpdf.hpp
@@ -67,45 +67,46 @@ return_type_t<T_y, T_shape, T_inv_scale> gamma_lpdf(const T_y& y,
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_inv_scale> beta_vec(beta);
+  size_t size_y = stan::math::size(y);
+  size_t size_alpha = stan::math::size(alpha);
+  size_t size_beta = stan::math::size(beta);
   size_t N = max_size(y, alpha, beta);
 
-  for (size_t n = 0; n < stan::math::size(y); n++) {
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    if (y_dbl < 0) {
+  for (size_t n = 0; n < size_y; n++) {
+    if (value_of(y_vec[n]) < 0) {
       return LOG_ZERO;
     }
   }
 
   VectorBuilder<include_summand<propto, T_y, T_shape>::value, T_partials_return,
                 T_y>
-      log_y(size(y));
+      log_y(size_y);
   if (include_summand<propto, T_y, T_shape>::value) {
-    for (size_t n = 0; n < stan::math::size(y); n++) {
-      if (value_of(y_vec[n]) > 0) {
-        log_y[n] = log(value_of(y_vec[n]));
-      }
+    for (size_t n = 0; n < size_y; n++) {
+      log_y[n] = log(value_of(y_vec[n]));
     }
   }
 
   VectorBuilder<include_summand<propto, T_shape>::value, T_partials_return,
                 T_shape>
-      lgamma_alpha(size(alpha));
+      lgamma_alpha(size_alpha);
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_alpha(size(alpha));
-  for (size_t n = 0; n < stan::math::size(alpha); n++) {
+      digamma_alpha(size_alpha);
+  for (size_t n = 0; n < size_alpha; n++) {
+    const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
     if (include_summand<propto, T_shape>::value) {
-      lgamma_alpha[n] = lgamma(value_of(alpha_vec[n]));
+      lgamma_alpha[n] = lgamma(alpha_dbl);
     }
     if (!is_constant_all<T_shape>::value) {
-      digamma_alpha[n] = digamma(value_of(alpha_vec[n]));
+      digamma_alpha[n] = digamma(alpha_dbl);
     }
   }
 
   VectorBuilder<include_summand<propto, T_shape, T_inv_scale>::value,
                 T_partials_return, T_inv_scale>
-      log_beta(size(beta));
+      log_beta(size_beta);
   if (include_summand<propto, T_shape, T_inv_scale>::value) {
-    for (size_t n = 0; n < stan::math::size(beta); n++) {
+    for (size_t n = 0; n < size_beta; n++) {
       log_beta[n] = log(value_of(beta_vec[n]));
     }
   }
diff --git a/stan/math/prim/prob/gamma_rng.hpp b/stan/math/prim/prob/gamma_rng.hpp
index 79fe6a11936..c06849b3bc1 100644
--- a/stan/math/prim/prob/gamma_rng.hpp
+++ b/stan/math/prim/prob/gamma_rng.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -47,8 +48,7 @@ inline typename VectorBuilder<true, double, T_shape, T_inv>::type gamma_rng(
   for (size_t n = 0; n < N; ++n) {
     // Convert rate (inverse scale) argument to scale for boost
     variate_generator<RNG&, gamma_distribution<> > gamma_rng(
-        rng, gamma_distribution<>(alpha_vec[n],
-                                  1 / static_cast<double>(beta_vec[n])));
+        rng, gamma_distribution<>(alpha_vec[n], inv(beta_vec[n])));
     output[n] = gamma_rng();
   }
 
diff --git a/stan/math/prim/prob/inv_chi_square_cdf.hpp b/stan/math/prim/prob/inv_chi_square_cdf.hpp
index fb279df3afd..820876ac450 100644
--- a/stan/math/prim/prob/inv_chi_square_cdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_cdf.hpp
@@ -8,9 +8,11 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/gamma_q.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/tgamma.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/operands_and_partials.hpp>
@@ -53,26 +55,36 @@ return_type_t<T_y, T_dof> inv_chi_square_cdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
+  for (size_t i = 0; i < size_y; i++) {
     if (value_of(y_vec[i]) == 0) {
       return ops_partials.build(0.0);
     }
   }
 
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+  VectorBuilder<true, T_partials_return, T_y> half_y_inv(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    half_y_inv[i] = 0.5 * inv(value_of(y_vec[i]));
+  }
 
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return nu_dbl = value_of(nu_vec[i]);
-      gamma_vec[i] = tgamma(0.5 * nu_dbl);
-      digamma_vec[i] = digamma(0.5 * nu_dbl);
+  VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
+  VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
+      tgamma_vec(size_nu);
+  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    half_nu[i] = half_nu_dbl;
+    if (!is_constant_all<T_y, T_dof>::value) {
+      tgamma_vec[i] = tgamma(half_nu[i]);
+    }
+    if (!is_constant_all<T_dof>::value) {
+      digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
 
@@ -83,36 +95,30 @@ return_type_t<T_y, T_dof> inv_chi_square_cdf(const T_y& y, const T_dof& nu) {
       continue;
     }
 
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return y_inv_dbl = 1.0 / y_dbl;
-    const T_partials_return nu_dbl = value_of(nu_vec[n]);
-
-    const T_partials_return Pn = gamma_q(0.5 * nu_dbl, 0.5 * y_inv_dbl);
-
+    const T_partials_return Pn = gamma_q(half_nu[n], half_y_inv[n]);
     P *= Pn;
 
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_[n]
-          += 0.5 * y_inv_dbl * y_inv_dbl * exp(-0.5 * y_inv_dbl)
-             * pow(0.5 * y_inv_dbl, 0.5 * nu_dbl - 1) / tgamma(0.5 * nu_dbl)
-             / Pn;
+          += 2 * square(half_y_inv[n]) * exp(-half_y_inv[n])
+             * pow(half_y_inv[n], half_nu[n] - 1) / (tgamma_vec[n] * Pn);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
           += 0.5
-             * grad_reg_inc_gamma(0.5 * nu_dbl, 0.5 * y_inv_dbl, gamma_vec[n],
+             * grad_reg_inc_gamma(half_nu[n], half_y_inv[n], tgamma_vec[n],
                                   digamma_vec[n])
              / Pn;
     }
   }
 
   if (!is_constant_all<T_y>::value) {
-    for (size_t n = 0; n < stan::math::size(y); ++n) {
+    for (size_t n = 0; n < size_y; ++n) {
       ops_partials.edge1_.partials_[n] *= P;
     }
   }
   if (!is_constant_all<T_dof>::value) {
-    for (size_t n = 0; n < stan::math::size(nu); ++n) {
+    for (size_t n = 0; n < size_nu; ++n) {
       ops_partials.edge2_.partials_[n] *= P;
     }
   }
diff --git a/stan/math/prim/prob/inv_chi_square_lccdf.hpp b/stan/math/prim/prob/inv_chi_square_lccdf.hpp
index 6500705a91a..ef401aecd02 100644
--- a/stan/math/prim/prob/inv_chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lccdf.hpp
@@ -8,10 +8,12 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/gamma_p.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/tgamma.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/operands_and_partials.hpp>
@@ -55,54 +57,56 @@ return_type_t<T_y, T_dof> inv_chi_square_lccdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (value_of(y_vec[i]) == 0) {
+  for (size_t i = 0; i < size_y; i++) {
+    const T_partials_return y_dbl = value_of(y_vec[i]);
+    if (y_dbl == 0) {
       return ops_partials.build(0.0);
     }
+    if (y_dbl == INFTY) {
+      return ops_partials.build(NEGATIVE_INFTY);
+    }
   }
 
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+  VectorBuilder<true, T_partials_return, T_y> half_y_inv(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    half_y_inv[i] = 0.5 * inv(value_of(y_vec[i]));
+  }
 
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return nu_dbl = value_of(nu_vec[i]);
-      gamma_vec[i] = tgamma(0.5 * nu_dbl);
-      digamma_vec[i] = digamma(0.5 * nu_dbl);
+  VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
+  VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
+      tgamma_vec(size_nu);
+  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    half_nu[i] = half_nu_dbl;
+    if (!is_constant_all<T_y, T_dof>::value) {
+      tgamma_vec[i] = tgamma(half_nu[i]);
+    }
+    if (!is_constant_all<T_dof>::value) {
+      digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
 
   for (size_t n = 0; n < N; n++) {
-    // Explicit results for extreme values
-    // The gradients are technically ill-defined, but treated as zero
-    if (value_of(y_vec[n]) == INFTY) {
-      return ops_partials.build(negative_infinity());
-    }
-
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return y_inv_dbl = 1.0 / y_dbl;
-    const T_partials_return nu_dbl = value_of(nu_vec[n]);
-
-    const T_partials_return Pn = gamma_p(0.5 * nu_dbl, 0.5 * y_inv_dbl);
-
+    const T_partials_return Pn = gamma_p(half_nu[n], half_y_inv[n]);
     P += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_[n]
-          -= 0.5 * y_inv_dbl * y_inv_dbl * exp(-0.5 * y_inv_dbl)
-             * pow(0.5 * y_inv_dbl, 0.5 * nu_dbl - 1) / tgamma(0.5 * nu_dbl)
-             / Pn;
+          -= 2 * square(half_y_inv[n]) * exp(-half_y_inv[n])
+             * pow(half_y_inv[n], half_nu[n] - 1) / (tgamma_vec[n] * Pn);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
           -= 0.5
-             * grad_reg_inc_gamma(0.5 * nu_dbl, 0.5 * y_inv_dbl, gamma_vec[n],
+             * grad_reg_inc_gamma(half_nu[n], half_y_inv[n], tgamma_vec[n],
                                   digamma_vec[n])
              / Pn;
     }
diff --git a/stan/math/prim/prob/inv_chi_square_lcdf.hpp b/stan/math/prim/prob/inv_chi_square_lcdf.hpp
index 51d1933be87..b0fb4d5d34e 100644
--- a/stan/math/prim/prob/inv_chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lcdf.hpp
@@ -8,10 +8,12 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/gamma_q.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/tgamma.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/operands_and_partials.hpp>
@@ -55,26 +57,36 @@ return_type_t<T_y, T_dof> inv_chi_square_lcdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
+  for (size_t i = 0; i < size_y; i++) {
     if (value_of(y_vec[i]) == 0) {
-      return ops_partials.build(negative_infinity());
+      return ops_partials.build(NEGATIVE_INFTY);
     }
   }
 
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+  VectorBuilder<true, T_partials_return, T_y> half_y_inv(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    half_y_inv[i] = 0.5 * inv(value_of(y_vec[i]));
+  }
 
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return nu_dbl = value_of(nu_vec[i]);
-      gamma_vec[i] = tgamma(0.5 * nu_dbl);
-      digamma_vec[i] = digamma(0.5 * nu_dbl);
+  VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
+  VectorBuilder<!is_constant_all<T_y, T_dof>::value, T_partials_return, T_dof>
+      tgamma_vec(size_nu);
+  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    half_nu[i] = half_nu_dbl;
+    if (!is_constant_all<T_y, T_dof>::value) {
+      tgamma_vec[i] = tgamma(half_nu[i]);
+    }
+    if (!is_constant_all<T_dof>::value) {
+      digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
 
@@ -85,24 +97,18 @@ return_type_t<T_y, T_dof> inv_chi_square_lcdf(const T_y& y, const T_dof& nu) {
       continue;
     }
 
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return y_inv_dbl = 1.0 / y_dbl;
-    const T_partials_return nu_dbl = value_of(nu_vec[n]);
-
-    const T_partials_return Pn = gamma_q(0.5 * nu_dbl, 0.5 * y_inv_dbl);
-
+    const T_partials_return Pn = gamma_q(half_nu[n], half_y_inv[n]);
     P += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_[n]
-          += 0.5 * y_inv_dbl * y_inv_dbl * exp(-0.5 * y_inv_dbl)
-             * pow(0.5 * y_inv_dbl, 0.5 * nu_dbl - 1) / tgamma(0.5 * nu_dbl)
-             / Pn;
+          += 2 * square(half_y_inv[n]) * exp(-half_y_inv[n])
+             * pow(half_y_inv[n], half_nu[n] - 1) / (tgamma_vec[n] * Pn);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
           += 0.5
-             * grad_reg_inc_gamma(0.5 * nu_dbl, 0.5 * y_inv_dbl, gamma_vec[n],
+             * grad_reg_inc_gamma(half_nu[n], half_y_inv[n], tgamma_vec[n],
                                   digamma_vec[n])
              / Pn;
     }
diff --git a/stan/math/prim/prob/inv_chi_square_lpdf.hpp b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
index 6f00dcaf0cc..3fa7c990c7c 100644
--- a/stan/math/prim/prob/inv_chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/inv_chi_square_lpdf.hpp
@@ -6,6 +6,7 @@
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/gamma_q.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
@@ -35,13 +36,13 @@ namespace math {
  * @tparam T_dof type of degrees of freedom
  * @param y A scalar variable.
  * @param nu Degrees of freedom.
- * @throw std::domain_error if nu is not greater than or equal to 0
- * @throw std::domain_error if y is not greater than or equal to 0.
+ * @throw std::domain_error if y is nan.
+ * @throw std::domain_error if nu is not positive.
  */
 template <bool propto, typename T_y, typename T_dof>
 return_type_t<T_y, T_dof> inv_chi_square_lpdf(const T_y& y, const T_dof& nu) {
-  using std::log;
   using T_partials_return = partials_return_t<T_y, T_dof>;
+  using std::log;
   static const char* function = "inv_chi_square_lpdf";
   check_positive_finite(function, "Degrees of freedom parameter", nu);
   check_not_nan(function, "Random variable", y);
@@ -57,54 +58,54 @@ return_type_t<T_y, T_dof> inv_chi_square_lpdf(const T_y& y, const T_dof& nu) {
 
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu);
 
-  for (size_t n = 0; n < stan::math::size(y); n++) {
+  for (size_t n = 0; n < size_y; n++) {
     if (value_of(y_vec[n]) <= 0) {
       return LOG_ZERO;
     }
   }
 
+  VectorBuilder<include_summand<propto, T_y>::value, T_partials_return, T_y>
+      inv_y(size_y);
   VectorBuilder<include_summand<propto, T_y, T_dof>::value, T_partials_return,
                 T_y>
-      log_y(size(y));
-  for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (include_summand<propto, T_y, T_dof>::value) {
-      log_y[i] = log(value_of(y_vec[i]));
-    }
-  }
-
-  VectorBuilder<include_summand<propto, T_y>::value, T_partials_return, T_y>
-      inv_y(size(y));
-  for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (include_summand<propto, T_y>::value) {
-      inv_y[i] = 1.0 / value_of(y_vec[i]);
+      log_y(size_y);
+  if (include_summand<propto, T_y, T_dof>::value) {
+    for (size_t i = 0; i < size_y; i++) {
+      const T_partials_return y_dbl = value_of(y_vec[i]);
+      log_y[i] = log(y_dbl);
+      if (include_summand<propto, T_y>::value) {
+        inv_y[i] = inv(y_dbl);
+      }
     }
   }
 
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      lgamma_half_nu(size(nu));
+      lgamma_half_nu(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_half_nu_over_two(size(nu));
-  for (size_t i = 0; i < stan::math::size(nu); i++) {
+      digamma_half_nu(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
     T_partials_return half_nu = 0.5 * value_of(nu_vec[i]);
     if (include_summand<propto, T_dof>::value) {
       lgamma_half_nu[i] = lgamma(half_nu);
     }
     if (!is_constant_all<T_dof>::value) {
-      digamma_half_nu_over_two[i] = digamma(half_nu) * 0.5;
+      digamma_half_nu[i] = digamma(half_nu);
     }
   }
 
   for (size_t n = 0; n < N; n++) {
     const T_partials_return nu_dbl = value_of(nu_vec[n]);
-    const T_partials_return half_nu = 0.5 * nu_dbl;
+    const T_partials_return half_nu_p1 = 0.5 * nu_dbl + 1.0;
 
     if (include_summand<propto, T_dof>::value) {
       logp -= nu_dbl * HALF_LOG_TWO + lgamma_half_nu[n];
     }
     if (include_summand<propto, T_y, T_dof>::value) {
-      logp -= (half_nu + 1.0) * log_y[n];
+      logp -= half_nu_p1 * log_y[n];
     }
     if (include_summand<propto, T_y>::value) {
       logp -= 0.5 * inv_y[n];
@@ -112,11 +113,11 @@ return_type_t<T_y, T_dof> inv_chi_square_lpdf(const T_y& y, const T_dof& nu) {
 
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_[n]
-          += -(half_nu + 1.0) * inv_y[n] + 0.5 * inv_y[n] * inv_y[n];
+          += (-half_nu_p1 + 0.5 * inv_y[n]) * inv_y[n];
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
-          -= HALF_LOG_TWO + digamma_half_nu_over_two[n] + 0.5 * log_y[n];
+          -= 0.5 * (LOG_TWO + digamma_half_nu[n] + log_y[n]);
     }
   }
   return ops_partials.build(logp);
diff --git a/stan/math/prim/prob/inv_gamma_cdf.hpp b/stan/math/prim/prob/inv_gamma_cdf.hpp
index a9cfcf04fb8..bacce4e3c3b 100644
--- a/stan/math/prim/prob/inv_gamma_cdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_cdf.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/gamma_q.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
@@ -60,26 +61,35 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_cdf(const T_y& y,
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_scale> beta_vec(beta);
+  size_t size_y = stan::math::size(y);
+  size_t size_alpha = stan::math::size(alpha);
   size_t N = max_size(y, alpha, beta);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
+  for (size_t i = 0; i < size_y; i++) {
     if (value_of(y_vec[i]) == 0) {
       return ops_partials.build(0.0);
     }
   }
 
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+  VectorBuilder<true, T_partials_return, T_y> inv_y(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    inv_y[i] = inv(value_of(y_vec[i]));
+  }
 
-  if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < stan::math::size(alpha); i++) {
+  VectorBuilder<!is_constant_all<T_y, T_shape, T_scale>::value,
+                T_partials_return, T_shape>
+      tgamma_vec(size_alpha);
+  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
+      digamma_vec(size_alpha);
+  if (!is_constant_all<T_y, T_shape, T_scale>::value) {
+    for (size_t i = 0; i < size_alpha; i++) {
       const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      tgamma_vec[i] = tgamma(alpha_dbl);
+      if (!is_constant_all<T_shape>::value) {
+        digamma_vec[i] = digamma(alpha_dbl);
+      }
     }
   }
 
@@ -90,42 +100,37 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_cdf(const T_y& y,
       continue;
     }
 
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return y_inv_dbl = 1.0 / y_dbl;
     const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
-    const T_partials_return beta_dbl = value_of(beta_vec[n]);
-
-    const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_inv_dbl);
-
+    const T_partials_return beta_over_y = value_of(beta_vec[n]) * inv_y[n];
+    const T_partials_return Pn = gamma_q(alpha_dbl, beta_over_y);
+    const T_partials_return rep_deriv
+        = is_constant_all<T_y, T_scale>::value
+              ? 0
+              : inv_y[n] * exp(-beta_over_y) * pow(beta_over_y, alpha_dbl - 1)
+                    / (tgamma_vec[n] * Pn);
     P *= Pn;
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n]
-          += beta_dbl * y_inv_dbl * y_inv_dbl * exp(-beta_dbl * y_inv_dbl)
-             * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) / tgamma(alpha_dbl)
-             / Pn;
+      ops_partials.edge1_.partials_[n] += beta_over_y * rep_deriv;
     }
     if (!is_constant_all<T_shape>::value) {
       ops_partials.edge2_.partials_[n]
-          += grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_inv_dbl, gamma_vec[n],
+          += grad_reg_inc_gamma(alpha_dbl, beta_over_y, tgamma_vec[n],
                                 digamma_vec[n])
              / Pn;
     }
     if (!is_constant_all<T_scale>::value) {
-      ops_partials.edge3_.partials_[n]
-          += -y_inv_dbl * exp(-beta_dbl * y_inv_dbl)
-             * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) / tgamma(alpha_dbl)
-             / Pn;
+      ops_partials.edge3_.partials_[n] -= rep_deriv;
     }
   }
 
   if (!is_constant_all<T_y>::value) {
-    for (size_t n = 0; n < stan::math::size(y); ++n) {
+    for (size_t n = 0; n < size_y; ++n) {
       ops_partials.edge1_.partials_[n] *= P;
     }
   }
   if (!is_constant_all<T_shape>::value) {
-    for (size_t n = 0; n < stan::math::size(alpha); ++n) {
+    for (size_t n = 0; n < size_alpha; ++n) {
       ops_partials.edge2_.partials_[n] *= P;
     }
   }
diff --git a/stan/math/prim/prob/inv_gamma_lccdf.hpp b/stan/math/prim/prob/inv_gamma_lccdf.hpp
index 014494e52d6..b4dd9156a3b 100644
--- a/stan/math/prim/prob/inv_gamma_lccdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lccdf.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/gamma_p.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
@@ -46,62 +47,65 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lccdf(const T_y& y,
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_scale> beta_vec(beta);
+  size_t size_y = stan::math::size(y);
+  size_t size_alpha = stan::math::size(alpha);
   size_t N = max_size(y, alpha, beta);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (value_of(y_vec[i]) == 0) {
+  for (size_t i = 0; i < size_y; i++) {
+    const T_partials_return y_dbl = value_of(y_vec[i]);
+    if (y_dbl == 0) {
       return ops_partials.build(0.0);
     }
+    if (y_dbl == INFTY) {
+      return ops_partials.build(NEGATIVE_INFTY);
+    }
   }
 
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+  VectorBuilder<true, T_partials_return, T_y> inv_y(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    inv_y[i] = inv(value_of(y_vec[i]));
+  }
 
-  if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < stan::math::size(alpha); i++) {
+  VectorBuilder<!is_constant_all<T_y, T_shape, T_scale>::value,
+                T_partials_return, T_shape>
+      tgamma_vec(size_alpha);
+  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
+      digamma_vec(size_alpha);
+  if (!is_constant_all<T_y, T_shape, T_scale>::value) {
+    for (size_t i = 0; i < size_alpha; i++) {
       const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      tgamma_vec[i] = tgamma(alpha_dbl);
+      if (!is_constant_all<T_shape>::value) {
+        digamma_vec[i] = digamma(alpha_dbl);
+      }
     }
   }
 
   for (size_t n = 0; n < N; n++) {
-    // Explicit results for extreme values
-    // The gradients are technically ill-defined, but treated as zero
-    if (value_of(y_vec[n]) == INFTY) {
-      return ops_partials.build(negative_infinity());
-    }
-
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return y_inv_dbl = 1.0 / y_dbl;
     const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
-    const T_partials_return beta_dbl = value_of(beta_vec[n]);
-
-    const T_partials_return Pn = gamma_p(alpha_dbl, beta_dbl * y_inv_dbl);
+    const T_partials_return beta_over_y = value_of(beta_vec[n]) * inv_y[n];
+    const T_partials_return Pn = gamma_p(alpha_dbl, beta_over_y);
+    const T_partials_return rep_deriv
+        = is_constant_all<T_y, T_scale>::value
+              ? 0
+              : inv_y[n] * exp(-beta_over_y) * pow(beta_over_y, alpha_dbl - 1)
+                    / (tgamma_vec[n] * Pn);
 
     P += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n]
-          -= beta_dbl * y_inv_dbl * y_inv_dbl * exp(-beta_dbl * y_inv_dbl)
-             * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) / tgamma(alpha_dbl)
-             / Pn;
+      ops_partials.edge1_.partials_[n] -= beta_over_y * rep_deriv;
     }
     if (!is_constant_all<T_shape>::value) {
       ops_partials.edge2_.partials_[n]
-          -= grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_inv_dbl, gamma_vec[n],
+          -= grad_reg_inc_gamma(alpha_dbl, beta_over_y, tgamma_vec[n],
                                 digamma_vec[n])
              / Pn;
     }
     if (!is_constant_all<T_scale>::value) {
-      ops_partials.edge3_.partials_[n]
-          += y_inv_dbl * exp(-beta_dbl * y_inv_dbl)
-             * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) / tgamma(alpha_dbl)
-             / Pn;
+      ops_partials.edge3_.partials_[n] += rep_deriv;
     }
   }
   return ops_partials.build(P);
diff --git a/stan/math/prim/prob/inv_gamma_lcdf.hpp b/stan/math/prim/prob/inv_gamma_lcdf.hpp
index e2108c178f8..7071de97790 100644
--- a/stan/math/prim/prob/inv_gamma_lcdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lcdf.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/gamma_q.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
@@ -46,26 +47,35 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lcdf(const T_y& y,
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_scale> beta_vec(beta);
+  size_t size_y = stan::math::size(y);
+  size_t size_alpha = stan::math::size(alpha);
   size_t N = max_size(y, alpha, beta);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
+  for (size_t i = 0; i < size_y; i++) {
     if (value_of(y_vec[i]) == 0) {
-      return ops_partials.build(negative_infinity());
+      return ops_partials.build(NEGATIVE_INFTY);
     }
   }
 
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      gamma_vec(size(alpha));
-  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_vec(size(alpha));
+  VectorBuilder<true, T_partials_return, T_y> inv_y(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    inv_y[i] = inv(value_of(y_vec[i]));
+  }
 
-  if (!is_constant_all<T_shape>::value) {
-    for (size_t i = 0; i < stan::math::size(alpha); i++) {
+  VectorBuilder<!is_constant_all<T_y, T_shape, T_scale>::value,
+                T_partials_return, T_shape>
+      tgamma_vec(size_alpha);
+  VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
+      digamma_vec(size_alpha);
+  if (!is_constant_all<T_y, T_shape, T_scale>::value) {
+    for (size_t i = 0; i < size_alpha; i++) {
       const T_partials_return alpha_dbl = value_of(alpha_vec[i]);
-      gamma_vec[i] = tgamma(alpha_dbl);
-      digamma_vec[i] = digamma(alpha_dbl);
+      tgamma_vec[i] = tgamma(alpha_dbl);
+      if (!is_constant_all<T_shape>::value) {
+        digamma_vec[i] = digamma(alpha_dbl);
+      }
     }
   }
 
@@ -76,32 +86,28 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lcdf(const T_y& y,
       continue;
     }
 
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return y_inv_dbl = 1.0 / y_dbl;
     const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
-    const T_partials_return beta_dbl = value_of(beta_vec[n]);
-
-    const T_partials_return Pn = gamma_q(alpha_dbl, beta_dbl * y_inv_dbl);
+    const T_partials_return beta_over_y = value_of(beta_vec[n]) * inv_y[n];
+    const T_partials_return Pn = gamma_q(alpha_dbl, beta_over_y);
+    const T_partials_return rep_deriv
+        = is_constant_all<T_y, T_scale>::value
+              ? 0
+              : inv_y[n] * exp(-beta_over_y) * pow(beta_over_y, alpha_dbl - 1)
+                    / (tgamma_vec[n] * Pn);
 
     P += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
-      ops_partials.edge1_.partials_[n]
-          += beta_dbl * y_inv_dbl * y_inv_dbl * exp(-beta_dbl * y_inv_dbl)
-             * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) / tgamma(alpha_dbl)
-             / Pn;
+      ops_partials.edge1_.partials_[n] += beta_over_y * rep_deriv;
     }
     if (!is_constant_all<T_shape>::value) {
       ops_partials.edge2_.partials_[n]
-          += grad_reg_inc_gamma(alpha_dbl, beta_dbl * y_inv_dbl, gamma_vec[n],
+          += grad_reg_inc_gamma(alpha_dbl, beta_over_y, tgamma_vec[n],
                                 digamma_vec[n])
              / Pn;
     }
     if (!is_constant_all<T_scale>::value) {
-      ops_partials.edge3_.partials_[n]
-          += -y_inv_dbl * exp(-beta_dbl * y_inv_dbl)
-             * pow(beta_dbl * y_inv_dbl, alpha_dbl - 1) / tgamma(alpha_dbl)
-             / Pn;
+      ops_partials.edge3_.partials_[n] -= rep_deriv;
     }
   }
   return ops_partials.build(P);
diff --git a/stan/math/prim/prob/inv_gamma_lpdf.hpp b/stan/math/prim/prob/inv_gamma_lpdf.hpp
index 8e6b71a8195..1b68dd8bed0 100644
--- a/stan/math/prim/prob/inv_gamma_lpdf.hpp
+++ b/stan/math/prim/prob/inv_gamma_lpdf.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
@@ -23,15 +24,15 @@ namespace math {
  * Shape and scale parameters must be greater than 0.
  * y must be greater than 0.
  *
+ * @tparam T_y type of scalar
+ * @tparam T_shape type of shape
+ * @tparam T_scale type of scale
  * @param y A scalar variable.
  * @param alpha Shape parameter.
  * @param beta Scale parameter.
  * @throw std::domain_error if alpha is not greater than 0.
  * @throw std::domain_error if beta is not greater than 0.
  * @throw std::domain_error if y is not greater than 0.
- * @tparam T_y Type of scalar.
- * @tparam T_shape Type of shape.
- * @tparam T_scale Type of scale.
  */
 template <bool propto, typename T_y, typename T_shape, typename T_scale>
 return_type_t<T_y, T_shape, T_scale> inv_gamma_lpdf(const T_y& y,
@@ -59,51 +60,53 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lpdf(const T_y& y,
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_shape> alpha_vec(alpha);
   scalar_seq_view<T_scale> beta_vec(beta);
+  size_t size_y = stan::math::size(y);
+  size_t size_alpha = stan::math::size(alpha);
+  size_t size_beta = stan::math::size(beta);
   size_t N = max_size(y, alpha, beta);
 
-  for (size_t n = 0; n < stan::math::size(y); n++) {
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    if (y_dbl <= 0) {
+  for (size_t n = 0; n < size_y; n++) {
+    if (value_of(y_vec[n]) <= 0) {
       return LOG_ZERO;
     }
   }
 
   VectorBuilder<include_summand<propto, T_y, T_shape>::value, T_partials_return,
                 T_y>
-      log_y(size(y));
+      log_y(size_y);
   VectorBuilder<include_summand<propto, T_y, T_scale>::value, T_partials_return,
                 T_y>
-      inv_y(size(y));
-  for (size_t n = 0; n < stan::math::size(y); n++) {
+      inv_y(size_y);
+  for (size_t n = 0; n < size_y; n++) {
+    const T_partials_return y_dbl = value_of(y_vec[n]);
     if (include_summand<propto, T_y, T_shape>::value) {
-      if (value_of(y_vec[n]) > 0) {
-        log_y[n] = log(value_of(y_vec[n]));
-      }
+      log_y[n] = log(y_dbl);
     }
     if (include_summand<propto, T_y, T_scale>::value) {
-      inv_y[n] = 1.0 / value_of(y_vec[n]);
+      inv_y[n] = inv(y_dbl);
     }
   }
 
   VectorBuilder<include_summand<propto, T_shape>::value, T_partials_return,
                 T_shape>
-      lgamma_alpha(size(alpha));
+      lgamma_alpha(size_alpha);
   VectorBuilder<!is_constant_all<T_shape>::value, T_partials_return, T_shape>
-      digamma_alpha(size(alpha));
-  for (size_t n = 0; n < stan::math::size(alpha); n++) {
+      digamma_alpha(size_alpha);
+  for (size_t n = 0; n < size_alpha; n++) {
+    const T_partials_return alpha_dbl = value_of(alpha_vec[n]);
     if (include_summand<propto, T_shape>::value) {
-      lgamma_alpha[n] = lgamma(value_of(alpha_vec[n]));
+      lgamma_alpha[n] = lgamma(alpha_dbl);
     }
     if (!is_constant_all<T_shape>::value) {
-      digamma_alpha[n] = digamma(value_of(alpha_vec[n]));
+      digamma_alpha[n] = digamma(alpha_dbl);
     }
   }
 
   VectorBuilder<include_summand<propto, T_shape, T_scale>::value,
                 T_partials_return, T_scale>
-      log_beta(size(beta));
+      log_beta(size_beta);
   if (include_summand<propto, T_shape, T_scale>::value) {
-    for (size_t n = 0; n < stan::math::size(beta); n++) {
+    for (size_t n = 0; n < size_beta; n++) {
       log_beta[n] = log(value_of(beta_vec[n]));
     }
   }
@@ -127,7 +130,7 @@ return_type_t<T_y, T_shape, T_scale> inv_gamma_lpdf(const T_y& y,
 
     if (!is_constant_all<scalar_type_t<T_y>>::value) {
       ops_partials.edge1_.partials_[n]
-          += -(alpha_dbl + 1) * inv_y[n] + beta_dbl * inv_y[n] * inv_y[n];
+          -= (alpha_dbl + 1 - beta_dbl * inv_y[n]) * inv_y[n];
     }
     if (!is_constant_all<scalar_type_t<T_shape>>::value) {
       ops_partials.edge2_.partials_[n]
diff --git a/stan/math/prim/prob/inv_gamma_rng.hpp b/stan/math/prim/prob/inv_gamma_rng.hpp
index 9adf2a3932f..4fc040787fc 100644
--- a/stan/math/prim/prob/inv_gamma_rng.hpp
+++ b/stan/math/prim/prob/inv_gamma_rng.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -17,8 +18,8 @@ namespace math {
  * alpha and beta can each be a scalar or a one-dimensional container. Any
  * non-scalar inputs must be the same size.
  *
- * @tparam T_shape Type of shape parameter
- * @tparam T_scale Type of scale parameter
+ * @tparam T_shape type of shape parameter
+ * @tparam T_scale type of scale parameter
  * @tparam RNG type of random number generator
  * @param alpha (Sequence of) positive shape parameter(s)
  * @param beta (Sequence of) positive scale parameter(s)
@@ -45,10 +46,9 @@ inv_gamma_rng(const T_shape& alpha, const T_scale& beta, RNG& rng) {
   VectorBuilder<true, double, T_shape, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, gamma_distribution<> > gamma_rng(
-        rng, gamma_distribution<>(alpha_vec[n],
-                                  1 / static_cast<double>(beta_vec[n])));
-    output[n] = 1 / gamma_rng();
+    variate_generator<RNG&, gamma_distribution<>> gamma_rng(
+        rng, gamma_distribution<>(alpha_vec[n], inv(beta_vec[n])));
+    output[n] = inv(gamma_rng());
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
index 4eb742443fc..673558babd3 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_cdf.hpp
@@ -8,9 +8,11 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/gamma_q.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/tgamma.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/operands_and_partials.hpp>
@@ -25,6 +27,7 @@ namespace math {
  *
  * @tparam T_y type of scalar
  * @tparam T_dof type of degrees of freedom
+ * @tparam T_scale type of scale
  * @param y A scalar variable.
  * @param nu Degrees of freedom.
  * @param s Scale parameter.
@@ -58,25 +61,30 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_cdf(const T_y& y,
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
   scalar_seq_view<T_scale> s_vec(s);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu, s);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
+  for (size_t i = 0; i < size_y; i++) {
     if (value_of(y_vec[i]) == 0) {
       return ops_partials.build(0.0);
     }
   }
 
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+  VectorBuilder<true, T_partials_return, T_y> inv_y(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    inv_y[i] = inv(value_of(y_vec[i]));
+  }
 
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
-      gamma_vec[i] = tgamma(half_nu_dbl);
+  VectorBuilder<true, T_partials_return, T_dof> tgamma_vec(size_nu);
+  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    tgamma_vec[i] = tgamma(half_nu_dbl);
+    if (!is_constant_all<T_dof>::value) {
       digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
@@ -88,38 +96,34 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_cdf(const T_y& y,
       continue;
     }
 
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return y_inv_dbl = 1.0 / y_dbl;
     const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[n]);
     const T_partials_return s_dbl = value_of(s_vec[n]);
-    const T_partials_return half_s2_overx_dbl = 0.5 * s_dbl * s_dbl * y_inv_dbl;
-    const T_partials_return half_nu_s2_overx_dbl
-        = 2.0 * half_nu_dbl * half_s2_overx_dbl;
-
-    const T_partials_return Pn = gamma_q(half_nu_dbl, half_nu_s2_overx_dbl);
+    const T_partials_return s2_over_y = square(s_dbl) * inv_y[n];
+    const T_partials_return half_nu_s2_over_y = half_nu_dbl * s2_over_y;
+    const T_partials_return Pn = gamma_q(half_nu_dbl, half_nu_s2_over_y);
     const T_partials_return gamma_p_deriv
-        = exp(-half_nu_s2_overx_dbl)
-          * pow(half_nu_s2_overx_dbl, half_nu_dbl - 1) / tgamma(half_nu_dbl);
+        = exp(-half_nu_s2_over_y) * pow(half_nu_s2_over_y, half_nu_dbl - 1)
+          / tgamma_vec[n];
 
     P *= Pn;
 
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_[n]
-          += half_nu_s2_overx_dbl * y_inv_dbl * gamma_p_deriv / Pn;
+          += half_nu_s2_over_y * inv_y[n] * gamma_p_deriv / Pn;
     }
 
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
-          += (0.5
-                  * grad_reg_inc_gamma(half_nu_dbl, half_nu_s2_overx_dbl,
-                                       gamma_vec[n], digamma_vec[n])
-              - half_s2_overx_dbl * gamma_p_deriv)
+          += 0.5
+             * (grad_reg_inc_gamma(half_nu_dbl, half_nu_s2_over_y,
+                                   tgamma_vec[n], digamma_vec[n])
+                - s2_over_y * gamma_p_deriv)
              / Pn;
     }
 
     if (!is_constant_all<T_scale>::value) {
       ops_partials.edge3_.partials_[n]
-          += -2.0 * half_nu_dbl * s_dbl * y_inv_dbl * gamma_p_deriv / Pn;
+          += -2.0 * half_nu_dbl * s_dbl * inv_y[n] * gamma_p_deriv / Pn;
     }
   }
 
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
index e96fb4963ae..89562a8a9a6 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lccdf.hpp
@@ -8,10 +8,12 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/gamma_p.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/tgamma.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/operands_and_partials.hpp>
@@ -46,66 +48,65 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lccdf(
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
   scalar_seq_view<T_scale> s_vec(s);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu, s);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (value_of(y_vec[i]) == 0) {
+  for (size_t i = 0; i < size_y; i++) {
+    const T_partials_return y_dbl = value_of(y_vec[i]);
+    if (y_dbl == 0) {
       return ops_partials.build(0.0);
     }
+    if (y_dbl == INFTY) {
+      return ops_partials.build(NEGATIVE_INFTY);
+    }
   }
 
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+  VectorBuilder<true, T_partials_return, T_y> inv_y(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    inv_y[i] = inv(value_of(y_vec[i]));
+  }
 
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
-      gamma_vec[i] = tgamma(half_nu_dbl);
+  VectorBuilder<true, T_partials_return, T_dof> tgamma_vec(size_nu);
+  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    tgamma_vec[i] = tgamma(half_nu_dbl);
+    if (!is_constant_all<T_dof>::value) {
       digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
 
   for (size_t n = 0; n < N; n++) {
-    // Explicit results for extreme values
-    // The gradients are technically ill-defined, but treated as zero
-    if (value_of(y_vec[n]) == INFTY) {
-      return ops_partials.build(negative_infinity());
-    }
-
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return y_inv_dbl = 1.0 / y_dbl;
     const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[n]);
     const T_partials_return s_dbl = value_of(s_vec[n]);
-    const T_partials_return half_s2_overx_dbl = 0.5 * s_dbl * s_dbl * y_inv_dbl;
-    const T_partials_return half_nu_s2_overx_dbl
-        = 2.0 * half_nu_dbl * half_s2_overx_dbl;
-
-    const T_partials_return Pn = gamma_p(half_nu_dbl, half_nu_s2_overx_dbl);
+    const T_partials_return s2_over_y = square(s_dbl) * inv_y[n];
+    const T_partials_return half_nu_s2_over_y = half_nu_dbl * s2_over_y;
+    const T_partials_return Pn = gamma_p(half_nu_dbl, half_nu_s2_over_y);
     const T_partials_return gamma_p_deriv
-        = exp(-half_nu_s2_overx_dbl)
-          * pow(half_nu_s2_overx_dbl, half_nu_dbl - 1) / tgamma(half_nu_dbl);
+        = exp(-half_nu_s2_over_y) * pow(half_nu_s2_over_y, half_nu_dbl - 1)
+          / tgamma_vec[n];
 
     P += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_[n]
-          -= half_nu_s2_overx_dbl * y_inv_dbl * gamma_p_deriv / Pn;
+          -= half_nu_s2_over_y * inv_y[n] * gamma_p_deriv / Pn;
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
-          -= (0.5
-                  * grad_reg_inc_gamma(half_nu_dbl, half_nu_s2_overx_dbl,
-                                       gamma_vec[n], digamma_vec[n])
-              - half_s2_overx_dbl * gamma_p_deriv)
+          -= 0.5
+             * (grad_reg_inc_gamma(half_nu_dbl, half_nu_s2_over_y,
+                                   tgamma_vec[n], digamma_vec[n])
+                - s2_over_y * gamma_p_deriv)
              / Pn;
     }
     if (!is_constant_all<T_scale>::value) {
       ops_partials.edge3_.partials_[n]
-          += 2.0 * half_nu_dbl * s_dbl * y_inv_dbl * gamma_p_deriv / Pn;
+          += 2.0 * half_nu_dbl * s_dbl * inv_y[n] * gamma_p_deriv / Pn;
     }
   }
   return ops_partials.build(P);
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
index 1143cacdf29..454718e3eaf 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lcdf.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/exp.hpp>
 #include <stan/math/prim/fun/gamma_q.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
@@ -46,25 +47,30 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lcdf(
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
   scalar_seq_view<T_scale> s_vec(s);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
   size_t N = max_size(y, nu, s);
 
   // Explicit return for extreme values
   // The gradients are technically ill-defined, but treated as zero
-  for (size_t i = 0; i < stan::math::size(y); i++) {
+  for (size_t i = 0; i < size_y; i++) {
     if (value_of(y_vec[i]) == 0) {
-      return ops_partials.build(negative_infinity());
+      return ops_partials.build(NEGATIVE_INFTY);
     }
   }
 
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      gamma_vec(size(nu));
-  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_vec(size(nu));
+  VectorBuilder<true, T_partials_return, T_y> inv_y(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    inv_y[i] = inv(value_of(y_vec[i]));
+  }
 
-  if (!is_constant_all<T_dof>::value) {
-    for (size_t i = 0; i < stan::math::size(nu); i++) {
-      const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
-      gamma_vec[i] = tgamma(half_nu_dbl);
+  VectorBuilder<true, T_partials_return, T_dof> tgamma_vec(size_nu);
+  VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
+      digamma_vec(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[i]);
+    tgamma_vec[i] = tgamma(half_nu_dbl);
+    if (!is_constant_all<T_dof>::value) {
       digamma_vec[i] = digamma(half_nu_dbl);
     }
   }
@@ -76,36 +82,32 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lcdf(
       continue;
     }
 
-    const T_partials_return y_dbl = value_of(y_vec[n]);
-    const T_partials_return y_inv_dbl = 1.0 / y_dbl;
     const T_partials_return half_nu_dbl = 0.5 * value_of(nu_vec[n]);
     const T_partials_return s_dbl = value_of(s_vec[n]);
-    const T_partials_return half_s2_overx_dbl = 0.5 * s_dbl * s_dbl * y_inv_dbl;
-    const T_partials_return half_nu_s2_overx_dbl
-        = 2.0 * half_nu_dbl * half_s2_overx_dbl;
-
-    const T_partials_return Pn = gamma_q(half_nu_dbl, half_nu_s2_overx_dbl);
+    const T_partials_return s2_over_y = square(s_dbl) * inv_y[n];
+    const T_partials_return half_nu_s2_over_y = half_nu_dbl * s2_over_y;
+    const T_partials_return Pn = gamma_q(half_nu_dbl, half_nu_s2_over_y);
     const T_partials_return gamma_p_deriv
-        = exp(-half_nu_s2_overx_dbl)
-          * pow(half_nu_s2_overx_dbl, half_nu_dbl - 1) / tgamma(half_nu_dbl);
+        = exp(-half_nu_s2_over_y) * pow(half_nu_s2_over_y, half_nu_dbl - 1)
+          / tgamma_vec[n];
 
     P += log(Pn);
 
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_[n]
-          += half_nu_s2_overx_dbl * y_inv_dbl * gamma_p_deriv / Pn;
+          += half_nu_s2_over_y * inv_y[n] * gamma_p_deriv / Pn;
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
-          += (0.5
-                  * grad_reg_inc_gamma(half_nu_dbl, half_nu_s2_overx_dbl,
-                                       gamma_vec[n], digamma_vec[n])
-              - half_s2_overx_dbl * gamma_p_deriv)
+          += 0.5
+             * (grad_reg_inc_gamma(half_nu_dbl, half_nu_s2_over_y,
+                                   tgamma_vec[n], digamma_vec[n])
+                - s2_over_y * gamma_p_deriv)
              / Pn;
     }
     if (!is_constant_all<T_scale>::value) {
       ops_partials.edge3_.partials_[n]
-          += -2.0 * half_nu_dbl * s_dbl * y_inv_dbl * gamma_p_deriv / Pn;
+          += -2.0 * half_nu_dbl * s_dbl * inv_y[n] * gamma_p_deriv / Pn;
     }
   }
   return ops_partials.build(P);
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp b/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
index 430da2df82e..a4ea41e6110 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_lpdf.hpp
@@ -6,11 +6,13 @@
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/digamma.hpp>
 #include <stan/math/prim/fun/grad_reg_inc_gamma.hpp>
+#include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/lgamma.hpp>
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/operands_and_partials.hpp>
 #include <cmath>
@@ -32,12 +34,12 @@ namespace math {
  *
  * @tparam T_y type of scalar
  * @tparam T_dof type of degrees of freedom
+ * @tparam T_scale type of scale parameter
  * @param y A scalar variable.
  * @param nu Degrees of freedom.
  * @param s Scale parameter.
- * @throw std::domain_error if nu is not greater than 0
- * @throw std::domain_error if s is not greater than 0.
- * @throw std::domain_error if y is not greater than 0.
+ * @throw std::domain_error if y is nan.
+ * @throw std::domain_error if nu or s is not greater than 0.
  */
 template <bool propto, typename T_y, typename T_dof, typename T_scale>
 return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
@@ -65,6 +67,10 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
   scalar_seq_view<T_y> y_vec(y);
   scalar_seq_view<T_dof> nu_vec(nu);
   scalar_seq_view<T_scale> s_vec(s);
+  size_t size_y = stan::math::size(y);
+  size_t size_nu = stan::math::size(nu);
+  size_t size_s = stan::math::size(s);
+  size_t size_y_s = max_size(y, s);
   size_t N = max_size(y, nu, s);
 
   for (size_t n = 0; n < N; n++) {
@@ -73,57 +79,50 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
     }
   }
 
-  VectorBuilder<include_summand<propto, T_dof, T_y, T_scale>::value,
-                T_partials_return, T_dof>
-      half_nu(size(nu));
-  for (size_t i = 0; i < stan::math::size(nu); i++) {
-    if (include_summand<propto, T_dof, T_y, T_scale>::value) {
-      half_nu[i] = 0.5 * value_of(nu_vec[i]);
-    }
+  VectorBuilder<true, T_partials_return, T_dof> half_nu(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
+    half_nu[i] = 0.5 * value_of(nu_vec[i]);
   }
 
-  VectorBuilder<include_summand<propto, T_dof, T_y>::value, T_partials_return,
+  VectorBuilder<true, T_partials_return, T_y> inv_y(size_y);
+  VectorBuilder<include_summand<propto, T_y, T_dof>::value, T_partials_return,
                 T_y>
-      log_y(size(y));
-  for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (include_summand<propto, T_dof, T_y>::value) {
-      log_y[i] = log(value_of(y_vec[i]));
-    }
-  }
-
-  VectorBuilder<include_summand<propto, T_dof, T_y, T_scale>::value,
-                T_partials_return, T_y>
-      inv_y(size(y));
-  for (size_t i = 0; i < stan::math::size(y); i++) {
-    if (include_summand<propto, T_dof, T_y, T_scale>::value) {
-      inv_y[i] = 1.0 / value_of(y_vec[i]);
+      log_y(size_y);
+  for (size_t i = 0; i < size_y; i++) {
+    T_partials_return y_dbl = value_of(y_vec[i]);
+    inv_y[i] = inv(y_dbl);
+    if (include_summand<propto, T_y, T_dof>::value) {
+      log_y[i] = log(y_dbl);
     }
   }
 
   VectorBuilder<include_summand<propto, T_dof, T_scale>::value,
                 T_partials_return, T_scale>
-      log_s(size(s));
-  for (size_t i = 0; i < stan::math::size(s); i++) {
-    if (include_summand<propto, T_dof, T_scale>::value) {
+      log_s(size_s);
+  if (include_summand<propto, T_dof, T_scale>::value) {
+    for (size_t i = 0; i < size_s; i++) {
       log_s[i] = log(value_of(s_vec[i]));
     }
   }
 
+  VectorBuilder<true, T_partials_return, T_y, T_scale> s2_over_y(size_y_s);
+  for (size_t i = 0; i < size_y_s; i++) {
+    s2_over_y[i] = square(value_of(s_vec[i])) * inv_y[i];
+  }
+
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      log_half_nu(size(nu));
+      log_half_nu(size_nu);
   VectorBuilder<include_summand<propto, T_dof>::value, T_partials_return, T_dof>
-      lgamma_half_nu(size(nu));
+      lgamma_half_nu(size_nu);
   VectorBuilder<!is_constant_all<T_dof>::value, T_partials_return, T_dof>
-      digamma_half_nu_over_two(size(nu));
-  for (size_t i = 0; i < stan::math::size(nu); i++) {
-    if (include_summand<propto, T_dof>::value) {
-      lgamma_half_nu[i] = lgamma(half_nu[i]);
-    }
+      digamma_half_nu(size_nu);
+  for (size_t i = 0; i < size_nu; i++) {
     if (include_summand<propto, T_dof>::value) {
       log_half_nu[i] = log(half_nu[i]);
+      lgamma_half_nu[i] = lgamma(half_nu[i]);
     }
     if (!is_constant_all<T_dof>::value) {
-      digamma_half_nu_over_two[i] = digamma(half_nu[i]) * 0.5;
+      digamma_half_nu[i] = digamma(half_nu[i]);
     }
   }
 
@@ -140,22 +139,23 @@ return_type_t<T_y, T_dof, T_scale> scaled_inv_chi_square_lpdf(
       logp -= (half_nu[n] + 1.0) * log_y[n];
     }
     if (include_summand<propto, T_dof, T_y, T_scale>::value) {
-      logp -= half_nu[n] * s_dbl * s_dbl * inv_y[n];
+      logp -= half_nu[n] * s2_over_y[n];
     }
 
     if (!is_constant_all<T_y>::value) {
       ops_partials.edge1_.partials_[n]
-          += -(half_nu[n] + 1.0) * inv_y[n]
-             + half_nu[n] * s_dbl * s_dbl * inv_y[n] * inv_y[n];
+          -= inv_y[n] * (half_nu[n] + 1.0 - half_nu[n] * s2_over_y[n]);
     }
     if (!is_constant_all<T_dof>::value) {
       ops_partials.edge2_.partials_[n]
-          += 0.5 * log_half_nu[n] + 0.5 - digamma_half_nu_over_two[n] + log_s[n]
-             - 0.5 * log_y[n] - 0.5 * s_dbl * s_dbl * inv_y[n];
+          += log_s[n]
+             + 0.5
+                   * (log_half_nu[n] + 1 - digamma_half_nu[n] - log_y[n]
+                      - s2_over_y[n]);
     }
     if (!is_constant_all<T_scale>::value) {
       ops_partials.edge3_.partials_[n]
-          += nu_dbl / s_dbl - nu_dbl * inv_y[n] * s_dbl;
+          += nu_dbl * (inv(s_dbl) - inv_y[n] * s_dbl);
     }
   }
   return ops_partials.build(logp);
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp b/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp
index 444e26e24a7..776cf759c5c 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp
@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/square.hpp>
 #include <boost/random/chi_squared_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 
@@ -48,7 +49,7 @@ scaled_inv_chi_square_rng(const T_deg& nu, const T_scale& s, RNG& rng) {
   for (size_t n = 0; n < N; ++n) {
     variate_generator<RNG&, chi_squared_distribution<> > chi_square_rng(
         rng, chi_squared_distribution<>(nu_vec[n]));
-    output[n] = nu_vec[n] * s_vec[n] * s_vec[n] / chi_square_rng();
+    output[n] = nu_vec[n] * square(s_vec[n]) / chi_square_rng();
   }
 
   return output.data();