Implement Prod.R_Op

[lamblin]

I think it can be implemented by taking the L_Op with an output grad of 1, to get the Jacobian in a compact form, and then doing an elemwise mul with the evaluation point, then summing over the same axes as the original prod.

Using L_Op would make it possible to benefit from the special-casing of 0 that we already to with ProdWithoutZeros.

The mul-sum could also be replaced by a dot or tensordot, but it could be complex to set up depending on the reduction axes.

Reported in #5792.