Use ArrayD as the return value of operator overloading functions

[SparrowLii]

In the current operator overloading function( which in impl_ops.rs ), we always use the dimension of the first operand as the dimension of the return value. This is obviously inconsistent with the facts, and this is the fundamental reason why we cannot achieve co-broadcasting.
I think the greatest significance of operator overloading and co_broadcast is to simplify the user's use, for which we can appropriately sacrifice a little efficiency. In this case, we should not ask users to ensure that the dimensions of the operands are the same. ( Otherwise they can use zip_mut_with directly).
So co_broadcast is what we need to implement in operator overloading. And because we cannot determine the dimension of its return value, using ArrayD is an inevitable trend.
( Or, we can implement separate operator overloading for each Dim (Ix0 ~ IxDyn) through procedural macros, but this will increase the amount of code dozens of times )

[jturner314]

We can implement co-broadcasting without resorting to IxDyn. I worked on co-broadcasting a few years ago with a couple of different approaches but didn't get a chance to finish it up. Basically, to determine the output dimension type, you can introduce a trait with two type inputs and an associated type output. (For a similar trait, consider std::ops::Add, which defines the Output type for two parameter types.) Then, you can have a function which computes the shape when two shapes are co-broadcasted together. Once the new shape is determined, you can broadcast each of the arrays with .broadcast(). (It would be nice to support the mutable case as well; this requires some additional checks for aliasing.) Then, the desired operation can be performed.

See also #11.

[bluss]

Dynamic dimensional arrays are so inefficient that it's embarassing