Modified order of operations in _approximate_mode() changes function behavior

### Describe the bug

In issue [#20774](https://github.com/scikit-learn/scikit-learn/issues/20774) and subsequent [PR #20904](https://github.com/scikit-learn/scikit-learn/pull/20904), the order of operations in the calculation of `continuous` was changed from
```
continuous = n_draws * class_counts / class_counts.sum()
```
to
```
continuous = class_counts / class_counts.sum() * n_draws
```

For cases where the class count values are low and number of classes is high, deficiency in float precision results in the `continuous` array being substantively different in these two arrangements and changes both the output of both this function and the subsequent splitting operations for which this function is used.  So while this change may be "fine" in the long-run, it was definitely not an entirely transparent change.

### Steps/Code to Reproduce

```
from sklearn.utils import check_random_state
import numpy as np

def _approximate_mode_0242(class_counts, n_draws, rng):
    rng = check_random_state(rng)
    continuous = n_draws * class_counts / class_counts.sum()
    floored = np.floor(continuous)
    need_to_add = int(n_draws - floored.sum())
    if need_to_add > 0:
        remainder = continuous - floored
        values = np.sort(np.unique(remainder))[::-1]
        for value in values:
            inds, = np.where(remainder == value)
            add_now = min(len(inds), need_to_add)
            inds = rng.choice(inds, size=add_now, replace=False)
            floored[inds] += 1
            need_to_add -= add_now
            if need_to_add == 0:
                break
    return floored.astype(int)


def _approximate_mode_102(class_counts, n_draws, rng):
    rng = check_random_state(rng)
    continuous = class_counts / class_counts.sum() * n_draws  # note the order of operations
    floored = np.floor(continuous)
    need_to_add = int(n_draws - floored.sum())
    if need_to_add > 0:
        remainder = continuous - floored
        values = np.sort(np.unique(remainder))[::-1]
        for value in values:
            inds, = np.where(remainder == value)
            add_now = min(len(inds), need_to_add)
            inds = rng.choice(inds, size=add_now, replace=False)
            floored[inds] += 1
            need_to_add -= add_now
            if need_to_add == 0:
                break
    return floored.astype(int)


class_counts = np.array([
    119,6,6,6,6,6,6,6,6,6,6,6,6,
    6,6,6,6,6,5,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,6,6,6,6,3,
    3,3,3,3,3,3,3,3,3,3,3,3,3,
    3,3,3,3,3,3,6,3,3,3,3,3,3,
    2,3,3,3,3,3,3,3,3,3,3,3,3,
    3,3,3,3,3,3,3,3,3,3,3,3,3,
    3,3,3,3,3,3,3,3,3,3,3,3,3,
    3,3,3,3,3,3,3,6,6,6,6,6,6,
    6,3,3,3,3,6,3,3,3,3,3,3,3,
    3,3,3,3,3,3,3,3,3,3,3,6,6,
    6,3,3,3,3,3,3,3,6,6,6,6,3,
    3,3,3,3,6,3,3,3,3,3,6,6,6,
    6,6,3,3,6,3,3,6,6,6,6,6,6,
    3,6,5,6,3,3,3,3,3,6,6,6,6,
    6,6,6,6,6,3,6,6,3,3,6,6,6,
    6,6,6,6,6,6,6,6,6,6,3,6,6,
    6,6,6,6,3,3,3,6,6,6,6,6,6,
    3,3,3,3,6,6,5,6,6,6,6,6,5,
    6,6,6,6,6,6,2,3,3,3,3,3,3,
    3,3,3,3,3,3,3,3,3,3,3,3,3,
    3,3,3,3,3,3,3,3,3,3,3,3,3,
    3,3,3,3,3,3,3,3,3,3,3,3])
n_draws = 682

test1 = _approximate_mode_0242(class_counts, n_draws, 123)
test2 = _approximate_mode_102(class_counts, n_draws, 123)
print(np.unique(test1 == test2, return_counts=True))
```
```
# 0.24.2 version
continuous = n_draws * class_counts / class_counts.sum()
floored = np.floor(continuous)
remainder = continuous - floored
remainder[:5]

>> array([0.5, 0. , 0. , 0. , 0. ])
```
```
# 1.0.2 version
continuous = class_counts / class_counts.sum() * n_draws
floored = np.floor(continuous)
remainder = continuous - floored
remainder[:5]

>> array([5.0000000e-01, 4.4408921e-16, 4.4408921e-16, 4.4408921e-16, 4.4408921e-16])
```
Notably, it's worth considering whether ^ IS an appropriate outcome, given the functions reliance on `np.unique()`.  The reason is that we should be sampling simultaneously from all of the stuff whose remainder (in exact arithmetic) would be 0.5, whereas in the new version we are sampling first from a subset, then another subset, then another subset. The loss of precision is deterministic, meaning that there could be systematic biases in which classes get sampled based on their subtly different remainders.

### Expected Results

```
test1 = _approximate_mode_0242(class_counts, n_draws, 123)
test2 = _approximate_mode_102(class_counts, n_draws, 123)
print(np.unique(test1 == test2, return_counts=True))

>> (array([True]), array([ 298]))
```

### Actual Results

```
test1 = _approximate_mode_0242(class_counts, n_draws, 123)
test2 = _approximate_mode_102(class_counts, n_draws, 123)
print(np.unique(test1 == test2, return_counts=True))

>> (array([False,  True]), array([ 76, 222]))
```

### Versions

System:
    python: 3.9.7 (default, Sep  8 2021, 18:24:04)  [Clang 12.0.5 (clang-1205.0.22.11)]
   machine: macOS-11.3-x86_64-i386-64bit

Python dependencies:
          pip: 22.0.4
   setuptools: 57.4.0
      sklearn: 0.24.2
        numpy: 1.22.3
        scipy: 1.8.0
       Cython: None
       pandas: 1.1.5
   matplotlib: 3.4.3
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True

and 

System:
    python: 3.9.7 (default, Sep  8 2021, 18:24:04)  [Clang 12.0.5 (clang-1205.0.22.11)]
   machine: macOS-11.3-x86_64-i386-64bit

Python dependencies:
          pip: 21.2.3
   setuptools: 57.4.0
      sklearn: 1.0.2
        numpy: 1.22.3
        scipy: 1.8.0
       Cython: None
       pandas: 1.1.5
   matplotlib: 3.4.3
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True

COMMENT:
I agree the change was not transparent and we should update our `Change Models` section to state that `train_test_split` and `StraitfiedKFold` may have different splits.

In this case, I think `_approximate_mode`'s docstring of not being off by one still holds:

https://github.com/scikit-learn/scikit-learn/blob/7a26816ff47f07337e544118fcee12be09ce660e/sklearn/utils/__init__.py#L1040-L1042

In summary, I am okay with keeping the current behavior and update our docs to communicate the change.

Labels: Bug, module:model_selection