From da42a48375aadce2d642f37d86a5dcc0439cf9f0 Mon Sep 17 00:00:00 2001
From: Steve Bronder <sbronder@flatironinstitute.org>
Date: Fri, 13 Feb 2026 15:28:02 -0500
Subject: [PATCH 1/6] update to allow finite_diff to handle higher order
 autodiff

---
 stan/math/fwd/core/fvar.hpp                  |   3 +-
 stan/math/fwd/fun/tangent_of.hpp             |  65 +++++
 stan/math/fwd/functor/finite_diff.hpp        |  13 +-
 stan/math/fwd/meta/is_fvar.hpp               |   3 -
 stan/math/prim/functor/make_holder_tuple.hpp |   2 +-
 stan/math/prim/functor/map.hpp               |  24 ++
 stan/math/prim/meta/is_fvar.hpp              |  67 +++++
 stan/math/prim/meta/is_matrix_cl.hpp         |   6 +
 stan/math/rev/fun/adjoint_of.hpp             |  96 ++++---
 test/unit/math/fwd/fun/tangent_of_test.cpp   | 219 ++++++++++++++
 test/unit/math/mix/functor/map_test.cpp      | 282 +++++++++++++++++++
 test/unit/math/prim/functor/map_if_test.cpp  |   4 +-
 test/unit/math/rev/fun/adjoint_of_test.cpp   | 166 +++++++++++
 test/unit/math/rev/functor/map_if_test.cpp   |  41 +++
 14 files changed, 934 insertions(+), 57 deletions(-)
 create mode 100644 stan/math/fwd/fun/tangent_of.hpp
 create mode 100644 stan/math/prim/functor/map.hpp
 create mode 100644 test/unit/math/fwd/fun/tangent_of_test.cpp
 create mode 100644 test/unit/math/mix/functor/map_test.cpp
 create mode 100644 test/unit/math/rev/fun/adjoint_of_test.cpp
 create mode 100644 test/unit/math/rev/functor/map_if_test.cpp

diff --git a/stan/math/fwd/core/fvar.hpp b/stan/math/fwd/core/fvar.hpp
index 187d6a5700b..a2dca06eef2 100644
--- a/stan/math/fwd/core/fvar.hpp
+++ b/stan/math/fwd/core/fvar.hpp
@@ -65,7 +65,8 @@ struct fvar {
    *
    * @return tangent of this variable
    */
-  Scalar d() const { return d_; }
+  const Scalar& d() const { return d_; }
+  Scalar& d() { return d_; }
 
   /**
    * Construct a forward variable with zero value and tangent.
diff --git a/stan/math/fwd/fun/tangent_of.hpp b/stan/math/fwd/fun/tangent_of.hpp
new file mode 100644
index 00000000000..33a26ce19b9
--- /dev/null
+++ b/stan/math/fwd/fun/tangent_of.hpp
@@ -0,0 +1,65 @@
+#ifndef STAN_MATH_REV_FUN_TANGENT_OF_HPP
+#define STAN_MATH_REV_FUN_TANGENT_OF_HPP
+
+#include <stan/math/rev/core.hpp>
+#include <stan/math/prim/functor/map.hpp>
+
+namespace stan::math {
+
+/**
+ * Returns a reference to a variable's tangent.
+ * @note In the case of a `std::vector`, this will return back an Eigen map of a vector of the tangents.
+ * @tparam T Any object contains a scalar type of `fvar`
+ * @param x a fvar
+ * @return reference to `x`'s tangent
+ */
+template <typename T>
+inline decltype(auto) tangent_of(T&& x) noexcept {
+  if constexpr (!is_all_fvar_scalar_v<T>) {
+    static_assert(sizeof(std::decay_t<T>*) == 0, "Non-fvar types do not have an tangent.");
+  } else {
+    return map([](auto&& x_) -> decltype(auto) {
+      using arg_t = decltype(x_);
+      using arg_decay_t = std::decay_t<decltype(x_)>;
+      if constexpr (is_tuple_v<arg_decay_t>) {
+        return stan::math::apply([](auto&&... args) -> decltype(auto) {
+          return make_holder_tuple(tangent_of(std::forward<decltype(args)>(args))...);
+        }, std::forward<arg_t>(x_));
+      } else if constexpr (is_fvar_v<arg_decay_t>) {
+        return std::forward<arg_t>(x_).d();
+      } else if constexpr (is_eigen_v<arg_decay_t>) {
+        return make_holder([](auto&& x_hold_) -> decltype(auto) {
+          return std::forward<decltype(x_hold_)>(x_hold_).d();
+        }, std::forward<arg_t>(x_));
+      } else if constexpr (is_std_vector_v<arg_decay_t>) {
+        if constexpr (is_fvar_v<value_type_t<arg_decay_t>>) {
+          return make_holder([](auto&& x_hold_) -> decltype(auto) {
+            constexpr bool const_x = std::is_const_v<std::remove_reference_t<decltype(x_hold_)>>;
+            using base_map_t = Eigen::Matrix<value_type_t<arg_decay_t>, -1, 1>;
+            using map_t = std::conditional_t<const_x, const base_map_t, base_map_t>;
+            return Eigen::Map<map_t>(x_hold_.data(), x_hold_.size()).d();
+          }, std::forward<arg_t>(x_));
+        } else {
+          if constexpr (std::is_rvalue_reference_v<arg_t&&>) {
+              std::vector<decltype(tangent_of(std::move(x[0])))> tangents;
+              tangents.reserve(x_.size());
+              for (auto&& elem : x_) {
+                tangents.push_back(tangent_of(std::move(elem)));
+              }
+              return tangents;
+          } else {
+              std::vector<decltype(tangent_of(x[0]))> tangents;
+              tangents.reserve(x_.size());
+              for (auto&& elem : x_) {
+                tangents.push_back(tangent_of(elem));
+              }
+              return tangents;
+          }
+        }
+      }
+    }, std::forward<T>(x));
+  }
+}
+}  // namespace stan::math
+
+#endif  // TANGENT_OF_HPP
diff --git a/stan/math/fwd/functor/finite_diff.hpp b/stan/math/fwd/functor/finite_diff.hpp
index 6def7b0bc65..be873d5be5e 100644
--- a/stan/math/fwd/functor/finite_diff.hpp
+++ b/stan/math/fwd/functor/finite_diff.hpp
@@ -1,11 +1,12 @@
 #ifndef STAN_MATH_FWD_FUNCTOR_FINITE_DIFF_HPP
 #define STAN_MATH_FWD_FUNCTOR_FINITE_DIFF_HPP
 
-#include <stan/math/prim/meta.hpp>
+#include <stan/math/fwd/fun/tangent_of.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/fun/value_of.hpp>
+#include <stan/math/fwd/fun/sum.hpp>
 #include <stan/math/prim/functor/apply_scalar_binary.hpp>
 #include <stan/math/prim/functor/finite_diff_gradient_auto.hpp>
