Add functions to write hessian vector products into an array

The functions

1. stan::math::internal::finite_diff_hessian_times_vector_auto()
2. stan::math::hessian_times_vector()

accept the hessian vector product (hvp or Hv, respectively) as an Eigen::VectorXd& out parameter. Similar to issue #2739, it would benefit BridgeStan to allow signatures that effectively write into double* hvp such that BridgeStan could offer cheaper hessian vector products.

I propose adding functions with signatures

1. stan::math::internal::finite_diff_hessian_times_vector_auto()

template <typename F, typename EigVec, typename InputIt, 
  require_eigen_vector_vt<std::is_arithmetic, EigVec>* = nullptr>
void finite_diff_hessian_times_vector_auto(const F& f, const EigVec& x,
                                           const EigVec& v, double& fx, InputIt first_hvp, InputIt last_hvp)

2. stan::math::hessian_times_vector()

template <typename F, typename EigVec, typename InputIt, 
  require_eigen_vector_vt<std::is_arithmetic, EigVec>* = nullptr>
void hessian_times_vector(const F& f,
                          const EigVec& x,
                          const EigVec& v,
                          double& fx,
                          InputIt first_Hv, InputIt last_Hv) 

In fact, stan::math::hessian could benefit from a similar added function, but I'm having a hard time wrapping my head around the design of the signature since the out parameter here is a matrix. @WardBrian, would you mind lending an opinion?

Current Version:

v4.7.0

Start with stan/math/rev/functor/finite_diff_hessian_times_vector_auto.hpp and stan/math/mix/functor/hessian_times_vector.hpp, then compare the existing APIs with the BridgeStan call in src/model_rng.cpp and the approach from issue #2739. Confirm how iterator output ranges should coexist with the current Eigen::VectorXd outputs, and resolve the open design question for hessian before documenting completion criteria.