Hi,

I would like to do something like the outer product for general operations. Basically, I f I have a function

```
f(x1, x2, [...], xn) = [...]
```

and I provide it with n arrays `a1`

… `an`

as input, then I would like to apply `f`

on all combination of input elements, such that the resulting array has the shape:

```
outer(f, a1, a2, [...], an) |> size == (size(a1)..., size(a2)..., [...], size(an)...)
```

*How can I do that in a clean and performant way?*

*Example:*

```
julia> f(a, b, c) = a+b+c
f (generic function with 1 method)
julia> x = 1:300; y = 1:400; z = 1:500;
# performant but super hard to read
julia> f.(x, reshape(y, (ones(Int64, ndims(x))..., size(y)...)), reshape(z, (ones(Int64, ndims(x) + ndims(y))..., size(z)...))) |> size == (size(x)..., size(y)..., size(z)...)
true
julia> @benchmark f.(x, reshape(y, (ones(Int64, ndims(x))..., size(y)...)), reshape(z, (ones(Int64, ndims(x) + ndims(y))..., size(z)...)))
BenchmarkTools.Trial:
memory estimate: 457.76 MiB
allocs estimate: 13
--------------
minimum time: 141.622 ms (0.29% GC)
median time: 150.833 ms (9.09% GC)
mean time: 158.361 ms (10.34% GC)
maximum time: 193.617 ms (22.94% GC)
--------------
samples: 32
evals/sample: 1
# easy to read but slow
julia> broadcast( u -> f(u...), Base.Iterators.product(x, y, z)) |> size == (size(x)..., size(y)..., size(z)...)
true
julia> @benchmark broadcast( u -> f(u...), Base.Iterators.product(x, y, z))
BenchmarkTools.Trial:
memory estimate: 1.79 GiB
allocs estimate: 6
--------------
minimum time: 731.471 ms (1.30% GC)
median time: 799.173 ms (13.66% GC)
mean time: 792.291 ms (11.79% GC)
maximum time: 820.479 ms (11.84% GC)
--------------
samples: 7
evals/sample: 1
```

The first version seems performant, the second one seems clean. Is there a standard function that I am overlooking?