Fix RNG segfaults when given size-0 input

stan/math/prim/prob/beta_binomial_rng.hpp

@@ -4,6 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/prob/binomial_rng.hpp>
 #include <stan/math/prim/prob/beta_rng.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 
 namespace stan {
 namespace math {
@@ -39,6 +40,10 @@ beta_binomial_rng(const T_N &N, const T_shape1 &alpha, const T_shape2 &beta,
   check_consistent_sizes(function, "First prior sample size parameter", alpha,
                          "Second prior sample size parameter", beta);
 
+  if (size_zero(N, alpha, beta)) {
+    return {};
+  }
+
   T_N_ref N_ref = N;
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;

stan/math/prim/prob/beta_proportion_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/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -39,6 +40,11 @@ beta_proportion_rng(const T_loc &mu, const T_prec &kappa, RNG &rng) {
   static constexpr const char *function = "beta_proportion_rng";
   check_consistent_sizes(function, "Location parameter", mu,
                          "Precision parameter", kappa);
+
+  if (size_zero(mu, kappa)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_kappa_ref kappa_ref = kappa;
   check_positive(function, "Location parameter", mu_ref);

stan/math/prim/prob/beta_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
@@ -43,6 +44,10 @@ inline typename VectorBuilder<true, double, T_shape1, T_shape2>::type beta_rng(
   static constexpr const char *function = "beta_rng";
   check_consistent_sizes(function, "First shape parameter", alpha,
                          "Second shape Parameter", beta);
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
+
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
   check_positive_finite(function, "First shape parameter", alpha_ref);

stan/math/prim/prob/binomial_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/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/binomial_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -38,6 +39,10 @@ inline typename VectorBuilder<true, int, T_N, T_theta>::type binomial_rng(
   static constexpr const char* function = "binomial_rng";
   check_consistent_sizes(function, "Population size parameter", N,
                          "Probability Parameter", theta);
+  if (size_zero(N, theta)) {
+    return {};
+  }
+
   T_N_ref N_ref = N;
   T_theta_ref theta_ref = theta;
   check_nonnegative(function, "Population size parameter", N_ref);

stan/math/prim/prob/cauchy_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/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/cauchy_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -39,6 +40,10 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type cauchy_rng(
   using T_sigma_ref = ref_type_t<T_scale>;
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);

stan/math/prim/prob/discrete_range_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/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/uniform_int_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -40,6 +41,9 @@ discrete_range_rng(const T_lower& lower, const T_upper& upper, RNG& rng) {
   check_consistent_sizes(function, "Lower bound parameter", lower,
                          "Upper bound parameter", upper);
   check_greater_or_equal(function, "Upper bound parameter", upper, lower);
+  if (size_zero(upper, lower)) {
+    return {};
+  }
 
   scalar_seq_view<T_lower> lower_vec(lower);
   scalar_seq_view<T_upper> upper_vec(upper);

stan/math/prim/prob/double_exponential_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/log1m.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -41,6 +42,10 @@ double_exponential_rng(const T_loc& mu, const T_scale& sigma, RNG& rng) {
   static constexpr const char* function = "double_exponential_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);

stan/math/prim/prob/exp_mod_normal_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/prob/exponential_rng.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/prob/normal_rng.hpp>
 #include <boost/random/normal_distribution.hpp>
@@ -45,7 +46,9 @@ exp_mod_normal_rng(const T_loc& mu, const T_scale& sigma,
   using T_lambda_ref = ref_type_t<T_inv_scale>;
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma, "Inv_scale Parameter", lambda);
-
+  // if (size_zero(mu, sigma, lambda)) {
+  //   return {};
+  // }
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   T_lambda_ref lambda_ref = lambda;

stan/math/prim/prob/frechet_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/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/weibull_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -38,6 +39,10 @@ inline typename VectorBuilder<true, double, T_shape, T_scale>::type frechet_rng(
   static constexpr const char* function = "frechet_rng";
   check_consistent_sizes(function, "Shape parameter", alpha, "Scale Parameter",
                          sigma);
+  if (size_zero(alpha, sigma)) {
+    return {};
+  }
+
   T_alpha_ref alpha_ref = alpha;
   T_sigma_ref sigma_ref = sigma;
   check_positive_finite(function, "Shape parameter", alpha_ref);

stan/math/prim/prob/gamma_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 
@@ -39,6 +40,10 @@ inline typename VectorBuilder<true, double, T_shape, T_inv>::type gamma_rng(
   static constexpr const char* function = "gamma_rng";
   check_consistent_sizes(function, "Shape parameter", alpha,
                          "Inverse scale Parameter", beta);
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
+
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
   check_positive_finite(function, "Shape parameter", alpha_ref);

stan/math/prim/prob/gumbel_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/uniform_01.hpp>
 #include <boost/random/variate_generator.hpp>
 #include <cmath>
@@ -40,6 +41,9 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type gumbel_rng(
   static constexpr const char* function = "gumbel_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          beta);
+  if (size_zero(mu, beta)) {
+    return {};
+  }
   T_mu_ref mu_ref = mu;
   T_beta_ref beta_ref = beta;
   check_finite(function, "Location parameter", mu_ref);

stan/math/prim/prob/inv_gamma_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 
@@ -39,6 +40,10 @@ inv_gamma_rng(const T_shape& alpha, const T_scale& beta, RNG& rng) {
   static constexpr const char* function = "inv_gamma_rng";
   check_consistent_sizes(function, "Shape parameter", alpha, "Scale Parameter",
                          beta);
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
+
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
   check_positive_finite(function, "Shape parameter", alpha_ref);

stan/math/prim/prob/logistic_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/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/exponential_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -40,6 +41,10 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type logistic_rng(
   static constexpr const char* function = "logistic_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);

stan/math/prim/prob/loglogistic_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/uniform_01.hpp>
 #include <boost/random/variate_generator.hpp>
 
@@ -40,6 +41,9 @@ loglogistic_rng(const T_scale& alpha, const T_shape& beta, RNG& rng) {
   static constexpr const char* function = "loglogistic_rng";
   check_consistent_sizes(function, "Scale parameter", alpha, "Shape Parameter",
                          beta);
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
   check_positive_finite(function, "Scale parameter", alpha_ref);

stan/math/prim/prob/lognormal_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/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/lognormal_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
@@ -39,6 +40,9 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type lognormal_rng(
   static constexpr const char* function = "lognormal_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);

stan/math/prim/prob/neg_binomial_2_log_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/poisson_distribution.hpp>
@@ -43,6 +44,9 @@ neg_binomial_2_log_rng(const T_loc& eta, const T_inv& phi, RNG& rng) {
   static constexpr const char* function = "neg_binomial_2_log_rng";
   check_consistent_sizes(function, "Log-location parameter", eta,
                          "Inverse dispersion parameter", phi);
+  if (size_zero(eta, phi)) {
+    return {};
+  }
   T_eta_ref eta_ref = eta;
   T_phi_ref phi_ref = phi;
   check_finite(function, "Log-location parameter", eta_ref);

stan/math/prim/prob/neg_binomial_2_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/poisson_distribution.hpp>
@@ -42,6 +43,9 @@ neg_binomial_2_rng(const T_loc& mu, const T_prec& phi, RNG& rng) {
   static constexpr const char* function = "neg_binomial_2_rng";
   check_consistent_sizes(function, "Location parameter", mu,
                          "Precision parameter", phi);
+  if (size_zero(mu, phi)) {
+    return {};
+  }
   T_mu_ref mu_ref = mu;
   T_phi_ref phi_ref = phi;
   check_positive_finite(function, "Location parameter", mu_ref);

stan/math/prim/prob/neg_binomial_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <boost/random/gamma_distribution.hpp>
 #include <boost/random/poisson_distribution.hpp>
@@ -44,6 +45,11 @@ inline typename VectorBuilder<true, int, T_shape, T_inv>::type neg_binomial_rng(
                          "Inverse scale Parameter", beta);
   T_alpha_ref alpha_ref = alpha;
   T_beta_ref beta_ref = beta;
+
+  if (size_zero(alpha, beta)) {
+    return {};
+  }
+
   check_positive_finite(function, "Shape parameter", alpha_ref);
   check_positive_finite(function, "Inverse scale parameter", beta_ref);
 

stan/math/prim/prob/normal_rng.hpp

@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
+#include <stan/math/prim/fun/size_zero.hpp>
 #include <boost/random/normal_distribution.hpp>
 #include <boost/random/variate_generator.hpp>
 
@@ -39,6 +40,10 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type normal_rng(
   static constexpr const char* function = "normal_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma);
+  if (size_zero(mu, sigma)) {
+    return {};
+  }
+
   T_mu_ref mu_ref = mu;
   T_sigma_ref sigma_ref = sigma;
   check_finite(function, "Location parameter", mu_ref);

stan/math/prim/prob/pareto_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/size_zero.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <boost/random/exponential_distribution.hpp>
@@ -39,6 +40,9 @@ inline typename VectorBuilder<true, double, T_shape, T_scale>::type pareto_rng(
   static constexpr const char* function = "pareto_rng";
   check_consistent_sizes(function, "Scale Parameter", y_min, "Shape parameter",
                          alpha);
+  if (size_zero(y_min, alpha)) {
+    return {};
+  }
   const auto& y_min_ref = to_ref(y_min);
   const auto& alpha_ref = to_ref(alpha);
   check_positive_finite(function, "Scale parameter", y_min_ref);

stan/math/prim/prob/pareto_type_2_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/size_zero.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/prob/exponential_rng.hpp>
@@ -46,6 +47,9 @@ pareto_type_2_rng(const T_loc& mu, const T_scale& lambda, const T_shape& alpha,
   static constexpr const char* function = "pareto_type_2_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          lambda, "Shape Parameter", alpha);
+  if (size_zero(mu, lambda, alpha)) {
+    return {};
+  }
   const auto& mu_ref = to_ref(mu);
   const auto& lambda_ref = to_ref(lambda);
   const auto& alpha_ref = to_ref(alpha);