From 2867e17ec6ff56b928a59b20582b37491e7cf92f Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Sun, 9 Jan 2022 21:13:05 +0800
Subject: [PATCH 1/7] Initial 1

---
 stan/math/fwd.hpp                             |  1 +
 stan/math/fwd/prob.hpp                        |  6 +++
 stan/math/fwd/prob/exponential_qf.hpp         | 42 +++++++++++++++++++
 stan/math/mix.hpp                             |  2 +
 stan/math/prim/prob.hpp                       |  1 +
 stan/math/prim/prob/exponential_qf.hpp        | 29 +++++++++++++
 stan/math/rev.hpp                             |  1 +
 stan/math/rev/prob.hpp                        |  6 +++
 stan/math/rev/prob/exponential_qf.hpp         | 40 ++++++++++++++++++
 .../math/mix/prob/exponential_qf_test.cpp     | 17 ++++++++
 .../math/prim/prob/exponential_qf_test.cpp    | 28 +++++++++++++
 11 files changed, 173 insertions(+)
 create mode 100644 stan/math/fwd/prob.hpp
 create mode 100644 stan/math/fwd/prob/exponential_qf.hpp
 create mode 100644 stan/math/prim/prob/exponential_qf.hpp
 create mode 100644 stan/math/rev/prob.hpp
 create mode 100644 stan/math/rev/prob/exponential_qf.hpp
 create mode 100644 test/unit/math/mix/prob/exponential_qf_test.cpp
 create mode 100644 test/unit/math/prim/prob/exponential_qf_test.cpp

diff --git a/stan/math/fwd.hpp b/stan/math/fwd.hpp
index 7f30955b0f2..3bab5be267f 100644
--- a/stan/math/fwd.hpp
+++ b/stan/math/fwd.hpp
@@ -11,6 +11,7 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun.hpp>
 #include <stan/math/fwd/functor.hpp>
+#include <stan/math/fwd/prob.hpp>
 
 #include <stan/math/prim.hpp>
 
diff --git a/stan/math/fwd/prob.hpp b/stan/math/fwd/prob.hpp
new file mode 100644
index 00000000000..725acd8bb23
--- /dev/null
+++ b/stan/math/fwd/prob.hpp
@@ -0,0 +1,6 @@
+#ifndef STAN_MATH_FWD_PROB_HPP
+#define STAN_MATH_FWD_PROB_HPP
+
+#include <stan/math/fwd/prob/exponential_qf.hpp>
+
+#endif
diff --git a/stan/math/fwd/prob/exponential_qf.hpp b/stan/math/fwd/prob/exponential_qf.hpp
new file mode 100644
index 00000000000..60667add484
--- /dev/null
+++ b/stan/math/fwd/prob/exponential_qf.hpp
@@ -0,0 +1,42 @@
+#ifndef STAN_MATH_FWD_PROB_EXPONENTIAL_QF_HPP
+#define STAN_MATH_FWD_PROB_EXPONENTIAL_QF_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/prob/exponential_qf.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename Tp, typename Tbeta,
+          require_any_st_fvar<Tp, Tbeta>* = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  using vector_t = decltype(exponential_qf(value_of(p),
+                                           value_of(beta)));
+  using fvar_t = return_type_t<Tp, Tbeta>;
+
+  auto p_val = value_of(p);
+  auto beta_val = value_of(beta);
+  promote_scalar_t<fvar_t, vector_t>
+    ret(exponential_qf(p_val, beta_val));
+  ret.d();// = vector_t(0.0);
+  if(is_fvar<Tp>::value) {
+    as_array_or_scalar(ret.d()) +=
+      as_array_or_scalar(forward_as<promote_scalar_t<fvar_t, Tp>>(p).d())
+      * inv(as_array_or_scalar(beta_val)
+            - as_array_or_scalar(beta_val)
+            * as_array_or_scalar(p_val));
+  }
+  if(is_fvar<Tbeta>::value) {
+    as_array_or_scalar(ret.d()) +=
+      as_array_or_scalar(forward_as<promote_scalar_t<fvar_t, Tbeta>>(beta).d())
+      * log1m(as_array_or_scalar(p_val)) * as_array_or_scalar(beta_val);
+  }
+
+  return ret;
+}
+
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/mix.hpp b/stan/math/mix.hpp
index 38b6a5e9c0e..63651abd76c 100644
--- a/stan/math/mix.hpp
+++ b/stan/math/mix.hpp
@@ -13,11 +13,13 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun.hpp>
 #include <stan/math/fwd/functor.hpp>
+#include <stan/math/fwd/prob.hpp>
 
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/fun.hpp>
 #include <stan/math/rev/functor.hpp>
