Corner case in `dirichlet_lpdf`: error when both alpha - 1 = theta = 0

In the code there is currently:

  const auto& theta_log
      = to_ref_if<!is_constant_all<T_prior_size>::value>(theta_dbl.log());

  if (include_summand<propto, T_prob, T_prior_size>::value) {
    lp += (theta_log * alpha_m_1).sum();
  }

I understand that by the convention (and also by Stan documentation) for Dirichlet distribution, we usually require theta > 0. But when you think about it, for the corner case when alpha -1 == theta == 0, it is reasonable to define the lpdf to be 0 by the same logic why we define pow(0,0)=1. This is relevant in my application, where I may have on input an information that certain index of theta is a forbidden state and may thus be represented as hard constrained to be equal to 0, but then the Dirichlet prior crashes no matter what. This would be very inefficient for me to workaround by deleting the offending array element as I need to do this over a matrix with many rows each of which may have different constraints, so currently I needed to make my custom implementation of dirichlet_lpdf just to deal with that. I thus propose we change the code to something like:

  if (include_summand<propto, T_prob, T_prior_size>::value) {
    lp += lmultiply(alpha_m_1, theta).sum();
  }

The partials_vec<1> may need a corresponding adjustment.

Start in stan/math/prim/prob/dirichlet_lpdf.hpp at the linked theta_log calculation, then inspect the corresponding partials_vec<1> handling. Determine whether the alpha - 1 = theta = 0 case should contribute zero without producing an error, and ensure the associated autodiff behavior is consistent.