Staged bsxfun and other broadcast operations

This patch adds an implementation of `bsxfun(op::Function, As::Array...)`
that generates inner loop functions depending on `op`, the number of arguments
`As`, and the dimension of the result.

The machinery is also used to provide
 * broadcasting operators `.+ .- .* ./` for `Array` arguments
 * `bsxfun!(op::Function, As::Array...)`, which stores the result in
   `As[1]` **Edit:** Removed this one, pending a rewrite with separate output argument.
 * `get_bsxfun(op::Function)` and `get_bsxfun!(op::Function)` to instantiate a
   broadcasting operation from a given per-element operation
 * `bsx_getindex(A::AbstractArray, inds::Array...)`,
   a pointwise broadcasting indexing operation as described in #2591,
   and a matching `bsx_setindex!(A::AbstractArray, X::Array, inds::Array...)`

It should be straightforward to extend the functions to work with `StridedArray` arguments instead of just `Array`, but I'm uncertain of the proper interface to work with subarrays using manual linear indexing.

The code generates inner loop functions such as

    function bsx_plus_inner!(shape::NTuple{2,Int}, arrays::NTuple{2,Array},
                             strides::Matrix{Int}, result::Array)
        @assert (size(strides)==(2,2))
        (n1,n2) = shape
        if (n1==0); return; end
        if (n2==0); return; end
        (A1,A2) = arrays
        k1 = 1
        k2 = 1
        (s1_1,s2_1,s1_2,s2_2) = strides
        k = 1
        for i2 = 1:n2
            for i1 = 1:n1
                begin 
                    result[k] = +(A1[k1],A2[k2])
                    k += 1
                end
                k1 += s1_1
                k2 += s2_1
            end
            k1 += s1_2
            k2 += s2_2
        end
    end

The example above is for `bsxfun(+, A, B)` with two-dimensional result.
A linear index is updated using strides for each array, which should allow to accomodate
both broadcasting and subarrays.

On my machine, the timing example from #2991 gives

    elapsed time: 0.806272502 seconds  # subcol_forloop(a, b)
    elapsed time: 3.1062779 seconds    # subcol_bsxfun(a, b)
    elapsed time: 0.936735113 seconds  # .- from this pull request

i.e. the generated `.-` is only 16% slower than the handcoded example in #2991.

It's a somewhat complex piece of code (that generates code that generates code), and I think that it warrants some discussion. I could polish this forever, but I think that I need some reactions before moving forward. What do you guys think? Would something like this be useful? What parts of it?

**Edit:** Renamed `bsxfun` and `bsx` to `broadcast`.