Add an inverse-CDF path to categorical sampling

[deescondrillas]

Proposed changes

Fixes #3847 for the single-distribution case.

categorical samples with Gumbel-max, drawing one perturbation per category per sample, so
N categories and M samples build an N x M intermediate to return M indices. At
N = M = 1e6 that is ~4 TB.

When logits holds a single distribution and more than one sample is requested, sample by
inverting the cumulative weights instead: memory goes from O(N * M) to O(N + M) and work
from O(N * M) to O(N + M log N). It is composed from existing ops, so CPU, Metal and CUDA
are covered with no new kernels. searchsorted needs a 1-D sequence, so batched logits
and the single-sample overload keep the Gumbel path.

Weights accumulate in float32, so a category holding less than roughly 1e-7 of the
total mass may be drawn with a distorted probability. This is noted in the docstring.

Behaviour changes

• For a fixed key the affected overloads return different, identically distributed samples.
• All -inf logits now draw uniformly rather than always returning index 0.
• Several +inf entries now split the samples rather than always returning the first.

Benchmarks

benchmarks/python/categorical_bench.py, Metal:

N	M	inverse-CDF	Gumbel
10⁴	10⁴	0.28 ms / 0.6 MB	97.8 ms / 400 MB
5·10⁴	10⁴	0.53 ms / 1.0 MB	516 ms / 2000 MB
10⁶	10⁶	3.2 ms / 25 MB	needs ~8 TB

CPU is 1233x and 4616x on the first two rows.

C++ 25/25 and Python 21/21 on Metal and DEVICE=cpu. CUDA was not built or tested (Apple
silicon); the change adds no backend code.

Checklist

Put an x in the boxes that apply.

• I have read the CONTRIBUTING document
• I have run pre-commit run --all-files to format my code / installed pre-commit prior to committing changes
• I have added tests that prove my fix is effective or that my feature works
• I have updated the necessary documentation (if needed)