+#include <stan/math/rev/prob.hpp>
 
 #include <stan/math/prim.hpp>
 
diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp
index 33584411ec3..77cf4ed8ee0 100644
--- a/stan/math/prim/prob.hpp
+++ b/stan/math/prim/prob.hpp
@@ -106,6 +106,7 @@
 #include <stan/math/prim/prob/exponential_lcdf.hpp>
 #include <stan/math/prim/prob/exponential_log.hpp>
 #include <stan/math/prim/prob/exponential_lpdf.hpp>
+#include <stan/math/prim/prob/exponential_qf.hpp>
 #include <stan/math/prim/prob/exponential_rng.hpp>
 #include <stan/math/prim/prob/frechet_ccdf_log.hpp>
 #include <stan/math/prim/prob/frechet_cdf.hpp>
diff --git a/stan/math/prim/prob/exponential_qf.hpp b/stan/math/prim/prob/exponential_qf.hpp
new file mode 100644
index 00000000000..e2e39124bf1
--- /dev/null
+++ b/stan/math/prim/prob/exponential_qf.hpp
@@ -0,0 +1,29 @@
+#ifndef STAN_MATH_PRIM_PROB_EXPONENTIAL_QF_HPP
+#define STAN_MATH_PRIM_PROB_EXPONENTIAL_QF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+template <typename Tp, typename Tbeta,
+          require_all_st_arithmetic<Tp, Tbeta>* = nullptr,
+          require_any_eigen_vector_t<Tp, Tbeta>* = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  return (-log1p(-as_array_or_scalar(p))
+            / as_array_or_scalar(beta)).matrix();
+}
+
+template <typename Tp, typename Tbeta,
+          require_all_arithmetic_t<Tp, Tbeta> * = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  return -log1p(-p) / beta;
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev.hpp b/stan/math/rev.hpp
index 866deb448d3..20d5a477846 100644
--- a/stan/math/rev.hpp
+++ b/stan/math/rev.hpp
@@ -9,6 +9,7 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/fun.hpp>
 #include <stan/math/rev/functor.hpp>
+#include <stan/math/rev/prob.hpp>
 
 #include <stan/math/prim.hpp>
 
diff --git a/stan/math/rev/prob.hpp b/stan/math/rev/prob.hpp
new file mode 100644
index 00000000000..48156af85ed
--- /dev/null
+++ b/stan/math/rev/prob.hpp
@@ -0,0 +1,6 @@
+#ifndef STAN_MATH_REV_PROB_HPP
+#define STAN_MATH_REV_PROB_HPP
+
+#include <stan/math/rev/prob/exponential_qf.hpp>
+
+#endif
diff --git a/stan/math/rev/prob/exponential_qf.hpp b/stan/math/rev/prob/exponential_qf.hpp
new file mode 100644
index 00000000000..136fce77b2f
--- /dev/null
+++ b/stan/math/rev/prob/exponential_qf.hpp
@@ -0,0 +1,40 @@
+#ifndef STAN_MATH_REV_PROB_EXPONENTIAL_QF_HPP
+#define STAN_MATH_REV_PROB_EXPONENTIAL_QF_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/prim/prob/exponential_qf.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename Tp, typename Tbeta,
+          require_any_st_var<Tp, Tbeta>* = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  using inner_ret_type = decltype(exponential_qf(value_of(p), value_of(beta)));
+  using ret_type = return_var_matrix_t<inner_ret_type, Tp, Tbeta>;
+
+  arena_t<Tp> arena_p = p;
+  arena_t<Tbeta> arena_beta = beta;
+  arena_t<ret_type> ret = exponential_qf(value_of(arena_p),
+                                         value_of(arena_beta));
+    reverse_pass_callback([ret, arena_p, arena_beta]() mutable {
+      decltype(auto) p_val = as_array_or_scalar(value_of(arena_p));
+      decltype(auto) beta_val = as_array_or_scalar(value_of(arena_beta));
+      decltype(auto) ret_adj = as_array_or_scalar(ret.adj());
+      if (!is_constant<Tp>::value) {
+        as_array_or_scalar(forward_as<promote_scalar_t<var, Tp>>(arena_p).adj()) += ret_adj
+          * inv(beta_val - beta_val * p_val);
+      }
+      if (!is_constant<Tbeta>::value) {
+        as_array_or_scalar(forward_as<promote_scalar_t<var, Tbeta>>(arena_beta).adj()) += ret_adj
+          * log1m(p_val) / square(beta_val);
+      }
+    });
+    return ret_type(ret);
+}
+
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/test/unit/math/mix/prob/exponential_qf_test.cpp b/test/unit/math/mix/prob/exponential_qf_test.cpp
new file mode 100644
index 00000000000..8fe52f4749e
--- /dev/null
+++ b/test/unit/math/mix/prob/exponential_qf_test.cpp
@@ -0,0 +1,17 @@
+#include <stan/math/mix.hpp>
+#include <test/unit/math/test_ad.hpp>
+
+TEST(mathMixProb, exponential_qf) {
+  auto f = [](const auto& p, const auto& beta) {
+    return stan::math::exponential_qf(p, beta);
+  };
+
+  Eigen::VectorXd p(2);
+  p << 0.2, 0.6;
+
+  Eigen::VectorXd beta(2);
+  beta << 16, 2;
+
+  stan::test::expect_ad(f, p, beta);
+  stan::test::expect_ad(f, 0.1, 8);
+}
diff --git a/test/unit/math/prim/prob/exponential_qf_test.cpp b/test/unit/math/prim/prob/exponential_qf_test.cpp
new file mode 100644
index 00000000000..42075295159
--- /dev/null
+++ b/test/unit/math/prim/prob/exponential_qf_test.cpp
@@ -0,0 +1,28 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <test/unit/util.hpp>
+
+TEST(ProbExponential, quantile_values) {
+  using stan::math::exponential_qf;
+
+  EXPECT_FLOAT_EQ(exponential_qf(0.1, 8), 0.0131700644572);
+  EXPECT_FLOAT_EQ(exponential_qf(0.4, 7), 0.0729750891094);
+  EXPECT_FLOAT_EQ(exponential_qf(0.1, 7), 0.0150515022368);
+
+  Eigen::VectorXd p(2);
+  p << 0.2, 0.6;
+
+  Eigen::VectorXd beta(2);
+  beta << 16, 2;
+
+  Eigen::VectorXd res(2);
+  res << 0.0318776501877, 0.1308986759820;
+
+  EXPECT_MATRIX_FLOAT_EQ(exponential_qf(p, 7), res);
+
+  res << 0.0222921839962, 0.1783374719694;
+  EXPECT_MATRIX_FLOAT_EQ(exponential_qf(0.3, beta), res);
+
+  res << 0.0139464719571, 0.4581453659371;
+  EXPECT_MATRIX_FLOAT_EQ(exponential_qf(p, beta), res);
+}

