update bessel so that reverse mode is chosen over the perfect forward specialization

Author: SteveBronder

stan/math/prim/fun/bessel_second_kind.hpp

@@ -52,7 +52,8 @@ inline T2 bessel_second_kind(int v, const T2 z) {
  * @param b Second input
  * @return Bessel second kind function applied to the two inputs.
  */
-template <typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr>
+template <typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr,
+ require_not_var_matrix_t<T2>* = nullptr>
 inline auto bessel_second_kind(T1&& a, T2&& b) {
   return apply_scalar_binary(
       [](auto&& c, auto&& d) {

stan/math/rev/fun/bessel_second_kind.hpp

@@ -8,31 +8,26 @@
 namespace stan {
 namespace math {
 
-inline var bessel_second_kind(int v, const var& a) {
-  double ret_val = bessel_second_kind(v, a.val());
-  auto precomp_bessel
-      = v * ret_val / a.val() - bessel_second_kind(v + 1, a.val());
-  return make_callback_var(ret_val, [precomp_bessel, a](auto& vi) mutable {
-    a.adj() += vi.adj() * precomp_bessel;
-  });
-}
-
-/**
- * Overload with `var_value<Matrix>` for `int`, `std::vector<int>`, and
- * `std::vector<std::vector<int>>`
- */
-template <typename T1, typename T2, require_st_integral<T1>* = nullptr,
-          require_eigen_t<T2>* = nullptr>
-inline auto bessel_second_kind(const T1& v, const var_value<T2>& a) {
-  auto ret_val = bessel_second_kind(v, a.val()).array().eval();
-  auto v_map = as_array_or_scalar(v);
-  auto precomp_bessel
-      = to_arena(v_map * ret_val / a.val().array()
-                 - bessel_second_kind(v_map + 1, a.val().array()));
-  return make_callback_var(
-      ret_val.matrix(), [precomp_bessel, a](const auto& vi) mutable {
-        a.adj().array() += vi.adj().array() * precomp_bessel;
-      });
+template <typename T1, typename T2, require_st_integral<T1>* = nullptr, require_var_t<T2>* = nullptr>
+inline auto bessel_second_kind(T1&& v, T2&& a) {
+  if constexpr (is_stan_scalar_v<T2>) {
+    double ret_val = bessel_second_kind(v, a.val());
+    auto precomp_bessel
+        = v * ret_val / a.val() - bessel_second_kind(v + 1, a.val());
+    return make_callback_var(ret_val, [precomp_bessel, a](auto& vi) mutable {
+      a.adj() += vi.adj() * precomp_bessel;
+    });
+  } else {
+    auto ret_val = bessel_second_kind(v, a.val()).array().eval();
+    auto v_map = as_array_or_scalar(v);
+    auto precomp_bessel
+        = to_arena(v_map * ret_val / a.val().array()
+                  - bessel_second_kind(v_map + 1, a.val().array()));
+    return make_callback_var(
+        ret_val.matrix(), [precomp_bessel, a](const auto& vi) mutable {
+          a.adj().array() += vi.adj().array() * precomp_bessel;
+        });
+  }
 }
 
 }  // namespace math

test/unit/math/mix/fun/bessel_second_kind_test.cpp

@@ -25,10 +25,12 @@ TEST(mathMixScalFun, besselSecondKind_vec) {
   Eigen::VectorXd in2(2);
   in2 << 0.5, 3.4;
   stan::test::expect_ad_vectorized_binary(f, std_in1, in2);
+  stan::test::expect_ad_matvar(f, std_in1, in2);
 
   std::vector<std::vector<int>> std_std_in1{std_in1, std_in1};
   Eigen::MatrixXd mat_in2 = in2.replicate(1, 2);
   stan::test::expect_ad_vectorized_binary(f, std_std_in1, mat_in2);
+  stan::test::expect_ad_matvar(f, std_std_in1, mat_in2);
 }
 
 TEST(mathMixScalFun, besselSecondKind_matvec) {