Generate gufunc input shapes with a gufunc= arg to mutually_broadcastable_shapes

[Zac-HD]

This is the tracking issue for a long-awaited feature: generating a collection of array shapes for a Numpy 'generalised ufunc'. It's fairly niche, but very useful to anyone in that niche.

As background, a strategy for gufunc_shapes was first proposed in #1760/#1784 and later released as the hypothesis-gufunc extension. Around the same time, #1829 proposed our boradcastable_shapes strategy (for single shapes), and we later decided (#1969) to release mutually_broadcastable_shapes without gufunc support in the first instance. That was added by #2153/#2173, so we're finally ready for gufunc support!

---

The implementation isn't too bad, but there are a bunch of little details we'll need to handle carefully.

• The gufunc argument is keyword-only, defaults to not_set, and exactly one of it or num_shapes must be passed.
• If a gufunc callable is passed, use its .signature attribute. IMO we should reject None signatures, since the ufunc is not generalised and it simply gives us num_shapes.
• Parse the gufunc signature string using numpy.lib.function_base._parse_gufunc_signature something which is actually NEP-20 compliant, i.e. our own regex-based parser.
• When drawing, we sort the size names (m and n in '(m,n),(n)->(m)' for matrix-vector multiplication), then draw an allowed size for each, and construct our initial tuples for each shape - the "core dimensions".
• Also handle fixed-size dimensions and optional dimensions and broadcastable dimensions (NEP rejected). Fiddly but no worse than our existing logic for other dimensions.
• min_dims and max_dims are interpreted (and documented!) as applying to the "loop dimensions" only. In the example above, max_dims=1 would mean that we could generate a three-dimensional output (loop,m,n) and the result shape would be (loop,m).
• The exception is that we (silently) stop adding to a shape when it reaches 32 dimensions even if this means that it will have fewer than min_dims dimensions. We should also adjust the biased_coin probability in these cases to avoid squashing our uniform distribution against the upper bound.
  If for example the cap were two dimensions, it would be impossible to use any loop dimensions or have a 2D result from matrix-vector multiplication without supporting this (or ugly per-dimension size limits). It's very unlikely to arise in practice, but if Hypothesis doesn't handle edge cases who will?
  I believe that optional dimensions (eg in matmul (m?,n),(n,p?)->(m?,p?)) do count towards the limit of 32 dimensions even if unused, but we should check this and generate them if not.

We'll also want to make sure it's well-tested and thoroughly documented before merging, but I think that's everything we need in for a viable implementation 😄

[Zac-HD]

For @rsokl - parsing signatures is mostly done, and that's the hard part... master...Zac-HD:ufunc-shapes

[rsokl]

You read my mind! I was going to ask you to ping me if you had been cooking up the parser 😄