From 074ee72501d314b9aa0e4bd163d4d48d137d8a06 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Mon, 10 Jan 2022 11:39:50 +0800
Subject: [PATCH 2/7] Fix fwd implementation

---
 stan/math/fwd/prob/exponential_qf.hpp | 58 +++++++++++++++++++--------
 1 file changed, 42 insertions(+), 16 deletions(-)

diff --git a/stan/math/fwd/prob/exponential_qf.hpp b/stan/math/fwd/prob/exponential_qf.hpp
index 60667add484..9a354fe58b0 100644
--- a/stan/math/fwd/prob/exponential_qf.hpp
+++ b/stan/math/fwd/prob/exponential_qf.hpp
@@ -9,33 +9,59 @@ namespace stan {
 namespace math {
 
 template <typename Tp, typename Tbeta,
-          require_any_st_fvar<Tp, Tbeta>* = nullptr>
+          require_any_st_fvar<Tp, Tbeta>* = nullptr,
+          require_all_stan_scalar_t<Tp, Tbeta>* = nullptr>
 inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
-  using vector_t = decltype(exponential_qf(value_of(p),
-                                           value_of(beta)));
   using fvar_t = return_type_t<Tp, Tbeta>;
 
   auto p_val = value_of(p);
   auto beta_val = value_of(beta);
-  promote_scalar_t<fvar_t, vector_t>
-    ret(exponential_qf(p_val, beta_val));
-  ret.d();// = vector_t(0.0);
-  if(is_fvar<Tp>::value) {
-    as_array_or_scalar(ret.d()) +=
-      as_array_or_scalar(forward_as<promote_scalar_t<fvar_t, Tp>>(p).d())
-      * inv(as_array_or_scalar(beta_val)
-            - as_array_or_scalar(beta_val)
-            * as_array_or_scalar(p_val));
+  fvar_t ret(exponential_qf(p_val, beta_val));
+
+  if(!is_constant<Tp>::value) {
+    ret.d_ += forward_as<fvar_t>(p).d_
+      * inv(beta_val - beta_val * p_val);
   }
-  if(is_fvar<Tbeta>::value) {
-    as_array_or_scalar(ret.d()) +=
-      as_array_or_scalar(forward_as<promote_scalar_t<fvar_t, Tbeta>>(beta).d())
-      * log1m(as_array_or_scalar(p_val)) * as_array_or_scalar(beta_val);
+  if(!is_constant<Tbeta>::value) {
+    ret.d_ += forward_as<fvar_t>(beta).d_
+      * log1m(p_val) / square(beta_val);
   }
 
   return ret;
 }
 
+template <typename Tp, typename Tbeta,
+          require_any_st_fvar<Tp, Tbeta>* = nullptr,
+          require_any_eigen_vector_t<Tp, Tbeta>* = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  using vector_t
+    = plain_type_t<decltype(exponential_qf(value_of(p), value_of(beta)))>;
+  using fvar_t = return_type_t<Tp, Tbeta>;
+
+  auto p_val = value_of(p);
+  auto beta_val = value_of(beta);
+  vector_t ret = exponential_qf(p_val, beta_val);
+
+  auto p_array = as_array_or_scalar(p_val);
+  auto beta_array = as_array_or_scalar(beta_val);
+  vector_t d_ = vector_t::Zero(ret.rows(), ret.cols());
+
+  if(!is_constant<Tp>::value) {
+    as_array_or_scalar(d_)
+      += as_array_or_scalar(forward_as<promote_scalar_t<fvar_t, Tp>>(p).d())
+      * inv(beta_array - beta_array * p_array);
+  }
+  if(!is_constant<Tbeta>::value) {
+    as_array_or_scalar(d_)
+      += as_array_or_scalar(forward_as<promote_scalar_t<fvar_t,
+                                                        Tbeta>>(beta).d())
+      * log1m(p_array) / square(beta_array);
+  }
+
+  return ret.binaryExpr(d_, [&](const auto& val, const auto& deriv) {
+    return fvar_t(val, deriv);
+  }).eval();
+}
 
 }  // namespace math
 }  // namespace stan

