reduce_log_sum_exp function, aka reduce_sum for mixtures

I have a situation where I'd like use something like reduce_sum, but it's for marginalizing out a discrete variable, so the reduction is log-sum-exp rather than just sum.

The general pattern for marginalizing out a discrete variable is this:

vector[K] lp;
for (z in 1:K)
  lp[z] = log p(z) + log p(y | z, x, theta);
target += log_sum_exp(lp);

It'd be nice to have a reduce_log_sum_exp function that works just like reduce_sum, only with log_sum_exp rather than sum as the final reduction before adding to target.

The log-sum-exp gradient itself can be done analytically as part of the reduction, which will save a whole bunch of storage relative to just using map_rect and doing the log-sum-exp on the result.

Start by tracing the existing reduce_sum and map_rect entry points in stanc3, then inspect how log_sum_exp is represented and differentiated. No files or tests are named in the issue, so identify the relevant compiler paths and add coverage showing that reduce_log_sum_exp marginalizes the discrete-variable terms and contributes correctly to target.