title: Compute x**k with x, k being arrays of arbitrary dimensionality
body: <p>I have two numpy arrays: One array <code>x</code> with shape <code>(n, a0, a1, ...)</code> and one array <code>k</code> with shape <code>(n, b0, b1, ...)</code>. I would like to compute and array of exponentials such that the output has dimension <code>(a0, a1, ..., b0, b1, ...)</code> and </p>
  
  
  
  <pre><code>out[i0, i1, ..., j0, j1, ...] == prod(x[:, i0, i1, ...] ** k[:, j0, j1, ...])
  </code></pre>
  
  <p>If there is only one <code>a_i</code> and one <code>b_j</code>, broadcasting does the trick via</p>
  
  
  
  <pre><code>import numpy
  
  x = numpy.random.rand(2, 31)
  k = numpy.random.randint(1, 10, size=(2, 101))
  
  out = numpy.prod(x[..., None]**k[:, None], axis=0)
  </code></pre>
  
  <p>If <code>x</code> has a few dimensions more, more <code>None</code>s have to be added:</p>
  
  
  
  <pre><code>x = numpy.random.rand(2, 31, 32, 33)
  k = numpy.random.randint(1, 10, size=(2, 101))
  
  out = numpy.prod(x[..., None]**k[:, None, None, None], axis=0)
  </code></pre>
  
  <p>If <code>x</code> has a few dimensions more, more <code>None</code>s have to be added at other places:</p>
  
  
  
  <pre><code>x = numpy.random.rand(2, 31)
  k = numpy.random.randint(1, 10, size=(2, 51, 51))
  
  out = numpy.prod(x[..., None, None]**k[:, None], axis=0)
  </code></pre>
  
  <p>How to make the computation of <code>out</code> generic with respect to the dimensionality of the input arrays?</p>