From fbb16e638eac107e6dbe453bd536bb44ddc50d41 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Mon, 10 Jan 2022 11:56:44 +0800
Subject: [PATCH 3/7] Update doc

---
 stan/math/prim/prob/exponential_qf.hpp        | 52 ++++++++++++++++---
 stan/math/rev/prob/exponential_qf.hpp         | 28 ++++++++--
 .../math/prim/prob/exponential_qf_test.cpp    | 25 +++++++++
 3 files changed, 94 insertions(+), 11 deletions(-)

diff --git a/stan/math/prim/prob/exponential_qf.hpp b/stan/math/prim/prob/exponential_qf.hpp
index e2e39124bf1..8abe68b63c5 100644
--- a/stan/math/prim/prob/exponential_qf.hpp
+++ b/stan/math/prim/prob/exponential_qf.hpp
@@ -9,18 +9,58 @@
 namespace stan {
 namespace math {
 
+/** \ingroup prob_dists
+ * The quantile of an exponential density for p with the specified
+ * inverse scale parameter.
+ * Inverse scale parameter must be greater than 0.
+ * p must be bounded by 0 and 1.
+ *
+ \f[
+   q = -\log(1 - p) / \beta
+ \f]
+ *
+ * @tparam Tp type of probability input
+ * @tparam Tbeta type of inverse scale
+ * @param p A scalar variable.
+ * @param beta Inverse scale parameter.
+ * @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.
+ */
 template <typename Tp, typename Tbeta,
-          require_all_st_arithmetic<Tp, Tbeta>* = nullptr,
-          require_any_eigen_vector_t<Tp, Tbeta>* = nullptr>
+          require_all_arithmetic_t<Tp, Tbeta> * = nullptr>
 inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
-  return (-log1p(-as_array_or_scalar(p))
-            / as_array_or_scalar(beta)).matrix();
+  static const char* function = "exponential_qf";
+  check_positive_finite(function, "Inverse scale parameter", beta);
+  check_bounded(function, "Probability parameter", p, 0, 1);
+  return -log1p(-p) / beta;
 }
 
+/** \ingroup prob_dists
+ * The quantile of an exponential density for p with the specified
+ * inverse scale parameter.
+ * Inverse scale parameter must be greater than 0.
+ * p must be bounded by 0 and 1.
+ * 
+ * Specialisation for use where any input is an Eigen vector
+ *
+ * @tparam Tp type of probability input
+ * @tparam Tbeta type of inverse scale
+ * @param p A vector or scalar of probabilities.
+ * @param beta Vector or scalar of inverse scale parameter.
+ * @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.
+ */
 template <typename Tp, typename Tbeta,
-          require_all_arithmetic_t<Tp, Tbeta> * = nullptr>
+          require_all_st_arithmetic<Tp, Tbeta>* = nullptr,
+          require_any_eigen_vector_t<Tp, Tbeta>* = nullptr>
 inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
-  return -log1p(-p) / beta;
+  static const char* function = "exponential_qf";
+  ref_type_t<Tp> p_ref = p;
+  ref_type_t<Tbeta> beta_ref = beta;
+  check_positive_finite(function, "Inverse scale parameter", beta_ref);
+  check_bounded(function, "Probability parameter", p_ref, 0, 1);
+  return (-log1p(-as_array_or_scalar(p_ref))
+            / as_array_or_scalar(beta_ref)).matrix();
 }
 
 }  // namespace math
