ENH: add an axis argument to np.bincount

This adds an `axis` argument to `np.bincount`, which can now take multidimensional
arrays and do the counting over an arbitrary number of axes. `axis` defaults
to `None`, which does the counting over all the axes, i.e. over the flattened
array.

The shape of the output is computed by removing from `list` (the input array)
all dimensions in `axis`, and appending a dimension of length
`max(max(list) + 1, minlength)` to the end. `out[..., n]` will hold the number
of occurrences of `n` at the given position over the selected axes.

If a `weights` argument is provided, its shape is broadcasted with `list`
before removing the axes. In this case, `axis` refers to the axes of `list`
*before* broadcasting, and `out[..., n]` will hold the sum of `weights` over
the selected axes, at all positions where `list` takes the value `n`.

The general case is handled with nested iterators, but shortcuts without
having to set up an iterator are provided for 1D cases, with no performance
loss against the previous version.

As a plus, this PR also solves #823, by providing specialized functions for
all integer types to find the max. There are also specialized functions for
all integer types for counting and doing weighted counting.