How to treat `out` arguments and in-place operations?

The Array API does not specify `out` arguments, but numpy ufuncs do.

A question is to what extent we should support it. For inplace operands like `__imul__`, at the moment we return a new `Quantity`, but with `value` shared with the old one (and a new unit), so
```
t = q1
q1 *= q2
t is q1
# False
t.value is q1.value
# True
t == q1
# depends on unit of q2; True if dimensionless, False otherwise.
```

We could make it similar for functions, i.e.,
```
r = np.multiply(q1, q2, out=q1)
r is q1
# False
r.value is q1.value
# True
r == q1
# depends on unit of q2
```
But this does break expectations one might have from numpy, where `r is out` always holds.

The alternative would be to be like `astropy`, where we allow the unit to be changed. But that breaks the immutability of `Quantity`.

p.s. The above also assumes that the unit of the `out` argument should be ignored. This is probably best (and what we do in astropy), but one could argue that since the data type of `out` is not ignored - the result is cast to it - perhaps the unit shouldn't either? Though if we do that, then in-place multiplication becomes impossible (as would, e.g., `np.sin(q, out=q)`), which seems strange. So, probably fine with changing unit, the question is just whether a new instance is returned.

Begin with Quantity.__imul__ and the NumPy multiply(out=...) behavior described in the issue; compare the current identity, value-sharing, and unit semantics. The issue names no files or tests, so first establish the relevant implementation and test entry points, then define and verify a single consistent policy for out and in-place operations.

Written by the indexing model from the issue text.