diff --git a/stan/math/rev/prob/exponential_qf.hpp b/stan/math/rev/prob/exponential_qf.hpp
index 136fce77b2f..2f2d1fecbc6 100644
--- a/stan/math/rev/prob/exponential_qf.hpp
+++ b/stan/math/rev/prob/exponential_qf.hpp
@@ -8,6 +8,24 @@
 namespace stan {
 namespace math {
 
+/** \ingroup prob_dists
+ * The quantile of an exponential density for p with the specified
+ * inverse scale parameter.
+ * Inverse scale parameter must be greater than 0.
+ * p must be bounded by 0 and 1.
+ *
+ \f[
+   \frac{\partial }{\partial p} = \frac{1}{\beta - \beta p} \\
+   \frac{\partial }{\partial \beta} = \frac{\log(1 - p)}{\beta^2}
+ \f]
+ *
+ * @tparam Tp type of probability input
+ * @tparam Tbeta type of inverse scale
+ * @param p A scalar variable.
+ * @param beta Inverse scale parameter.
+ * @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.
+ */
 template <typename Tp, typename Tbeta,
           require_any_st_var<Tp, Tbeta>* = nullptr>
 inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
@@ -23,18 +41,18 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
       decltype(auto) beta_val = as_array_or_scalar(value_of(arena_beta));
       decltype(auto) ret_adj = as_array_or_scalar(ret.adj());
       if (!is_constant<Tp>::value) {
-        as_array_or_scalar(forward_as<promote_scalar_t<var, Tp>>(arena_p).adj()) += ret_adj
-          * inv(beta_val - beta_val * p_val);
+        as_array_or_scalar(forward_as<promote_scalar_t<var, Tp>>(arena_p).adj())
+          += ret_adj * inv(beta_val - beta_val * p_val);
       }
       if (!is_constant<Tbeta>::value) {
-        as_array_or_scalar(forward_as<promote_scalar_t<var, Tbeta>>(arena_beta).adj()) += ret_adj
-          * log1m(p_val) / square(beta_val);
+        as_array_or_scalar(
+          forward_as<promote_scalar_t<var, Tbeta>>(arena_beta).adj())
+          += ret_adj * log1m(p_val) / square(beta_val);
       }
     });
     return ret_type(ret);
 }
 
-
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/test/unit/math/prim/prob/exponential_qf_test.cpp b/test/unit/math/prim/prob/exponential_qf_test.cpp
index 42075295159..dc7d1a84bd1 100644
--- a/test/unit/math/prim/prob/exponential_qf_test.cpp
+++ b/test/unit/math/prim/prob/exponential_qf_test.cpp
@@ -26,3 +26,28 @@ TEST(ProbExponential, quantile_values) {
   res << 0.0139464719571, 0.4581453659371;
   EXPECT_MATRIX_FLOAT_EQ(exponential_qf(p, beta), res);
 }
+
+TEST(ProbExponential, quantile_throws) {
+  using stan::math::exponential_qf;
+
+  EXPECT_THROW(exponential_qf(2, 8), std::domain_error);
+  EXPECT_THROW(exponential_qf(0.4, -7), std::domain_error);
+
+  Eigen::VectorXd p(2);
+  p << 0.2, 0.6;
+
+  Eigen::VectorXd p_invalid(2);
+  p_invalid << 0.2, 2.2;
+
+  Eigen::VectorXd beta(2);
+  beta << 16, 2;
+
+  Eigen::VectorXd beta_invalid(2);
+  beta_invalid << 16, -8;
+
+  EXPECT_THROW(exponential_qf(p_invalid, 8), std::domain_error);
+  EXPECT_THROW(exponential_qf(p_invalid, beta), std::domain_error);
+  EXPECT_THROW(exponential_qf(0.4, beta_invalid), std::domain_error);
+  EXPECT_THROW(exponential_qf(p, beta_invalid), std::domain_error);
+  EXPECT_THROW(exponential_qf(p_invalid, beta_invalid), std::domain_error);
+}
\ No newline at end of file

From 580fc48bddbe6bf669ea2749591313810449fbe9 Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Mon, 10 Jan 2022 12:01:09 +0800
Subject: [PATCH 4/7] fix cpplint

