<snip>

Wow, really nice work!  It would be great if that could make into NumPy
:-)  Regarding your comment on numexpr being faster, I'm not sure (your
new_iterator branch does not work for me; it gives me an error like:
AttributeError: 'module' object has no attribute 'newiter'),

but my
guess is that your approach seems actually faster:

>>> a = np.random.random((50,50,50,10))

>>> b = np.random.random((50,50,1,10))

>>> c = np.random.random((50,50,50,1))

>>> timeit 3*a+b-(a/c)

10 loops, best of 3: 67.5 ms per loop
>>> import numexpr as ne
>>> ne.evaluate("3*a+b-(a/c)
>>> timeit ne.evaluate("3*a+b-(a/c)")
10 loops, best of 3: 42.8 ms per loop

i.e. numexpr is not able to achieve the 2x speedup mark that you are
getting with ``luf`` (using a Core2 @ 3 GHz here).