Use linearly indexed array math only when ndims match

[jishnub]

This was broken in #59961, as map deals with trailing singleton axes differently from broadcasting:

julia> map(+, ones(1), ones(1,1)) |> size
(1,)

julia> broadcast(+, ones(1), ones(1,1)) |> size
(1, 1)

This PR limits the new method to the case where the ndims match, in which case there are no trailing axes and the two are equivalent. The alternate approach suggested in #59961 (comment) is to reshape the arrays, but this adds overhead that nullifies the performance improvement for small arrays.

[adienes]

sometimes 1d is a special case; could you add one or two more test cases comparing 2d to 3d? (and maybe similarly some for 0d which is also often special cased)

[adienes]

seems kind of a shame to have performance depend on the presence of otherwise inconsequential trailing dimensions but I don't have a better proposal