---
 stan/math/fwd/prob/exponential_qf.hpp            | 8 ++++----
 test/unit/math/prim/prob/exponential_qf_test.cpp | 2 +-
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/stan/math/fwd/prob/exponential_qf.hpp b/stan/math/fwd/prob/exponential_qf.hpp
index 9a354fe58b0..4f698bb21cb 100644
--- a/stan/math/fwd/prob/exponential_qf.hpp
+++ b/stan/math/fwd/prob/exponential_qf.hpp
@@ -18,11 +18,11 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   auto beta_val = value_of(beta);
   fvar_t ret(exponential_qf(p_val, beta_val));
 
-  if(!is_constant<Tp>::value) {
+  if (!is_constant<Tp>::value) {
     ret.d_ += forward_as<fvar_t>(p).d_
       * inv(beta_val - beta_val * p_val);
   }
-  if(!is_constant<Tbeta>::value) {
+  if (!is_constant<Tbeta>::value) {
     ret.d_ += forward_as<fvar_t>(beta).d_
       * log1m(p_val) / square(beta_val);
   }
@@ -46,12 +46,12 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   auto beta_array = as_array_or_scalar(beta_val);
   vector_t d_ = vector_t::Zero(ret.rows(), ret.cols());
 
-  if(!is_constant<Tp>::value) {
+  if (!is_constant<Tp>::value) {
     as_array_or_scalar(d_)
       += as_array_or_scalar(forward_as<promote_scalar_t<fvar_t, Tp>>(p).d())
       * inv(beta_array - beta_array * p_array);
   }
-  if(!is_constant<Tbeta>::value) {
+  if (!is_constant<Tbeta>::value) {
     as_array_or_scalar(d_)
       += as_array_or_scalar(forward_as<promote_scalar_t<fvar_t,
                                                         Tbeta>>(beta).d())
diff --git a/test/unit/math/prim/prob/exponential_qf_test.cpp b/test/unit/math/prim/prob/exponential_qf_test.cpp
index dc7d1a84bd1..2ceca7d76c0 100644
--- a/test/unit/math/prim/prob/exponential_qf_test.cpp
+++ b/test/unit/math/prim/prob/exponential_qf_test.cpp
@@ -50,4 +50,4 @@ TEST(ProbExponential, quantile_throws) {
   EXPECT_THROW(exponential_qf(0.4, beta_invalid), std::domain_error);
   EXPECT_THROW(exponential_qf(p, beta_invalid), std::domain_error);
   EXPECT_THROW(exponential_qf(p_invalid, beta_invalid), std::domain_error);
-}
\ No newline at end of file
+}

From c4dd5281987d013af99e07c98bf130b7175b0f33 Mon Sep 17 00:00:00 2001
From: Stan Jenkins <mc.stanislaw@gmail.com>
Date: Mon, 10 Jan 2022 04:13:49 +0000
Subject: [PATCH 5/7] [Jenkins] auto-formatting by clang-format version
 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final)

---
 stan/math/fwd/prob/exponential_qf.hpp  | 25 ++++++++++-----------
 stan/math/prim/prob/exponential_qf.hpp |  8 +++----
 stan/math/rev/prob/exponential_qf.hpp  | 31 +++++++++++++-------------
 3 files changed, 31 insertions(+), 33 deletions(-)

diff --git a/stan/math/fwd/prob/exponential_qf.hpp b/stan/math/fwd/prob/exponential_qf.hpp
index 4f698bb21cb..c7f305d27f4 100644
--- a/stan/math/fwd/prob/exponential_qf.hpp
+++ b/stan/math/fwd/prob/exponential_qf.hpp
@@ -19,12 +19,10 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   fvar_t ret(exponential_qf(p_val, beta_val));
 
   if (!is_constant<Tp>::value) {
-    ret.d_ += forward_as<fvar_t>(p).d_
-      * inv(beta_val - beta_val * p_val);
+    ret.d_ += forward_as<fvar_t>(p).d_ * inv(beta_val - beta_val * p_val);
   }
   if (!is_constant<Tbeta>::value) {
-    ret.d_ += forward_as<fvar_t>(beta).d_
-      * log1m(p_val) / square(beta_val);
+    ret.d_ += forward_as<fvar_t>(beta).d_ * log1m(p_val) / square(beta_val);
   }
 
   return ret;
@@ -35,7 +33,7 @@ template <typename Tp, typename Tbeta,
           require_any_eigen_vector_t<Tp, Tbeta>* = nullptr>
 inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   using vector_t
-    = plain_type_t<decltype(exponential_qf(value_of(p), value_of(beta)))>;
+      = plain_type_t<decltype(exponential_qf(value_of(p), value_of(beta)))>;
   using fvar_t = return_type_t<Tp, Tbeta>;
 
   auto p_val = value_of(p);
@@ -48,19 +46,20 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
 
   if (!is_constant<Tp>::value) {
     as_array_or_scalar(d_)
-      += as_array_or_scalar(forward_as<promote_scalar_t<fvar_t, Tp>>(p).d())
-      * inv(beta_array - beta_array * p_array);
+        += as_array_or_scalar(forward_as<promote_scalar_t<fvar_t, Tp>>(p).d())
+           * inv(beta_array - beta_array * p_array);
   }
   if (!is_constant<Tbeta>::value) {
     as_array_or_scalar(d_)
-      += as_array_or_scalar(forward_as<promote_scalar_t<fvar_t,
-                                                        Tbeta>>(beta).d())
-      * log1m(p_array) / square(beta_array);
+        += as_array_or_scalar(
+               forward_as<promote_scalar_t<fvar_t, Tbeta>>(beta).d())
+           * log1m(p_array) / square(beta_array);
   }
 
-  return ret.binaryExpr(d_, [&](const auto& val, const auto& deriv) {
-    return fvar_t(val, deriv);
-  }).eval();
+  return ret
+      .binaryExpr(d_, [&](const auto& val,
+                          const auto& deriv) { return fvar_t(val, deriv); })
+      .eval();
 }
 
 }  // namespace math
diff --git a/stan/math/prim/prob/exponential_qf.hpp b/stan/math/prim/prob/exponential_qf.hpp
index 8abe68b63c5..ae541f5f80f 100644
--- a/stan/math/prim/prob/exponential_qf.hpp
+++ b/stan/math/prim/prob/exponential_qf.hpp
@@ -27,7 +27,7 @@ namespace math {
  * @throw std::domain_error if y is not greater than or equal to 0.
  */
 template <typename Tp, typename Tbeta,
-          require_all_arithmetic_t<Tp, Tbeta> * = nullptr>
+          require_all_arithmetic_t<Tp, Tbeta>* = nullptr>
 inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   static const char* function = "exponential_qf";
   check_positive_finite(function, "Inverse scale parameter", beta);
@@ -40,7 +40,7 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
  * inverse scale parameter.
  * Inverse scale parameter must be greater than 0.
  * p must be bounded by 0 and 1.
- * 
+ *
  * Specialisation for use where any input is an Eigen vector
  *
  * @tparam Tp type of probability input
@@ -59,8 +59,8 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   ref_type_t<Tbeta> beta_ref = beta;
   check_positive_finite(function, "Inverse scale parameter", beta_ref);
   check_bounded(function, "Probability parameter", p_ref, 0, 1);
-  return (-log1p(-as_array_or_scalar(p_ref))
-            / as_array_or_scalar(beta_ref)).matrix();
+  return (-log1p(-as_array_or_scalar(p_ref)) / as_array_or_scalar(beta_ref))
+      .matrix();
 }
 
 }  // namespace math
