vectorizing a method that operates on an array

I have a pybind11ified function that operates on a vector of dimension 1,

#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

namespace py = pybind11;

void distill(py::array_t<double, py::array::c_style | py::array::forcecast> p) {
  auto r = p.mutable_unchecked<1>();
  for (ssize_t i = 1; i < r.shape(0); i++) {
    r(i) = 3.14 * (r(i) + r(i-1));  // whatever
  }
}

PYBIND11_MODULE(_accupy, m) {
  m.def("distill", &distill);
}

I would now like to change the function to have it operate on the first dimension of numpy arrays of arbitrary dimensionality. The first thing that comes to mind is adding an inner loop

for (ssize_t j = 0; j < r.shape(1); j++) {
}

to the above and wrapping the method in Python code à la

def distill(p):
    q = p.reshape(p.shape[0], numpy.prod(p.shape[1:]))
    out = _mymod.distill(q)
    return out.reshape(p.shape)

This however wouldn't work with arrays of one dimension anymore. Also, instead of for (ssize_t j = 0; j < r.shape(1); j++) {} it's probably better to use BLAS's multiply-add functions there; a rabbit hole I'd rather not descend into.

Is there a canonical way for vectorizing pybind11 methods that operate on an array?

Begin with the py::array_t and mutable_unchecked<1>() usage in the issue's distill example, then compare it with the reshape-based Python wrapper. Determine a canonical approach that handles both one-dimensional and arbitrary-dimensional NumPy arrays and defines the expected result for each shape; no test or source file is named in the issue.