-#include <stan/math/prim/fun/value_of.hpp>
-#include <stan/math/prim/fun/sum.hpp>
 #include <stan/math/prim/fun/serializer.hpp>
 
 namespace stan {
@@ -43,9 +44,9 @@ inline constexpr double aggregate_tangent(const FuncTangent& tangent,
  */
 template <typename FuncTangent, typename InputArg,
           require_st_fvar<InputArg>* = nullptr>
-inline auto aggregate_tangent(const FuncTangent& tangent, const InputArg& arg) {
-  return sum(apply_scalar_binary([](auto&& x, auto&& y) { return x * y.d_; },
-                                 tangent, arg));
+inline auto aggregate_tangent(FuncTangent&& tangent, InputArg&& arg) {
+  return sum(apply_scalar_binary([](auto&& x, auto&& y) { return x * y; },
+                                 std::forward<FuncTangent>(tangent), tangent_of(std::forward<InputArg>(arg))));
 }
 }  // namespace internal
 
diff --git a/stan/math/fwd/meta/is_fvar.hpp b/stan/math/fwd/meta/is_fvar.hpp
index c5d67de8894..e208d08bc71 100644
--- a/stan/math/fwd/meta/is_fvar.hpp
+++ b/stan/math/fwd/meta/is_fvar.hpp
@@ -21,8 +21,5 @@ struct is_fvar<T,
                std::enable_if_t<internal::is_fvar_impl<std::decay_t<T>>::value>>
     : std::true_type {};
 
-template <typename T>
-inline constexpr bool is_fvar_v = is_fvar<T>::value;
-
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/functor/make_holder_tuple.hpp b/stan/math/prim/functor/make_holder_tuple.hpp
index f436ccfcc28..ef9ff8fbd57 100644
--- a/stan/math/prim/functor/make_holder_tuple.hpp
+++ b/stan/math/prim/functor/make_holder_tuple.hpp
@@ -28,7 +28,7 @@ namespace internal {
 template <typename T>
 struct deduce_cvr {
   using type
-      = std::conditional_t<std::is_rvalue_reference_v<T>, std::decay_t<T>, T&&>;
+      = std::conditional_t<std::is_rvalue_reference_v<T&&>, std::decay_t<T>, T&&>;
 };
 
 template <typename T>
diff --git a/stan/math/prim/functor/map.hpp b/stan/math/prim/functor/map.hpp
new file mode 100644
index 00000000000..4717f1c56bf
--- /dev/null
+++ b/stan/math/prim/functor/map.hpp
@@ -0,0 +1,24 @@
+#ifndef STAN_MATH_PRIM_FUNCTOR_MAP_HPP
+#define STAN_MATH_PRIM_FUNCTOR_MAP_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <type_traits>
+
+namespace stan::math {
+template <typename F, typename T>
+inline decltype(auto) map(F&& f, T&& arg) {
+  if constexpr (is_tuple_v<T>) {
+    return stan::math::apply(
+        [](auto&& f, auto&&... args) {
+          return make_holder_tuple(
+              std::forward<decltype(f)>(f)(
+              std::forward<decltype(args)>(args))...);
+        },
+        std::forward<T>(arg), std::forward<F>(f));
+  } else {
+    return std::forward<F>(f)(std::forward<T>(arg));
+  }
+}
+}
+
+#endif  // STAN_MATH_PRIM_FUNCTOR_MAP_HPP
diff --git a/stan/math/prim/meta/is_fvar.hpp b/stan/math/prim/meta/is_fvar.hpp
index 0a26f93c1b5..4e8687f2bbf 100644
--- a/stan/math/prim/meta/is_fvar.hpp
+++ b/stan/math/prim/meta/is_fvar.hpp
@@ -14,6 +14,9 @@ namespace stan {
 template <typename T, typename = void>
 struct is_fvar : std::false_type {};
 
+template <typename T>
+inline constexpr bool is_fvar_v = is_fvar<T>::value;
+
 /** \ingroup type_trait
  * Specialization for pointers returns the underlying value the pointer is
  * pointing to.
@@ -23,6 +26,70 @@ struct value_type<T, std::enable_if_t<is_fvar<std::decay_t<T>>::value>> {
   using type = typename std::decay_t<T>::Scalar;
 };
 
+template <typename T>
+using has_fvar_scalar_type = is_fvar<scalar_type_t<T>>;
+
+
+namespace internal {
+
+template <typename... Types>
+struct is_any_fvar_scalar_impl {
+  static constexpr bool value
+      = (has_fvar_scalar_type<std::decay_t<Types>>::value || ...);
+};
+
+template <typename... Types>
+struct is_any_fvar_scalar_impl<std::tuple<Types...>> {
+  static constexpr bool value
+      = (is_any_fvar_scalar_impl<std::decay_t<Types>>::value || ...);
+};
+
+template <typename T, typename... VecArgs>
+struct is_any_fvar_scalar_impl<std::vector<T, VecArgs...>> {
+  static constexpr bool value = is_any_fvar_scalar_impl<std::decay_t<T>>::value;
+};
+
+}  // namespace internal
+
+template <typename... Types>
+struct is_any_fvar_scalar
+    : std::disjunction<
+          internal::is_any_fvar_scalar_impl<std::decay_t<Types>>...> {};
+
+template <typename... Types>
+constexpr bool is_any_fvar_scalar_v
+    = is_any_fvar_scalar<std::decay_t<Types>...>::value;
+
+namespace internal {
+template <typename... Types>
+struct is_all_fvar_scalar_impl {
+  static constexpr bool value
+      = (has_fvar_scalar_type<std::decay_t<Types>>::value && ...);
+};
+template <typename... Types>
+struct is_all_fvar_scalar_impl<std::tuple<Types...>> {
+  static constexpr bool value
+      = (is_all_fvar_scalar_impl<std::decay_t<Types>>::value && ...);
+};
+
+template <typename... Types, typename... VecArgs>
+struct is_all_fvar_scalar_impl<std::vector<std::tuple<Types...>, VecArgs...>> {
+  static constexpr bool value
+      = (is_all_fvar_scalar_impl<std::decay_t<Types>>::value && ...);
+};
+
+}  // namespace internal
+
+template <typename... Types>
+struct is_all_fvar_scalar
+    : std::disjunction<
+          internal::is_all_fvar_scalar_impl<std::decay_t<Types>>...> {};
+
+template <typename... Types>
+constexpr bool is_all_fvar_scalar_v
+    = is_all_fvar_scalar<std::decay_t<Types>...>::value;
+
+
 /*! \ingroup require_stan_scalar_real */
 /*! \defgroup fvar_types fvar  */
 /*! \addtogroup fvar_types */
diff --git a/stan/math/prim/meta/is_matrix_cl.hpp b/stan/math/prim/meta/is_matrix_cl.hpp
index a61e9817d59..7375cbb5857 100644
--- a/stan/math/prim/meta/is_matrix_cl.hpp
+++ b/stan/math/prim/meta/is_matrix_cl.hpp
@@ -25,6 +25,12 @@ template <typename T>
 struct is_matrix_cl
     : public std::is_base_of<math::matrix_cl_base, std::decay_t<T>> {};
 
+
+/** \ingroup matrix_cl_group
+ * Checks if the decayed type of T is a matrix_cl.
+ */
+template <typename T>
+inline constexpr bool is_matrix_cl_v = is_matrix_cl<std::decay_t<T>>::value;
 /*! \ingroup matrix_cl_group */
 /*! \defgroup matrix_cl_types matrix_cl  */
 /*! \addtogroup matrix_cl_types */
diff --git a/stan/math/rev/fun/adjoint_of.hpp b/stan/math/rev/fun/adjoint_of.hpp
index 7eb135f3167..67a83234af6 100644
--- a/stan/math/rev/fun/adjoint_of.hpp
+++ b/stan/math/rev/fun/adjoint_of.hpp
@@ -2,58 +2,64 @@
 #define STAN_MATH_REV_FUN_ADJOINT_OF_HPP
 
 #include <stan/math/rev/core.hpp>
+#include <stan/math/prim/functor/map.hpp>
 
-namespace stan {
-namespace math {
-
-namespace internal {
-struct nonexisting_adjoint {
-  template <typename T>
-  nonexisting_adjoint operator+(const T&) {
-    return *this;
-  }
-  template <typename T>
-  nonexisting_adjoint operator+=(T) const {
-    throw std::runtime_error(
-        "internal::nonexisting_adjoint::operator+= should never be called! "
-        "Please file a bug report. rev/fun/adjoint_of.hpp line "
-        + std::to_string(__LINE__));
-  }
-  template <typename T>
-  nonexisting_adjoint operator-=(T) const {
-    throw std::runtime_error(
-        "internal::nonexisting_adjoint::operator-= should never be called! "
-        "Please file a bug report. rev/fun/adjoint_of.hpp line "
-        + std::to_string(__LINE__));
-  }
-};
-}  // namespace internal
+namespace stan::math {
 
 /**
  * Returns a reference to a variable's adjoint.
- *
+ * @note In the case of a `std::vector`, this will return back an Eigen map of a vector of the adjoints.
+ * @tparam T Any object contains a scalar type of `var`
  * @param x a var
  * @return reference to `x`'s adjoint
  */
-template <typename T, require_var_t<T>* = nullptr>
-inline auto& adjoint_of(const T& x) noexcept {
-  return x.adj();
-}
-
-/**
- * Returns a reference to a variable's adjoint. If the input object is not var,
- * it does not have an adjoint and this returns a dummy object. It defines
- * operators += and -=, but they should not actually be called.
- *
- * @param x any non-var object
- * @return a dummy adjoint
- */
-template <typename T, require_not_var_t<T>* = nullptr>
-internal::nonexisting_adjoint adjoint_of(const T& x) {
-  return {};
+template <typename T>
+inline decltype(auto) adjoint_of(T&& x) noexcept {
+  if constexpr (!is_all_var_scalar_v<T>) {
+    static_assert(sizeof(std::decay_t<T>*) == 0, "Non-var types do not have an adjoint.");
+  } else {
+    return map([](auto&& x_) -> decltype(auto) {
+      using arg_t = decltype(x_);
+      using arg_decay_t = std::decay_t<decltype(x_)>;
+      if constexpr (is_tuple_v<arg_decay_t>) {
+        return stan::math::apply([](auto&&... args) -> decltype(auto) {
+          return make_holder_tuple(adjoint_of(std::forward<decltype(args)>(args))...);
+        }, std::forward<arg_t>(x_));
+      } else if constexpr (is_var_v<arg_decay_t>) {
+        return std::forward<arg_t>(x_).adj();
+      } else if constexpr (is_eigen_v<arg_decay_t>) {
+        return make_holder([](auto&& x_hold_) -> decltype(auto) {
+          return std::forward<decltype(x_hold_)>(x_hold_).adj();
+        }, std::forward<arg_t>(x_));
+      } else if constexpr (is_std_vector_v<arg_decay_t>) {
+        if constexpr (is_var_v<value_type_t<arg_decay_t>>) {
+          return make_holder([](auto&& x_hold_) -> decltype(auto) {
+            constexpr bool const_x = std::is_const_v<std::remove_reference_t<decltype(x_hold_)>>;
+            using base_map_t = Eigen::Matrix<var, -1, 1>;
+            using map_t = std::conditional_t<const_x, const base_map_t, base_map_t>;
+            return Eigen::Map<map_t>(x_hold_.data(), x_hold_.size()).adj();
+          }, std::forward<arg_t>(x_));
+        } else {
+          if constexpr (std::is_rvalue_reference_v<arg_t&&>) {
+              std::vector<decltype(adjoint_of(std::move(x[0])))> adjoints;
+              adjoints.reserve(x_.size());
+              for (auto&& elem : x_) {
+                adjoints.push_back(adjoint_of(std::move(elem)));
+              }
+              return adjoints;
+          } else {
+              std::vector<decltype(adjoint_of(x[0]))> adjoints;
+              adjoints.reserve(x_.size());
+              for (auto&& elem : x_) {
+                adjoints.push_back(adjoint_of(elem));
+              }
+              return adjoints;
+          }
+        }
+      }
+    }, std::forward<T>(x));
+  }
 }
-
-}  // namespace math
-}  // namespace stan
+}  // namespace stan::math
 
 #endif  // ADJOINT_OF_HPP
diff --git a/test/unit/math/fwd/fun/tangent_of_test.cpp b/test/unit/math/fwd/fun/tangent_of_test.cpp
new file mode 100644
index 00000000000..9eddb112277
--- /dev/null
+++ b/test/unit/math/fwd/fun/tangent_of_test.cpp
@@ -0,0 +1,219 @@
+#include <stan/math/fwd.hpp>
+#include <stan/math/rev.hpp>
+#include <test/unit/math/prim/util.hpp>
+#include <test/unit/math/rev/util.hpp>
+#include <test/unit/pretty_print_types.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+TEST_F(AgradRev, FwdFun_tangent_of_fvar_var_scalar_lvalue) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+
+  fvar_var x(stan::math::var(2.0), stan::math::var(0.0));
+  auto&& x_tan = stan::math::tangent_of(x);
+  EXPECT_SAME_TYPE(stan::math::var&, decltype(x_tan));
+  x_tan = 1.75;
+  EXPECT_FLOAT_EQ(x.d_.val(), 1.75);
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_fvar_var_scalar_rvalue) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+
+  auto&& x_tan
+      = stan::math::tangent_of(fvar_var(stan::math::var(2.0), stan::math::var(0.0)));
+  EXPECT_SAME_TYPE(stan::math::var&, decltype(x_tan));
+
+  x_tan = 3.25;
+  EXPECT_FLOAT_EQ(x_tan.val(), 3.25);
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_matrix_fvar_var_lvalue_type_and_values) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+  using matrix_fvar_var
+      = Eigen::Matrix<fvar_var, Eigen::Dynamic, Eigen::Dynamic>;
+
+  matrix_fvar_var x(2, 2);
+  for (Eigen::Index i = 0; i < x.rows(); ++i) {
+    for (Eigen::Index j = 0; j < x.cols(); ++j) {
+      x(i, j) = fvar_var(stan::math::var(0.0), stan::math::var(0.0));
+    }
+  }
+  decltype(auto) x_tan = stan::math::tangent_of(x);
+  EXPECT_SAME_TYPE(decltype(x.d()), decltype(x_tan));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 1.0, 2.0, 3.0, 4.0;
+  x_tan = expected;
+  for (Eigen::Index i = 0; i < x.rows(); ++i) {
+    for (Eigen::Index j = 0; j < x.cols(); ++j) {
+      EXPECT_FLOAT_EQ(x(i, j).d_.val(), expected(i, j));
+    }
+  }
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_matrix_fvar_var_lvalue) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+  using matrix_fvar_var
+      = Eigen::Matrix<fvar_var, Eigen::Dynamic, Eigen::Dynamic>;
+
+  matrix_fvar_var x(2, 3);
+  for (Eigen::Index i = 0; i < x.rows(); ++i) {
+    for (Eigen::Index j = 0; j < x.cols(); ++j) {
+      x(i, j) = fvar_var(stan::math::var(0.0), stan::math::var(0.0));
+    }
+  }
+  auto&& x_tan = stan::math::tangent_of(x);
+
+  Eigen::MatrixXd expected(2, 3);
+  expected << 0.5, 1.5, 2.5, 3.5, 4.5, 5.5;
+  x_tan = expected;
+
+  for (Eigen::Index i = 0; i < x.rows(); ++i) {
+    for (Eigen::Index j = 0; j < x.cols(); ++j) {
+      EXPECT_FLOAT_EQ(x(i, j).d_.val(), expected(i, j));
+    }
+  }
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_accepts_eigen_matrix_fvar_var_expression_rvalue) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+  using matrix_fvar_var
+      = Eigen::Matrix<fvar_var, Eigen::Dynamic, Eigen::Dynamic>;
+
+  matrix_fvar_var a(2, 2);
+  matrix_fvar_var b(2, 2);
+  for (Eigen::Index i = 0; i < a.rows(); ++i) {
+    for (Eigen::Index j = 0; j < a.cols(); ++j) {
+      const double aval = static_cast<double>(i + j + 1);
+      const double bval = static_cast<double>(2 * i + j + 0.5);
+      a(i, j) = fvar_var(stan::math::var(0.0), stan::math::var(aval));
+      b(i, j) = fvar_var(stan::math::var(0.0), stan::math::var(bval));
+    }
+  }
+  auto a_plus_b = a + b;
+  decltype(auto) expr_tan = stan::math::tangent_of(a_plus_b);
+  Eigen::MatrixXd expected(2, 2);
+  expected << 1.5, 3.5, 4.5, 6.5;
+  EXPECT_MATRIX_EQ(stan::math::value_of(expr_tan), expected);
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_std_vector_fvar_var_lvalue) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+
+  std::vector<fvar_var> x(3);
+  x[0] = fvar_var(stan::math::var(1.0), stan::math::var(0.0));
+  x[1] = fvar_var(stan::math::var(2.0), stan::math::var(0.0));
+  x[2] = fvar_var(stan::math::var(3.0), stan::math::var(0.0));
+  const auto x_size_before = x.size();
+
+  auto&& x_tan = stan::math::tangent_of(x);
+  EXPECT_EQ(x_tan.size(), x_size_before);
+  x_tan[0] = 2.0;
+  x_tan[1] = 4.0;
+  x_tan[2] = 6.0;
+
+  EXPECT_FLOAT_EQ(x[0].d_.val(), 2.0);
+  EXPECT_FLOAT_EQ(x[1].d_.val(), 4.0);
+  EXPECT_FLOAT_EQ(x[2].d_.val(), 6.0);
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_std_vector_fvar_var_rvalue) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+
+  auto x_tan
+      = stan::math::tangent_of(std::vector<fvar_var>{
+          fvar_var(stan::math::var(1.0), stan::math::var(0.0)),
+          fvar_var(stan::math::var(2.0), stan::math::var(0.0))});
+
+  EXPECT_EQ(x_tan.size(), 2);
+  x_tan[0] = 5.0;
+  x_tan[1] = 7.0;
+  EXPECT_FLOAT_EQ(x_tan[0].val(), 5.0);
+  EXPECT_FLOAT_EQ(x_tan[1].val(), 7.0);
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_nested_std_vector_fvar_var_lvalue) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+
+  std::vector<std::vector<fvar_var>> x{
+      {fvar_var(stan::math::var(1.0), stan::math::var(0.0)),
+       fvar_var(stan::math::var(2.0), stan::math::var(0.0))},
+      {fvar_var(stan::math::var(3.0), stan::math::var(0.0))}};
+
+  const auto outer_size_before = x.size();
+  const auto first_inner_size_before = x[0].size();
+  const auto second_inner_size_before = x[1].size();
+  auto&& x_tan = stan::math::tangent_of(x);
+  EXPECT_EQ(x_tan.size(), outer_size_before);
+  EXPECT_EQ(x_tan[0].size(), first_inner_size_before);
+  EXPECT_EQ(x_tan[1].size(), second_inner_size_before);
+
+  x_tan[0][0] = 1.5;
+  x_tan[0][1] = 2.5;
+  x_tan[1][0] = 3.5;
+
+  EXPECT_FLOAT_EQ(x[0][0].d_.val(), 1.5);
+  EXPECT_FLOAT_EQ(x[0][1].d_.val(), 2.5);
+  EXPECT_FLOAT_EQ(x[1][0].d_.val(), 3.5);
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_nested_std_vector_fvar_var_rvalue) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+
+  auto x_tan = stan::math::tangent_of(std::vector<std::vector<fvar_var>>{
+      {fvar_var(stan::math::var(1.0), stan::math::var(0.0)),
+       fvar_var(stan::math::var(2.0), stan::math::var(0.0))},
+      {fvar_var(stan::math::var(3.0), stan::math::var(0.0)),
+       fvar_var(stan::math::var(4.0), stan::math::var(0.0)),
+       fvar_var(stan::math::var(5.0), stan::math::var(0.0))}});
+  ASSERT_EQ(x_tan.size(), 2);
+  EXPECT_EQ(x_tan[0].size(), 2);
+  EXPECT_EQ(x_tan[1].size(), 3);
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_tuple_mixed_inputs) {
+  using fvar_var = stan::math::fvar<stan::math::var>;
+  using stan::test::is_lvalue_ref_element_v;
+
+  fvar_var scalar(stan::math::var(1.0), stan::math::var(0.0));
+  Eigen::Matrix<fvar_var, Eigen::Dynamic, Eigen::Dynamic> mat(2, 2);
+  for (Eigen::Index i = 0; i < mat.rows(); ++i) {
+    for (Eigen::Index j = 0; j < mat.cols(); ++j) {
+      mat(i, j) = fvar_var(stan::math::var(0.0), stan::math::var(0.0));
+    }
+  }
+  std::vector<fvar_var> vec{
+      fvar_var(stan::math::var(0.0), stan::math::var(0.0)),
+      fvar_var(stan::math::var(0.0), stan::math::var(0.0)),
+      fvar_var(stan::math::var(0.0), stan::math::var(0.0))};
+
+  auto out = stan::math::tangent_of(std::forward_as_tuple(scalar, mat, vec));
+  using out_t = decltype(out);
+  EXPECT_EQ(std::tuple_size_v<out_t>, 3);
+  EXPECT_TRUE((is_lvalue_ref_element_v<0, out_t>));
+
+  std::get<0>(out) = 7.0;
+  Eigen::MatrixXd mat_expected = Eigen::MatrixXd::Constant(2, 2, 3.0);
+  std::get<1>(out) = mat_expected;
+  for (int i = 0; i < 3; ++i) {
+    std::get<2>(out)[i] = 4.0;
+  }
+
+  EXPECT_FLOAT_EQ(scalar.d_.val(), 7.0);
+  for (Eigen::Index i = 0; i < mat.rows(); ++i) {
+    for (Eigen::Index j = 0; j < mat.cols(); ++j) {
+      EXPECT_FLOAT_EQ(mat(i, j).d_.val(), mat_expected(i, j));
+    }
+  }
+  EXPECT_FLOAT_EQ(vec[0].d_.val(), 4.0);
+  EXPECT_FLOAT_EQ(vec[1].d_.val(), 4.0);
+  EXPECT_FLOAT_EQ(vec[2].d_.val(), 4.0);
+}
+
+TEST_F(AgradRev, FwdFun_tangent_of_empty_tuple) {
+  auto out = stan::math::tangent_of(std::tuple<>{});
+  EXPECT_EQ(std::tuple_size_v<decltype(out)>, 0);
+}
diff --git a/test/unit/math/mix/functor/map_test.cpp b/test/unit/math/mix/functor/map_test.cpp
new file mode 100644
index 00000000000..e37acce083a
--- /dev/null
+++ b/test/unit/math/mix/functor/map_test.cpp
@@ -0,0 +1,282 @@
+#include <stan/math/mix.hpp>
+#include <test/unit/math/prim/util.hpp>
+#include <test/unit/pretty_print_types.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+#include <numeric>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+namespace {
+
+struct mixed_output_functor {
+  int operator()(double x) const { return static_cast<int>(x * 2.0); }
+
+  double operator()(const stan::math::var& x) const { return x.val(); }
+
+  Eigen::Index operator()(const Eigen::VectorXd& x) const { return x.size(); }
+
+  double operator()(const Eigen::MatrixXd& x) const { return x.sum(); }
+
+  std::size_t operator()(const std::vector<int>& x) const { return x.size(); }
+
+  std::vector<int> operator()(const std::vector<std::vector<int>>& x) const {
+    std::vector<int> out;
+    out.reserve(x.size());
+    for (const auto& inner : x) {
+      out.push_back(std::accumulate(inner.begin(), inner.end(), 0));
+    }
+    return out;
+  }
+};
+
+TEST(MathMixFunctor, map_scalar_support_and_return_type_changes) {
+  int i = 4;
+  auto&& i_ref = stan::math::map(
+      [](auto& x) -> auto& { return x; }, i);
+  EXPECT_SAME_TYPE(int&, decltype(i_ref));
+  i_ref = 11;
+  EXPECT_EQ(i, 11);
+
+  const int ci = 5;
+  auto&& ci_ref = stan::math::map(
+      [](const auto& x) -> const auto& { return x; }, ci);
+  EXPECT_SAME_TYPE(const int&, decltype(ci_ref));
+
+  auto as_bool = stan::math::map([](double x) { return x > 0.0; }, -1.5);
+  EXPECT_TYPE(bool, as_bool);
+  EXPECT_FALSE(as_bool);
+
+  stan::math::var ad = 2.0;
+  auto ad_out = stan::math::map([](const auto& x) { return x * 3.0; }, ad);
+  EXPECT_TYPE(stan::math::var, ad_out);
+  EXPECT_FLOAT_EQ(ad_out.val(), 6.0);
+
+  stan::math::fvar<stan::math::var> fvad(2.0, 1.0);
+  auto fvad_out
+      = stan::math::map([](const auto& x) { return x * 4.0; }, fvad);
+  EXPECT_TYPE(stan::math::fvar<stan::math::var>, fvad_out);
+  EXPECT_FLOAT_EQ(fvad_out.val_.val(), 8.0);
+  EXPECT_FLOAT_EQ(fvad_out.d_.val(), 4.0);
+}
+
+TEST(MathMixFunctor, map_non_tuple_accepts_rvalue_scalar) {
+  auto out = stan::math::map([](double&& x) { return x + 1.0; }, 2.5);
+  EXPECT_TYPE(double, out);
+  EXPECT_FLOAT_EQ(out, 3.5);
+}
+
+TEST(MathMixFunctor, map_non_tuple_eigen_vector) {
+  Eigen::VectorXd eigen_vec(3);
+  eigen_vec << 1.0, 2.0, 3.0;
+  auto vec_out = stan::math::map(
+                      [](const Eigen::VectorXd& x) {
+                        return (2.0 * x.array()).matrix();
+                      },
+                      eigen_vec)
+                     .eval();
+  EXPECT_SAME_TYPE(Eigen::VectorXd, std::decay_t<decltype(vec_out)>);
+  Eigen::VectorXd expected_vec(3);
+  expected_vec << 2.0, 4.0, 6.0;
+  EXPECT_MATRIX_EQ(vec_out, expected_vec);
+}
+
+TEST(MathMixFunctor, map_non_tuple_eigen_matrix) {
+  Eigen::MatrixXd eigen_mat(2, 2);
+  eigen_mat << 1.0, 2.0, 3.0, 4.0;
+  auto mat_out
+      = stan::math::map([](const Eigen::MatrixXd& x) { return (x * 2.0); },
+                        eigen_mat)
+            .eval();
+  EXPECT_SAME_TYPE(Eigen::MatrixXd, std::decay_t<decltype(mat_out)>);
+  EXPECT_MATRIX_EQ(mat_out, (eigen_mat * 2.0).eval());
+}
+
+TEST(MathMixFunctor, map_non_tuple_accepts_eigen_vector_expression) {
+  Eigen::VectorXd lhs(3);
+  lhs << 1.0, 2.0, 3.0;
+  Eigen::VectorXd rhs(3);
+  rhs << 4.0, 5.0, 6.0;
+
+  auto vec_out
+      = stan::math::map([](const auto& x) { return x.eval(); }, lhs + rhs)
+            .eval();
+  EXPECT_SAME_TYPE(Eigen::VectorXd, std::decay_t<decltype(vec_out)>);
+  EXPECT_MATRIX_EQ(vec_out, (lhs + rhs).eval());
+}
+
+TEST(MathMixFunctor, map_non_tuple_accepts_eigen_matrix_expression) {
+  Eigen::MatrixXd lhs(2, 2);
+  lhs << 1.0, 2.0, 3.0, 4.0;
+  Eigen::MatrixXd rhs(2, 2);
+  rhs << -1.0, 0.5, 1.5, -2.0;
+
+  auto mat_out
+      = stan::math::map([](const auto& x) { return x.eval(); }, lhs * rhs)
+            .eval();
+  EXPECT_SAME_TYPE(Eigen::MatrixXd, std::decay_t<decltype(mat_out)>);
+  EXPECT_MATRIX_EQ(mat_out, (lhs * rhs).eval());
+}
+
+TEST(MathMixFunctor, map_non_tuple_std_vector) {
+  std::vector<int> std_vec{1, 2, 3};
+  auto std_vec_out = stan::math::map(
+      [](const std::vector<int>& x) {
+        std::vector<double> out;
+        out.reserve(x.size());
+        for (int v : x) {
+          out.push_back(v + 0.5);
+        }
+        return out;
+      },
+      std_vec);
+  EXPECT_SAME_TYPE(std::vector<double>, std::decay_t<decltype(std_vec_out)>);
+  ASSERT_EQ(std_vec_out.size(), std_vec.size());
+  EXPECT_FLOAT_EQ(std_vec_out[0], 1.5);
+  EXPECT_FLOAT_EQ(std_vec_out[1], 2.5);
+  EXPECT_FLOAT_EQ(std_vec_out[2], 3.5);
+}
+
+TEST(MathMixFunctor, map_non_tuple_accepts_rvalue_std_vector) {
+  auto std_vec_out = stan::math::map(
+      [](std::vector<int>&& x) {
+        x.push_back(4);
+        return std::move(x);
+      },
+      std::vector<int>{1, 2, 3});
+  EXPECT_SAME_TYPE(std::vector<int>, std::decay_t<decltype(std_vec_out)>);
+  std::vector<int> expected_std_vec{1, 2, 3, 4};
+  EXPECT_EQ(std_vec_out, expected_std_vec);
+}
+
+TEST(MathMixFunctor, map_non_tuple_nested_std_vector) {
+  std::vector<std::vector<int>> nested_vec{{1, 2, 3}, {4}, {}};
+  auto nested_out = stan::math::map(
+      [](const std::vector<std::vector<int>>& x) {
+        std::vector<int> out;
+        out.reserve(x.size());
+        for (const auto& inner : x) {
+          out.push_back(std::accumulate(inner.begin(), inner.end(), 0));
+        }
+        return out;
+      },
+      nested_vec);
+  EXPECT_SAME_TYPE(std::vector<int>, std::decay_t<decltype(nested_out)>);
+  std::vector<int> expected_nested{6, 4, 0};
+  EXPECT_EQ(nested_out, expected_nested);
+}
+
+TEST(MathMixFunctor, map_non_tuple_accepts_rvalue_nested_std_vector) {
+  auto nested_out = stan::math::map(
+      [](std::vector<std::vector<int>>&& x) {
+        x.push_back({7, 8});
+        return std::move(x);
+      },
+      std::vector<std::vector<int>>{{1, 2}, {3, 4}});
+  EXPECT_SAME_TYPE(std::vector<std::vector<int>>,
+                   std::decay_t<decltype(nested_out)>);
+  ASSERT_EQ(nested_out.size(), 3);
+  EXPECT_EQ(nested_out.back(), (std::vector<int>{7, 8}));
+}
+
+TEST(MathMixFunctor, map_tuple_mixed_inputs_respect_output_types) {
+  using stan::test::is_ref_element_v;
+  using stan::test::is_same_tuple_element_v;
+
+  double scalar = 1.5;
+  stan::math::var ad_scalar = 2.5;
+  Eigen::VectorXd eigen_vec(2);
+  eigen_vec << 3.0, 4.0;
+  Eigen::MatrixXd eigen_mat(2, 2);
+  eigen_mat << 1.0, 2.0, 3.0, 4.0;
+  std::vector<int> std_vec{5, 6, 7};
+  std::vector<std::vector<int>> nested_vec{{1, 2}, {3, 4, 5}};
+
+  auto out = stan::math::map(
+      mixed_output_functor{},
+      std::forward_as_tuple(scalar, ad_scalar, eigen_vec, eigen_mat, std_vec,
+                            nested_vec));
+  using out_t = decltype(out);
+  EXPECT_EQ(std::tuple_size_v<out_t>, 6);
+  EXPECT_TRUE((is_same_tuple_element_v<0, out_t, int>));
+  EXPECT_TRUE((is_same_tuple_element_v<1, out_t, double>));
+  EXPECT_TRUE((is_same_tuple_element_v<2, out_t, Eigen::Index>));
+  EXPECT_TRUE((is_same_tuple_element_v<3, out_t, double>));
+  EXPECT_TRUE((is_same_tuple_element_v<4, out_t, std::size_t>));
+  EXPECT_TRUE((is_same_tuple_element_v<5, out_t, std::vector<int>>));
+  EXPECT_FALSE((is_ref_element_v<0, out_t>));
+  EXPECT_FALSE((is_ref_element_v<1, out_t>));
+  EXPECT_FALSE((is_ref_element_v<2, out_t>));
+  EXPECT_FALSE((is_ref_element_v<3, out_t>));
+  EXPECT_FALSE((is_ref_element_v<4, out_t>));
+  EXPECT_FALSE((is_ref_element_v<5, out_t>));
+
+  EXPECT_EQ(std::get<0>(out), 3);
+  EXPECT_FLOAT_EQ(std::get<1>(out), 2.5);
+  EXPECT_EQ(std::get<2>(out), 2);
+  EXPECT_FLOAT_EQ(std::get<3>(out), 10.0);
+  EXPECT_EQ(std::get<4>(out), 3);
+  std::vector<int> expected_nested{3, 12};
+  EXPECT_EQ(std::get<5>(out), expected_nested);
+}
+
+TEST(MathMixFunctor, map_tuple_preserves_lvalue_references_when_returned) {
+  using stan::test::is_lvalue_ref_element_v;
+
+  int scalar = 2;
+  stan::math::var ad_scalar = 3.0;
+  Eigen::VectorXd eigen_vec(2);
+  eigen_vec << 1.0, 2.0;
+  std::vector<int> std_vec{5, 6};
+  std::vector<std::vector<int>> nested_vec{{1}, {2, 3}};
+
+  auto out = stan::math::map(
+      [](auto&& x) -> decltype(auto) { return std::forward<decltype(x)>(x); },
+      std::forward_as_tuple(scalar, ad_scalar, eigen_vec, std_vec, nested_vec));
+  using out_t = decltype(out);
+  EXPECT_TRUE((is_lvalue_ref_element_v<0, out_t>));
+  EXPECT_TRUE((is_lvalue_ref_element_v<1, out_t>));
+  EXPECT_TRUE((is_lvalue_ref_element_v<2, out_t>));
+  EXPECT_TRUE((is_lvalue_ref_element_v<3, out_t>));
+  EXPECT_TRUE((is_lvalue_ref_element_v<4, out_t>));
+
+  std::get<0>(out) = 11;
+  std::get<1>(out) = 7.0;
+  std::get<2>(out)(0) = -1.0;
+  std::get<3>(out)[1] = 42;
+  std::get<4>(out)[1].push_back(9);
+
+  EXPECT_EQ(scalar, 11);
+  EXPECT_FLOAT_EQ(ad_scalar.val(), 7.0);
+  EXPECT_FLOAT_EQ(eigen_vec(0), -1.0);
+  EXPECT_EQ(std_vec[1], 42);
+  EXPECT_EQ(nested_vec[1].back(), 9);
+}
+
+TEST(MathMixFunctor, map_tuple_rvalues_are_stored_by_value) {
+  using stan::test::is_ref_element_v;
+
+  auto out = stan::math::map(
+      [](auto&& x) -> decltype(auto) { return std::forward<decltype(x)>(x); },
+      std::make_tuple(4.0, std::vector<int>{1, 2, 3},
+                      Eigen::VectorXd::Constant(2, 5.0)));
+  using out_t = decltype(out);
+  EXPECT_FALSE((is_ref_element_v<0, out_t>));
+  EXPECT_FALSE((is_ref_element_v<1, out_t>));
+  EXPECT_FALSE((is_ref_element_v<2, out_t>));
+
+  EXPECT_FLOAT_EQ(std::get<0>(out), 4.0);
+  std::vector<int> expected_std_vec{1, 2, 3};
+  EXPECT_EQ(std::get<1>(out), expected_std_vec);
+  EXPECT_FLOAT_EQ(std::get<2>(out)(0), 5.0);
+  EXPECT_FLOAT_EQ(std::get<2>(out)(1), 5.0);
+}
+
+TEST(MathMixFunctor, map_empty_tuple_returns_empty_tuple) {
+  auto out = stan::math::map([](auto&& x) { return x; }, std::tuple<>{});
+  EXPECT_EQ(std::tuple_size_v<decltype(out)>, 0);
+}
+
+}  // namespace
diff --git a/test/unit/math/prim/functor/map_if_test.cpp b/test/unit/math/prim/functor/map_if_test.cpp
index a96640ba0f9..63b2a3d5a40 100644
--- a/test/unit/math/prim/functor/map_if_test.cpp
+++ b/test/unit/math/prim/functor/map_if_test.cpp
@@ -10,12 +10,14 @@ using is_floating_point_decay = std::is_floating_point<std::decay_t<T>>;
 TEST(MathPrim, map_if_base) {
   using stan::math::map_if;
   auto x
-      = map_if<is_floating_point_decay>([](auto&& x) noexcept { return x + 1; },
+      = map_if<is_floating_point_decay>([](auto&& x) noexcept { return static_cast<double>(x) + 1; },
                                         std::forward_as_tuple(1, 2.0, 3, 4.0));
   EXPECT_EQ(std::get<0>(x), 1);
   EXPECT_EQ(std::get<1>(x), 3.0);
   EXPECT_EQ(std::get<2>(x), 3);
   EXPECT_EQ(std::get<3>(x), 5.0);
+  static_assert(std::is_same_v<decltype(x), std::tuple<int, double, int, double>>,
+                "Should be tuple of int, double, int, double!");
 }
 
 TEST(MathPrim, map_if_eigen) {
diff --git a/test/unit/math/rev/fun/adjoint_of_test.cpp b/test/unit/math/rev/fun/adjoint_of_test.cpp
new file mode 100644
index 00000000000..d8865ed65a1
--- /dev/null
+++ b/test/unit/math/rev/fun/adjoint_of_test.cpp
@@ -0,0 +1,166 @@
+#include <stan/math/rev.hpp>
+#include <test/unit/math/prim/util.hpp>
+#include <test/unit/math/rev/util.hpp>
+#include <test/unit/pretty_print_types.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+TEST_F(AgradRev, RevFun_adjoint_of_var_scalar_lvalue) {
+  stan::math::var x = 2.0;
+  auto&& x_adj = stan::math::adjoint_of(x);
+  EXPECT_SAME_TYPE(double&, decltype(x_adj));
+  x_adj = 1.75;
+  EXPECT_FLOAT_EQ(x.adj(), 1.75);
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_var_scalar_rvalue) {
+  auto&& x_adj = stan::math::adjoint_of(stan::math::var(2.0));
+  EXPECT_SAME_TYPE(double&, decltype(x_adj));
+
+  x_adj = 3.25;
+  EXPECT_FLOAT_EQ(x_adj, 3.25);
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_var_value_matrix_lvalue) {
+  using stan::math::var_value;
+  var_value<Eigen::MatrixXd> x(Eigen::MatrixXd::Zero(2, 2));
+  auto&& x_adj = stan::math::adjoint_of(x);
+  EXPECT_SAME_TYPE(decltype(x.adj()), decltype(x_adj));
+
+  Eigen::MatrixXd expected(2, 2);
+  expected << 1.0, 2.0, 3.0, 4.0;
+  x_adj = expected;
+  EXPECT_MATRIX_EQ(x.adj(), expected);
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_matrix_var_lvalue) {
+  using matrix_v = Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic>;
+  matrix_v x = Eigen::MatrixXd::Zero(2, 3).template cast<stan::math::var>();
+  auto&& x_adj = stan::math::adjoint_of(x);
+
+  Eigen::MatrixXd expected(2, 3);
+  expected << 0.5, 1.5, 2.5, 3.5, 4.5, 5.5;
+  x_adj = expected;
+
+  for (Eigen::Index i = 0; i < x.rows(); ++i) {
+    for (Eigen::Index j = 0; j < x.cols(); ++j) {
+      EXPECT_FLOAT_EQ(x(i, j).adj(), expected(i, j));
+    }
+  }
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_accepts_eigen_var_value_expression_rvalue) {
+  using stan::math::set_zero_all_adjoints;
+  using stan::math::var_value;
+
+  var_value<Eigen::MatrixXd> a(Eigen::MatrixXd::Random(2, 2));
+  var_value<Eigen::MatrixXd> b(Eigen::MatrixXd::Random(2, 2));
+  auto a_plus_b = a + b;
+  decltype(auto) expr_adj = stan::math::adjoint_of(a_plus_b);
+  auto sum_var = stan::math::sum(a_plus_b);
+  stan::math::grad(sum_var.vi_);
+  Eigen::MatrixXd expected = Eigen::MatrixXd::Constant(2, 2, 1);
+  EXPECT_MATRIX_EQ(expr_adj, expected);
+  EXPECT_MATRIX_EQ(a.adj(), expected);
+  EXPECT_MATRIX_EQ(b.adj(), expected);
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_std_vector_var_lvalue) {
+  std::vector<stan::math::var> x(3);
+  x[0] = 1.0;
+  x[1] = 2.0;
+  x[2] = 3.0;
+  const auto x_size_before = x.size();
+
+  auto&& x_adj = stan::math::adjoint_of(x);
+  EXPECT_EQ(x.size(), x_size_before);
+
+  Eigen::VectorXd expected(3);
+  expected << 2.0, 4.0, 6.0;
+  x_adj = expected;
+
+  EXPECT_FLOAT_EQ(x[0].adj(), 2.0);
+  EXPECT_FLOAT_EQ(x[1].adj(), 4.0);
+  EXPECT_FLOAT_EQ(x[2].adj(), 6.0);
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_std_vector_var_rvalue) {
+  auto x_adj = stan::math::adjoint_of(
+      std::vector<stan::math::var>{stan::math::var(1.0), stan::math::var(2.0)});
+
+  EXPECT_EQ(x_adj.size(), 2);
+  Eigen::VectorXd expected(2);
+  expected << 5.0, 7.0;
+  x_adj = expected;
+  EXPECT_MATRIX_EQ(x_adj, expected);
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_nested_std_vector_var_lvalue) {
+  std::vector<std::vector<stan::math::var>> x{
+      {stan::math::var(1.0), stan::math::var(2.0)},
+      {stan::math::var(3.0)}};
+
+  const auto outer_size_before = x.size();
+  const auto first_inner_size_before = x[0].size();
+  const auto second_inner_size_before = x[1].size();
+  auto&& x_adj = stan::math::adjoint_of(x);
+  EXPECT_EQ(x.size(), outer_size_before);
+  EXPECT_EQ(x[0].size(), first_inner_size_before);
+  EXPECT_EQ(x[1].size(), second_inner_size_before);
+
+  Eigen::VectorXd expected0(2);
+  expected0 << 1.5, 2.5;
+  Eigen::VectorXd expected1(1);
+  expected1 << 3.5;
+  x_adj[0] = expected0;
+  x_adj[1] = expected1;
+
+  EXPECT_FLOAT_EQ(x[0][0].adj(), 1.5);
+  EXPECT_FLOAT_EQ(x[0][1].adj(), 2.5);
+  EXPECT_FLOAT_EQ(x[1][0].adj(), 3.5);
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_nested_std_vector_var_rvalue) {
+  auto x_adj = stan::math::adjoint_of(
+      std::vector<std::vector<stan::math::var>>{
+          {stan::math::var(1.0), stan::math::var(2.0)},
+          {stan::math::var(3.0), stan::math::var(4.0), stan::math::var(5.0)}});
+  ASSERT_EQ(x_adj.size(), 2);
+  EXPECT_EQ(x_adj[0].size(), 2);
+  EXPECT_EQ(x_adj[1].size(), 3);
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_tuple_mixed_inputs) {
+  using stan::test::is_lvalue_ref_element_v;
+
+  stan::math::var scalar = 1.0;
+  Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> mat
+      = Eigen::MatrixXd::Zero(2, 2).template cast<stan::math::var>();
+  std::vector<stan::math::var> vec{stan::math::var(0.0), stan::math::var(0.0),
+                                   stan::math::var(0.0)};
+
+  auto out = stan::math::adjoint_of(std::forward_as_tuple(scalar, mat, vec));
+  using out_t = decltype(out);
+  EXPECT_EQ(std::tuple_size_v<out_t>, 3);
+  EXPECT_TRUE((is_lvalue_ref_element_v<0, out_t>));
+
+  std::get<0>(out) = 7.0;
+  Eigen::MatrixXd mat_expected = Eigen::MatrixXd::Constant(2, 2, 3.0);
+  std::get<1>(out) = mat_expected;
+  Eigen::VectorXd vec_expected = Eigen::VectorXd::Constant(3, 4.0);
+  std::get<2>(out) = vec_expected;
+
+  EXPECT_FLOAT_EQ(scalar.adj(), 7.0);
+  EXPECT_MATRIX_EQ(mat.adj(), mat_expected);
+  EXPECT_FLOAT_EQ(vec[0].adj(), 4.0);
+  EXPECT_FLOAT_EQ(vec[1].adj(), 4.0);
+  EXPECT_FLOAT_EQ(vec[2].adj(), 4.0);
+}
+
+TEST_F(AgradRev, RevFun_adjoint_of_empty_tuple) {
+  auto out = stan::math::adjoint_of(std::tuple<>{});
+  EXPECT_EQ(std::tuple_size_v<decltype(out)>, 0);
+}
diff --git a/test/unit/math/rev/functor/map_if_test.cpp b/test/unit/math/rev/functor/map_if_test.cpp
new file mode 100644
index 00000000000..ca32dc0ff72
--- /dev/null
+++ b/test/unit/math/rev/functor/map_if_test.cpp
@@ -0,0 +1,41 @@
+#include <stan/math/rev.hpp>
+#include <test/unit/util.hpp>
+#include <test/unit/math/prim/util.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+#include <stdexcept>
+namespace {
+
+TEST(MathPrim, map_if_base) {
+  using stan::is_any_var_scalar;
+  using stan::math::map_if;
+  using stan::math::var;
+  auto x
+      = map_if<is_any_var_scalar>([](auto&& x) noexcept { return x + 1; },
+                                        std::forward_as_tuple(1, var(2.0), std::vector<double>{3.0, 4.0}, std::vector<int>{4, 3, 2, 1}));
+  static_assert(std::is_same_v<decltype(x), std::tuple<int, var, std::vector<double>, std::vector<int>>>,
+                "Should be tuple of int, var, vector<double>, vector<int>!");
+}
+
+TEST(MathPrim, map_if_eigen) {
+  using stan::math::map_if;
+  using stan::test::is_ref_element_v;
+  using stan::test::is_same_tuple_element_v;
+  Eigen::MatrixXd a = Eigen::MatrixXd::Random(3, 3);
+  std::vector<double> b{1, 2, 3};
+  auto x = map_if<stan::is_eigen>(
+      [](auto&& x) { return (x * x).eval(); },
+      std::forward_as_tuple(a, b, a * a, std::vector<int>{1, 2, 3}));
+  EXPECT_MATRIX_EQ(std::get<0>(x), (a * a).eval());
+  static_assert(is_same_tuple_element_v<0, decltype(x), Eigen::MatrixXd>,
+                "1st should be MatrixXd!");
+  static_assert(!is_ref_element_v<0, decltype(x)>, "0th should be an lvalue!");
+  for (int i = 0; i < b.size(); i++) {
+    EXPECT_EQ(std::get<1>(x)[i], b[i]);
+  }
+  static_assert(is_ref_element_v<1, decltype(x)>,
+                "1st should be an reference!");
+  static_assert(!is_ref_element_v<2, decltype(x)>, "2nd should be an lvalue!");
+  static_assert(!is_ref_element_v<3, decltype(x)>, "3rd should be an lvalue!");
+}
+}  // namespace

From 0f97d9a162fbff7f7140b08836f225a733faaadc Mon Sep 17 00:00:00 2001
From: Steve Bronder <sbronder@flatironinstitute.org>
Date: Thu, 19 Feb 2026 17:09:05 -0500
Subject: [PATCH 2/6] add tuple version of finite_diff_gradient_auto

---
 stan/math/fwd/functor/finite_diff.hpp         |  91 +++-
 .../functor/conditional_copy_and_promote.hpp  | 102 +---
 stan/math/mix/functor/laplace_likelihood.hpp  |   2 +-
 .../mix/functor/laplace_marginal_density.hpp  |   2 +-
 .../functor/finite_diff_gradient_auto.hpp     | 443 +++++++++++++++++-
 stan/math/prim/meta.hpp                       |   1 +
 .../meta/conditional_copy_and_promote.hpp     | 114 +++++
 .../math/rev/core/filter_var_scalar_types.hpp |   6 +-
 .../math/fwd/functor/finite_diff_test.cpp     |  66 +++
 .../finite_diff_gradient_auto_test.cpp        |  38 ++
 10 files changed, 737 insertions(+), 128 deletions(-)
 create mode 100644 stan/math/prim/meta/conditional_copy_and_promote.hpp
 create mode 100644 test/unit/math/fwd/functor/finite_diff_test.cpp

diff --git a/stan/math/fwd/functor/finite_diff.hpp b/stan/math/fwd/functor/finite_diff.hpp
index be873d5be5e..a0b56a5be2d 100644
--- a/stan/math/fwd/functor/finite_diff.hpp
+++ b/stan/math/fwd/functor/finite_diff.hpp
@@ -4,10 +4,10 @@
 #include <stan/math/fwd/fun/tangent_of.hpp>
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun/value_of.hpp>
-#include <stan/math/fwd/fun/sum.hpp>
-#include <stan/math/prim/functor/apply_scalar_binary.hpp>
+#include <stan/math/prim/fun/as_array_or_scalar.hpp>
 #include <stan/math/prim/functor/finite_diff_gradient_auto.hpp>
-#include <stan/math/prim/fun/serializer.hpp>
+#include <stan/math/prim/functor/for_each.hpp>
+#include <tuple>
 
 namespace stan {
 namespace math {
@@ -45,11 +45,59 @@ inline constexpr double aggregate_tangent(const FuncTangent& tangent,
 template <typename FuncTangent, typename InputArg,
           require_st_fvar<InputArg>* = nullptr>
 inline auto aggregate_tangent(FuncTangent&& tangent, InputArg&& arg) {
-  return sum(apply_scalar_binary([](auto&& x, auto&& y) { return x * y; },
-                                 std::forward<FuncTangent>(tangent), tangent_of(std::forward<InputArg>(arg))));
+  auto tangent_arr = as_array_or_scalar(std::forward<FuncTangent>(tangent));
+  auto arg_tangent_arr
+      = as_array_or_scalar(tangent_of(std::forward<InputArg>(arg)));
+  if constexpr (is_stan_scalar_v<std::decay_t<decltype(tangent_arr)>>) {
+    return tangent_arr * arg_tangent_arr;
+  } else {
+    using RetType = return_type_t<std::decay_t<decltype(tangent_arr(0))>,
+                                  std::decay_t<decltype(arg_tangent_arr(0))>>;
+    RetType rtn = 0;
+    for (Eigen::Index i = 0; i < tangent_arr.size(); ++i) {
+      rtn += tangent_arr(i) * arg_tangent_arr(i);
+    }
+    return rtn;
+  }
 }
 }  // namespace internal
 
+template <typename SeqT, SeqT... s, SeqT... t>
+inline constexpr std::integer_sequence<SeqT, s..., t...> concat_sequences(
+    std::integer_sequence<SeqT, s...>, std::integer_sequence<SeqT, t...>) {
+  return {};
+}
+
+template <typename SeqT, SeqT... s, SeqT... t, class... R>
+inline constexpr auto concat_sequences(std::integer_sequence<SeqT, s...>,
+                                std::integer_sequence<SeqT, t...>, R...) {
+  return concat_sequences(std::integer_sequence<SeqT, s..., t...>{}, R{}...);
+}
+
+template <template <std::size_t...> class Predicate, class SeqT, SeqT a>
+constexpr auto filter_single_seq(std::integer_sequence<SeqT, a>) {
+  if constexpr (Predicate<a>::value)
+    return std::integer_sequence<SeqT, a>{};
+  else
+    return std::integer_sequence<SeqT>{};
+}
+
+template <template <std::size_t...> class Predicate, class SeqT, SeqT... b>
+constexpr auto filter_integer_seq(std::integer_sequence<SeqT, b...>) {
+  if constexpr (sizeof...(b) > 0)  // non empty sequence
+    return concat_sequences(
+        filter_single_seq<Predicate>(std::integer_sequence<SeqT, b>{})...);
+  else  // empty sequence case
+    return std::integer_sequence<SeqT>{};
+}
+
+
+template <typename SeqT>
+constexpr std::integer_sequence<SeqT> concat_sequences(
+    std::integer_sequence<SeqT>) {
+  return {};
+}
+
 /**
  * Construct an fvar<T> where the tangent is calculated by finite-differencing.
  * Finite-differencing is only perfomed where the scalar type to be evaluated is
@@ -70,26 +118,23 @@ template <typename F, typename... TArgs,
 inline auto finite_diff(const F& func, const TArgs&... args) {
   using FvarT = return_type_t<TArgs...>;
   using FvarInnerT = typename FvarT::Scalar;
-
-  std::vector<FvarInnerT> serialised_args
-      = serialize<FvarInnerT>(value_of(args)...);
-
-  auto serial_functor = [&](auto&& v) {
-    auto v_deserializer = to_deserializer(v);
-    return func(v_deserializer.read(args)...);
-  };
-
+  auto val_args_tuple = std::forward_as_tuple(value_of(args)...);
+  auto autodiff_args = filter_ad_scalar_types(std::forward_as_tuple(args...));
   FvarInnerT rtn_value;
-  std::vector<FvarInnerT> grad;
-  finite_diff_gradient_auto(serial_functor, serialised_args, rtn_value, grad);
-
+  auto grads = zeroed_container(std::forward_as_tuple(value_of(args)...));
+  constexpr auto autodiff_idxs = filter_integer_seq<internal::contains_autodiff>(std::index_sequence_for<TArgs...>{});
+  constexpr auto grad_idxs = stan::math::apply([](auto&&... grads) {
+    return std::index_sequence_for<decltype(grads)...>{};
+  }, grads);
+  finite_diff_gradient_auto(func, rtn_value, val_args_tuple, grads, autodiff_idxs, grad_idxs);
   FvarInnerT rtn_grad = 0;
-  auto grad_deserializer = to_deserializer(grad);
-  // Use a fold-expression to aggregate tangents for input arguments
-  static_cast<void>(
-      std::initializer_list<int>{(rtn_grad += internal::aggregate_tangent(
-                                      grad_deserializer.read(args), args),
-                                  0)...});
+  stan::math::for_each(
+      [&rtn_grad](auto&& grad_i, auto&& arg_i) {
+        rtn_grad += internal::aggregate_tangent(
+            std::forward<decltype(grad_i)>(grad_i),
+            std::forward<decltype(arg_i)>(arg_i));
+      },
+      grads, autodiff_args);
 
   return FvarT(rtn_value, rtn_grad);
 }
diff --git a/stan/math/mix/functor/conditional_copy_and_promote.hpp b/stan/math/mix/functor/conditional_copy_and_promote.hpp
index cc9410bfa02..0c9ea3f2d77 100644
--- a/stan/math/mix/functor/conditional_copy_and_promote.hpp
+++ b/stan/math/mix/functor/conditional_copy_and_promote.hpp
@@ -1,106 +1,6 @@
 #ifndef STAN_MATH_MIX_FUNCTOR_CONDITIONAL_COPY_AND_PROMOTE_HPP
 #define STAN_MATH_MIX_FUNCTOR_CONDITIONAL_COPY_AND_PROMOTE_HPP
 
-#include <stan/math/mix/functor/hessian_block_diag.hpp>
-#include <stan/math/prim/functor.hpp>
-#include <stan/math/prim/fun.hpp>
-
-namespace stan::math::internal {
-
-/**
- * Decide if object should be deep or shallow copied when
- * using @ref conditional_copy_and_promote .
- */
-enum class COPY_TYPE : uint8_t { SHALLOW = 0, DEEP = 1 };
-
-/**
- * Conditional copy and promote a type's scalar type to a `PromotedType`.
- * @tparam Filter type trait with a static constexpr bool member `value`
- *  that is true if the type should be promoted. Otherwise, the type is
- *  left unchanged.
- * @tparam PromotedType type to promote the scalar to.
- * @tparam CopyType type of copy to perform.
- * @tparam Args variadic arguments.
- * @param args variadic arguments to conditionally copy and promote.
- * @return a tuple where each element is either a reference to the original
- * argument or a promoted copy of the argument.
- */
-template <template <typename...> class Filter,
-          typename PromotedType = stan::math::var,
-          COPY_TYPE CopyType = COPY_TYPE::DEEP, typename... Args>
-inline auto conditional_copy_and_promote(Args&&... args) {
-  return map_if<Filter>(
-      [](auto&& arg) {
-        if constexpr (is_tuple_v<decltype(arg)>) {
-          return stan::math::apply(
-              [](auto&&... inner_args) {
-                return make_holder_tuple(
-                    conditional_copy_and_promote<Filter, PromotedType,
-                                                 CopyType>(
-                        std::forward<decltype(inner_args)>(inner_args))...);
-              },
-              std::forward<decltype(arg)>(arg));
-        } else if constexpr (is_std_vector_v<decltype(arg)>) {
-          std::vector<decltype(conditional_copy_and_promote<
-                               Filter, PromotedType, CopyType>(arg[0]))>
-              ret;
-          for (std::size_t i = 0; i < arg.size(); ++i) {
-            ret.push_back(
-                conditional_copy_and_promote<Filter, PromotedType, CopyType>(
-                    arg[i]));
-          }
-          return ret;
-        } else {
-          if constexpr (CopyType == COPY_TYPE::DEEP) {
-            return stan::math::eval(promote_scalar<PromotedType>(
-                value_of_rec(std::forward<decltype(arg)>(arg))));
-          } else if (CopyType == COPY_TYPE::SHALLOW) {
-            if constexpr (std::is_same_v<PromotedType,
-                                         scalar_type_t<decltype(arg)>>) {
-              return std::forward<decltype(arg)>(arg);
-            } else {
-              return stan::math::eval(promote_scalar<PromotedType>(
-                  std::forward<decltype(arg)>(arg)));
-            }
-          }
-        }
-      },
-      std::forward<Args>(args)...);
-}
-
-/**
- * Conditional deep copy types with a `var` scalar type to `PromotedType`.
- * @tparam PromotedType type to promote the scalar to.
- * @tparam Args variadic arguments.
- * @param args variadic arguments to conditionally copy and promote.
- * @return a tuple where each element is either a reference to the original
- * argument or a promoted copy of the argument.
- */
-template <typename PromotedType, typename... Args>
-inline auto deep_copy_vargs(Args&&... args) {
-  return conditional_copy_and_promote<is_any_var_scalar, PromotedType,
-                                      COPY_TYPE::DEEP>(
-      std::forward<Args>(args)...);
-}
-
-/**
- * Conditional shallow copy types with a `var` scalar type to `PromotedType`.
- * @note This function is useful whenever you are inside of nested autodiff
- *  and want to allow the input arguments from an outer autodiff to be used
- *  in an inner autodiff without making a hard copy of the input arguments.
- * @tparam PromotedType type to promote the scalar to.
- * @tparam Args variadic arguments.
- * @param args variadic arguments to conditionally copy and promote.
- * @return a tuple where each element is either a reference to the original
- * argument or a promoted copy of the argument.
- */
-template <typename PromotedType, typename... Args>
-inline auto shallow_copy_vargs(Args&&... args) {
-  return conditional_copy_and_promote<is_any_var_scalar, PromotedType,
-                                      COPY_TYPE::SHALLOW>(
-      std::forward<Args>(args)...);
-}
-
-}  // namespace stan::math::internal
+#include <stan/math/prim/meta/conditional_copy_and_promote.hpp>
 
 #endif
diff --git a/stan/math/mix/functor/laplace_likelihood.hpp b/stan/math/mix/functor/laplace_likelihood.hpp
index 88c36163b8d..14b1e5e55c0 100644
--- a/stan/math/mix/functor/laplace_likelihood.hpp
+++ b/stan/math/mix/functor/laplace_likelihood.hpp
@@ -2,7 +2,7 @@
 #define STAN_MATH_MIX_FUNCTOR_LAPLACE_LIKELIHOOD_HPP
 
 #include <stan/math/mix/functor/hessian_block_diag.hpp>
-#include <stan/math/mix/functor/conditional_copy_and_promote.hpp>
+#include <stan/math/prim/meta/conditional_copy_and_promote.hpp>
 #include <stan/math/prim/functor.hpp>
 #include <stan/math/prim/fun.hpp>
 
diff --git a/stan/math/mix/functor/laplace_marginal_density.hpp b/stan/math/mix/functor/laplace_marginal_density.hpp
index fe00cf6dccf..c4c900004a9 100644
--- a/stan/math/mix/functor/laplace_marginal_density.hpp
+++ b/stan/math/mix/functor/laplace_marginal_density.hpp
@@ -3,7 +3,7 @@
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/mix/functor/laplace_likelihood.hpp>
 #include <stan/math/mix/functor/laplace_marginal_density_estimator.hpp>
-#include <stan/math/mix/functor/conditional_copy_and_promote.hpp>
+#include <stan/math/prim/meta/conditional_copy_and_promote.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun.hpp>
diff --git a/stan/math/prim/functor/finite_diff_gradient_auto.hpp b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
index 8c9951ae37b..c4d0653d136 100644
--- a/stan/math/prim/functor/finite_diff_gradient_auto.hpp
+++ b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
@@ -3,10 +3,331 @@
 
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/finite_diff_stepsize.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+#include <stan/math/prim/meta/conditional_copy_and_promote.hpp>
+#include <stan/math/prim/functor/apply.hpp>
+#include <stan/math/prim/functor/filter_map.hpp>
+#include <stan/math/prim/functor/for_each.hpp>
+#include <stan/math/prim/functor/make_holder_tuple.hpp>
 #include <cmath>
+#include <functional>
+#include <tuple>
+#include <utility>
+#include <vector>
 
 namespace stan {
 namespace math {
+namespace internal {
+
+/**
+ * Type-dependent false helper for `static_assert` in unreachable branches.
+ */
+template <typename...>
+struct dependent_false : std::false_type {};
+
+/**
+ * Trait indicating whether a type contains at least one autodiff scalar.
+ *
+ * Scalars are checked directly, while container specializations recurse into
+ * contained value types.
+ */
+template <typename T, typename = void>
+struct contains_autodiff : bool_constant<is_autodiff_v<std::decay_t<T>>> {};
+
+template <typename... Args>
+struct contains_autodiff<std::tuple<Args...>, void>
+    : bool_constant<(contains_autodiff<std::decay_t<Args>>::value || ...)> {};
+
+template <typename T, typename... VecArgs>
+struct contains_autodiff<std::vector<T, VecArgs...>, void>
+    : bool_constant<contains_autodiff<std::decay_t<T>>::value> {};
+
+template <typename T>
+inline constexpr bool contains_autodiff_v
+    = contains_autodiff<std::decay_t<T>>::value;
+
+template <typename T>
+struct contains_autodiff_decayed : contains_autodiff<std::decay_t<T>> {};
+
+template <typename T, typename = void>
+struct count_autodiff_args;
+
+template <typename... Args>
+struct count_autodiff_args<std::tuple<Args...>, void>
+    : std::integral_constant<
+          std::size_t,
+          (static_cast<std::size_t>(contains_autodiff_v<Args>) + ... + 0)> {};
+
+template <typename TupleArgs>
+inline constexpr std::size_t count_autodiff_args_v
+    = count_autodiff_args<std::decay_t<TupleArgs>>::value;
+
+template <typename T, template <class...> class Cond, std::size_t pos>
+struct indices_of_if_impl {
+  static constexpr std::size_t value = (Cond<T>::value ? pos : -1);
+};
+
+/**
+ * Wrapper that carries a top-level argument by reference and exposes a
+ * potentially different mask type for filtering.
+ *
+ * @tparam MaskT type used for filter decisions
+ * @tparam ArgT referenced argument type
+ */
+template <typename MaskT, typename ArgT>
+struct masked_arg_ref {
+  using arg_t = std::remove_reference_t<ArgT>;
+  std::reference_wrapper<arg_t> arg_;
+
+  explicit masked_arg_ref(arg_t& arg) : arg_(arg) {}
+  inline arg_t& get() const { return arg_.get(); }
+};
+
+/**
+ * Specialization so `filter_map<contains_autodiff>` can filter wrapped
+ * top-level arguments based on the explicit mask type.
+ *
+ * @tparam MaskT mask type used for autodiff detection
+ * @tparam ArgT referenced argument type
+ */
+template <typename MaskT, typename ArgT>
+struct contains_autodiff<masked_arg_ref<MaskT, ArgT>, void>
+    : contains_autodiff<std::decay_t<MaskT>> {};
+
+/**
+ * Trait to detect wrapped top-level arguments.
+ */
+template <typename T>
+struct is_masked_arg_ref : std::false_type {};
+
+template <typename MaskT, typename ArgT>
+struct is_masked_arg_ref<masked_arg_ref<MaskT, ArgT>> : std::true_type {};
+
+/**
+ * Internal implementation for creating a tuple of masked top-level argument
+ * wrappers.
+ *
+ * @tparam TupleMask tuple type used for autodiff mask decisions
+ * @tparam TupleArgs tuple instance to filter
+ * @tparam Is compile-time tuple indices
+ * @param args tuple instance to filter
+ * @return tuple of wrapped top-level arguments carrying explicit mask types
+ */
+template <typename TupleMask, typename TupleArgs, std::size_t... Is>
+inline auto make_masked_arg_ref_tuple_impl(
+    TupleArgs&& args, std::index_sequence<Is...>) {
+  return make_holder_tuple(
+      masked_arg_ref<std::tuple_element_t<Is, std::decay_t<TupleMask>>,
+                     decltype(std::get<Is>(std::forward<TupleArgs>(args)))>(
+          std::get<Is>(std::forward<TupleArgs>(args)))...);
+}
+
+/**
+ * Create a tuple of masked top-level argument wrappers.
+ *
+ * @tparam TupleMask tuple type used for autodiff mask decisions
+ * @tparam TupleArgs tuple instance to wrap
+ * @param args tuple instance to wrap
+ * @return tuple of wrapped top-level arguments carrying explicit mask types
+ */
+template <typename TupleMask, typename TupleArgs>
+inline auto make_masked_arg_ref_tuple(TupleArgs&& args) {
+  return make_masked_arg_ref_tuple_impl<TupleMask>(
+      std::forward<TupleArgs>(args),
+      std::make_index_sequence<std::tuple_size<std::decay_t<TupleMask>>::value>{
+      });
+}
+
+/**
+ * Filter a tuple and return a tuple containing references to the elements
+ * whose types satisfy `contains_autodiff`.
+ *
+ * This mirrors the `filter_var_scalar_types` pattern in rev.
+ *
+ * @tparam T input tuple or nested tuple-like structure
+ * @param t input to filter
+ * @return tuple of references to autodiff-containing top-level arguments
+ */
+template <typename T>
+inline constexpr decltype(auto) filter_autodiff_types(T&& t) {
+  return stan::math::filter_map<contains_autodiff>(
+      [](auto&& arg) -> decltype(auto) {
+        using arg_t = std::decay_t<decltype(arg)>;
+        if constexpr (is_tuple_v<arg_t>) {
+          return filter_autodiff_types(std::forward<decltype(arg)>(arg));
+        } else if constexpr (is_masked_arg_ref<arg_t>::value) {
+          return std::forward<decltype(arg)>(arg).get();
+        } else {
+          return std::forward<decltype(arg)>(arg);
+        }
+      },
+      std::forward<T>(t));
+}
+
+/**
+ * Filter top-level arguments using an explicit tuple mask.
+ *
+ * @tparam TupleMask tuple type used for autodiff mask decisions
+ * @tparam TupleArgs tuple instance to filter
+ * @param args tuple instance to filter
+ * @return tuple of references to autodiff-containing top-level arguments
+ */
+template <typename TupleMask, typename TupleArgs>
+inline constexpr decltype(auto) filter_autodiff_types_with_mask(
+    TupleArgs&& args) {
+  return filter_autodiff_types(
+      make_masked_arg_ref_tuple<TupleMask>(std::forward<TupleArgs>(args)));
+}
+
+/**
+ * Builds a zero-initialized gradient container matching `arg` shape for
+ * autodiff leaves, while preserving non-autodiff leaves by value.
+ *
+ * @tparam ScalarT scalar type
+ * @tparam Arg argument/container type used as the shape template
+ * @param arg shape template argument
+ * @return gradient container aligned with `arg`
+ */
+template <typename Arg>
+inline auto zeroed_container(const Arg& arg) {
+  using ArgDec = std::decay_t<Arg>;
+  using PartialType = partials_type_t<scalar_type_t<ArgDec>>;
+  if constexpr (is_tuple_v<ArgDec>) {
+    return stan::math::apply(
+        [](const auto&... inner_args) {
+          return std::make_tuple(zeroed_container(inner_args)...);
+        }, arg);
+  } else if constexpr (is_std_vector_v<ArgDec>) {
+    std::vector<PartialType> ret;
+    ret.reserve(arg.size());
+    for (const auto& el : arg) {
+      ret.push_back(zeroed_container(el));
+    }
+    return ret;
+  } else if constexpr (is_eigen_v<ArgDec>) {
+      using eigen_t = promote_scalar_t<PartialType, plain_type_t<ArgDec>>;
+      return eigen_t::Zero(arg.rows(), arg.cols()).eval();
+  } else if constexpr (is_stan_scalar_v<ArgDec>) {
+      return PartialType(0);
+  } else {
+    static_assert(
+        dependent_false<ArgDec>::value,
+        "Unsupported gradient container in finite_diff_gradient_auto.");
+  }
+}
+
+/**
+ * Extracts a double-valued coordinate for finite-difference step-size
+ * computation from scalar-like autodiff stacks.
+ *
+ * @tparam T scalar type
+ * @param x scalar value
+ * @return value used to compute finite-difference step size
+ */
+template <typename T>
+inline double finite_diff_stepsize_value(const T& x) {
+  if constexpr (is_var<std::decay_t<T>>::value) {
+    return x.val();
+  } else if constexpr (is_fvar<std::decay_t<T>>::value) {
+    return finite_diff_stepsize_value(x.val_);
+  } else {
+    return value_of_rec(x);
+  }
+}
+
+/**
+ * Recursively traverses `arg_mask`, `arg_work`, and `grad` with matching
+ * structure and invokes `fn` on each autodiff scalar coordinate pair
+ * `(arg_work_coord, grad_coord)`.
+ *
+ * Non-autodiff scalar leaves in `arg_mask` are skipped.
+ *
+ * @tparam ArgMask mask argument type (uses original scalar categories)
+ * @tparam ArgWork mutable working argument type
+ * @tparam Grad gradient container type
+ * @tparam F callable with signature compatible with `(arg_work_coord, grad_coord)`
+ * @param arg_mask mask argument used to decide perturbable coordinates
+ * @param arg_work mutable working argument coordinates
+ * @param grad gradient output coordinates
+ * @param fn callback applied at each perturbable coordinate
+ */
+template <typename ArgMask, typename Grad, typename F>
+inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, Grad&& grad, F&& fn) {
+  using arg_mask_t = std::decay_t<ArgMask>;
+  if constexpr (is_stan_scalar_v<arg_mask_t>) {
+    if constexpr (is_autodiff_v<arg_mask_t>) {
+      fn(arg_mask, grad);
+    }
+  } else if constexpr (is_eigen_v<arg_mask_t>) {
+    for (Eigen::Index i = 0; i < arg_mask.size(); ++i) {
+      for_each_coordinate_pair_mut(arg_mask.coeffRef(i),
+                                   grad.coeffRef(i), fn);
+    }
+  } else if constexpr (is_std_vector_v<arg_mask_t>) {
+    for (std::size_t i = 0; i < arg_mask.size(); ++i) {
+      for_each_coordinate_pair_mut(arg_mask[i], grad[i], fn);
+    }
+  } else if constexpr (is_tuple_v<arg_mask_t>) {
+    static_assert(std::tuple_size<std::decay_t<ArgMask>>::value
+                      == std::tuple_size<std::decay_t<Grad>>::value,
+                  "Tuple size mismatch in finite_diff_gradient_auto traversal");
+    stan::math::for_each(
+        [&fn](auto&& arg_mask_i, auto&& grad_i) {
+          for_each_coordinate_pair_mut(arg_mask_i, grad_i, fn);
+        },
+        std::forward<ArgMask>(arg_mask),
+        std::forward<Grad>(grad));
+  } else {
+    static_assert(dependent_false<arg_mask_t>::value,
+                  "Unsupported container in finite_diff_gradient_auto.");
+  }
+}
+
+/**
+ * Applies the six-point finite-difference stencil to every perturbable
+ * coordinate in one autodiff top-level argument and writes the corresponding
+ * gradient coordinate.
+ *
+ * Coordinates are always restored to their original value after each stencil.
+ *
+ * @tparam ScalarT scalar type of function value and gradients
+ * @tparam F callable type
+ * @tparam ArgMask mask argument type (original categories)
+ * @tparam ArgWork mutable working argument type
+ * @tparam NthGrad gradient container type for this top-level argument
+ * @tparam TupleArgsWork mutable tuple of all working arguments
+ * @param f function being finite-differenced
+ * @param arg_mask mask argument used to choose perturbable coordinates
+ * @param arg_work mutable working argument coordinates
+ * @param grad gradient output container for this top-level argument
+ * @param args_work mutable tuple of all working arguments passed to `f`
+ */
+template <typename ScalarT, typename F, typename ArgMask,
+          typename NthGrad, typename TupleArgsWork>
+inline void finite_diff_gradient_auto_impl(const F& f, ArgMask&& arg_mask,
+                                           NthGrad&& grad,
+                                           TupleArgsWork&& args_work) {
+  static constexpr int h_scale[6] = {3, 2, 1, -3, -2, -1};
+  static constexpr int mults[6] = {1, -9, 45, -1, 9, -45};
+
+  for_each_coordinate_pair_mut(
+      std::forward<ArgMask>(arg_mask), std::forward<NthGrad>(grad),
+      [&f, &args_work](auto& arg_coord, auto& grad_coord) {
+        const auto orig = arg_coord;
+        const double h = finite_diff_stepsize(finite_diff_stepsize_value(orig));
+        ScalarT delta_f = 0;
+        for (int j = 0; j < 6; ++j) {
+          arg_coord = orig + h * h_scale[j];
+          delta_f += stan::math::apply(
+                         [&f](auto&&... args_i) { return f(args_i...); },
+                         args_work)
+                     * mults[j];
+        }
+        arg_coord = orig;
+        grad_coord = delta_f / (60 * h);
+      });
+}
+}  // namespace internal
 
 /**
  * Calculate the value and the gradient of the specified function
@@ -47,7 +368,8 @@ namespace math {
  * @param[out] grad_fx gradient of function at argument
  */
 template <typename F, typename VectorT, typename GradVectorT,
-          typename ScalarT = return_type_t<VectorT>>
+          typename ScalarT = return_type_t<VectorT>,
+          require_vector_t<VectorT>* = nullptr>
 inline void finite_diff_gradient_auto(const F& f, VectorT&& x, ScalarT& fx,
                                       GradVectorT& grad_fx) {
   using EigT = Eigen::Matrix<ScalarT, -1, 1>;
@@ -72,6 +394,125 @@ inline void finite_diff_gradient_auto(const F& f, VectorT&& x, ScalarT& fx,
   });
 }
 
+/**
+ * Filter a tuple and return a tuple with references to the types with an autodiff
+ * scalar type.
+ * @tparam T Possibly a tuple, std::vector, Eigen type, or scalar
+ * @param[in] t Input to filter
+ * @return Filtered input with only var scalar types
+ */
+template <typename T>
+inline constexpr decltype(auto) filter_ad_scalar_types(T&& t) {
+  return stan::math::filter_map<internal::contains_autodiff>(
+      [](auto&& arg) -> decltype(auto) {
+        if constexpr (is_tuple_v<std::decay_t<decltype(arg)>>) {
+          return filter_ad_scalar_types(std::forward<decltype(arg)>(arg));
+        } else {
+          return std::forward<decltype(arg)>(arg);
+        }
+      },
+      std::forward<T>(t));
+}
+
+/**
+ * Calculate the function value and gradients for heterogeneous tuple arguments
+ * while preserving non-autodiff argument scalar types at the call boundary.
+ *
+ * Pattern used by this overload:
+ * 1. Build `args_ad` by filtering top-level autodiff-containing entries from
+ *    the original `args` tuple.
+ * 2. Build one mutable working tuple `args_work` where autodiff-containing
+ *    entries are deep-copied and promoted to `ScalarT`, and non-autodiff
+ *    entries are forwarded unchanged.
+ * 3. Build `args_work_ad` by applying the same top-level filter shape as
+ *    `args_ad` (via the `TupleArgs` mask) so both filtered tuples align.
+ * 4. Zip `args_ad`, `args_work_ad`, and `grads` and run
+ *    `finite_diff_gradient_auto_impl` for each top-level autodiff entry.
+ *
+ * The split between original mask (`args_ad`) and mutable working values
+ * (`args_work_ad`) ensures coordinate traversal only perturbs locations that
+ * correspond to autodiff coordinates in the original argument pack, while all
+ * function evaluations run against the promoted mutable working tuple.
+ *
+ * @tparam F callable type
+ * @tparam ScalarT scalar type used for function value and gradient storage
+ * @tparam TupleArgs tuple of input arguments
+ * @tparam TupleGrads tuple of gradient containers aligned to autodiff
+ *         top-level arguments
+ * @param[in] f function to differentiate
+ * @param[out] fx function value at `args`
+ * @param[in] args tuple of function arguments
+ * @param[out] grads compact tuple of gradient containers
+ */
+template <typename F, typename ScalarT, typename TupleArgs, typename TupleGrads,
+          require_tuple_t<TupleArgs>* = nullptr,
+          require_tuple_t<TupleGrads>* = nullptr,
+          require_t<bool_constant<
+              (internal::count_autodiff_args_v<TupleArgs> > 0)>>* = nullptr>
+inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
+                                      TupleGrads&& grads) {
+  fx = stan::math::apply([&f](auto&&... args_i) { return f(args_i...); },
+                         args);
+  stan::math::for_each([&f, &args, &grads](auto&& arg, auto&& grad) {
+    if (!std::is_integral_v<scalar_type_t<decltype(grad)>> && !std::is_integral_v<scalar_type_t<decltype(arg)>>) {
+      internal::finite_diff_gradient_auto_impl<ScalarT>(
+          f, arg, grad, args);
+    }
+  }, args, grads);
+ }
+
+
+/**
+ * Calculate the function value and gradients for heterogeneous tuple arguments
+ * while preserving non-autodiff argument scalar types at the call boundary.
+ *
+ * Pattern used by this overload:
+ * 1. Build `args_ad` by filtering top-level autodiff-containing entries from
+ *    the original `args` tuple.
+ * 2. Build one mutable working tuple `args_work` where autodiff-containing
+ *    entries are deep-copied and promoted to `ScalarT`, and non-autodiff
+ *    entries are forwarded unchanged.
+ * 3. Build `args_work_ad` by applying the same top-level filter shape as
+ *    `args_ad` (via the `TupleArgs` mask) so both filtered tuples align.
+ * 4. Zip `args_ad`, `args_work_ad`, and `grads` and run
+ *    `finite_diff_gradient_auto_impl` for each top-level autodiff entry.
+ *
+ * The split between original mask (`args_ad`) and mutable working values
+ * (`args_work_ad`) ensures coordinate traversal only perturbs locations that
+ * correspond to autodiff coordinates in the original argument pack, while all
+ * function evaluations run against the promoted mutable working tuple.
+ *
+ * @tparam F callable type
+ * @tparam ScalarT scalar type used for function value and gradient storage
+ * @tparam TupleArgs tuple of input arguments
+ * @tparam TupleGrads tuple of gradient containers aligned to autodiff
+ *         top-level arguments
+ * @param[in] f function to differentiate
+ * @param[out] fx function value at `args`
+ * @param[in] args tuple of function arguments
+ * @param[out] grads compact tuple of gradient containers
+ */
+template <typename F, typename ScalarT, typename TupleArgs, typename TupleGrads,
+          require_tuple_t<TupleArgs>* = nullptr,
+          require_tuple_t<TupleGrads>* = nullptr,
+          std::size_t... Idxs,
+          std::size_t... GradIdxs,
+          require_t<bool_constant<
+              (internal::count_autodiff_args_v<TupleArgs> > 0)>>* = nullptr>
+inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
+                                      TupleGrads&& grads,
+                                      std::index_sequence<Idxs...> /* idxs */,
+                                      std::index_sequence<GradIdxs...> /* grad_idxs */) {
+  fx = stan::math::apply([&f](auto&&... args_i) { return f(args_i...); },
+                         args);
+   using Swallow = int[];
+  static_cast<void>(Swallow{(
+    static_cast<void>(internal::finite_diff_gradient_auto_impl<ScalarT>(
+        f, std::get<Idxs>(args),
+        std::get<GradIdxs>(grads), args)), 0)...});
+
+
+ }
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/meta.hpp b/stan/math/prim/meta.hpp
index 5407f128133..af2966e2351 100644
--- a/stan/math/prim/meta.hpp
+++ b/stan/math/prim/meta.hpp
@@ -136,5 +136,6 @@
 #include <stan/math/prim/meta/void_t.hpp>
 #include <stan/math/prim/meta/StdVectorBuilder.hpp>
 #include <stan/math/prim/meta/VectorBuilder.hpp>