diff --git a/stan/math/rev/prob/exponential_qf.hpp b/stan/math/rev/prob/exponential_qf.hpp
index 2f2d1fecbc6..ff2d591ee00 100644
--- a/stan/math/rev/prob/exponential_qf.hpp
+++ b/stan/math/rev/prob/exponential_qf.hpp
@@ -26,31 +26,30 @@ namespace math {
  * @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.
  */
-template <typename Tp, typename Tbeta,
-          require_any_st_var<Tp, Tbeta>* = nullptr>
+template <typename Tp, typename Tbeta, require_any_st_var<Tp, Tbeta>* = nullptr>
 inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   using inner_ret_type = decltype(exponential_qf(value_of(p), value_of(beta)));
   using ret_type = return_var_matrix_t<inner_ret_type, Tp, Tbeta>;
 
   arena_t<Tp> arena_p = p;
   arena_t<Tbeta> arena_beta = beta;
-  arena_t<ret_type> ret = exponential_qf(value_of(arena_p),
-                                         value_of(arena_beta));
-    reverse_pass_callback([ret, arena_p, arena_beta]() mutable {
-      decltype(auto) p_val = as_array_or_scalar(value_of(arena_p));
-      decltype(auto) beta_val = as_array_or_scalar(value_of(arena_beta));
-      decltype(auto) ret_adj = as_array_or_scalar(ret.adj());
-      if (!is_constant<Tp>::value) {
-        as_array_or_scalar(forward_as<promote_scalar_t<var, Tp>>(arena_p).adj())
+  arena_t<ret_type> ret
+      = exponential_qf(value_of(arena_p), value_of(arena_beta));
+  reverse_pass_callback([ret, arena_p, arena_beta]() mutable {
+    decltype(auto) p_val = as_array_or_scalar(value_of(arena_p));
+    decltype(auto) beta_val = as_array_or_scalar(value_of(arena_beta));
+    decltype(auto) ret_adj = as_array_or_scalar(ret.adj());
+    if (!is_constant<Tp>::value) {
+      as_array_or_scalar(forward_as<promote_scalar_t<var, Tp>>(arena_p).adj())
           += ret_adj * inv(beta_val - beta_val * p_val);
-      }
-      if (!is_constant<Tbeta>::value) {
-        as_array_or_scalar(
+    }
+    if (!is_constant<Tbeta>::value) {
+      as_array_or_scalar(
           forward_as<promote_scalar_t<var, Tbeta>>(arena_beta).adj())
           += ret_adj * log1m(p_val) / square(beta_val);
-      }
-    });
-    return ret_type(ret);
+    }
+  });
+  return ret_type(ret);
 }
 
 }  // namespace math

