Provide an exact internal_method for da.percentile #12686

Dask has no exact percentile for a chunked 1-D array, and da.median and da.quantile raise for axis=None on chunked arrays. da.percentile is approximate by design (#3099). There are existing algorithms that compute these statistics exactly out of core. In particular, multi-pass radix selection, which is what CUDA torch.kthvalue and the radix variants of RAFT's select_k use (Alabi et al., ACM JEA 2012), does it in at most four passes over the input for 64-bit dtypes, two for 32-bit and one for 16-bit, with one 65,536-bin int64 histogram per tracked rank per task. Each pass is a map over the chunks plus a sum, and the next pass depends only on the summed histogram, so chunks are re-read rather than retained. It should be possible to match np.percentile(a.compute(), q, method=method), with NaN propagating as in NumPy.

I'd propose this as an opt-in internal_method="exact" for da.percentile, leaving the default alone (#6566 is about changing it). It would be 1-D only like the existing methods; quantile(axis=None) on N-D input is the same computation on raveled chunks and could follow. If there's interest, I would be happy to put together a PR.