+#include <stan/math/prim/meta/conditional_copy_and_promote.hpp>
 
 #endif
diff --git a/stan/math/prim/meta/conditional_copy_and_promote.hpp b/stan/math/prim/meta/conditional_copy_and_promote.hpp
new file mode 100644
index 00000000000..b895a98eb61
--- /dev/null
+++ b/stan/math/prim/meta/conditional_copy_and_promote.hpp
@@ -0,0 +1,114 @@
+#ifndef STAN_MATH_PRIM_META_CONDITIONAL_COPY_AND_PROMOTE_HPP
+#define STAN_MATH_PRIM_META_CONDITIONAL_COPY_AND_PROMOTE_HPP
+
+#include <stan/math/prim/fun/eval.hpp>
+#include <stan/math/prim/fun/promote_scalar.hpp>
+#include <stan/math/prim/fun/value_of.hpp>
+#include <stan/math/prim/functor/map_if.hpp>
+#include <cstdint>
+#include <vector>
+
+namespace stan {
+namespace math {
+namespace internal {
+
+/**
+ * Decide if object should be deep or shallow copied when
+ * using @ref conditional_copy_and_promote .
+ */
+enum class COPY_TYPE : uint8_t { SHALLOW = 0, DEEP = 1 };
+
+/**
+ * Conditional copy and promote a type's scalar type to a `PromotedType`.
+ * @tparam Filter type trait with a static constexpr bool member `value`
+ *  that is true if the type should be promoted. Otherwise, the type is
+ *  left unchanged.
+ * @tparam PromotedType type to promote the scalar to.
+ * @tparam CopyType type of copy to perform.
+ * @tparam Args variadic arguments.
+ * @param args variadic arguments to conditionally copy and promote.
+ * @return a tuple where each element is either a reference to the original
+ * argument or a promoted copy of the argument.
+ */
+template <template <typename...> class Filter, typename PromotedType = double,
+          COPY_TYPE CopyType = COPY_TYPE::DEEP, typename... Args>
+inline auto conditional_copy_and_promote(Args&&... args) {
+  return map_if<Filter>(
+      [](auto&& arg) {
+        using arg_t = std::decay_t<decltype(arg)>;
+        if constexpr (is_tuple_v<arg_t>) {
+          return stan::math::apply(
+              [](auto&&... inner_args) {
+                return make_holder_tuple(
+                    conditional_copy_and_promote<Filter, PromotedType,
+                                                 CopyType>(
+                        std::forward<decltype(inner_args)>(inner_args))...);
+              },
+              std::forward<decltype(arg)>(arg));
+        } else if constexpr (is_std_vector_v<arg_t>) {
+          using elem_t = std::decay_t<decltype(conditional_copy_and_promote<
+              Filter, PromotedType, CopyType>(
+              std::declval<const value_type_t<arg_t>&>()))>;
+          std::vector<elem_t> ret;
+          ret.reserve(arg.size());
+          for (std::size_t i = 0; i < arg.size(); ++i) {
+            ret.push_back(
+                conditional_copy_and_promote<Filter, PromotedType, CopyType>(
+                    arg[i]));
+          }
+          return ret;
+        } else {
+          if constexpr (CopyType == COPY_TYPE::DEEP) {
+            return stan::math::eval(promote_scalar<PromotedType>(
+                value_of(std::forward<decltype(arg)>(arg))));
+          } else {
+            if constexpr (std::is_same_v<PromotedType, scalar_type_t<arg_t>>) {
+              return std::forward<decltype(arg)>(arg);
+            } else {
+              return stan::math::eval(promote_scalar<PromotedType>(
+                  std::forward<decltype(arg)>(arg)));
+            }
+          }
+        }
+      },
+      std::forward<Args>(args)...);
+}
+
+/**
+ * Conditional deep copy types with a `var` scalar type to `PromotedType`.
+ * @tparam PromotedType type to promote the scalar to.
+ * @tparam Args variadic arguments.
+ * @param args variadic arguments to conditionally copy and promote.
+ * @return a tuple where each element is either a reference to the original
+ * argument or a promoted copy of the argument.
+ */
+template <typename PromotedType, typename... Args>
+inline auto deep_copy_vargs(Args&&... args) {
+  return conditional_copy_and_promote<is_any_var_scalar, PromotedType,
+                                      COPY_TYPE::DEEP>(
+      std::forward<Args>(args)...);
+}
+
+/**
+ * Conditional shallow copy types with a `var` scalar type to `PromotedType`.
+ * @note This function is useful whenever you are inside of nested autodiff
+ *  and want to allow the input arguments from an outer autodiff to be used
+ *  in an inner autodiff without making a hard copy of the input arguments.
+ * @tparam PromotedType type to promote the scalar to.
+ * @tparam Args variadic arguments.
+ * @param args variadic arguments to conditionally copy and promote.
+ * @return a tuple where each element is either a reference to the original
+ * argument or a promoted copy of the argument.
+ */
+template <typename PromotedType, typename... Args>
+inline auto shallow_copy_vargs(Args&&... args) {
+  return conditional_copy_and_promote<is_any_var_scalar, PromotedType,
+                                      COPY_TYPE::SHALLOW>(
+      std::forward<Args>(args)...);
+}
+
+}  // namespace internal
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev/core/filter_var_scalar_types.hpp b/stan/math/rev/core/filter_var_scalar_types.hpp
index 4b56cf80cf7..f18df646d7b 100644
--- a/stan/math/rev/core/filter_var_scalar_types.hpp
+++ b/stan/math/rev/core/filter_var_scalar_types.hpp
@@ -18,7 +18,11 @@ template <typename T>
 inline constexpr decltype(auto) filter_var_scalar_types(T&& t) {
   return stan::math::filter_map<is_any_var_scalar>(
       [](auto&& arg) -> decltype(auto) {
-        return std::forward<decltype(arg)>(arg);
+        if constexpr (is_tuple_v<std::decay_t<decltype(arg)>>) {
+          return filter_var_scalar_types(std::forward<decltype(arg)>(arg));
+        } else {
+          return std::forward<decltype(arg)>(arg);
+        }
       },
       std::forward<T>(t));
 }
diff --git a/test/unit/math/fwd/functor/finite_diff_test.cpp b/test/unit/math/fwd/functor/finite_diff_test.cpp
new file mode 100644
index 00000000000..9cdb4bf1ec3
--- /dev/null
+++ b/test/unit/math/fwd/functor/finite_diff_test.cpp
@@ -0,0 +1,66 @@
+#include <stan/math/fwd.hpp>
+#include <gtest/gtest.h>
+#include <vector>
+
+struct finite_diff_typed_args_functor {
+  template <typename T_theta>
+  auto operator()(const std::vector<T_theta>& theta,
+                  const std::vector<int>& x_i) const {
+    return theta[0] * x_i[0] + theta[1] * x_i[1];
+  }
+};
+
+TEST(FwdFunctor, finiteDiffTypedIntArgs) {
+  using stan::math::fvar;
+  std::vector<fvar<double>> theta(2);
+  theta[0] = fvar<double>(1.2, 0.5);
+  theta[1] = fvar<double>(-0.7, -1.0);
+  const std::vector<int> x_i{3, 4};
+
+  auto out = stan::math::finite_diff(finite_diff_typed_args_functor{}, theta,
+                                     x_i);
+
+  EXPECT_NEAR(0.8, out.val_, 1e-10);
+  EXPECT_NEAR(-2.5, out.d_, 1e-8);
+}
+
+TEST(FwdFunctor, finiteDiffTypedIntArgsFvarVar) {
+  using stan::math::fvar;
+  using stan::math::var;
+  std::vector<fvar<var>> theta(2);
+  theta[0] = fvar<var>(1.2, 0.5);
+  theta[1] = fvar<var>(-0.7, -1.0);
+  const std::vector<int> x_i{3, 4};
+
+  auto out = stan::math::finite_diff(finite_diff_typed_args_functor{}, theta,
+                                     x_i);
+
+  EXPECT_NEAR(0.8, out.val_.val(), 1e-10);
+  EXPECT_NEAR(-2.5, out.d_.val(), 1e-8);
+}
+
+struct integrate_1d_finite_diff_functor {
+  template <typename T_x, typename T_xc, typename T_theta>
+  auto operator()(const T_x& x, const T_xc& /* xc */, const std::vector<T_theta>& theta,
+                  const std::vector<double>& x_r,
+                  const std::vector<int>& x_i, std::ostream* /* msgs */) const {
+    return theta[0] * x + x_r[0] + x_i[0];
+  }
+};
+
+TEST(FwdFunctor, finiteDiffIntegrate1dMixedArgs) {
+  using stan::math::fvar;
+  const double relative_tolerance = 1e-8;
+  std::ostream* msgs = nullptr;
+  std::vector<fvar<double>> theta(1);
+  theta[0] = fvar<double>(2.0, 1.0);
+  const std::vector<double> x_r{1.5};
+  const std::vector<int> x_i{2};
+
+  auto out = stan::math::integrate_1d(integrate_1d_finite_diff_functor{}, 0.0,
+                                      1.0, theta, x_r, x_i, msgs,
+                                      relative_tolerance);
+
+  EXPECT_NEAR(4.5, out.val_, 1e-8);
+  EXPECT_NEAR(0.5, out.d_, 1e-6);
+}
diff --git a/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp b/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp
index eeae792742a..9d135ce0a13 100644
--- a/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp
+++ b/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp
@@ -1,6 +1,7 @@
 #include <stan/math/mix.hpp>
 #include <test/unit/math/expect_near_rel.hpp>
 #include <gtest/gtest.h>
+#include <tuple>
 #include <vector>
 
 template <typename F>
@@ -48,3 +49,40 @@ TEST(MathMixMatFunctor, FiniteDiffGradientAuto) {
   w << 1, 2, 3, 4, 5;
   expect_match_autodiff(log_fun, w);
 }
+
+TEST(MathMixMatFunctor, FiniteDiffGradientAutoTupleNested) {
+  using stan::math::fvar;
+  using fvar_stdvec = std::vector<fvar<double>>;
+  fvar_stdvec a(2);
+  a[0] = fvar<double>(1.3, 0.0);
+  a[1] = fvar<double>(-0.7, 0.0);
+  using fvar_mat = Eigen::Matrix<fvar<double>, -1, 1>;
+  fvar_mat b{{fvar<double>(2.0, 0.0), fvar<double>(-3.0, 0.0)}};
+  const std::vector<int> x_i{5, 6};
+  Eigen::VectorXd c{{4.0}};
+  auto args = std::make_tuple(a,
+    std::make_tuple(b, x_i),
+    c);
+  fvar_mat b_grad = fvar_mat::Zero(b.size());
+  auto grads = std::make_tuple(
+      fvar_stdvec(a.size(), 0.0),
+      std::make_tuple(b_grad, std::vector<int>(x_i.size(), 0)),
+      Eigen::VectorXd{{0.0}});
+
+  auto f = [](const auto& arg_a, const auto& arg_b_tuple, const auto& arg_c) {
+    const auto& arg_b = std::get<0>(arg_b_tuple);
+    const auto& arg_i = std::get<1>(arg_b_tuple);
+    return arg_a[0] * arg_a[1] + arg_b(0) * arg_b(1) + arg_c(0) * arg_i[0];
+  };
+
+  fvar<double> fx = 0;
+  stan::math::finite_diff_gradient_auto(f, fx, args, grads);
+
+  EXPECT_NEAR(13.09, fx.val(), 1e-8);
+  EXPECT_NEAR(a[1].val_, std::get<0>(grads)[0].val_, 1e-8);
+  EXPECT_NEAR(a[0].val_, std::get<0>(grads)[1].val_, 1e-8);
+  EXPECT_NEAR(b(1).val_, std::get<0>(std::get<1>(grads))(0).val_, 1e-8);
+  EXPECT_NEAR(b(0).val_, std::get<0>(std::get<1>(grads))(1).val_, 1e-8);
+  EXPECT_EQ(0, std::get<1>(std::get<1>(grads))[0]);
+  EXPECT_EQ(0, std::get<1>(std::get<1>(grads))[1]);
+}

From 8c49b5ad830f38eb63eba9aa5b059cd206925457 Mon Sep 17 00:00:00 2001
From: Steve Bronder <sbronder@flatironinstitute.org>
Date: Fri, 20 Feb 2026 12:01:49 -0500
Subject: [PATCH 3/6] fixup finite_diff to use the finite_diff_gradient_auto
 with tuples

---
 stan/math/fwd/functor/finite_diff.hpp         |  34 ++--
 .../functor/finite_diff_gradient_auto.hpp     | 177 ++++++------------
 2 files changed, 71 insertions(+), 140 deletions(-)

diff --git a/stan/math/fwd/functor/finite_diff.hpp b/stan/math/fwd/functor/finite_diff.hpp
index a0b56a5be2d..27e38dc61e5 100644
--- a/stan/math/fwd/functor/finite_diff.hpp
+++ b/stan/math/fwd/functor/finite_diff.hpp
@@ -68,36 +68,37 @@ inline constexpr std::integer_sequence<SeqT, s..., t...> concat_sequences(
   return {};
 }
 
+template <typename SeqT>
+constexpr std::integer_sequence<SeqT> concat_sequences(
+    std::integer_sequence<SeqT>) {
+  return {};
+}
+
 template <typename SeqT, SeqT... s, SeqT... t, class... R>
 inline constexpr auto concat_sequences(std::integer_sequence<SeqT, s...>,
-                                std::integer_sequence<SeqT, t...>, R...) {
+                                       std::integer_sequence<SeqT, t...>,
+                                       R...) {
   return concat_sequences(std::integer_sequence<SeqT, s..., t...>{}, R{}...);
 }
 
-template <template <std::size_t...> class Predicate, class SeqT, SeqT a>
+template <template <class...> class Predicate, typename Arg, class SeqT, SeqT a>
 constexpr auto filter_single_seq(std::integer_sequence<SeqT, a>) {
-  if constexpr (Predicate<a>::value)
+  if constexpr (Predicate<Arg>::value)
     return std::integer_sequence<SeqT, a>{};
   else
     return std::integer_sequence<SeqT>{};
 }
 
-template <template <std::size_t...> class Predicate, class SeqT, SeqT... b>
+template <template <class...> class Predicate, typename... Args, class SeqT,
+          SeqT... b>
 constexpr auto filter_integer_seq(std::integer_sequence<SeqT, b...>) {
   if constexpr (sizeof...(b) > 0)  // non empty sequence
-    return concat_sequences(
-        filter_single_seq<Predicate>(std::integer_sequence<SeqT, b>{})...);
+    return concat_sequences(filter_single_seq<Predicate, Args>(
+        std::integer_sequence<SeqT, b>{})...);
   else  // empty sequence case
     return std::integer_sequence<SeqT>{};
 }
 
-
-template <typename SeqT>
-constexpr std::integer_sequence<SeqT> concat_sequences(
-    std::integer_sequence<SeqT>) {
-  return {};
-}
-
 /**
  * Construct an fvar<T> where the tangent is calculated by finite-differencing.
  * Finite-differencing is only perfomed where the scalar type to be evaluated is
@@ -118,15 +119,14 @@ template <typename F, typename... TArgs,
 inline auto finite_diff(const F& func, const TArgs&... args) {
   using FvarT = return_type_t<TArgs...>;
   using FvarInnerT = typename FvarT::Scalar;
-  auto val_args_tuple = std::forward_as_tuple(value_of(args)...);
   auto autodiff_args = filter_ad_scalar_types(std::forward_as_tuple(args...));
   FvarInnerT rtn_value;
-  auto grads = zeroed_container(std::forward_as_tuple(value_of(args)...));
-  constexpr auto autodiff_idxs = filter_integer_seq<internal::contains_autodiff>(std::index_sequence_for<TArgs...>{});
+  auto grads = internal::zeroed_container(std::forward_as_tuple(args...));
+  constexpr auto autodiff_idxs = filter_integer_seq<internal::contains_autodiff, TArgs...>(std::index_sequence_for<TArgs...>{});
   constexpr auto grad_idxs = stan::math::apply([](auto&&... grads) {
     return std::index_sequence_for<decltype(grads)...>{};
   }, grads);
-  finite_diff_gradient_auto(func, rtn_value, val_args_tuple, grads, autodiff_idxs, grad_idxs);
+  finite_diff_gradient_auto(func, rtn_value, std::forward_as_tuple(value_of(args)...), grads, autodiff_idxs, grad_idxs);
   FvarInnerT rtn_grad = 0;
   stan::math::for_each(
       [&rtn_grad](auto&& grad_i, auto&& arg_i) {
diff --git a/stan/math/prim/functor/finite_diff_gradient_auto.hpp b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
index c4d0653d136..d51a0f60961 100644
--- a/stan/math/prim/functor/finite_diff_gradient_auto.hpp
+++ b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
@@ -31,33 +31,35 @@ struct dependent_false : std::false_type {};
  * Scalars are checked directly, while container specializations recurse into
  * contained value types.
  */
-template <typename T, typename = void>
-struct contains_autodiff : bool_constant<is_autodiff_v<std::decay_t<T>>> {};
+template <typename T>
+struct contains_autodiff_impl : bool_constant<is_autodiff_v<std::decay_t<T>>> {};
 
 template <typename... Args>
-struct contains_autodiff<std::tuple<Args...>, void>
-    : bool_constant<(contains_autodiff<std::decay_t<Args>>::value || ...)> {};
+struct contains_autodiff_impl<std::tuple<Args...>>
+    : bool_constant<(contains_autodiff_impl<std::decay_t<Args>>::value || ...)> {};
 
 template <typename T, typename... VecArgs>
-struct contains_autodiff<std::vector<T, VecArgs...>, void>
-    : bool_constant<contains_autodiff<std::decay_t<T>>::value> {};
+struct contains_autodiff_impl<std::vector<T, VecArgs...>>
+    : bool_constant<contains_autodiff_impl<std::decay_t<T>>::value> {};
 
 template <typename T>
-inline constexpr bool contains_autodiff_v
-    = contains_autodiff<std::decay_t<T>>::value;
+struct contains_autodiff : contains_autodiff_impl<std::decay_t<T>> {};
 
 template <typename T>
-struct contains_autodiff_decayed : contains_autodiff<std::decay_t<T>> {};
+inline constexpr bool contains_autodiff_v = contains_autodiff<T>::value;
 
-template <typename T, typename = void>
-struct count_autodiff_args;
+template <typename T>
+struct count_autodiff_args_impl;
 
 template <typename... Args>
-struct count_autodiff_args<std::tuple<Args...>, void>
+struct count_autodiff_args_impl<std::tuple<Args...>>
     : std::integral_constant<
           std::size_t,
           (static_cast<std::size_t>(contains_autodiff_v<Args>) + ... + 0)> {};
 
+template <typename T>
+struct count_autodiff_args : count_autodiff_args_impl<std::decay_t<T>> {};
+
 template <typename TupleArgs>
 inline constexpr std::size_t count_autodiff_args_v
     = count_autodiff_args<std::decay_t<TupleArgs>>::value;
@@ -67,76 +69,6 @@ struct indices_of_if_impl {
   static constexpr std::size_t value = (Cond<T>::value ? pos : -1);
 };
 
-/**
- * Wrapper that carries a top-level argument by reference and exposes a
- * potentially different mask type for filtering.
- *
- * @tparam MaskT type used for filter decisions
- * @tparam ArgT referenced argument type
- */
-template <typename MaskT, typename ArgT>
-struct masked_arg_ref {
-  using arg_t = std::remove_reference_t<ArgT>;
-  std::reference_wrapper<arg_t> arg_;
-
-  explicit masked_arg_ref(arg_t& arg) : arg_(arg) {}
-  inline arg_t& get() const { return arg_.get(); }
-};
-
-/**
- * Specialization so `filter_map<contains_autodiff>` can filter wrapped
- * top-level arguments based on the explicit mask type.
- *
- * @tparam MaskT mask type used for autodiff detection
- * @tparam ArgT referenced argument type
- */
-template <typename MaskT, typename ArgT>
-struct contains_autodiff<masked_arg_ref<MaskT, ArgT>, void>
-    : contains_autodiff<std::decay_t<MaskT>> {};
-
-/**
- * Trait to detect wrapped top-level arguments.
- */
-template <typename T>
-struct is_masked_arg_ref : std::false_type {};
-
-template <typename MaskT, typename ArgT>
-struct is_masked_arg_ref<masked_arg_ref<MaskT, ArgT>> : std::true_type {};
-
-/**
- * Internal implementation for creating a tuple of masked top-level argument
- * wrappers.
- *
- * @tparam TupleMask tuple type used for autodiff mask decisions
- * @tparam TupleArgs tuple instance to filter
- * @tparam Is compile-time tuple indices
- * @param args tuple instance to filter
- * @return tuple of wrapped top-level arguments carrying explicit mask types
- */
-template <typename TupleMask, typename TupleArgs, std::size_t... Is>
-inline auto make_masked_arg_ref_tuple_impl(
-    TupleArgs&& args, std::index_sequence<Is...>) {
-  return make_holder_tuple(
-      masked_arg_ref<std::tuple_element_t<Is, std::decay_t<TupleMask>>,
-                     decltype(std::get<Is>(std::forward<TupleArgs>(args)))>(
-          std::get<Is>(std::forward<TupleArgs>(args)))...);
-}
-
-/**
- * Create a tuple of masked top-level argument wrappers.
- *
- * @tparam TupleMask tuple type used for autodiff mask decisions
- * @tparam TupleArgs tuple instance to wrap
- * @param args tuple instance to wrap
- * @return tuple of wrapped top-level arguments carrying explicit mask types
- */
-template <typename TupleMask, typename TupleArgs>
-inline auto make_masked_arg_ref_tuple(TupleArgs&& args) {
-  return make_masked_arg_ref_tuple_impl<TupleMask>(
-      std::forward<TupleArgs>(args),
-      std::make_index_sequence<std::tuple_size<std::decay_t<TupleMask>>::value>{
-      });
-}
 
 /**
  * Filter a tuple and return a tuple containing references to the elements
@@ -155,8 +87,6 @@ inline constexpr decltype(auto) filter_autodiff_types(T&& t) {
         using arg_t = std::decay_t<decltype(arg)>;
         if constexpr (is_tuple_v<arg_t>) {
           return filter_autodiff_types(std::forward<decltype(arg)>(arg));
-        } else if constexpr (is_masked_arg_ref<arg_t>::value) {
-          return std::forward<decltype(arg)>(arg).get();
         } else {
           return std::forward<decltype(arg)>(arg);
         }
@@ -165,50 +95,31 @@ inline constexpr decltype(auto) filter_autodiff_types(T&& t) {
 }
 
 /**
- * Filter top-level arguments using an explicit tuple mask.
- *
- * @tparam TupleMask tuple type used for autodiff mask decisions
- * @tparam TupleArgs tuple instance to filter
- * @param args tuple instance to filter
- * @return tuple of references to autodiff-containing top-level arguments
- */
-template <typename TupleMask, typename TupleArgs>
-inline constexpr decltype(auto) filter_autodiff_types_with_mask(
-    TupleArgs&& args) {
-  return filter_autodiff_types(
-      make_masked_arg_ref_tuple<TupleMask>(std::forward<TupleArgs>(args)));
-}
-
-/**
- * Builds a zero-initialized gradient container matching `arg` shape for
- * autodiff leaves, while preserving non-autodiff leaves by value.
+ * Build a zero-initialized gradient object for a single autodiff-containing
+ * leaf/container argument.
  *
- * @tparam ScalarT scalar type
- * @tparam Arg argument/container type used as the shape template
- * @param arg shape template argument
- * @return gradient container aligned with `arg`
+ * @tparam Arg leaf/container type
+ * @param arg argument value used for shape
+ * @return zero-initialized gradient container matching `arg` shape
  */
 template <typename Arg>
-inline auto zeroed_container(const Arg& arg) {
+inline auto zeroed_container_leaf(const Arg& arg) {
   using ArgDec = std::decay_t<Arg>;
   using PartialType = partials_type_t<scalar_type_t<ArgDec>>;
-  if constexpr (is_tuple_v<ArgDec>) {
-    return stan::math::apply(
-        [](const auto&... inner_args) {
-          return std::make_tuple(zeroed_container(inner_args)...);
-        }, arg);
-  } else if constexpr (is_std_vector_v<ArgDec>) {
-    std::vector<PartialType> ret;
+  if constexpr (is_std_vector_v<ArgDec>) {
+    using ElemType = std::decay_t<decltype(
+        zeroed_container_leaf(std::declval<const value_type_t<ArgDec>&>()))>;
+    std::vector<ElemType> ret;
     ret.reserve(arg.size());
     for (const auto& el : arg) {
-      ret.push_back(zeroed_container(el));
+      ret.push_back(zeroed_container_leaf(el));
     }
     return ret;
   } else if constexpr (is_eigen_v<ArgDec>) {
-      using eigen_t = promote_scalar_t<PartialType, plain_type_t<ArgDec>>;
-      return eigen_t::Zero(arg.rows(), arg.cols()).eval();
+    using eigen_t = promote_scalar_t<PartialType, plain_type_t<ArgDec>>;
+    return eigen_t::Zero(arg.rows(), arg.cols()).eval();
   } else if constexpr (is_stan_scalar_v<ArgDec>) {
-      return PartialType(0);
+    return PartialType(0);
   } else {
     static_assert(
         dependent_false<ArgDec>::value,
@@ -216,6 +127,30 @@ inline auto zeroed_container(const Arg& arg) {
   }
 }
 
+/**
+ * Build a compact tuple of zero-initialized gradient containers for
+ * autodiff-containing entries in `args`.
+ *
+ * Non-autodiff top-level entries are filtered out.
+ *
+ * @tparam TupleArgs tuple of input arguments
+ * @param args tuple of input arguments
+ * @return tuple of zeroed gradient containers for autodiff-containing entries
+ */
+template <typename TupleArgs, require_tuple_t<TupleArgs>* = nullptr>
+inline auto zeroed_container(TupleArgs&& args) {
+  return stan::math::filter_map<contains_autodiff>(
+      [](auto&& arg) {
+        using Arg = std::decay_t<decltype(arg)>;
+        if constexpr (is_tuple_v<Arg>) {
+          return zeroed_container(std::forward<decltype(arg)>(arg));
+        } else {
+          return zeroed_container_leaf(arg);
+        }
+      },
+      std::forward<TupleArgs>(args));
+}
+
 /**
  * Extracts a double-valued coordinate for finite-difference step-size
  * computation from scalar-like autodiff stacks.
@@ -255,9 +190,7 @@ template <typename ArgMask, typename Grad, typename F>
 inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, Grad&& grad, F&& fn) {
   using arg_mask_t = std::decay_t<ArgMask>;
   if constexpr (is_stan_scalar_v<arg_mask_t>) {
-    if constexpr (is_autodiff_v<arg_mask_t>) {
       fn(arg_mask, grad);
-    }
   } else if constexpr (is_eigen_v<arg_mask_t>) {
     for (Eigen::Index i = 0; i < arg_mask.size(); ++i) {
       for_each_coordinate_pair_mut(arg_mask.coeffRef(i),
@@ -454,7 +387,7 @@ inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
   fx = stan::math::apply([&f](auto&&... args_i) { return f(args_i...); },
                          args);
   stan::math::for_each([&f, &args, &grads](auto&& arg, auto&& grad) {
-    if (!std::is_integral_v<scalar_type_t<decltype(grad)>> && !std::is_integral_v<scalar_type_t<decltype(arg)>>) {
+    if constexpr (!std::is_integral_v<scalar_type_t<decltype(grad)>> && !std::is_integral_v<scalar_type_t<decltype(arg)>>) {
       internal::finite_diff_gradient_auto_impl<ScalarT>(
           f, arg, grad, args);
     }
@@ -496,9 +429,7 @@ template <typename F, typename ScalarT, typename TupleArgs, typename TupleGrads,
           require_tuple_t<TupleArgs>* = nullptr,
           require_tuple_t<TupleGrads>* = nullptr,
           std::size_t... Idxs,
-          std::size_t... GradIdxs,
-          require_t<bool_constant<
-              (internal::count_autodiff_args_v<TupleArgs> > 0)>>* = nullptr>
+          std::size_t... GradIdxs>
 inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
                                       TupleGrads&& grads,
                                       std::index_sequence<Idxs...> /* idxs */,

From 44bc539706cab1755e254a156232c842c9c4a6c3 Mon Sep 17 00:00:00 2001
From: Steve Bronder <sbronder@flatironinstitute.org>
Date: Fri, 20 Feb 2026 12:43:41 -0500
Subject: [PATCH 4/6] cleanup and add tests for meta helpers

---
 stan/math/fwd/functor/finite_diff.hpp         |  58 +--
 stan/math/prim/fun.hpp                        |   1 +
 stan/math/prim/fun/zeroed_container.hpp       |  85 +++++
 .../functor/finite_diff_gradient_auto.hpp     | 349 ++++++------------
 stan/math/prim/meta.hpp                       |   3 +
 stan/math/prim/meta/contains_autodiff.hpp     |  37 ++
 stan/math/prim/meta/filter_index_sequence.hpp |  76 ++++
 .../math/prim/meta/filtered_tuple_indices.hpp |  47 +++
 .../math/prim/fun/zeroed_container_test.cpp   |  56 +++
 .../math/prim/meta/contains_autodiff_test.cpp |  31 ++
 .../prim/meta/filter_index_sequence_test.cpp  |  35 ++
 .../prim/meta/filtered_tuple_indices_test.cpp |  37 ++
 12 files changed, 542 insertions(+), 273 deletions(-)
 create mode 100644 stan/math/prim/fun/zeroed_container.hpp
 create mode 100644 stan/math/prim/meta/contains_autodiff.hpp
 create mode 100644 stan/math/prim/meta/filter_index_sequence.hpp
 create mode 100644 stan/math/prim/meta/filtered_tuple_indices.hpp
 create mode 100644 test/unit/math/prim/fun/zeroed_container_test.cpp
 create mode 100644 test/unit/math/prim/meta/contains_autodiff_test.cpp
 create mode 100644 test/unit/math/prim/meta/filter_index_sequence_test.cpp
 create mode 100644 test/unit/math/prim/meta/filtered_tuple_indices_test.cpp

diff --git a/stan/math/fwd/functor/finite_diff.hpp b/stan/math/fwd/functor/finite_diff.hpp
index 27e38dc61e5..483c82324c0 100644
--- a/stan/math/fwd/functor/finite_diff.hpp
+++ b/stan/math/fwd/functor/finite_diff.hpp
@@ -5,8 +5,11 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun/value_of.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
+#include <stan/math/prim/fun/zeroed_container.hpp>
 #include <stan/math/prim/functor/finite_diff_gradient_auto.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
+#include <stan/math/prim/meta/contains_autodiff.hpp>
+#include <stan/math/prim/meta/filtered_tuple_indices.hpp>
 #include <tuple>
 
 namespace stan {
@@ -28,7 +31,7 @@ template <typename FuncTangent, typename InputArg,
 inline constexpr double aggregate_tangent(const FuncTangent& tangent,
                                           const InputArg& arg) {
   return 0;
-}
+}  // namespace internal
 
 /**
  * Helper function for aggregating tangents if the respective input argument
@@ -60,43 +63,6 @@ inline auto aggregate_tangent(FuncTangent&& tangent, InputArg&& arg) {
     return rtn;
   }
 }
-}  // namespace internal
-
-template <typename SeqT, SeqT... s, SeqT... t>
-inline constexpr std::integer_sequence<SeqT, s..., t...> concat_sequences(
-    std::integer_sequence<SeqT, s...>, std::integer_sequence<SeqT, t...>) {
-  return {};
-}
-
-template <typename SeqT>
-constexpr std::integer_sequence<SeqT> concat_sequences(
-    std::integer_sequence<SeqT>) {
-  return {};
-}
-
-template <typename SeqT, SeqT... s, SeqT... t, class... R>
-inline constexpr auto concat_sequences(std::integer_sequence<SeqT, s...>,
-                                       std::integer_sequence<SeqT, t...>,
-                                       R...) {
-  return concat_sequences(std::integer_sequence<SeqT, s..., t...>{}, R{}...);
-}
-
-template <template <class...> class Predicate, typename Arg, class SeqT, SeqT a>
-constexpr auto filter_single_seq(std::integer_sequence<SeqT, a>) {
-  if constexpr (Predicate<Arg>::value)
-    return std::integer_sequence<SeqT, a>{};
-  else
-    return std::integer_sequence<SeqT>{};
-}
-
-template <template <class...> class Predicate, typename... Args, class SeqT,
-          SeqT... b>
-constexpr auto filter_integer_seq(std::integer_sequence<SeqT, b...>) {
-  if constexpr (sizeof...(b) > 0)  // non empty sequence
-    return concat_sequences(filter_single_seq<Predicate, Args>(
-        std::integer_sequence<SeqT, b>{})...);
-  else  // empty sequence case
-    return std::integer_sequence<SeqT>{};
 }
 
 /**
@@ -121,12 +87,16 @@ inline auto finite_diff(const F& func, const TArgs&... args) {
   using FvarInnerT = typename FvarT::Scalar;
   auto autodiff_args = filter_ad_scalar_types(std::forward_as_tuple(args...));
   FvarInnerT rtn_value;
-  auto grads = internal::zeroed_container(std::forward_as_tuple(args...));
-  constexpr auto autodiff_idxs = filter_integer_seq<internal::contains_autodiff, TArgs...>(std::index_sequence_for<TArgs...>{});
-  constexpr auto grad_idxs = stan::math::apply([](auto&&... grads) {
-    return std::index_sequence_for<decltype(grads)...>{};
-  }, grads);
-  finite_diff_gradient_auto(func, rtn_value, std::forward_as_tuple(value_of(args)...), grads, autodiff_idxs, grad_idxs);
+  auto grads = stan::math::zeroed_container<contains_autodiff>(
+      std::forward_as_tuple(args...));
+  using args_tuple_t = std::tuple<TArgs...>;
+  using autodiff_idxs_t
+      = filtered_tuple_indices_t<contains_autodiff, args_tuple_t>;
+  using grads_t = std::decay_t<decltype(grads)>;
+  using grad_idxs_t = std::make_index_sequence<std::tuple_size<grads_t>::value>;
+  finite_diff_gradient_auto(func, rtn_value, std::forward_as_tuple(args...),
+                            std::forward_as_tuple(value_of(args)...), grads,
+                            autodiff_idxs_t{}, grad_idxs_t{});
   FvarInnerT rtn_grad = 0;
   stan::math::for_each(
       [&rtn_grad](auto&& grad_i, auto&& arg_i) {
diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp
index bb88f8eddf1..518c87efa41 100644
--- a/stan/math/prim/fun.hpp
+++ b/stan/math/prim/fun.hpp
@@ -344,6 +344,7 @@
 #include <stan/math/prim/fun/variance.hpp>
 #include <stan/math/prim/fun/welford_covar_estimator.hpp>
 #include <stan/math/prim/fun/welford_var_estimator.hpp>
+#include <stan/math/prim/fun/zeroed_container.hpp>
 #include <stan/math/prim/fun/zeros_array.hpp>
 #include <stan/math/prim/fun/zeros_int_array.hpp>
 #include <stan/math/prim/fun/zeros_row_vector.hpp>
diff --git a/stan/math/prim/fun/zeroed_container.hpp b/stan/math/prim/fun/zeroed_container.hpp
new file mode 100644
index 00000000000..adceb94574c
--- /dev/null
+++ b/stan/math/prim/fun/zeroed_container.hpp
@@ -0,0 +1,85 @@
+#ifndef STAN_MATH_PRIM_FUN_ZEROED_CONTAINER_HPP
+#define STAN_MATH_PRIM_FUN_ZEROED_CONTAINER_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/functor/filter_map.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+namespace stan {
+namespace math {
+namespace internal {
+
+/**
+ * Type-dependent false helper for `static_assert` in unreachable branches.
+ */
+template <typename...>
+struct zeroed_container_dependent_false : std::false_type {};
+
+/**
+ * Build a zero-initialized container matching the shape of `arg`.
+ *
+ * Scalars produce zero partials scalars, Eigen produces zero matrices of the
+ * same dimensions, and std::vector recurses by element.
+ */
+template <typename Arg>
+inline auto zeroed_container_leaf(const Arg& arg) {
+  using arg_t = std::decay_t<Arg>;
+  using partial_t = partials_type_t<scalar_type_t<arg_t>>;
+  if constexpr (is_std_vector_v<arg_t>) {
+    using elem_t = std::decay_t<decltype(
+        zeroed_container_leaf(std::declval<const value_type_t<arg_t>&>()))>;
+    std::vector<elem_t> ret;
+    ret.reserve(arg.size());
+    for (const auto& element : arg) {
+      ret.push_back(zeroed_container_leaf(element));
+    }
+    return ret;
+  } else if constexpr (is_eigen_v<arg_t>) {
+    using eig_t = promote_scalar_t<partial_t, plain_type_t<arg_t>>;
+    return eig_t::Zero(arg.rows(), arg.cols()).eval();
+  } else if constexpr (is_stan_scalar_v<arg_t>) {
+    return partial_t(0);
+  } else {
+    static_assert(
+        zeroed_container_dependent_false<arg_t>::value,
+        "Unsupported container in zeroed_container.");
+  }
+}
+
+}  // namespace internal
+
+/**
+ * Build a compact tuple of zero-initialized containers for entries in `args`
+ * whose types satisfy `Predicate`.
+ *
+ * Tuple entries that do not satisfy `Predicate` are filtered out. Matching
+ * tuple entries recurse so nested tuple output remains compact as well.
+ *
+ * @tparam Predicate unary type-trait predicate
+ * @tparam TupleArgs tuple input type
+ * @param args tuple of arguments used for shape and filtering
+ * @return tuple of zero-initialized containers for filtered entries
+ */
+template <template <class...> class Predicate, typename TupleArgs,
+          require_tuple_t<TupleArgs>* = nullptr>
+inline auto zeroed_container(TupleArgs&& args) {
+  return stan::math::filter_map<Predicate>(
+      [](auto&& arg) {
+        using arg_t = std::decay_t<decltype(arg)>;
+        if constexpr (is_tuple_v<arg_t>) {
+          return zeroed_container<Predicate>(std::forward<decltype(arg)>(arg));
+        } else {
+          return internal::zeroed_container_leaf(arg);
+        }
+      },
+      std::forward<TupleArgs>(args));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/functor/finite_diff_gradient_auto.hpp b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
index d51a0f60961..5b6dd7774ec 100644
--- a/stan/math/prim/functor/finite_diff_gradient_auto.hpp
+++ b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
@@ -4,14 +4,15 @@
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/finite_diff_stepsize.hpp>
 #include <stan/math/prim/fun/value_of_rec.hpp>
-#include <stan/math/prim/meta/conditional_copy_and_promote.hpp>
 #include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/functor/filter_map.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
 #include <stan/math/prim/functor/make_holder_tuple.hpp>
+#include <stan/math/prim/meta/contains_autodiff.hpp>
+#include <stan/math/prim/meta/filtered_tuple_indices.hpp>
 #include <cmath>
-#include <functional>
 #include <tuple>
+#include <type_traits>
 #include <utility>
 #include <vector>
 
@@ -25,51 +26,6 @@ namespace internal {
 template <typename...>
 struct dependent_false : std::false_type {};
 
-/**
- * Trait indicating whether a type contains at least one autodiff scalar.
- *
- * Scalars are checked directly, while container specializations recurse into
- * contained value types.
- */
-template <typename T>
-struct contains_autodiff_impl : bool_constant<is_autodiff_v<std::decay_t<T>>> {};
-
-template <typename... Args>
-struct contains_autodiff_impl<std::tuple<Args...>>
-    : bool_constant<(contains_autodiff_impl<std::decay_t<Args>>::value || ...)> {};
-
-template <typename T, typename... VecArgs>
-struct contains_autodiff_impl<std::vector<T, VecArgs...>>
-    : bool_constant<contains_autodiff_impl<std::decay_t<T>>::value> {};
-
-template <typename T>
-struct contains_autodiff : contains_autodiff_impl<std::decay_t<T>> {};
-
-template <typename T>
-inline constexpr bool contains_autodiff_v = contains_autodiff<T>::value;
-
-template <typename T>
-struct count_autodiff_args_impl;
-
-template <typename... Args>
-struct count_autodiff_args_impl<std::tuple<Args...>>
-    : std::integral_constant<
-          std::size_t,
-          (static_cast<std::size_t>(contains_autodiff_v<Args>) + ... + 0)> {};
-
-template <typename T>
-struct count_autodiff_args : count_autodiff_args_impl<std::decay_t<T>> {};
-
-template <typename TupleArgs>
-inline constexpr std::size_t count_autodiff_args_v
-    = count_autodiff_args<std::decay_t<TupleArgs>>::value;
-
-template <typename T, template <class...> class Cond, std::size_t pos>
-struct indices_of_if_impl {
-  static constexpr std::size_t value = (Cond<T>::value ? pos : -1);
-};
-
-
 /**
  * Filter a tuple and return a tuple containing references to the elements
  * whose types satisfy `contains_autodiff`.
@@ -78,7 +34,7 @@ struct indices_of_if_impl {
  *
  * @tparam T input tuple or nested tuple-like structure
  * @param t input to filter
- * @return tuple of references to autodiff-containing top-level arguments
+ * @return tuple of references to autodiff-containing entries
  */
 template <typename T>
 inline constexpr decltype(auto) filter_autodiff_types(T&& t) {
@@ -94,65 +50,13 @@ inline constexpr decltype(auto) filter_autodiff_types(T&& t) {
       std::forward<T>(t));
 }
 
-/**
- * Build a zero-initialized gradient object for a single autodiff-containing
- * leaf/container argument.
- *
- * @tparam Arg leaf/container type
- * @param arg argument value used for shape
- * @return zero-initialized gradient container matching `arg` shape
- */
-template <typename Arg>
-inline auto zeroed_container_leaf(const Arg& arg) {
-  using ArgDec = std::decay_t<Arg>;
-  using PartialType = partials_type_t<scalar_type_t<ArgDec>>;
-  if constexpr (is_std_vector_v<ArgDec>) {
-    using ElemType = std::decay_t<decltype(
-        zeroed_container_leaf(std::declval<const value_type_t<ArgDec>&>()))>;
-    std::vector<ElemType> ret;
-    ret.reserve(arg.size());
-    for (const auto& el : arg) {
-      ret.push_back(zeroed_container_leaf(el));
-    }
-    return ret;
-  } else if constexpr (is_eigen_v<ArgDec>) {
-    using eigen_t = promote_scalar_t<PartialType, plain_type_t<ArgDec>>;
-    return eigen_t::Zero(arg.rows(), arg.cols()).eval();
-  } else if constexpr (is_stan_scalar_v<ArgDec>) {
-    return PartialType(0);
-  } else {
-    static_assert(
-        dependent_false<ArgDec>::value,
-        "Unsupported gradient container in finite_diff_gradient_auto.");
-  }
+template <typename Tuple, std::size_t... Is>
+inline auto tuple_subset(Tuple&& tuple, std::index_sequence<Is...>) {
+  return make_holder_tuple(std::get<Is>(std::forward<Tuple>(tuple))...);
 }
 
 /**
- * Build a compact tuple of zero-initialized gradient containers for
- * autodiff-containing entries in `args`.
- *
- * Non-autodiff top-level entries are filtered out.
- *
- * @tparam TupleArgs tuple of input arguments
- * @param args tuple of input arguments
- * @return tuple of zeroed gradient containers for autodiff-containing entries
- */
-template <typename TupleArgs, require_tuple_t<TupleArgs>* = nullptr>
-inline auto zeroed_container(TupleArgs&& args) {
-  return stan::math::filter_map<contains_autodiff>(
-      [](auto&& arg) {
-        using Arg = std::decay_t<decltype(arg)>;
-        if constexpr (is_tuple_v<Arg>) {
-          return zeroed_container(std::forward<decltype(arg)>(arg));
-        } else {
-          return zeroed_container_leaf(arg);
-        }
-      },
-      std::forward<TupleArgs>(args));
-}
-
-/**
- * Extracts a double-valued coordinate for finite-difference step-size
+ * Extract a double-valued coordinate used for finite-difference step size
  * computation from scalar-like autodiff stacks.
  *
  * @tparam T scalar type
@@ -171,45 +75,72 @@ inline double finite_diff_stepsize_value(const T& x) {
 }
 
 /**
- * Recursively traverses `arg_mask`, `arg_work`, and `grad` with matching
- * structure and invokes `fn` on each autodiff scalar coordinate pair
- * `(arg_work_coord, grad_coord)`.
+ * Recursively traverses `arg_mask`, `arg_work`, and `grad` and invokes `fn` on
+ * each perturbable coordinate pair `(arg_work_coord, grad_coord)`.
  *
- * Non-autodiff scalar leaves in `arg_mask` are skipped.
+ * `arg_mask` controls whether a scalar coordinate is perturbable, while
+ * `arg_work` is the mutable coordinate that gets perturbed.
  *
- * @tparam ArgMask mask argument type (uses original scalar categories)
- * @tparam ArgWork mutable working argument type
- * @tparam Grad gradient container type
- * @tparam F callable with signature compatible with `(arg_work_coord, grad_coord)`
- * @param arg_mask mask argument used to decide perturbable coordinates
- * @param arg_work mutable working argument coordinates
+ * @tparam ArgMask traversal mask type
+ * @tparam ArgWork mutable argument type
+ * @tparam Grad gradient output type
+ * @tparam F callable compatible with `(arg_work_coord, grad_coord)`
+ * @param arg_mask mask selecting perturbable coordinates
+ * @param arg_work mutable argument coordinates
  * @param grad gradient output coordinates
- * @param fn callback applied at each perturbable coordinate
+ * @param fn callback applied for each perturbable coordinate
  */
-template <typename ArgMask, typename Grad, typename F>
-inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, Grad&& grad, F&& fn) {
+template <typename ArgMask, typename ArgWork, typename Grad, typename F>
+inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, ArgWork&& arg_work,
+                                         Grad&& grad, F&& fn) {
   using arg_mask_t = std::decay_t<ArgMask>;
   if constexpr (is_stan_scalar_v<arg_mask_t>) {
-      fn(arg_mask, grad);
+    if constexpr (contains_autodiff_v<arg_mask_t>) {
+      fn(arg_work, grad);
+    }
   } else if constexpr (is_eigen_v<arg_mask_t>) {
     for (Eigen::Index i = 0; i < arg_mask.size(); ++i) {
-      for_each_coordinate_pair_mut(arg_mask.coeffRef(i),
+      for_each_coordinate_pair_mut(arg_mask.coeffRef(i), arg_work.coeffRef(i),
                                    grad.coeffRef(i), fn);
     }
   } else if constexpr (is_std_vector_v<arg_mask_t>) {
     for (std::size_t i = 0; i < arg_mask.size(); ++i) {
-      for_each_coordinate_pair_mut(arg_mask[i], grad[i], fn);
+      for_each_coordinate_pair_mut(arg_mask[i], arg_work[i], grad[i], fn);
     }
   } else if constexpr (is_tuple_v<arg_mask_t>) {
-    static_assert(std::tuple_size<std::decay_t<ArgMask>>::value
-                      == std::tuple_size<std::decay_t<Grad>>::value,
-                  "Tuple size mismatch in finite_diff_gradient_auto traversal");
-    stan::math::for_each(
-        [&fn](auto&& arg_mask_i, auto&& grad_i) {
-          for_each_coordinate_pair_mut(arg_mask_i, grad_i, fn);
-        },
-        std::forward<ArgMask>(arg_mask),
-        std::forward<Grad>(grad));
+    constexpr auto arg_mask_size = std::tuple_size<std::decay_t<ArgMask>>::value;
+    constexpr auto arg_work_size = std::tuple_size<std::decay_t<ArgWork>>::value;
+    static_assert(arg_mask_size == arg_work_size,
+                  "Tuple size mismatch in finite_diff_gradient_auto traversal "
+                  "between arg_mask and arg_work.");
+    constexpr auto grad_size = std::tuple_size<std::decay_t<Grad>>::value;
+    if constexpr (arg_mask_size == grad_size) {
+      stan::math::for_each(
+          [&fn](auto&& arg_mask_i, auto&& arg_work_i, auto&& grad_i) {
+            for_each_coordinate_pair_mut(arg_mask_i, arg_work_i, grad_i, fn);
+          },
+          std::forward<ArgMask>(arg_mask), std::forward<ArgWork>(arg_work),
+          std::forward<Grad>(grad));
+    } else {
+      using selected_idxs_t
+          = filtered_tuple_indices_t<contains_autodiff, std::decay_t<ArgMask>>;
+      auto arg_mask_ad
+          = tuple_subset(std::forward<ArgMask>(arg_mask), selected_idxs_t{});
+      auto arg_work_ad
+          = tuple_subset(std::forward<ArgWork>(arg_work), selected_idxs_t{});
+      static_assert(
+          std::tuple_size<std::decay_t<decltype(arg_mask_ad)>>::value
+                  == std::tuple_size<std::decay_t<Grad>>::value
+              && std::tuple_size<std::decay_t<decltype(arg_work_ad)>>::value
+                     == std::tuple_size<std::decay_t<Grad>>::value,
+          "Tuple size mismatch in finite_diff_gradient_auto traversal");
+      stan::math::for_each(
+          [&fn](auto&& arg_mask_i, auto&& arg_work_i, auto&& grad_i) {
+            for_each_coordinate_pair_mut(arg_mask_i, arg_work_i, grad_i, fn);
+          },
+          std::move(arg_mask_ad), std::move(arg_work_ad),
+          std::forward<Grad>(grad));
+    }
   } else {
     static_assert(dependent_false<arg_mask_t>::value,
                   "Unsupported container in finite_diff_gradient_auto.");
@@ -217,35 +148,35 @@ inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, Grad&& grad, F&& fn
 }
 
 /**
- * Applies the six-point finite-difference stencil to every perturbable
- * coordinate in one autodiff top-level argument and writes the corresponding
- * gradient coordinate.
+ * Apply the six-point finite-difference stencil to every perturbable
+ * coordinate in one top-level argument and write gradient coordinates.
  *
- * Coordinates are always restored to their original value after each stencil.
+ * Coordinates in `arg_work` are restored to their original value after each
+ * stencil evaluation.
  *
  * @tparam ScalarT scalar type of function value and gradients
  * @tparam F callable type
- * @tparam ArgMask mask argument type (original categories)
+ * @tparam ArgMask mask argument type
  * @tparam ArgWork mutable working argument type
  * @tparam NthGrad gradient container type for this top-level argument
- * @tparam TupleArgsWork mutable tuple of all working arguments
+ * @tparam TupleArgsWork mutable tuple of all working arguments passed to `f`
  * @param f function being finite-differenced
  * @param arg_mask mask argument used to choose perturbable coordinates
  * @param arg_work mutable working argument coordinates
  * @param grad gradient output container for this top-level argument
  * @param args_work mutable tuple of all working arguments passed to `f`
  */
-template <typename ScalarT, typename F, typename ArgMask,
+template <typename ScalarT, typename F, typename ArgMask, typename ArgWork,
           typename NthGrad, typename TupleArgsWork>
 inline void finite_diff_gradient_auto_impl(const F& f, ArgMask&& arg_mask,
-                                           NthGrad&& grad,
+                                           ArgWork&& arg_work, NthGrad&& grad,
                                            TupleArgsWork&& args_work) {
   static constexpr int h_scale[6] = {3, 2, 1, -3, -2, -1};
   static constexpr int mults[6] = {1, -9, 45, -1, 9, -45};
 
   for_each_coordinate_pair_mut(
-      std::forward<ArgMask>(arg_mask), std::forward<NthGrad>(grad),
-      [&f, &args_work](auto& arg_coord, auto& grad_coord) {
+      std::forward<ArgMask>(arg_mask), std::forward<ArgWork>(arg_work),
+      std::forward<NthGrad>(grad), [&f, &args_work](auto& arg_coord, auto& grad_coord) {
         const auto orig = arg_coord;
         const double h = finite_diff_stepsize(finite_diff_stepsize_value(orig));
         ScalarT delta_f = 0;
@@ -260,6 +191,7 @@ inline void finite_diff_gradient_auto_impl(const F& f, ArgMask&& arg_mask,
         grad_coord = delta_f / (60 * h);
       });
 }
+
 }  // namespace internal
 
 /**
@@ -317,10 +249,9 @@ inline void finite_diff_gradient_auto(const F& f, VectorT&& x, ScalarT& fx,
     double h = finite_diff_stepsize(value_of_rec(x[i]));
     ScalarT delta_f = 0;
     for (int j = 0; j < 6; ++j) {
-      auto x_temp
-          = EigT::NullaryExpr(x.size(), [&x, &i, &h, &j](Eigen::Index k) {
-              return k == i ? x[i] + h * h_scale[j] : x[k];
-            });
+      auto x_temp = EigT::NullaryExpr(x.size(), [&x, &i, &h, &j](Eigen::Index k) {
+        return k == i ? x[i] + h * h_scale[j] : x[k];
+      });
       delta_f += f(std::move(x_temp)) * mults[j];
     }
     return delta_f / (60 * h);
@@ -328,122 +259,82 @@ inline void finite_diff_gradient_auto(const F& f, VectorT&& x, ScalarT& fx,
 }
 
 /**
- * Filter a tuple and return a tuple with references to the types with an autodiff
+ * Filter a tuple and return a tuple with references to entries with autodiff
  * scalar type.
+ *
  * @tparam T Possibly a tuple, std::vector, Eigen type, or scalar
  * @param[in] t Input to filter
- * @return Filtered input with only var scalar types
+ * @return Filtered input with only autodiff scalar entries
  */
 template <typename T>
 inline constexpr decltype(auto) filter_ad_scalar_types(T&& t) {
-  return stan::math::filter_map<internal::contains_autodiff>(
-      [](auto&& arg) -> decltype(auto) {
-        if constexpr (is_tuple_v<std::decay_t<decltype(arg)>>) {
-          return filter_ad_scalar_types(std::forward<decltype(arg)>(arg));
-        } else {
-          return std::forward<decltype(arg)>(arg);
-        }
-      },
-      std::forward<T>(t));
+  return internal::filter_autodiff_types(std::forward<T>(t));
 }
 
 /**
  * Calculate the function value and gradients for heterogeneous tuple arguments
- * while preserving non-autodiff argument scalar types at the call boundary.
- *
- * Pattern used by this overload:
- * 1. Build `args_ad` by filtering top-level autodiff-containing entries from
- *    the original `args` tuple.
- * 2. Build one mutable working tuple `args_work` where autodiff-containing
- *    entries are deep-copied and promoted to `ScalarT`, and non-autodiff
- *    entries are forwarded unchanged.
- * 3. Build `args_work_ad` by applying the same top-level filter shape as
- *    `args_ad` (via the `TupleArgs` mask) so both filtered tuples align.
- * 4. Zip `args_ad`, `args_work_ad`, and `grads` and run
- *    `finite_diff_gradient_auto_impl` for each top-level autodiff entry.
+ * while preserving argument scalar types at the call boundary.
  *
- * The split between original mask (`args_ad`) and mutable working values
- * (`args_work_ad`) ensures coordinate traversal only perturbs locations that
- * correspond to autodiff coordinates in the original argument pack, while all
- * function evaluations run against the promoted mutable working tuple.
+ * This overload expects `grads` to mirror the top-level structure of `args`.
+ * Non-autodiff top-level entries are skipped.
  *
  * @tparam F callable type
  * @tparam ScalarT scalar type used for function value and gradient storage
  * @tparam TupleArgs tuple of input arguments
- * @tparam TupleGrads tuple of gradient containers aligned to autodiff
- *         top-level arguments
+ * @tparam TupleGrads tuple of gradient containers
  * @param[in] f function to differentiate
  * @param[out] fx function value at `args`
  * @param[in] args tuple of function arguments
- * @param[out] grads compact tuple of gradient containers
+ * @param[out] grads tuple of gradient containers
  */
 template <typename F, typename ScalarT, typename TupleArgs, typename TupleGrads,
           require_tuple_t<TupleArgs>* = nullptr,
-          require_tuple_t<TupleGrads>* = nullptr,
-          require_t<bool_constant<
-              (internal::count_autodiff_args_v<TupleArgs> > 0)>>* = nullptr>
+          require_tuple_t<TupleGrads>* = nullptr>
 inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
                                       TupleGrads&& grads) {
-  fx = stan::math::apply([&f](auto&&... args_i) { return f(args_i...); },
-                         args);
-  stan::math::for_each([&f, &args, &grads](auto&& arg, auto&& grad) {
-    if constexpr (!std::is_integral_v<scalar_type_t<decltype(grad)>> && !std::is_integral_v<scalar_type_t<decltype(arg)>>) {
-      internal::finite_diff_gradient_auto_impl<ScalarT>(
-          f, arg, grad, args);
-    }
-  }, args, grads);
- }
-
+  fx = stan::math::apply([&f](auto&&... args_i) { return f(args_i...); }, args);
+  stan::math::for_each(
+      [&f, &args](auto&& arg_i, auto&& grad_i) {
+        if constexpr (contains_autodiff_v<std::decay_t<decltype(arg_i)>>) {
+          internal::finite_diff_gradient_auto_impl<ScalarT>(f, arg_i, arg_i,
+                                                            grad_i, args);
+        }
+      },
+      args, grads);
+}
 
 /**
- * Calculate the function value and gradients for heterogeneous tuple arguments
- * while preserving non-autodiff argument scalar types at the call boundary.
- *
- * Pattern used by this overload:
- * 1. Build `args_ad` by filtering top-level autodiff-containing entries from
- *    the original `args` tuple.
- * 2. Build one mutable working tuple `args_work` where autodiff-containing
- *    entries are deep-copied and promoted to `ScalarT`, and non-autodiff
- *    entries are forwarded unchanged.
- * 3. Build `args_work_ad` by applying the same top-level filter shape as
- *    `args_ad` (via the `TupleArgs` mask) so both filtered tuples align.
- * 4. Zip `args_ad`, `args_work_ad`, and `grads` and run
- *    `finite_diff_gradient_auto_impl` for each top-level autodiff entry.
- *
- * The split between original mask (`args_ad`) and mutable working values
- * (`args_work_ad`) ensures coordinate traversal only perturbs locations that
- * correspond to autodiff coordinates in the original argument pack, while all
- * function evaluations run against the promoted mutable working tuple.
+ * Indexed tuple overload used for compact gradient tuples. `Idxs` select
+ * autodiff-containing top-level entries from both `args_mask` and `args_work`,
+ * while `GradIdxs` selects matching gradient entries from `grads`.
  *
  * @tparam F callable type
- * @tparam ScalarT scalar type used for function value and gradient storage
- * @tparam TupleArgs tuple of input arguments
- * @tparam TupleGrads tuple of gradient containers aligned to autodiff
- *         top-level arguments
- * @param[in] f function to differentiate
- * @param[out] fx function value at `args`
- * @param[in] args tuple of function arguments
- * @param[out] grads compact tuple of gradient containers
+ * @tparam ScalarT scalar type used for function value and gradients
+ * @tparam TupleArgsMask tuple used only as perturbation mask
+ * @tparam TupleArgsWork mutable tuple used for function evaluation
+ * @tparam TupleGrads compact tuple of gradient containers
  */
-template <typename F, typename ScalarT, typename TupleArgs, typename TupleGrads,
-          require_tuple_t<TupleArgs>* = nullptr,
-          require_tuple_t<TupleGrads>* = nullptr,
-          std::size_t... Idxs,
-          std::size_t... GradIdxs>
-inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
-                                      TupleGrads&& grads,
-                                      std::index_sequence<Idxs...> /* idxs */,
-                                      std::index_sequence<GradIdxs...> /* grad_idxs */) {
+template <typename F, typename ScalarT, typename TupleArgsMask,
+          typename TupleArgsWork, typename TupleGrads, std::size_t... Idxs,
+          std::size_t... GradIdxs, require_tuple_t<TupleArgsMask>* = nullptr,
+          require_tuple_t<TupleArgsWork>* = nullptr,
+          require_tuple_t<TupleGrads>* = nullptr>
+inline void finite_diff_gradient_auto(
+    const F& f, ScalarT& fx, TupleArgsMask&& args_mask, TupleArgsWork&& args_work,
+    TupleGrads&& grads, std::index_sequence<Idxs...> /* idxs */,
+    std::index_sequence<GradIdxs...> /* grad_idxs */) {
+  static_assert(sizeof...(Idxs) == sizeof...(GradIdxs),
+                "Index sequence size mismatch in finite_diff_gradient_auto.");
   fx = stan::math::apply([&f](auto&&... args_i) { return f(args_i...); },
-                         args);
-   using Swallow = int[];
-  static_cast<void>(Swallow{(
-    static_cast<void>(internal::finite_diff_gradient_auto_impl<ScalarT>(
-        f, std::get<Idxs>(args),
-        std::get<GradIdxs>(grads), args)), 0)...});
-
+                         args_work);
+  using Swallow = int[];
+  static_cast<void>(Swallow{
+      (static_cast<void>(internal::finite_diff_gradient_auto_impl<ScalarT>(
+           f, std::get<Idxs>(args_mask), std::get<Idxs>(args_work),
+           std::get<GradIdxs>(grads), args_work)),
+       0)...});
+}
 
- }
 }  // namespace math
 }  // namespace stan
 #endif
diff --git a/stan/math/prim/meta.hpp b/stan/math/prim/meta.hpp
index af2966e2351..bdb692c46bc 100644
--- a/stan/math/prim/meta.hpp
+++ b/stan/math/prim/meta.hpp
@@ -71,9 +71,12 @@
 #include <stan/math/prim/meta/ad_promotable.hpp>
 #include <stan/math/prim/meta/append_return_type.hpp>
 #include <stan/math/prim/meta/base_type.hpp>
+#include <stan/math/prim/meta/contains_autodiff.hpp>
 #include <stan/math/prim/meta/contains_std_vector.hpp>
 #include <stan/math/prim/meta/contains_tuple.hpp>
 #include <stan/math/prim/meta/error_index.hpp>
+#include <stan/math/prim/meta/filter_index_sequence.hpp>
+#include <stan/math/prim/meta/filtered_tuple_indices.hpp>
 #include <stan/math/prim/meta/forward_as.hpp>
 #include <stan/math/prim/meta/holder.hpp>
 #include <stan/math/prim/meta/include_summand.hpp>
diff --git a/stan/math/prim/meta/contains_autodiff.hpp b/stan/math/prim/meta/contains_autodiff.hpp
new file mode 100644
index 00000000000..40b3447caa1
--- /dev/null
+++ b/stan/math/prim/meta/contains_autodiff.hpp
@@ -0,0 +1,37 @@
+#ifndef STAN_MATH_PRIM_META_CONTAINS_AUTODIFF_HPP
+#define STAN_MATH_PRIM_META_CONTAINS_AUTODIFF_HPP
+
+#include <stan/math/prim/meta/bool_constant.hpp>
+#include <stan/math/prim/meta/is_autodiff.hpp>
+#include <tuple>
+#include <type_traits>
+#include <vector>
+
+namespace stan {
+
+/**
+ * Trait indicating whether a type contains at least one autodiff scalar.
+ *
+ * Scalars are checked directly. Tuple and std::vector specializations recurse
+ * into their contained element types.
+ *
+ * @tparam T type to inspect
+ */
+template <typename T>
+struct contains_autodiff : bool_constant<is_autodiff_v<std::decay_t<T>>> {};
+
+template <typename... Args>
+struct contains_autodiff<std::tuple<Args...>>
+    : bool_constant<(contains_autodiff<std::decay_t<Args>>::value || ...)> {};
+
+template <typename T, typename... VecArgs>
+struct contains_autodiff<std::vector<T, VecArgs...>>
+    : bool_constant<contains_autodiff<std::decay_t<T>>::value> {};
+
+template <typename T>
+inline constexpr bool contains_autodiff_v
+    = contains_autodiff<std::decay_t<T>>::value;
+
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/meta/filter_index_sequence.hpp b/stan/math/prim/meta/filter_index_sequence.hpp
new file mode 100644
index 00000000000..569f406c29b
--- /dev/null
+++ b/stan/math/prim/meta/filter_index_sequence.hpp
@@ -0,0 +1,76 @@
+#ifndef STAN_MATH_PRIM_META_FILTER_INDEX_SEQUENCE_HPP
+#define STAN_MATH_PRIM_META_FILTER_INDEX_SEQUENCE_HPP
+
+#include <utility>
+
+namespace stan {
+namespace math {
+
+/**
+ * Concatenates integer sequences.
+ */
+template <typename SeqT, SeqT... Is, SeqT... Js>
+inline constexpr std::integer_sequence<SeqT, Is..., Js...>
+concat_integer_sequences(std::integer_sequence<SeqT, Is...>,
+                         std::integer_sequence<SeqT, Js...>) {
+  return {};
+}
+
+template <typename SeqT>
+inline constexpr std::integer_sequence<SeqT> concat_integer_sequences(
+    std::integer_sequence<SeqT>) {
+  return {};
+}
+
+template <typename SeqT, SeqT... Is, SeqT... Js, class... Rest>
+inline constexpr auto concat_integer_sequences(
+    std::integer_sequence<SeqT, Is...>, std::integer_sequence<SeqT, Js...>,
+    Rest...) {
+  return concat_integer_sequences(std::integer_sequence<SeqT, Is..., Js...>{},
+                                  Rest{}...);
+}
+
+/**
+ * Returns a singleton index sequence if `Predicate<Arg>::value` is true,
+ * otherwise an empty sequence.
+ */
+template <template <class...> class Predicate, typename Arg, class SeqT,
+          SeqT I>
+inline constexpr auto filter_single_index_sequence(
+    std::integer_sequence<SeqT, I>) {
+  if constexpr (Predicate<Arg>::value) {
+    return std::integer_sequence<SeqT, I>{};
+  } else {
+    return std::integer_sequence<SeqT>{};
+  }
+}
+
+/**
+ * Filters an index sequence by applying `Predicate` to matching `Args...`.
+ *
+ * `Args...` and the index sequence must have the same length.
+ */
+template <template <class...> class Predicate, typename... Args, class SeqT,
+          SeqT... Is>
+inline constexpr auto filter_index_sequence(std::integer_sequence<SeqT, Is...>) {
+  static_assert(sizeof...(Args) == sizeof...(Is),
+                "filter_index_sequence requires the same number of argument "
+                "types and indices.");
+  if constexpr (sizeof...(Is) > 0) {
+    return concat_integer_sequences(
+        filter_single_index_sequence<Predicate, Args>(
+            std::integer_sequence<SeqT, Is>{})...);
+  } else {
+    return std::integer_sequence<SeqT>{};
+  }
+}
+
+template <template <class...> class Predicate, typename Sequence,
+          typename... Args>
+using filter_index_sequence_t
+    = decltype(filter_index_sequence<Predicate, Args...>(Sequence{}));
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/meta/filtered_tuple_indices.hpp b/stan/math/prim/meta/filtered_tuple_indices.hpp
new file mode 100644
index 00000000000..7b6e57dbdd4
--- /dev/null
+++ b/stan/math/prim/meta/filtered_tuple_indices.hpp
@@ -0,0 +1,47 @@
+#ifndef STAN_MATH_PRIM_META_FILTERED_TUPLE_INDICES_HPP
+#define STAN_MATH_PRIM_META_FILTERED_TUPLE_INDICES_HPP
+
+#include <stan/math/prim/meta/filter_index_sequence.hpp>
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+namespace stan {
+namespace math {
+
+template <typename Sequence>
+struct index_sequence_size;
+
+template <typename T, T... Is>
+struct index_sequence_size<std::integer_sequence<T, Is...>>
+    : std::integral_constant<std::size_t, sizeof...(Is)> {};
+
+/**
+ * Produces a compact index sequence of top-level tuple elements whose
+ * types satisfy `Predicate`.
+ */
+template <template <class...> class Predicate, typename Tuple>
+struct filtered_tuple_indices;
+
+template <template <class...> class Predicate, typename... Args>
+struct filtered_tuple_indices<Predicate, std::tuple<Args...>> {
+  using type = filter_index_sequence_t<Predicate, std::index_sequence_for<Args...>,
+                                       std::decay_t<Args>...>;
+};
+
+template <template <class...> class Predicate, typename Tuple>
+using filtered_tuple_indices_t
+    = typename filtered_tuple_indices<Predicate, std::decay_t<Tuple>>::type;
+
+template <template <class...> class Predicate, typename Tuple>
+using filtered_tuple_size
+    = index_sequence_size<filtered_tuple_indices_t<Predicate, Tuple>>;
+
+template <template <class...> class Predicate, typename Tuple>
+inline constexpr std::size_t filtered_tuple_size_v
+    = filtered_tuple_size<Predicate, Tuple>::value;
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/test/unit/math/prim/fun/zeroed_container_test.cpp b/test/unit/math/prim/fun/zeroed_container_test.cpp
new file mode 100644
index 00000000000..e8c7d16e165
--- /dev/null
+++ b/test/unit/math/prim/fun/zeroed_container_test.cpp
@@ -0,0 +1,56 @@
+#include <stan/math/fwd/core/fvar.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/fun/zeroed_container.hpp>
+#include <stan/math/prim/meta/contains_autodiff.hpp>
+#include <gtest/gtest.h>
+#include <tuple>
+#include <type_traits>
+#include <vector>
+
+TEST(MathPrimFunZeroedContainer, autodiff_predicate_filters_and_zeros) {
+  using stan::math::fvar;
+  std::vector<fvar<double>> a{
+      fvar<double>(1.3, 0.0),
+      fvar<double>(-0.7, 0.0),
+  };
+  Eigen::Matrix<fvar<double>, -1, 1> b(2);
+  b << fvar<double>(2.0, 0.0), fvar<double>(-3.0, 0.0);
+  std::vector<int> x_i{5, 6};
+  Eigen::VectorXd c{{4.0}};
+
+  auto args = std::make_tuple(a, std::make_tuple(b, x_i), c);
+  auto grads = stan::math::zeroed_container<stan::contains_autodiff>(args);
+
+  using grads_t = std::decay_t<decltype(grads)>;
+  EXPECT_EQ(2, std::tuple_size<grads_t>::value);
+
+  const auto& a_grads = std::get<0>(grads);
+  EXPECT_EQ(a.size(), a_grads.size());
+  for (std::size_t i = 0; i < a_grads.size(); ++i) {
+    EXPECT_FLOAT_EQ(0.0, a_grads[i]);
+  }
+
+  const auto& b_grads = std::get<0>(std::get<1>(grads));
+  EXPECT_EQ(b.size(), b_grads.size());
+  for (Eigen::Index i = 0; i < b_grads.size(); ++i) {
+    EXPECT_FLOAT_EQ(0.0, b_grads(i));
+  }
+}
+
+TEST(MathPrimFunZeroedContainer, integral_predicate_filters_and_zeros) {
+  auto args = std::make_tuple(std::vector<int>{1, 2, 3}, 1.5,
+                              std::make_tuple(2.5, 7));
+  auto zeros = stan::math::zeroed_container<std::is_integral>(args);
+
+  using zeros_t = std::decay_t<decltype(zeros)>;
+  EXPECT_EQ(2, std::tuple_size<zeros_t>::value);
+
+  const auto& vec_zero = std::get<0>(zeros);
+  EXPECT_EQ(3, vec_zero.size());
+  EXPECT_EQ(0, vec_zero[0]);
+  EXPECT_EQ(0, vec_zero[1]);
+  EXPECT_EQ(0, vec_zero[2]);
+
+  const auto& nested_zero = std::get<0>(std::get<1>(zeros));
+  EXPECT_EQ(0, nested_zero);
+}
diff --git a/test/unit/math/prim/meta/contains_autodiff_test.cpp b/test/unit/math/prim/meta/contains_autodiff_test.cpp
new file mode 100644
index 00000000000..e7206302557
--- /dev/null
+++ b/test/unit/math/prim/meta/contains_autodiff_test.cpp
@@ -0,0 +1,31 @@
+#include <stan/math/fwd/core/fvar.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/meta/contains_autodiff.hpp>
+#include <gtest/gtest.h>
+#include <tuple>
+#include <vector>
+
+TEST(MathMetaContainsAutodiff, scalar_detection) {
+  using stan::math::fvar;
+  EXPECT_FALSE(stan::contains_autodiff_v<double>);
+  EXPECT_FALSE(stan::contains_autodiff_v<int>);
+  EXPECT_TRUE(stan::contains_autodiff_v<fvar<double>>);
+}
+
+TEST(MathMetaContainsAutodiff, container_detection) {
+  using stan::math::fvar;
+  EXPECT_FALSE((stan::contains_autodiff_v<std::vector<double>>));
+  EXPECT_TRUE((stan::contains_autodiff_v<std::vector<fvar<double>>>));
+  EXPECT_FALSE(
+      (stan::contains_autodiff_v<std::tuple<int, std::vector<double>>>));
+  EXPECT_TRUE((stan::contains_autodiff_v<
+               std::tuple<int, std::vector<fvar<double>>>>));
+}
+
+TEST(MathMetaContainsAutodiff, nested_tuple_detection) {
+  using stan::math::fvar;
+  using nested_ad_t = std::tuple<double, std::tuple<int, fvar<double>>>;
+  using nested_non_ad_t = std::tuple<double, std::tuple<int, double>>;
+  EXPECT_TRUE((stan::contains_autodiff_v<nested_ad_t>));
+  EXPECT_FALSE((stan::contains_autodiff_v<nested_non_ad_t>));
+}
diff --git a/test/unit/math/prim/meta/filter_index_sequence_test.cpp b/test/unit/math/prim/meta/filter_index_sequence_test.cpp
new file mode 100644
index 00000000000..8680bfcb520
--- /dev/null
+++ b/test/unit/math/prim/meta/filter_index_sequence_test.cpp
@@ -0,0 +1,35 @@
+#include <stan/math/fwd/core/fvar.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/meta/contains_autodiff.hpp>
+#include <stan/math/prim/meta/filter_index_sequence.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+TEST(MathMetaFilterIndexSequence, integral_filtering) {
+  using in_seq_t = std::index_sequence<0, 1, 2, 3>;
+  using expected_t = std::index_sequence<0, 2>;
+  using filtered_t = stan::math::filter_index_sequence_t<
+      std::is_integral, in_seq_t, int, double, long, float>;
+  EXPECT_SAME_TYPE(expected_t, filtered_t);
+}
+
+TEST(MathMetaFilterIndexSequence, empty_filtering) {
+  using in_seq_t = std::index_sequence<0, 1, 2>;
+  using expected_t = std::index_sequence<>;
+  using filtered_t = stan::math::filter_index_sequence_t<
+      std::is_integral, in_seq_t, double, float, long double>;
+  EXPECT_SAME_TYPE(expected_t, filtered_t);
+}
+
+TEST(MathMetaFilterIndexSequence, autodiff_filtering) {
+  using stan::math::fvar;
+  using in_seq_t = std::index_sequence<0, 1, 2>;
+  using expected_t = std::index_sequence<1>;
+  using filtered_t = stan::math::filter_index_sequence_t<
+      stan::contains_autodiff, in_seq_t, double, fvar<double>,
+      std::vector<int>>;
+  EXPECT_SAME_TYPE(expected_t, filtered_t);
+}
diff --git a/test/unit/math/prim/meta/filtered_tuple_indices_test.cpp b/test/unit/math/prim/meta/filtered_tuple_indices_test.cpp
new file mode 100644
index 00000000000..9776e142f0e
--- /dev/null
+++ b/test/unit/math/prim/meta/filtered_tuple_indices_test.cpp
@@ -0,0 +1,37 @@
+#include <stan/math/fwd/core/fvar.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/meta/contains_autodiff.hpp>
+#include <stan/math/prim/meta/filtered_tuple_indices.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+#include <tuple>
+#include <type_traits>
+#include <vector>
+
+TEST(MathMetaFilteredTupleIndices, autodiff_indices) {
+  using stan::math::fvar;
+  using tuple_t
+      = std::tuple<int, fvar<double>, std::vector<double>,
+                   std::tuple<fvar<double>, double>>;
+  using expected_t = std::index_sequence<1, 3>;
+  using idxs_t = stan::math::filtered_tuple_indices_t<stan::contains_autodiff,
+                                                      tuple_t>;
+  EXPECT_SAME_TYPE(expected_t, idxs_t);
+}
+
+TEST(MathMetaFilteredTupleIndices, integral_indices_and_size) {
+  using tuple_t = std::tuple<double, int, float, long>;
+  using expected_t = std::index_sequence<1, 3>;
+  using idxs_t
+      = stan::math::filtered_tuple_indices_t<std::is_integral, tuple_t>;
+  EXPECT_SAME_TYPE(expected_t, idxs_t);
+  EXPECT_EQ(2, (stan::math::filtered_tuple_size_v<std::is_integral, tuple_t>));
+}
+
+TEST(MathMetaFilteredTupleIndices, tuple_decay) {
+  using tuple_ref_t = const std::tuple<const int, double, const long>&;
+  using expected_t = std::index_sequence<0, 2>;
+  using idxs_t
+      = stan::math::filtered_tuple_indices_t<std::is_integral, tuple_ref_t>;
+  EXPECT_SAME_TYPE(expected_t, idxs_t);
+}

From 73172c16e90de3dce4b359b6659fb65d2094841c Mon Sep 17 00:00:00 2001
From: Steve Bronder <sbronder@flatironinstitute.org>
Date: Fri, 20 Feb 2026 14:50:33 -0500
Subject: [PATCH 5/6] clean up docs and tests

---
 stan/math/fwd/functor/finite_diff.hpp         |  53 +++++--
 stan/math/prim/fun.hpp                        |   2 +-
 ...ontainer.hpp => zeroed_filtered_tuple.hpp} |  38 +++--
 stan/math/prim/functor.hpp                    |   1 +
 stan/math/prim/functor/filter_types.hpp       |  38 +++++
 .../functor/finite_diff_gradient_auto.hpp     | 135 ++++++++++-------
 .../meta/conditional_copy_and_promote.hpp     |  12 +-
 stan/math/prim/meta/contains_autodiff.hpp     |  29 ++++
 stan/math/prim/meta/filter_index_sequence.hpp |  74 ++++++++-
 .../math/prim/meta/filtered_tuple_indices.hpp |  55 +++++++
 .../math/rev/core/filter_var_scalar_types.hpp |  13 +-
 .../math/{fwd => mix}/fun/tangent_of_test.cpp |   0
 .../finite_diff_gradient_auto_test.cpp        |  53 +++++++
 .../{fwd => mix}/functor/finite_diff_test.cpp |  49 +++++-
 .../meta/contains_autodiff_test.cpp           |  20 +++
 .../meta/filtered_tuple_indices_test.cpp      |   0
 .../math/prim/fun/zeroed_container_test.cpp   |  56 -------
 .../prim/fun/zeroed_filtered_tuple_test.cpp   | 143 ++++++++++++++++++
 .../math/prim/functor/filter_types_test.cpp   |  41 +++++
 19 files changed, 645 insertions(+), 167 deletions(-)
 rename stan/math/prim/fun/{zeroed_container.hpp => zeroed_filtered_tuple.hpp} (59%)
 create mode 100644 stan/math/prim/functor/filter_types.hpp
 rename test/unit/math/{fwd => mix}/fun/tangent_of_test.cpp (100%)
 rename test/unit/math/{fwd => mix}/functor/finite_diff_test.cpp (51%)
 rename test/unit/math/{prim => mix}/meta/contains_autodiff_test.cpp (62%)
 rename test/unit/math/{prim => mix}/meta/filtered_tuple_indices_test.cpp (100%)
 delete mode 100644 test/unit/math/prim/fun/zeroed_container_test.cpp
 create mode 100644 test/unit/math/prim/fun/zeroed_filtered_tuple_test.cpp
 create mode 100644 test/unit/math/prim/functor/filter_types_test.cpp

diff --git a/stan/math/fwd/functor/finite_diff.hpp b/stan/math/fwd/functor/finite_diff.hpp
index 483c82324c0..939708b6d65 100644
--- a/stan/math/fwd/functor/finite_diff.hpp
+++ b/stan/math/fwd/functor/finite_diff.hpp
@@ -2,12 +2,13 @@
 #define STAN_MATH_FWD_FUNCTOR_FINITE_DIFF_HPP
 
 #include <stan/math/fwd/fun/tangent_of.hpp>
-#include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun/value_of.hpp>
+#include <stan/math/fwd/meta.hpp>
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
-#include <stan/math/prim/fun/zeroed_container.hpp>
+#include <stan/math/prim/fun/zeroed_filtered_tuple.hpp>
 #include <stan/math/prim/functor/finite_diff_gradient_auto.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
+#include <stan/math/prim/functor/filter_types.hpp>
 #include <stan/math/prim/meta/contains_autodiff.hpp>
 #include <stan/math/prim/meta/filtered_tuple_indices.hpp>
 #include <tuple>
@@ -23,15 +24,18 @@ namespace internal {
  *
  * @tparam FuncTangent Type of tangent calculated by finite-differences
  * @tparam InputArg Type of the function input argument
- * @param tangent Calculated tangent
- * @param arg Input argument
+ * @param[in] tangent Calculated tangent
+ * @param[in] arg Input argument
+ * @return `0.0`, because non-fvar inputs do not contribute to the tangent
+ * aggregation.
+ * @throw None.
  */
 template <typename FuncTangent, typename InputArg,
           require_not_st_fvar<InputArg>* = nullptr>
 inline constexpr double aggregate_tangent(const FuncTangent& tangent,
                                           const InputArg& arg) {
   return 0;
-}  // namespace internal
+}
 
 /**
  * Helper function for aggregating tangents if the respective input argument
@@ -42,8 +46,10 @@ inline constexpr double aggregate_tangent(const FuncTangent& tangent,
  *
  * @tparam FuncTangent Type of tangent calculated by finite-differences
  * @tparam InputArg Type of the function input argument
- * @param tangent Calculated tangent
- * @param arg Input argument
+ * @param[in] tangent Calculated tangent
+ * @param[in] arg Input argument
+ * @return dot product of finite-difference tangent and input tangent.
+ * @throw None.
  */
 template <typename FuncTangent, typename InputArg,
           require_st_fvar<InputArg>* = nullptr>
@@ -63,7 +69,8 @@ inline auto aggregate_tangent(FuncTangent&& tangent, InputArg&& arg) {
     return rtn;
   }
 }
-}
+
+}  // namespace internal
 
 /**
  * Construct an fvar<T> where the tangent is calculated by finite-differencing.
@@ -77,23 +84,35 @@ inline auto aggregate_tangent(FuncTangent&& tangent, InputArg&& arg) {
  * @tparam TArgs Template parameter pack of the types passed in the `operator()`
  *                of the functor type `F`. Must contain at least on type whose
  *                scalar type is `fvar<T>`
+ *
+ * Internal pattern used by this overload:
+ * 1. Build `autodiff_args` with `filter_types<contains_autodiff>(args...)`.
+ * 2. Build compact zero-initialized `grads` with
+ *    `zeroed_filtered_tuple<contains_autodiff>(args...)`.
+ * 3. Use filtered top-level index sequences to dispatch tuple-native finite
+ *    differencing.
+ * 4. Aggregate tangent contributions by zipping compact `grads` and
+ *    `autodiff_args`.
+ *
  * @param func Functor for which fvar<T> support is needed
  * @param args Parameter pack of arguments to be passed to functor.
+ * @return Function value and tangent packed in the resulting `fvar`.
+ * @throw Any exception thrown by `func`.
  */
 template <typename F, typename... TArgs,
           require_any_st_fvar<TArgs...>* = nullptr>
 inline auto finite_diff(const F& func, const TArgs&... args) {
   using FvarT = return_type_t<TArgs...>;
   using FvarInnerT = typename FvarT::Scalar;
-  auto autodiff_args = filter_ad_scalar_types(std::forward_as_tuple(args...));
-  FvarInnerT rtn_value;
-  auto grads = stan::math::zeroed_container<contains_autodiff>(
-      std::forward_as_tuple(args...));
-  using args_tuple_t = std::tuple<TArgs...>;
+  auto autodiff_args
+      = stan::math::filter_types<contains_autodiff>(std::forward_as_tuple(args...));
   using autodiff_idxs_t
-      = filtered_tuple_indices_t<contains_autodiff, args_tuple_t>;
-  using grads_t = std::decay_t<decltype(grads)>;
-  using grad_idxs_t = std::make_index_sequence<std::tuple_size<grads_t>::value>;
+      = filtered_tuple_indices_t<contains_autodiff, std::tuple<TArgs...>>;
+  auto grads = stan::math::zeroed_filtered_tuple<contains_autodiff>(
+      std::forward_as_tuple(args...));
+  using grad_idxs_t = std::make_index_sequence<
+      std::tuple_size<std::decay_t<decltype(grads)>>::value>;
+  FvarInnerT rtn_value;
   finite_diff_gradient_auto(func, rtn_value, std::forward_as_tuple(args...),
                             std::forward_as_tuple(value_of(args)...), grads,
                             autodiff_idxs_t{}, grad_idxs_t{});
@@ -123,6 +142,8 @@ inline auto finite_diff(const F& func, const TArgs&... args) {
  *                scalar type is `fvar<T>`
  * @param func Functor
  * @param args... Parameter pack of arguments to be passed to functor.
+ * @return direct result of `func(args...)`.
+ * @throw Any exception thrown by `func`.
  */
 template <typename F, typename... TArgs,
           require_all_not_st_fvar<TArgs...>* = nullptr>
diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp
index 518c87efa41..42317813f3a 100644
--- a/stan/math/prim/fun.hpp
+++ b/stan/math/prim/fun.hpp
@@ -344,7 +344,7 @@
 #include <stan/math/prim/fun/variance.hpp>
 #include <stan/math/prim/fun/welford_covar_estimator.hpp>
 #include <stan/math/prim/fun/welford_var_estimator.hpp>
-#include <stan/math/prim/fun/zeroed_container.hpp>
+#include <stan/math/prim/fun/zeroed_filtered_tuple.hpp>
 #include <stan/math/prim/fun/zeros_array.hpp>
 #include <stan/math/prim/fun/zeros_int_array.hpp>
 #include <stan/math/prim/fun/zeros_row_vector.hpp>
diff --git a/stan/math/prim/fun/zeroed_container.hpp b/stan/math/prim/fun/zeroed_filtered_tuple.hpp
similarity index 59%
rename from stan/math/prim/fun/zeroed_container.hpp
rename to stan/math/prim/fun/zeroed_filtered_tuple.hpp
index adceb94574c..7e5b73a38d2 100644
--- a/stan/math/prim/fun/zeroed_container.hpp
+++ b/stan/math/prim/fun/zeroed_filtered_tuple.hpp
@@ -1,5 +1,5 @@
-#ifndef STAN_MATH_PRIM_FUN_ZEROED_CONTAINER_HPP
-#define STAN_MATH_PRIM_FUN_ZEROED_CONTAINER_HPP
+#ifndef STAN_MATH_PRIM_FUN_ZEROED_FILTERED_TUPLE_HPP
+#define STAN_MATH_PRIM_FUN_ZEROED_FILTERED_TUPLE_HPP
 
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/functor/filter_map.hpp>
@@ -15,27 +15,37 @@ namespace internal {
 
 /**
  * Type-dependent false helper for `static_assert` in unreachable branches.
+ *
+ * @tparam ... ignored dependent types used to delay `static_assert`
  */
 template <typename...>
-struct zeroed_container_dependent_false : std::false_type {};
+struct zeroed_filtered_tuple_dependent_false : std::false_type {};
 
 /**
  * Build a zero-initialized container matching the shape of `arg`.
  *
  * Scalars produce zero partials scalars, Eigen produces zero matrices of the
  * same dimensions, and std::vector recurses by element.
+ *
+ * @tparam Arg input container/scalar type
+ * @param[in] arg input value used only for shape
+ * @return zero-initialized container with the same structure and dimensions
+ * as `arg`
+ * @throw `std::bad_alloc` if allocation of nested `std::vector` storage fails.
+ * @throw None otherwise. Unsupported categories fail at compile time via
+ * `static_assert`.
  */
 template <typename Arg>
-inline auto zeroed_container_leaf(const Arg& arg) {
+inline auto zeroed_filtered_tuple_leaf(const Arg& arg) {
   using arg_t = std::decay_t<Arg>;
   using partial_t = partials_type_t<scalar_type_t<arg_t>>;
   if constexpr (is_std_vector_v<arg_t>) {
     using elem_t = std::decay_t<decltype(
-        zeroed_container_leaf(std::declval<const value_type_t<arg_t>&>()))>;
+        zeroed_filtered_tuple_leaf(std::declval<const value_type_t<arg_t>&>()))>;
     std::vector<elem_t> ret;
     ret.reserve(arg.size());
     for (const auto& element : arg) {
-      ret.push_back(zeroed_container_leaf(element));
+      ret.push_back(zeroed_filtered_tuple_leaf(element));
     }
     return ret;
   } else if constexpr (is_eigen_v<arg_t>) {
@@ -45,8 +55,8 @@ inline auto zeroed_container_leaf(const Arg& arg) {
     return partial_t(0);
   } else {
     static_assert(
-        zeroed_container_dependent_false<arg_t>::value,
-        "Unsupported container in zeroed_container.");
+        zeroed_filtered_tuple_dependent_false<arg_t>::value,
+        "Unsupported container in zeroed_filtered_tuple.");
   }
 }
 
@@ -61,19 +71,23 @@ inline auto zeroed_container_leaf(const Arg& arg) {
  *
  * @tparam Predicate unary type-trait predicate
  * @tparam TupleArgs tuple input type
- * @param args tuple of arguments used for shape and filtering
+ * @param[in] args tuple of arguments used for shape and filtering
  * @return tuple of zero-initialized containers for filtered entries
+ * @throw `std::bad_alloc` if nested vector allocation fails while constructing
+ * the returned containers.
+ * @throw None otherwise. Unsupported categories fail at compile time via
+ * `static_assert`.
  */
 template <template <class...> class Predicate, typename TupleArgs,
           require_tuple_t<TupleArgs>* = nullptr>
-inline auto zeroed_container(TupleArgs&& args) {
+inline auto zeroed_filtered_tuple(TupleArgs&& args) {
   return stan::math::filter_map<Predicate>(
       [](auto&& arg) {
         using arg_t = std::decay_t<decltype(arg)>;
         if constexpr (is_tuple_v<arg_t>) {
-          return zeroed_container<Predicate>(std::forward<decltype(arg)>(arg));
+          return zeroed_filtered_tuple<Predicate>(std::forward<decltype(arg)>(arg));
         } else {
-          return internal::zeroed_container_leaf(arg);
+          return internal::zeroed_filtered_tuple_leaf(arg);
         }
       },
       std::forward<TupleArgs>(args));
diff --git a/stan/math/prim/functor.hpp b/stan/math/prim/functor.hpp
index d735a5143ea..b1827f1d195 100644
--- a/stan/math/prim/functor.hpp
+++ b/stan/math/prim/functor.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/functor/apply_vector_unary.hpp>
 #include <stan/math/prim/functor/coupled_ode_system.hpp>
 #include <stan/math/prim/functor/filter_map.hpp>
+#include <stan/math/prim/functor/filter_types.hpp>
 #include <stan/math/prim/functor/finite_diff_gradient.hpp>
 #include <stan/math/prim/functor/finite_diff_gradient_auto.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
diff --git a/stan/math/prim/functor/filter_types.hpp b/stan/math/prim/functor/filter_types.hpp
new file mode 100644
index 00000000000..5c81ee65866
--- /dev/null
+++ b/stan/math/prim/functor/filter_types.hpp
@@ -0,0 +1,38 @@
+#ifndef STAN_MATH_PRIM_FUNCTOR_FILTER_TYPES_HPP
+#define STAN_MATH_PRIM_FUNCTOR_FILTER_TYPES_HPP
+
+#include <stan/math/prim/functor/filter_map.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <utility>
+
+namespace stan {
+namespace math {
+
+/**
+ * Filter a tuple and return a tuple containing references to entries whose
+ * types satisfy `Predicate`.
+ *
+ * This helper preserves nested tuple structure and delegates recursion to
+ * `filter_map`.
+ *
+ * @tparam Predicate unary type-trait predicate
+ * @tparam T tuple-like input type
+ * @param[in] t input tuple to filter
+ * @return tuple containing only filtered entries; returned elements preserve
+ * reference/value category from `t`
+ * @throw None.
+ */
+template <template <typename...> class Predicate, typename T,
+          require_tuple_t<T>* = nullptr>
+inline constexpr decltype(auto) filter_types(T&& t) {
+  return stan::math::filter_map<Predicate>(
+      [](auto&& arg) -> decltype(auto) {
+        return std::forward<decltype(arg)>(arg);
+      },
+      std::forward<T>(t));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/functor/finite_diff_gradient_auto.hpp b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
index 5b6dd7774ec..deaaa9746aa 100644
--- a/stan/math/prim/functor/finite_diff_gradient_auto.hpp
+++ b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
@@ -22,36 +22,28 @@ namespace internal {
 
 /**
  * Type-dependent false helper for `static_assert` in unreachable branches.
+ *
+ * @tparam ... ignored dependent types used to delay `static_assert`
  */
 template <typename...>
 struct dependent_false : std::false_type {};
 
 /**
- * Filter a tuple and return a tuple containing references to the elements
- * whose types satisfy `contains_autodiff`.
+ * Construct a holder tuple containing elements selected by `Is...`.
  *
- * This mirrors the `filter_var_scalar_types` pattern in rev.
+ * Returned elements preserve reference/value category through
+ * `make_holder_tuple`, so this helper can be used with tuples of references.
  *
- * @tparam T input tuple or nested tuple-like structure
- * @param t input to filter
- * @return tuple of references to autodiff-containing entries
+ * @tparam Tuple tuple-like input type
+ * @tparam Is selected index pack
+ * @param[in] tuple input tuple
+ * @param[in] selected_idxs index sequence selecting tuple entries
+ * @return holder tuple containing selected entries
+ * @throw None.
  */
-template <typename T>
-inline constexpr decltype(auto) filter_autodiff_types(T&& t) {
-  return stan::math::filter_map<contains_autodiff>(
-      [](auto&& arg) -> decltype(auto) {
-        using arg_t = std::decay_t<decltype(arg)>;
-        if constexpr (is_tuple_v<arg_t>) {
-          return filter_autodiff_types(std::forward<decltype(arg)>(arg));
-        } else {
-          return std::forward<decltype(arg)>(arg);
-        }
-      },
-      std::forward<T>(t));
-}
-
 template <typename Tuple, std::size_t... Is>
-inline auto tuple_subset(Tuple&& tuple, std::index_sequence<Is...>) {
+inline auto tuple_subset(Tuple&& tuple,
+                         [[maybe_unused]] std::index_sequence<Is...> selected_idxs) {
   return make_holder_tuple(std::get<Is>(std::forward<Tuple>(tuple))...);
 }
 
@@ -60,8 +52,9 @@ inline auto tuple_subset(Tuple&& tuple, std::index_sequence<Is...>) {
  * computation from scalar-like autodiff stacks.
  *
  * @tparam T scalar type
- * @param x scalar value
+ * @param[in] x scalar value
  * @return value used to compute finite-difference step size
+ * @throw None.
  */
 template <typename T>
 inline double finite_diff_stepsize_value(const T& x) {
@@ -81,18 +74,32 @@ inline double finite_diff_stepsize_value(const T& x) {
  * `arg_mask` controls whether a scalar coordinate is perturbable, while
  * `arg_work` is the mutable coordinate that gets perturbed.
  *
+ * Supported containers are:
+ * - Stan scalars
+ * - Eigen dense objects
+ * - `std::vector` (possibly nested)
+ * - `std::tuple` (including tuple-of-tuples)
+ *
+ * The `arg_mask` and `arg_work` structures must match exactly. The `grad`
+ * structure may either:
+ * - match top-level tuple arity with `arg_mask`, or
+ * - be compacted to only autodiff-containing top-level tuple entries.
+ *
  * @tparam ArgMask traversal mask type
  * @tparam ArgWork mutable argument type
  * @tparam Grad gradient output type
  * @tparam F callable compatible with `(arg_work_coord, grad_coord)`
- * @param arg_mask mask selecting perturbable coordinates
- * @param arg_work mutable argument coordinates
- * @param grad gradient output coordinates
- * @param fn callback applied for each perturbable coordinate
+ * @param[in] arg_mask mask selecting perturbable coordinates
+ * @param[in,out] arg_work mutable argument coordinates
+ * @param[out] grad gradient output coordinates
+ * @param[in] fn callback applied for each perturbable coordinate
+ * @return None.
+ * @throw None. Fails compilation with `static_assert` for unsupported
+ * container categories or mismatched tuple arity.
  */
 template <typename ArgMask, typename ArgWork, typename Grad, typename F>
-inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, ArgWork&& arg_work,
-                                         Grad&& grad, F&& fn) {
+inline void for_each_leaf_pair(ArgMask&& arg_mask, ArgWork&& arg_work,
+                               Grad&& grad, F&& fn) {
   using arg_mask_t = std::decay_t<ArgMask>;
   if constexpr (is_stan_scalar_v<arg_mask_t>) {
     if constexpr (contains_autodiff_v<arg_mask_t>) {
@@ -100,12 +107,12 @@ inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, ArgWork&& arg_work,
     }
   } else if constexpr (is_eigen_v<arg_mask_t>) {
     for (Eigen::Index i = 0; i < arg_mask.size(); ++i) {
-      for_each_coordinate_pair_mut(arg_mask.coeffRef(i), arg_work.coeffRef(i),
-                                   grad.coeffRef(i), fn);
+      for_each_leaf_pair(arg_mask.coeffRef(i), arg_work.coeffRef(i),
+                         grad.coeffRef(i), fn);
     }
   } else if constexpr (is_std_vector_v<arg_mask_t>) {
     for (std::size_t i = 0; i < arg_mask.size(); ++i) {
-      for_each_coordinate_pair_mut(arg_mask[i], arg_work[i], grad[i], fn);
+      for_each_leaf_pair(arg_mask[i], arg_work[i], grad[i], fn);
     }
   } else if constexpr (is_tuple_v<arg_mask_t>) {
     constexpr auto arg_mask_size = std::tuple_size<std::decay_t<ArgMask>>::value;
@@ -117,7 +124,7 @@ inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, ArgWork&& arg_work,
     if constexpr (arg_mask_size == grad_size) {
       stan::math::for_each(
           [&fn](auto&& arg_mask_i, auto&& arg_work_i, auto&& grad_i) {
-            for_each_coordinate_pair_mut(arg_mask_i, arg_work_i, grad_i, fn);
+            for_each_leaf_pair(arg_mask_i, arg_work_i, grad_i, fn);
           },
           std::forward<ArgMask>(arg_mask), std::forward<ArgWork>(arg_work),
           std::forward<Grad>(grad));
@@ -136,7 +143,7 @@ inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, ArgWork&& arg_work,
           "Tuple size mismatch in finite_diff_gradient_auto traversal");
       stan::math::for_each(
           [&fn](auto&& arg_mask_i, auto&& arg_work_i, auto&& grad_i) {
-            for_each_coordinate_pair_mut(arg_mask_i, arg_work_i, grad_i, fn);
+            for_each_leaf_pair(arg_mask_i, arg_work_i, grad_i, fn);
           },
           std::move(arg_mask_ad), std::move(arg_work_ad),
           std::forward<Grad>(grad));
@@ -160,11 +167,14 @@ inline void for_each_coordinate_pair_mut(ArgMask&& arg_mask, ArgWork&& arg_work,
  * @tparam ArgWork mutable working argument type
  * @tparam NthGrad gradient container type for this top-level argument
  * @tparam TupleArgsWork mutable tuple of all working arguments passed to `f`
- * @param f function being finite-differenced
- * @param arg_mask mask argument used to choose perturbable coordinates
- * @param arg_work mutable working argument coordinates
- * @param grad gradient output container for this top-level argument
- * @param args_work mutable tuple of all working arguments passed to `f`
+ * @param[in] f function being finite-differenced
+ * @param[in] arg_mask mask argument used to choose perturbable coordinates
+ * @param[in,out] arg_work mutable working argument coordinates
+ * @param[out] grad gradient output container for this top-level argument
+ * @param[in,out] args_work mutable tuple of all working arguments passed to
+ * `f`
+ * @return None.
+ * @throw Any exception thrown by `f`.
  */
 template <typename ScalarT, typename F, typename ArgMask, typename ArgWork,
           typename NthGrad, typename TupleArgsWork>
@@ -174,7 +184,7 @@ inline void finite_diff_gradient_auto_impl(const F& f, ArgMask&& arg_mask,
   static constexpr int h_scale[6] = {3, 2, 1, -3, -2, -1};
   static constexpr int mults[6] = {1, -9, 45, -1, 9, -45};
 
-  for_each_coordinate_pair_mut(
+  for_each_leaf_pair(
       std::forward<ArgMask>(arg_mask), std::forward<ArgWork>(arg_work),
       std::forward<NthGrad>(grad), [&f, &args_work](auto& arg_coord, auto& grad_coord) {
         const auto orig = arg_coord;
@@ -231,6 +241,8 @@ inline void finite_diff_gradient_auto_impl(const F& f, ArgMask&& arg_mask,
  * @param[in] x argument to function
  * @param[out] fx function applied to argument
  * @param[out] grad_fx gradient of function at argument
+ * @return None.
+ * @throw Any exception thrown by `f`.
  */
 template <typename F, typename VectorT, typename GradVectorT,
           typename ScalarT = return_type_t<VectorT>,
@@ -258,25 +270,19 @@ inline void finite_diff_gradient_auto(const F& f, VectorT&& x, ScalarT& fx,
   });
 }
 
-/**
- * Filter a tuple and return a tuple with references to entries with autodiff
- * scalar type.
- *
- * @tparam T Possibly a tuple, std::vector, Eigen type, or scalar
- * @param[in] t Input to filter
- * @return Filtered input with only autodiff scalar entries
- */
-template <typename T>
-inline constexpr decltype(auto) filter_ad_scalar_types(T&& t) {
-  return internal::filter_autodiff_types(std::forward<T>(t));
-}
-
 /**
  * Calculate the function value and gradients for heterogeneous tuple arguments
  * while preserving argument scalar types at the call boundary.
  *
  * This overload expects `grads` to mirror the top-level structure of `args`.
- * Non-autodiff top-level entries are skipped.
+ * Non-autodiff top-level entries are skipped in the finite-difference pass.
+ *
+ * Pattern used by this overload:
+ * 1. Evaluate `f(args...)` once for `fx`.
+ * 2. Zip top-level `args` and `grads`.
+ * 3. Use `contains_autodiff_v` to skip non-autodiff top-level entries.
+ * 4. Traverse nested leaves with `for_each_leaf_pair`, perturbing one
+ *    coordinate at a time and restoring it after each stencil evaluation.
  *
  * @tparam F callable type
  * @tparam ScalarT scalar type used for function value and gradient storage
@@ -286,6 +292,8 @@ inline constexpr decltype(auto) filter_ad_scalar_types(T&& t) {
  * @param[out] fx function value at `args`
  * @param[in] args tuple of function arguments
  * @param[out] grads tuple of gradient containers
+ * @return None.
+ * @throw Any exception thrown by `f`.
  */
 template <typename F, typename ScalarT, typename TupleArgs, typename TupleGrads,
           require_tuple_t<TupleArgs>* = nullptr,
@@ -306,13 +314,28 @@ inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
 /**
  * Indexed tuple overload used for compact gradient tuples. `Idxs` select
  * autodiff-containing top-level entries from both `args_mask` and `args_work`,
- * while `GradIdxs` selects matching gradient entries from `grads`.
+ * while `GradIdxs` selects matching compact gradient entries from `grads`.
+ *
+ * This overload is used by `stan::math::finite_diff` where:
+ * - `args_mask` preserves original autodiff type information
+ * - `args_work` is the value-only work tuple passed to function evaluations
+ * - `grads` stores only autodiff top-level gradient containers
  *
  * @tparam F callable type
  * @tparam ScalarT scalar type used for function value and gradients
  * @tparam TupleArgsMask tuple used only as perturbation mask
  * @tparam TupleArgsWork mutable tuple used for function evaluation
  * @tparam TupleGrads compact tuple of gradient containers
+ * @param[in] f function to differentiate
+ * @param[out] fx function value at `args_work`
+ * @param[in] args_mask tuple carrying autodiff mask information
+ * @param[in,out] args_work mutable value tuple used for finite-difference
+ * evaluations
+ * @param[out] grads compact tuple of gradient containers
+ * @param[in] idxs selected top-level indices in `args_mask`/`args_work`
+ * @param[in] grad_idxs selected top-level indices in `grads`
+ * @return None.
+ * @throw Any exception thrown by `f`.
  */
 template <typename F, typename ScalarT, typename TupleArgsMask,
           typename TupleArgsWork, typename TupleGrads, std::size_t... Idxs,
@@ -321,8 +344,8 @@ template <typename F, typename ScalarT, typename TupleArgsMask,
           require_tuple_t<TupleGrads>* = nullptr>
 inline void finite_diff_gradient_auto(
     const F& f, ScalarT& fx, TupleArgsMask&& args_mask, TupleArgsWork&& args_work,
-    TupleGrads&& grads, std::index_sequence<Idxs...> /* idxs */,
-    std::index_sequence<GradIdxs...> /* grad_idxs */) {
+    TupleGrads&& grads, [[maybe_unused]] std::index_sequence<Idxs...> idxs,
+    [[maybe_unused]] std::index_sequence<GradIdxs...> grad_idxs) {
   static_assert(sizeof...(Idxs) == sizeof...(GradIdxs),
                 "Index sequence size mismatch in finite_diff_gradient_auto.");
   fx = stan::math::apply([&f](auto&&... args_i) { return f(args_i...); },
diff --git a/stan/math/prim/meta/conditional_copy_and_promote.hpp b/stan/math/prim/meta/conditional_copy_and_promote.hpp
index b895a98eb61..6a81a6d3369 100644
--- a/stan/math/prim/meta/conditional_copy_and_promote.hpp
+++ b/stan/math/prim/meta/conditional_copy_and_promote.hpp
@@ -26,9 +26,11 @@ enum class COPY_TYPE : uint8_t { SHALLOW = 0, DEEP = 1 };
  * @tparam PromotedType type to promote the scalar to.
  * @tparam CopyType type of copy to perform.
  * @tparam Args variadic arguments.
- * @param args variadic arguments to conditionally copy and promote.
+ * @param[in] args variadic arguments to conditionally copy and promote.
  * @return a tuple where each element is either a reference to the original
  * argument or a promoted copy of the argument.
+ * @throw `std::bad_alloc` if allocation of nested `std::vector` storage fails.
+ * @throw Any exception thrown by element copy/promotion routines.
  */
 template <template <typename...> class Filter, typename PromotedType = double,
           COPY_TYPE CopyType = COPY_TYPE::DEEP, typename... Args>
@@ -78,9 +80,11 @@ inline auto conditional_copy_and_promote(Args&&... args) {
  * Conditional deep copy types with a `var` scalar type to `PromotedType`.
  * @tparam PromotedType type to promote the scalar to.
  * @tparam Args variadic arguments.
- * @param args variadic arguments to conditionally copy and promote.
+ * @param[in] args variadic arguments to conditionally copy and promote.
  * @return a tuple where each element is either a reference to the original
  * argument or a promoted copy of the argument.
+ * @throw `std::bad_alloc` if allocation of nested `std::vector` storage fails.
+ * @throw Any exception thrown by element copy/promotion routines.
  */
 template <typename PromotedType, typename... Args>
 inline auto deep_copy_vargs(Args&&... args) {
@@ -96,9 +100,11 @@ inline auto deep_copy_vargs(Args&&... args) {
  *  in an inner autodiff without making a hard copy of the input arguments.
  * @tparam PromotedType type to promote the scalar to.
  * @tparam Args variadic arguments.
- * @param args variadic arguments to conditionally copy and promote.
+ * @param[in] args variadic arguments to conditionally copy and promote.
  * @return a tuple where each element is either a reference to the original
  * argument or a promoted copy of the argument.
+ * @throw `std::bad_alloc` if allocation of nested `std::vector` storage fails.
+ * @throw Any exception thrown by element copy/promotion routines.
  */
 template <typename PromotedType, typename... Args>
 inline auto shallow_copy_vargs(Args&&... args) {
diff --git a/stan/math/prim/meta/contains_autodiff.hpp b/stan/math/prim/meta/contains_autodiff.hpp
index 40b3447caa1..f2830310b05 100644
--- a/stan/math/prim/meta/contains_autodiff.hpp
+++ b/stan/math/prim/meta/contains_autodiff.hpp
@@ -16,18 +16,47 @@ namespace stan {
  * into their contained element types.
  *
  * @tparam T type to inspect
+ * @return Inherits `bool_constant<true>` when `T` is an autodiff scalar;
+ * otherwise `bool_constant<false>`.
+ * @throw None.
  */
 template <typename T>
 struct contains_autodiff : bool_constant<is_autodiff_v<std::decay_t<T>>> {};
 
+/**
+ * Tuple specialization. Returns true if any top-level tuple element recursively
+ * contains an autodiff scalar.
+ *
+ * @tparam Args tuple element types
+ * @return Inherits `bool_constant<true>` when any tuple element recursively
+ * contains an autodiff scalar.
+ * @throw None.
+ */
 template <typename... Args>
 struct contains_autodiff<std::tuple<Args...>>
     : bool_constant<(contains_autodiff<std::decay_t<Args>>::value || ...)> {};
 
+/**
+ * `std::vector` specialization. Returns whether the element type recursively
+ * contains an autodiff scalar.
+ *
+ * @tparam T vector element type
+ * @tparam VecArgs allocator and container argument types
+ * @return Inherits `bool_constant<true>` when `T` recursively contains an
+ * autodiff scalar.
+ * @throw None.
+ */
 template <typename T, typename... VecArgs>
 struct contains_autodiff<std::vector<T, VecArgs...>>
     : bool_constant<contains_autodiff<std::decay_t<T>>::value> {};
 
+/**
+ * Convenience variable template for `contains_autodiff`.
+ *
+ * @tparam T type to inspect
+ * @return `true` when `T` recursively contains at least one autodiff scalar.
+ * @throw None.
+ */
 template <typename T>
 inline constexpr bool contains_autodiff_v
     = contains_autodiff<std::decay_t<T>>::value;
diff --git a/stan/math/prim/meta/filter_index_sequence.hpp b/stan/math/prim/meta/filter_index_sequence.hpp
index 569f406c29b..cfbfd32ef37 100644
--- a/stan/math/prim/meta/filter_index_sequence.hpp
+++ b/stan/math/prim/meta/filter_index_sequence.hpp
@@ -7,25 +7,55 @@ namespace stan {
 namespace math {
 
 /**
- * Concatenates integer sequences.
+ * Concatenate two integer sequences of the same value type.
+ *
+ * @tparam SeqT integral value type used by the sequences
+ * @tparam Is first sequence indices
+ * @tparam Js second sequence indices
+ * @param[in] first first integer sequence
+ * @param[in] second second integer sequence
+ * @return concatenated sequence containing `Is...` followed by `Js...`
+ * @throw None.
  */
 template <typename SeqT, SeqT... Is, SeqT... Js>
 inline constexpr std::integer_sequence<SeqT, Is..., Js...>
-concat_integer_sequences(std::integer_sequence<SeqT, Is...>,
-                         std::integer_sequence<SeqT, Js...>) {
+concat_integer_sequences([[maybe_unused]] std::integer_sequence<SeqT, Is...> first,
+                         [[maybe_unused]] std::integer_sequence<SeqT, Js...> second) {
   return {};
 }
 
+/**
+ * Identity overload for single-sequence concatenation.
+ *
+ * @tparam SeqT integral value type used by the sequence
+ * @param[in] seq input integer sequence
+ * @return `seq` unchanged
+ * @throw None.
+ */
 template <typename SeqT>
 inline constexpr std::integer_sequence<SeqT> concat_integer_sequences(
-    std::integer_sequence<SeqT>) {
+    [[maybe_unused]] std::integer_sequence<SeqT> seq) {
   return {};
 }
 
+/**
+ * Concatenate three or more integer sequences.
+ *
+ * @tparam SeqT integral value type used by the sequences
+ * @tparam Is indices in the first sequence
+ * @tparam Js indices in the second sequence
+ * @tparam Rest remaining sequence types
+ * @param[in] first first integer sequence
+ * @param[in] second second integer sequence
+ * @param[in] rest remaining integer sequences
+ * @return concatenated sequence preserving argument order
+ * @throw None.
+ */
 template <typename SeqT, SeqT... Is, SeqT... Js, class... Rest>
 inline constexpr auto concat_integer_sequences(
-    std::integer_sequence<SeqT, Is...>, std::integer_sequence<SeqT, Js...>,
-    Rest...) {
+    [[maybe_unused]] std::integer_sequence<SeqT, Is...> first,
+    [[maybe_unused]] std::integer_sequence<SeqT, Js...> second,
+    [[maybe_unused]] Rest... rest) {
   return concat_integer_sequences(std::integer_sequence<SeqT, Is..., Js...>{},
                                   Rest{}...);
 }
@@ -33,11 +63,20 @@ inline constexpr auto concat_integer_sequences(
 /**
  * Returns a singleton index sequence if `Predicate<Arg>::value` is true,
  * otherwise an empty sequence.
+ *
+ * @tparam Predicate unary predicate type trait
+ * @tparam Arg argument type evaluated by `Predicate`
+ * @tparam SeqT integral value type used by the sequence
+ * @tparam I index candidate
+ * @param[in] idx singleton input sequence containing `I`
+ * @return `std::integer_sequence<SeqT, I>` if `Predicate<Arg>::value`,
+ * otherwise `std::integer_sequence<SeqT>{}`
+ * @throw None.
  */
 template <template <class...> class Predicate, typename Arg, class SeqT,
           SeqT I>
 inline constexpr auto filter_single_index_sequence(
-    std::integer_sequence<SeqT, I>) {
+    [[maybe_unused]] std::integer_sequence<SeqT, I> idx) {
   if constexpr (Predicate<Arg>::value) {
     return std::integer_sequence<SeqT, I>{};
   } else {
@@ -49,10 +88,21 @@ inline constexpr auto filter_single_index_sequence(
  * Filters an index sequence by applying `Predicate` to matching `Args...`.
  *
  * `Args...` and the index sequence must have the same length.
+ *
+ * @tparam Predicate unary predicate type trait
+ * @tparam Args argument type pack aligned with `Is...`
+ * @tparam SeqT integral value type used by the sequence
+ * @tparam Is index pack to filter
+ * @param[in] idxs input index sequence
+ * @return compact integer sequence of indices where
+ * `Predicate<Args>::value` is true
+ * @throw None. Fails compilation with `static_assert` if `sizeof...(Args) !=
+ * sizeof...(Is)`.
  */
 template <template <class...> class Predicate, typename... Args, class SeqT,
           SeqT... Is>
-inline constexpr auto filter_index_sequence(std::integer_sequence<SeqT, Is...>) {
+inline constexpr auto filter_index_sequence(
+    [[maybe_unused]] std::integer_sequence<SeqT, Is...> idxs) {
   static_assert(sizeof...(Args) == sizeof...(Is),
                 "filter_index_sequence requires the same number of argument "
                 "types and indices.");
@@ -65,6 +115,14 @@ inline constexpr auto filter_index_sequence(std::integer_sequence<SeqT, Is...>)
   }
 }
 
+/**
+ * Alias returning the filtered sequence type produced by
+ * `filter_index_sequence`.
+ *
+ * @tparam Predicate unary predicate type trait
+ * @tparam Sequence input index sequence type
+ * @tparam Args argument type pack aligned with `Sequence`
+ */
 template <template <class...> class Predicate, typename Sequence,
           typename... Args>
 using filter_index_sequence_t
diff --git a/stan/math/prim/meta/filtered_tuple_indices.hpp b/stan/math/prim/meta/filtered_tuple_indices.hpp
index 7b6e57dbdd4..5f620b2fa46 100644
--- a/stan/math/prim/meta/filtered_tuple_indices.hpp
+++ b/stan/math/prim/meta/filtered_tuple_indices.hpp
@@ -9,9 +9,25 @@
 namespace stan {
 namespace math {
 
+/**
+ * Type trait returning the number of elements in an integer sequence.
+ *
+ * @tparam Sequence `std::integer_sequence` type
+ * @return Inherits `std::integral_constant<std::size_t, N>` where `N` is the
+ * sequence length.
+ * @throw None.
+ */
 template <typename Sequence>
 struct index_sequence_size;
 
+/**
+ * Specialization extracting the size of a concrete integer sequence.
+ *
+ * @tparam T integral value type used by the sequence
+ * @tparam Is sequence values
+ * @return Inherits `std::integral_constant<std::size_t, sizeof...(Is)>`.
+ * @throw None.
+ */
 template <typename T, T... Is>
 struct index_sequence_size<std::integer_sequence<T, Is...>>
     : std::integral_constant<std::size_t, sizeof...(Is)> {};
@@ -19,24 +35,63 @@ struct index_sequence_size<std::integer_sequence<T, Is...>>
 /**
  * Produces a compact index sequence of top-level tuple elements whose
  * types satisfy `Predicate`.
+ *
+ * @tparam Predicate unary predicate type trait
+ * @tparam Tuple tuple type to filter
+ * @return Provides nested `type` alias with filtered indices.
+ * @throw None.
  */
 template <template <class...> class Predicate, typename Tuple>
 struct filtered_tuple_indices;
 
+/**
+ * Tuple specialization that filters top-level tuple indices using `Predicate`.
+ *
+ * @tparam Predicate unary predicate type trait
+ * @tparam Args tuple element types
+ * @return Provides nested `type` alias equal to the filtered index sequence.
+ * @throw None.
+ */
 template <template <class...> class Predicate, typename... Args>
 struct filtered_tuple_indices<Predicate, std::tuple<Args...>> {
   using type = filter_index_sequence_t<Predicate, std::index_sequence_for<Args...>,
                                        std::decay_t<Args>...>;
 };
 
+/**
+ * Alias for the compact sequence of tuple indices whose element types satisfy
+ * `Predicate`.
+ *
+ * @tparam Predicate unary predicate type trait
+ * @tparam Tuple tuple type to filter
+ * @return Filtered tuple index-sequence type.
+ * @throw None.
+ */
 template <template <class...> class Predicate, typename Tuple>
 using filtered_tuple_indices_t
     = typename filtered_tuple_indices<Predicate, std::decay_t<Tuple>>::type;
 
+/**
+ * Type trait returning the number of filtered tuple indices.
+ *
+ * @tparam Predicate unary predicate type trait
+ * @tparam Tuple tuple type to filter
+ * @return `std::integral_constant<std::size_t, N>` where `N` is the filtered
+ * arity.
+ * @throw None.
+ */
 template <template <class...> class Predicate, typename Tuple>
 using filtered_tuple_size
     = index_sequence_size<filtered_tuple_indices_t<Predicate, Tuple>>;
 
+/**
+ * Convenience variable template for `filtered_tuple_size`.
+ *
+ * @tparam Predicate unary predicate type trait
+ * @tparam Tuple tuple type to filter
+ * @return Number of top-level tuple elements satisfying `Predicate`.
+ * @throw None.
+ */
 template <template <class...> class Predicate, typename Tuple>
 inline constexpr std::size_t filtered_tuple_size_v
     = filtered_tuple_size<Predicate, Tuple>::value;
diff --git a/stan/math/rev/core/filter_var_scalar_types.hpp b/stan/math/rev/core/filter_var_scalar_types.hpp
index f18df646d7b..fc2f3a8139a 100644
--- a/stan/math/rev/core/filter_var_scalar_types.hpp
+++ b/stan/math/rev/core/filter_var_scalar_types.hpp
@@ -2,7 +2,7 @@
 #define STAN_MATH_REV_CORE_FILTER_VAR_SCALAR_TYPES_HPP
 
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <stan/math/prim/functor/filter_map.hpp>
+#include <stan/math/prim/functor/filter_types.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <type_traits>
 
@@ -13,18 +13,11 @@ namespace stan::math::internal {
  * @tparam T Possibly a tuple, std::vector, Eigen type, or scalar
  * @param[in] t Input to filter
  * @return Filtered input with only var scalar types
+ * @throw None.
  */
 template <typename T>
 inline constexpr decltype(auto) filter_var_scalar_types(T&& t) {
-  return stan::math::filter_map<is_any_var_scalar>(
-      [](auto&& arg) -> decltype(auto) {
-        if constexpr (is_tuple_v<std::decay_t<decltype(arg)>>) {
-          return filter_var_scalar_types(std::forward<decltype(arg)>(arg));
-        } else {
-          return std::forward<decltype(arg)>(arg);
-        }
-      },
-      std::forward<T>(t));
+  return stan::math::filter_types<is_any_var_scalar>(std::forward<T>(t));
 }
 
 }  // namespace stan::math::internal
diff --git a/test/unit/math/fwd/fun/tangent_of_test.cpp b/test/unit/math/mix/fun/tangent_of_test.cpp
similarity index 100%
rename from test/unit/math/fwd/fun/tangent_of_test.cpp
rename to test/unit/math/mix/fun/tangent_of_test.cpp
diff --git a/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp b/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp
index 9d135ce0a13..5ded429fcbd 100644
--- a/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp
+++ b/test/unit/math/mix/functor/finite_diff_gradient_auto_test.cpp
@@ -86,3 +86,56 @@ TEST(MathMixMatFunctor, FiniteDiffGradientAutoTupleNested) {
   EXPECT_EQ(0, std::get<1>(std::get<1>(grads))[0]);
   EXPECT_EQ(0, std::get<1>(std::get<1>(grads))[1]);
 }
+
+TEST(MathMixMatFunctor, FiniteDiffGradientAutoTupleRestoresInputs) {
+  using stan::math::fvar;
+  using fvar_stdvec = std::vector<fvar<double>>;
+  fvar_stdvec a(2);
+  a[0] = fvar<double>(1.3, 0.0);
+  a[1] = fvar<double>(-0.7, 0.0);
+  using fvar_mat = Eigen::Matrix<fvar<double>, -1, 1>;
+  fvar_mat b{{fvar<double>(2.0, 0.0), fvar<double>(-3.0, 0.0)}};
+  std::vector<int> x_i{5, 6};
+  Eigen::VectorXd c{{4.0}};
+
+  auto args = std::make_tuple(a, std::make_tuple(b, x_i), c);
+  auto args_before = args;
+  fvar_mat b_grad = fvar_mat::Zero(b.size());
+  auto grads = std::make_tuple(
+      fvar_stdvec(a.size(), 0.0), std::make_tuple(b_grad, x_i),
+      Eigen::VectorXd::Zero(c.size()));
+
+  auto f = [](const auto& arg_a, const auto& arg_b_tuple, const auto& arg_c) {
+    const auto& arg_b = std::get<0>(arg_b_tuple);
+    const auto& arg_i = std::get<1>(arg_b_tuple);
+    return arg_a[0] * arg_a[1] + arg_b(0) * arg_b(1) + arg_c(0) * arg_i[0];
+  };
+
+  fvar<double> fx = 0;
+  stan::math::finite_diff_gradient_auto(f, fx, args, grads);
+
+  const auto& a_after = std::get<0>(args);
+  const auto& a_orig = std::get<0>(args_before);
+  for (std::size_t i = 0; i < a_after.size(); ++i) {
+    EXPECT_FLOAT_EQ(a_orig[i].val_, a_after[i].val_);
+    EXPECT_FLOAT_EQ(a_orig[i].d_, a_after[i].d_);
+  }
+
+  const auto& b_after = std::get<0>(std::get<1>(args));
+  const auto& b_orig = std::get<0>(std::get<1>(args_before));
+  for (Eigen::Index i = 0; i < b_after.size(); ++i) {
+    EXPECT_FLOAT_EQ(b_orig(i).val_, b_after(i).val_);
+    EXPECT_FLOAT_EQ(b_orig(i).d_, b_after(i).d_);
+  }
+
+  const auto& x_i_after = std::get<1>(std::get<1>(args));
+  const auto& x_i_orig = std::get<1>(std::get<1>(args_before));
+  EXPECT_EQ(x_i_orig, x_i_after);
+
+  const auto& c_after = std::get<2>(args);
+  const auto& c_orig = std::get<2>(args_before);
+  ASSERT_EQ(c_orig.size(), c_after.size());
+  for (Eigen::Index i = 0; i < c_after.size(); ++i) {
+    EXPECT_FLOAT_EQ(c_orig(i), c_after(i));
+  }
+}
diff --git a/test/unit/math/fwd/functor/finite_diff_test.cpp b/test/unit/math/mix/functor/finite_diff_test.cpp
similarity index 51%
rename from test/unit/math/fwd/functor/finite_diff_test.cpp
rename to test/unit/math/mix/functor/finite_diff_test.cpp
index 9cdb4bf1ec3..6e80a7b2e8c 100644
--- a/test/unit/math/fwd/functor/finite_diff_test.cpp
+++ b/test/unit/math/mix/functor/finite_diff_test.cpp
@@ -10,6 +10,14 @@ struct finite_diff_typed_args_functor {
   }
 };
 
+struct finite_diff_interleaved_args_functor {
+  template <typename T_a, typename T_b>
+  auto operator()(const std::vector<T_a>& a, const std::vector<int>& x_i,
+                  const std::vector<T_b>& b, const double c) const {
+    return a[0] * x_i[0] + b[0] * x_i[1] + c;
+  }
+};
+
 TEST(FwdFunctor, finiteDiffTypedIntArgs) {
   using stan::math::fvar;
   std::vector<fvar<double>> theta(2);
@@ -39,25 +47,56 @@ TEST(FwdFunctor, finiteDiffTypedIntArgsFvarVar) {
   EXPECT_NEAR(-2.5, out.d_.val(), 1e-8);
 }
 
-struct integrate_1d_finite_diff_functor {
+TEST(FwdFunctor, finiteDiffInterleavedTopLevelArgs) {
+  using stan::math::fvar;
+  std::vector<fvar<double>> a(1);
+  a[0] = fvar<double>(1.2, 0.7);
+  const std::vector<int> x_i{3, 4};
+  std::vector<fvar<double>> b(1);
+  b[0] = fvar<double>(-0.5, -1.3);
+  const double c = 1.1;
+
+  auto out
+      = stan::math::finite_diff(finite_diff_interleaved_args_functor{}, a, x_i,
+                                b, c);
+
+  EXPECT_NEAR(2.7, out.val_, 1e-10);
+  EXPECT_NEAR(-3.1, out.d_, 1e-8);
+}
+
+/**
+ * Regression for https://github.com/stan-dev/math/issues/3280
+ *
+ * Before the tuple-native finite-diff refactor, the forward-mode finite_diff
+ * path could serialize and deserialize arguments in a way that converted
+ * `x_i` from `std::vector<int>` into floating-point vectors. That then failed
+ * to match integrate_1d's required callback signature.
+ *
+ * This test verifies `finite_diff` + `integrate_1d` now works with mixed
+ * autodiff/real/int argument packs and produces the expected derivative.
+ */
+struct issue_3280_integrate_1d_functor {
   template <typename T_x, typename T_xc, typename T_theta>
-  auto operator()(const T_x& x, const T_xc& /* xc */, const std::vector<T_theta>& theta,
+  auto operator()(const T_x& x, const T_xc& /* xc */,
+                  const std::vector<T_theta>& theta,
                   const std::vector<double>& x_r,
-                  const std::vector<int>& x_i, std::ostream* /* msgs */) const {
+                  const std::vector<int>& x_i,
+                  std::ostream* /* msgs */) const {
     return theta[0] * x + x_r[0] + x_i[0];
   }
 };
 
-TEST(FwdFunctor, finiteDiffIntegrate1dMixedArgs) {
+TEST(FwdFunctor, finiteDiffIntegrate1dIssue3280Regression) {
   using stan::math::fvar;
   const double relative_tolerance = 1e-8;
   std::ostream* msgs = nullptr;
+
   std::vector<fvar<double>> theta(1);
   theta[0] = fvar<double>(2.0, 1.0);
   const std::vector<double> x_r{1.5};
   const std::vector<int> x_i{2};
 
-  auto out = stan::math::integrate_1d(integrate_1d_finite_diff_functor{}, 0.0,
+  auto out = stan::math::integrate_1d(issue_3280_integrate_1d_functor{}, 0.0,
                                       1.0, theta, x_r, x_i, msgs,
                                       relative_tolerance);
 
diff --git a/test/unit/math/prim/meta/contains_autodiff_test.cpp b/test/unit/math/mix/meta/contains_autodiff_test.cpp
similarity index 62%
rename from test/unit/math/prim/meta/contains_autodiff_test.cpp
rename to test/unit/math/mix/meta/contains_autodiff_test.cpp
index e7206302557..3c4b95e8a67 100644
--- a/test/unit/math/prim/meta/contains_autodiff_test.cpp
+++ b/test/unit/math/mix/meta/contains_autodiff_test.cpp
@@ -29,3 +29,23 @@ TEST(MathMetaContainsAutodiff, nested_tuple_detection) {
   EXPECT_TRUE((stan::contains_autodiff_v<nested_ad_t>));
   EXPECT_FALSE((stan::contains_autodiff_v<nested_non_ad_t>));
 }
+
+TEST(MathMetaContainsAutodiff, eigen_detection) {
+  using stan::math::fvar;
+  using ad_eigen_t = Eigen::Matrix<fvar<double>, -1, 1>;
+  using non_ad_eigen_t = Eigen::VectorXd;
+
+  EXPECT_TRUE((stan::contains_autodiff_v<ad_eigen_t>));
+  EXPECT_FALSE((stan::contains_autodiff_v<non_ad_eigen_t>));
+}
+
+TEST(MathMetaContainsAutodiff, nested_tuple_with_eigen_detection) {
+  using stan::math::fvar;
+  using nested_ad_eigen_t
+      = std::tuple<int, std::tuple<Eigen::Matrix<fvar<double>, -1, 1>, double>>;
+  using nested_non_ad_eigen_t
+      = std::tuple<int, std::tuple<Eigen::VectorXd, double>>;
+
+  EXPECT_TRUE((stan::contains_autodiff_v<nested_ad_eigen_t>));
+  EXPECT_FALSE((stan::contains_autodiff_v<nested_non_ad_eigen_t>));
+}
diff --git a/test/unit/math/prim/meta/filtered_tuple_indices_test.cpp b/test/unit/math/mix/meta/filtered_tuple_indices_test.cpp
similarity index 100%
rename from test/unit/math/prim/meta/filtered_tuple_indices_test.cpp
rename to test/unit/math/mix/meta/filtered_tuple_indices_test.cpp
diff --git a/test/unit/math/prim/fun/zeroed_container_test.cpp b/test/unit/math/prim/fun/zeroed_container_test.cpp
deleted file mode 100644
index e8c7d16e165..00000000000
--- a/test/unit/math/prim/fun/zeroed_container_test.cpp
+++ /dev/null
@@ -1,56 +0,0 @@
-#include <stan/math/fwd/core/fvar.hpp>
-#include <stan/math/fwd/meta.hpp>
-#include <stan/math/prim/fun/zeroed_container.hpp>
-#include <stan/math/prim/meta/contains_autodiff.hpp>
-#include <gtest/gtest.h>
-#include <tuple>
-#include <type_traits>
-#include <vector>
-
-TEST(MathPrimFunZeroedContainer, autodiff_predicate_filters_and_zeros) {
-  using stan::math::fvar;
-  std::vector<fvar<double>> a{
-      fvar<double>(1.3, 0.0),
-      fvar<double>(-0.7, 0.0),
-  };
-  Eigen::Matrix<fvar<double>, -1, 1> b(2);
-  b << fvar<double>(2.0, 0.0), fvar<double>(-3.0, 0.0);
-  std::vector<int> x_i{5, 6};
-  Eigen::VectorXd c{{4.0}};
-
-  auto args = std::make_tuple(a, std::make_tuple(b, x_i), c);
-  auto grads = stan::math::zeroed_container<stan::contains_autodiff>(args);
-
-  using grads_t = std::decay_t<decltype(grads)>;
-  EXPECT_EQ(2, std::tuple_size<grads_t>::value);
-
-  const auto& a_grads = std::get<0>(grads);
-  EXPECT_EQ(a.size(), a_grads.size());
-  for (std::size_t i = 0; i < a_grads.size(); ++i) {
-    EXPECT_FLOAT_EQ(0.0, a_grads[i]);
-  }
-
-  const auto& b_grads = std::get<0>(std::get<1>(grads));
-  EXPECT_EQ(b.size(), b_grads.size());
-  for (Eigen::Index i = 0; i < b_grads.size(); ++i) {
-    EXPECT_FLOAT_EQ(0.0, b_grads(i));
-  }
-}
-
-TEST(MathPrimFunZeroedContainer, integral_predicate_filters_and_zeros) {
-  auto args = std::make_tuple(std::vector<int>{1, 2, 3}, 1.5,
-                              std::make_tuple(2.5, 7));
-  auto zeros = stan::math::zeroed_container<std::is_integral>(args);
-
-  using zeros_t = std::decay_t<decltype(zeros)>;
-  EXPECT_EQ(2, std::tuple_size<zeros_t>::value);
-
-  const auto& vec_zero = std::get<0>(zeros);
-  EXPECT_EQ(3, vec_zero.size());
-  EXPECT_EQ(0, vec_zero[0]);
-  EXPECT_EQ(0, vec_zero[1]);
-  EXPECT_EQ(0, vec_zero[2]);
-
-  const auto& nested_zero = std::get<0>(std::get<1>(zeros));
-  EXPECT_EQ(0, nested_zero);
-}
diff --git a/test/unit/math/prim/fun/zeroed_filtered_tuple_test.cpp b/test/unit/math/prim/fun/zeroed_filtered_tuple_test.cpp
new file mode 100644
index 00000000000..c9a58488615
--- /dev/null
+++ b/test/unit/math/prim/fun/zeroed_filtered_tuple_test.cpp
@@ -0,0 +1,143 @@
+#include <stan/math/prim/fun/zeroed_filtered_tuple.hpp>
+#include <stan/math/prim/meta/contains_autodiff.hpp>
+#include <test/unit/math/prim/util.hpp>
+#include <gtest/gtest.h>
+#include <tuple>
+#include <type_traits>
+#include <vector>
+
+namespace {
+
+template <typename T>
+struct is_any_fp_scalar_impl {
+  static constexpr bool value
+      = std::is_floating_point_v<stan::scalar_type_t<std::decay_t<T>>>;
+};
+
+template <typename... Types>
+struct is_any_fp_scalar_impl<std::tuple<Types...>> {
+  static constexpr bool value
+      = (is_any_fp_scalar_impl<std::decay_t<Types>>::value || ...);
+};
+
+template <typename T, typename... VecArgs>
+struct is_any_fp_scalar_impl<std::vector<T, VecArgs...>> {
+  static constexpr bool value = is_any_fp_scalar_impl<std::decay_t<T>>::value;
+};
+
+template <typename T>
+struct is_any_fp_scalar {
+  static constexpr bool value = is_any_fp_scalar_impl<std::decay_t<T>>::value;
+};
+
+}  // namespace
+
+TEST(MathPrimFunZeroedContainer, floating_point_predicate_filters_and_zeros) {
+  std::vector<double> a{
+      1.3,
+      -0.7,
+  };
+  Eigen::Matrix<double, -1, 1> b(2);
+  b << 2.0, -3.0;
+  std::vector<int> x_i{5, 6};
+  Eigen::VectorXd c{{4.0}};
+
+  auto args = std::make_tuple(a, std::make_tuple(b, x_i), c);
+  auto grads = stan::math::zeroed_filtered_tuple<is_any_fp_scalar>(args);
+
+  using grads_t = std::decay_t<decltype(grads)>;
+  EXPECT_EQ(3, std::tuple_size<grads_t>::value);
+
+  const auto& a_grads = std::get<0>(grads);
+  EXPECT_EQ(a.size(), a_grads.size());
+  for (std::size_t i = 0; i < a_grads.size(); ++i) {
+    EXPECT_FLOAT_EQ(0.0, a_grads[i]);
+  }
+
+  const auto& b_grads = std::get<0>(std::get<1>(grads));
+  EXPECT_EQ(b.size(), b_grads.size());
+  for (Eigen::Index i = 0; i < b_grads.size(); ++i) {
+    EXPECT_FLOAT_EQ(0.0, b_grads(i));
+  }
+
+  const auto& c_grads = std::get<2>(grads);
+  EXPECT_EQ(c.size(), c_grads.size());
+  for (Eigen::Index i = 0; i < c_grads.size(); ++i) {
+    EXPECT_FLOAT_EQ(0.0, c_grads(i));
+  }
+}
+
+TEST(MathPrimFunZeroedContainer, integral_predicate_filters_and_zeros) {
+  auto args = std::make_tuple(std::vector<int>{1, 2, 3}, 1.5,
+                              std::make_tuple(2.5, 7));
+  auto zeros = stan::math::zeroed_filtered_tuple<std::is_integral>(args);
+
+  using zeros_t = std::decay_t<decltype(zeros)>;
+  EXPECT_EQ(2, std::tuple_size<zeros_t>::value);
+
+  const auto& vec_zero = std::get<0>(zeros);
+  EXPECT_EQ(3, vec_zero.size());
+  EXPECT_EQ(0, vec_zero[0]);
+  EXPECT_EQ(0, vec_zero[1]);
+  EXPECT_EQ(0, vec_zero[2]);
+
+  const auto& nested_zero = std::get<0>(std::get<1>(zeros));
+  EXPECT_EQ(0, nested_zero);
+}
+
+TEST(MathPrimFunZeroedContainer, nested_tuple_values_are_zeroed) {
+  auto args = std::make_tuple(
+      1.25, std::make_tuple(2.5, -3.0, 1, std::make_tuple(4.75)));
+  auto zeros = stan::math::zeroed_filtered_tuple<std::is_floating_point>(args);
+
+  using expected_t
+      = std::tuple<double, std::tuple<double, double, std::tuple<double>>>;
+  static_assert(std::is_same_v<std::decay_t<decltype(zeros)>, expected_t>,
+                "zeroed_filtered_tuple should preserve nested tuple structure.");
+
+  EXPECT_FLOAT_EQ(0.0, std::get<0>(zeros));
+  EXPECT_FLOAT_EQ(0.0, std::get<0>(std::get<1>(zeros)));
+  EXPECT_FLOAT_EQ(0.0, std::get<1>(std::get<1>(zeros)));
+  EXPECT_FLOAT_EQ(0.0, std::get<0>(std::get<2>(std::get<1>(zeros))));
+}
+
+TEST(MathPrimFunZeroedContainer, tuple_containing_tuple_is_zeroed) {
+  auto args = std::make_tuple(std::make_tuple(3.14), 3.14, std::make_tuple(1));
+  auto zeros = stan::math::zeroed_filtered_tuple<std::is_floating_point>(args);
+
+  using expected_t = std::tuple<std::tuple<double>, double>;
+  static_assert(
+      std::is_same_v<std::decay_t<decltype(zeros)>, expected_t>,
+      "zeroed_filtered_tuple should support tuple<tuple<...>> inputs.");
+
+  EXPECT_FLOAT_EQ(0.0, std::get<0>(std::get<0>(zeros)));
+}
+
+TEST(MathPrimFunZeroedContainer, empty_tuple_returns_empty_tuple) {
+  auto zeros = stan::math::zeroed_filtered_tuple<std::is_integral>(std::make_tuple());
+  static_assert(std::tuple_size<std::decay_t<decltype(zeros)>>::value == 0,
+                "Expected empty tuple output for empty tuple input.");
+}
+
+TEST(MathPrimFunZeroedContainer, no_match_returns_empty_tuple) {
+  auto args = std::make_tuple(1.0, std::make_tuple(2.0, 3.0));
+  auto zeros = stan::math::zeroed_filtered_tuple<std::is_integral>(args);
+  static_assert(std::tuple_size<std::decay_t<decltype(zeros)>>::value == 0,
+                "Expected empty tuple when predicate matches no entries.");
+}
+
+TEST(MathPrimFunZeroedContainer, vector_of_tuples_preserves_shape) {
+  using vec_tuple_t = std::vector<std::tuple<double, int>>;
+  vec_tuple_t vec{{1.2, 3}, {-0.4, 5}};
+  auto args = std::make_tuple(vec);
+  auto zeros = stan::math::zeroed_filtered_tuple<is_any_fp_scalar>(args);
+
+  using expected_t = std::tuple<std::vector<std::tuple<double>>>;
+  static_assert(std::is_same_v<std::decay_t<decltype(zeros)>, expected_t>,
+                "Expected compact tuple shape for vector-of-tuples.");
+
+  const auto& vec_zero = std::get<0>(zeros);
+  ASSERT_EQ(2, vec_zero.size());
+  EXPECT_FLOAT_EQ(0.0, std::get<0>(vec_zero[0]));
+  EXPECT_FLOAT_EQ(0.0, std::get<0>(vec_zero[1]));
+}
diff --git a/test/unit/math/prim/functor/filter_types_test.cpp b/test/unit/math/prim/functor/filter_types_test.cpp
new file mode 100644
index 00000000000..6c5b3877315
--- /dev/null
+++ b/test/unit/math/prim/functor/filter_types_test.cpp
@@ -0,0 +1,41 @@
+#include <stan/math/fwd/core/fvar.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/prim/functor/filter_types.hpp>
+#include <stan/math/prim/meta/contains_autodiff.hpp>
+#include <gtest/gtest.h>
+#include <tuple>
+#include <type_traits>
+#include <vector>
+
+TEST(MathPrimFunctorFilterTypes, autodiff_predicate_preserves_nested_tuples) {
+  using stan::math::fvar;
+  std::vector<fvar<double>> a{
+      fvar<double>(1.0, 0.0),
+      fvar<double>(2.0, 0.0),
+  };
+  std::vector<int> b{3, 4};
+  fvar<double> c(5.0, 0.0);
+  int x = 1;
+  auto nested = std::forward_as_tuple(b, c);
+  auto args = std::forward_as_tuple(x, a, nested);
+
+  auto filtered = stan::math::filter_types<stan::contains_autodiff>(args);
+  using expected_t = std::tuple<std::vector<fvar<double>>&, std::tuple<fvar<double>&>>;
+  static_assert(std::is_same_v<decltype(filtered), expected_t>,
+                "filter_types should preserve nested tuple structure.");
+
+  std::get<0>(filtered)[0].val_ = 7.0;
+  EXPECT_FLOAT_EQ(7.0, a[0].val_);
+  std::get<0>(std::get<1>(filtered)).val_ = -9.0;
+  EXPECT_FLOAT_EQ(-9.0, c.val_);
+}
+
+TEST(MathPrimFunctorFilterTypes, no_match_returns_empty_tuple) {
+  int x = 1;
+  std::tuple<double, double> y{2.0, 3.0};
+  std::vector<int> z{4, 5};
+  auto args = std::forward_as_tuple(x, y, z);
+  auto filtered = stan::math::filter_types<stan::contains_autodiff>(args);
+  static_assert(std::tuple_size<std::decay_t<decltype(filtered)>>::value == 0,
+                "Expected empty tuple when predicate matches no entries.");
+}

From 0367f36701cf08ef8b86a6528b5453a131b0cfd9 Mon Sep 17 00:00:00 2001
From: Steve Bronder <sbronder@flatironinstitute.org>
Date: Fri, 20 Feb 2026 15:31:19 -0500
Subject: [PATCH 6/6] update

---
 stan/math/fwd/functor/finite_diff.hpp         | 13 ++--
 stan/math/prim/fun/zeroed_filtered_tuple.hpp  | 10 +---
 .../functor/finite_diff_gradient_auto.hpp     | 59 +++++++++----------
 3 files changed, 32 insertions(+), 50 deletions(-)

diff --git a/stan/math/fwd/functor/finite_diff.hpp b/stan/math/fwd/functor/finite_diff.hpp
index 939708b6d65..1261b6802da 100644
--- a/stan/math/fwd/functor/finite_diff.hpp
+++ b/stan/math/fwd/functor/finite_diff.hpp
@@ -10,7 +10,6 @@
 #include <stan/math/prim/functor/for_each.hpp>
 #include <stan/math/prim/functor/filter_types.hpp>
 #include <stan/math/prim/meta/contains_autodiff.hpp>
-#include <stan/math/prim/meta/filtered_tuple_indices.hpp>
 #include <tuple>
 
 namespace stan {
@@ -89,8 +88,9 @@ inline auto aggregate_tangent(FuncTangent&& tangent, InputArg&& arg) {
  * 1. Build `autodiff_args` with `filter_types<contains_autodiff>(args...)`.
  * 2. Build compact zero-initialized `grads` with
  *    `zeroed_filtered_tuple<contains_autodiff>(args...)`.
- * 3. Use filtered top-level index sequences to dispatch tuple-native finite
- *    differencing.
+ * 3. Dispatch tuple-native finite differencing with:
+ *    - `args_mask = args...` (original autodiff typing)
+ *    - `args_work = value_of(args)...` (mutable finite-diff work tuple)
  * 4. Aggregate tangent contributions by zipping compact `grads` and
  *    `autodiff_args`.
  *
@@ -106,16 +106,11 @@ inline auto finite_diff(const F& func, const TArgs&... args) {
   using FvarInnerT = typename FvarT::Scalar;
   auto autodiff_args
       = stan::math::filter_types<contains_autodiff>(std::forward_as_tuple(args...));
-  using autodiff_idxs_t
-      = filtered_tuple_indices_t<contains_autodiff, std::tuple<TArgs...>>;
   auto grads = stan::math::zeroed_filtered_tuple<contains_autodiff>(
       std::forward_as_tuple(args...));
-  using grad_idxs_t = std::make_index_sequence<
-      std::tuple_size<std::decay_t<decltype(grads)>>::value>;
   FvarInnerT rtn_value;
   finite_diff_gradient_auto(func, rtn_value, std::forward_as_tuple(args...),
-                            std::forward_as_tuple(value_of(args)...), grads,
-                            autodiff_idxs_t{}, grad_idxs_t{});
+                            std::forward_as_tuple(value_of(args)...), grads);
   FvarInnerT rtn_grad = 0;
   stan::math::for_each(
       [&rtn_grad](auto&& grad_i, auto&& arg_i) {
diff --git a/stan/math/prim/fun/zeroed_filtered_tuple.hpp b/stan/math/prim/fun/zeroed_filtered_tuple.hpp
index 7e5b73a38d2..09f90c7722c 100644
--- a/stan/math/prim/fun/zeroed_filtered_tuple.hpp
+++ b/stan/math/prim/fun/zeroed_filtered_tuple.hpp
@@ -13,14 +13,6 @@ namespace stan {
 namespace math {
 namespace internal {
 
-/**
- * Type-dependent false helper for `static_assert` in unreachable branches.
- *
- * @tparam ... ignored dependent types used to delay `static_assert`
- */
-template <typename...>
-struct zeroed_filtered_tuple_dependent_false : std::false_type {};
-
 /**
  * Build a zero-initialized container matching the shape of `arg`.
  *
@@ -55,7 +47,7 @@ inline auto zeroed_filtered_tuple_leaf(const Arg& arg) {
     return partial_t(0);
   } else {
     static_assert(
-        zeroed_filtered_tuple_dependent_false<arg_t>::value,
+        sizeof(std::decay_t<Arg>*) == 0,
         "Unsupported container in zeroed_filtered_tuple.");
   }
 }
diff --git a/stan/math/prim/functor/finite_diff_gradient_auto.hpp b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
index deaaa9746aa..16bf42f3f51 100644
--- a/stan/math/prim/functor/finite_diff_gradient_auto.hpp
+++ b/stan/math/prim/functor/finite_diff_gradient_auto.hpp
@@ -20,14 +20,6 @@ namespace stan {
 namespace math {
 namespace internal {
 
-/**
- * Type-dependent false helper for `static_assert` in unreachable branches.
- *
- * @tparam ... ignored dependent types used to delay `static_assert`
- */
-template <typename...>
-struct dependent_false : std::false_type {};
-
 /**
  * Construct a holder tuple containing elements selected by `Is...`.
  *
@@ -67,6 +59,7 @@ inline double finite_diff_stepsize_value(const T& x) {
   }
 }
 
+
 /**
  * Recursively traverses `arg_mask`, `arg_work`, and `grad` and invokes `fn` on
  * each perturbable coordinate pair `(arg_work_coord, grad_coord)`.
@@ -145,11 +138,10 @@ inline void for_each_leaf_pair(ArgMask&& arg_mask, ArgWork&& arg_work,
           [&fn](auto&& arg_mask_i, auto&& arg_work_i, auto&& grad_i) {
             for_each_leaf_pair(arg_mask_i, arg_work_i, grad_i, fn);
           },
-          std::move(arg_mask_ad), std::move(arg_work_ad),
-          std::forward<Grad>(grad));
+          arg_mask_ad, arg_work_ad, std::forward<Grad>(grad));
     }
   } else {
-    static_assert(dependent_false<arg_mask_t>::value,
+    static_assert(sizeof(std::decay_t<arg_mask_t>*) == 0,
                   "Unsupported container in finite_diff_gradient_auto.");
   }
 }
@@ -275,13 +267,12 @@ inline void finite_diff_gradient_auto(const F& f, VectorT&& x, ScalarT& fx,
  * while preserving argument scalar types at the call boundary.
  *
  * This overload expects `grads` to mirror the top-level structure of `args`.
- * Non-autodiff top-level entries are skipped in the finite-difference pass.
+ * Non-autodiff leaves are skipped in the finite-difference pass.
  *
  * Pattern used by this overload:
  * 1. Evaluate `f(args...)` once for `fx`.
  * 2. Zip top-level `args` and `grads`.
- * 3. Use `contains_autodiff_v` to skip non-autodiff top-level entries.
- * 4. Traverse nested leaves with `for_each_leaf_pair`, perturbing one
+ * 3. Traverse nested leaves with `for_each_leaf_pair`, perturbing one
  *    coordinate at a time and restoring it after each stencil evaluation.
  *
  * @tparam F callable type
@@ -303,7 +294,7 @@ inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
   fx = stan::math::apply([&f](auto&&... args_i) { return f(args_i...); }, args);
   stan::math::for_each(
       [&f, &args](auto&& arg_i, auto&& grad_i) {
-        if constexpr (contains_autodiff_v<std::decay_t<decltype(arg_i)>>) {
+        if constexpr (!std::is_integral_v<scalar_type_t<decltype(arg_i)>>) {
           internal::finite_diff_gradient_auto_impl<ScalarT>(f, arg_i, arg_i,
                                                             grad_i, args);
         }
@@ -312,9 +303,7 @@ inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
 }
 
 /**
- * Indexed tuple overload used for compact gradient tuples. `Idxs` select
- * autodiff-containing top-level entries from both `args_mask` and `args_work`,
- * while `GradIdxs` selects matching compact gradient entries from `grads`.
+ * Tuple overload used for compact gradient tuples with an explicit mask tuple.
  *
  * This overload is used by `stan::math::finite_diff` where:
  * - `args_mask` preserves original autodiff type information
@@ -332,30 +321,36 @@ inline void finite_diff_gradient_auto(const F& f, ScalarT& fx, TupleArgs&& args,
  * @param[in,out] args_work mutable value tuple used for finite-difference
  * evaluations
  * @param[out] grads compact tuple of gradient containers
- * @param[in] idxs selected top-level indices in `args_mask`/`args_work`
- * @param[in] grad_idxs selected top-level indices in `grads`
  * @return None.
  * @throw Any exception thrown by `f`.
  */
 template <typename F, typename ScalarT, typename TupleArgsMask,
-          typename TupleArgsWork, typename TupleGrads, std::size_t... Idxs,
-          std::size_t... GradIdxs, require_tuple_t<TupleArgsMask>* = nullptr,
+          typename TupleArgsWork, typename TupleGrads,
+          require_tuple_t<TupleArgsMask>* = nullptr,
           require_tuple_t<TupleArgsWork>* = nullptr,
           require_tuple_t<TupleGrads>* = nullptr>
 inline void finite_diff_gradient_auto(
     const F& f, ScalarT& fx, TupleArgsMask&& args_mask, TupleArgsWork&& args_work,
-    TupleGrads&& grads, [[maybe_unused]] std::index_sequence<Idxs...> idxs,
-    [[maybe_unused]] std::index_sequence<GradIdxs...> grad_idxs) {
-  static_assert(sizeof...(Idxs) == sizeof...(GradIdxs),
-                "Index sequence size mismatch in finite_diff_gradient_auto.");
+    TupleGrads&& grads) {
   fx = stan::math::apply([&f](auto&&... args_i) { return f(args_i...); },
                          args_work);
-  using Swallow = int[];
-  static_cast<void>(Swallow{
-      (static_cast<void>(internal::finite_diff_gradient_auto_impl<ScalarT>(
-           f, std::get<Idxs>(args_mask), std::get<Idxs>(args_work),
-           std::get<GradIdxs>(grads), args_work)),
-       0)...});
+  using selected_idxs_t
+      = filtered_tuple_indices_t<contains_autodiff, std::decay_t<TupleArgsMask>>;
+  auto args_mask_ad = internal::tuple_subset(args_mask, selected_idxs_t{});
+  auto args_work_ad = internal::tuple_subset(args_work, selected_idxs_t{});
+  static_assert(std::tuple_size<std::decay_t<decltype(args_mask_ad)>>::value
+                        == std::tuple_size<std::decay_t<decltype(args_work_ad)>>::value
+                    && std::tuple_size<std::decay_t<decltype(args_mask_ad)>>::value
+                           == std::tuple_size<std::decay_t<TupleGrads>>::value,
+                "Tuple size mismatch in finite_diff_gradient_auto.");
+  stan::math::for_each(
+      [&f, &args_work](auto&& arg_mask_i, auto&& arg_work_i, auto&& grad_i) {
+        internal::finite_diff_gradient_auto_impl<ScalarT>(
+            f, std::forward<decltype(arg_mask_i)>(arg_mask_i),
+            std::forward<decltype(arg_work_i)>(arg_work_i),
+            std::forward<decltype(grad_i)>(grad_i), args_work);
+      },
+      args_mask_ad, args_work_ad, std::forward<TupleGrads>(grads));
 }
 
 }  // namespace math