From 855adee654376b32c4ccd441cad1807950f4392b Mon Sep 17 00:00:00 2001
From: Andrew Johnson <andrew.johnson@arjohnsonau.com>
Date: Mon, 10 Jan 2022 13:51:29 +0800
Subject: [PATCH 6/7] Eval for use-after-free?

---
 stan/math/prim/prob/exponential_qf.hpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/stan/math/prim/prob/exponential_qf.hpp b/stan/math/prim/prob/exponential_qf.hpp
index ae541f5f80f..2d400602706 100644
--- a/stan/math/prim/prob/exponential_qf.hpp
+++ b/stan/math/prim/prob/exponential_qf.hpp
@@ -27,7 +27,7 @@ namespace math {
  * @throw std::domain_error if y is not greater than or equal to 0.
  */
 template <typename Tp, typename Tbeta,
-          require_all_arithmetic_t<Tp, Tbeta>* = nullptr>
+          require_all_arithmetic_t<Tp, Tbeta> * = nullptr>
 inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   static const char* function = "exponential_qf";
   check_positive_finite(function, "Inverse scale parameter", beta);
@@ -40,7 +40,7 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
  * inverse scale parameter.
  * Inverse scale parameter must be greater than 0.
  * p must be bounded by 0 and 1.
- *
+ * 
  * Specialisation for use where any input is an Eigen vector
  *
  * @tparam Tp type of probability input
@@ -59,8 +59,8 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   ref_type_t<Tbeta> beta_ref = beta;
   check_positive_finite(function, "Inverse scale parameter", beta_ref);
   check_bounded(function, "Probability parameter", p_ref, 0, 1);
-  return (-log1p(-as_array_or_scalar(p_ref)) / as_array_or_scalar(beta_ref))
-      .matrix();
+  return (-log1p(-as_array_or_scalar(p_ref))
+            / as_array_or_scalar(beta_ref)).matrix().eval();
 }
 
 }  // namespace math

From cf33e691308a1bc2c42acce7c383a4e42bf5ae37 Mon Sep 17 00:00:00 2001
From: Stan Jenkins <mc.stanislaw@gmail.com>
Date: Mon, 10 Jan 2022 05:57:07 +0000
Subject: [PATCH 7/7] [Jenkins] auto-formatting by clang-format version
 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final)

---
 stan/math/prim/prob/exponential_qf.hpp | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/stan/math/prim/prob/exponential_qf.hpp b/stan/math/prim/prob/exponential_qf.hpp
index 2d400602706..7d4fa96e0fd 100644
--- a/stan/math/prim/prob/exponential_qf.hpp
+++ b/stan/math/prim/prob/exponential_qf.hpp
@@ -27,7 +27,7 @@ namespace math {
  * @throw std::domain_error if y is not greater than or equal to 0.
  */
 template <typename Tp, typename Tbeta,
-          require_all_arithmetic_t<Tp, Tbeta> * = nullptr>
+          require_all_arithmetic_t<Tp, Tbeta>* = nullptr>
 inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   static const char* function = "exponential_qf";
   check_positive_finite(function, "Inverse scale parameter", beta);
@@ -40,7 +40,7 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
  * inverse scale parameter.
  * Inverse scale parameter must be greater than 0.
  * p must be bounded by 0 and 1.
- * 
+ *
  * Specialisation for use where any input is an Eigen vector
  *
  * @tparam Tp type of probability input
@@ -59,8 +59,9 @@ inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
   ref_type_t<Tbeta> beta_ref = beta;
   check_positive_finite(function, "Inverse scale parameter", beta_ref);
   check_bounded(function, "Probability parameter", p_ref, 0, 1);
-  return (-log1p(-as_array_or_scalar(p_ref))
-            / as_array_or_scalar(beta_ref)).matrix().eval();
+  return (-log1p(-as_array_or_scalar(p_ref)) / as_array_or_scalar(beta_ref))
+      .matrix()
+      .eval();
 }
 
 }  // namespace math