Repository: numpy/numpy
File: numpy/_core/src/npysort/binsearch.cpp (C++)

Code before:
    All the computations at a given step are independent across different keys.
    Meaning that the CPU can execute multiple keys out-of-order in parallel.

    Invariant:
    - cmp(arr[i], key_val) == true  for all i < base
    - cmp(arr[i], key_val) == false for all i >= base + length

    The insertion index candidates are i in range [base, base+length] and
    on each iteration we shrink the range into either
        [base, ceil(length / 2)]
    or
        [base + floor(length / 2), ceil(length / 2)]

    Optimization: we unroll the first iteration for the following reasons:
        1. ret is not initialized with the bases, so we save |keys| writes
        by not having to intialize it with 0s.
        2. By assuming the initial base for every key is 0, we also save
        |keys| reads.
        3. In the first iteration, all elements are compared against the
        median. So we can store it in a variable and use it for all keys.

    This initial block replaces the initialization loop:

    for (npy_intp i = 0; i < key_len; ++i) {
        *(npy_intp *)(ret + i * ret_str) = 0;
    }

    Note that when arr_len = 1, then half is 0 so the following block
    initializes the array as with 0s.
    */
    npy_intp half = arr_len >> 1;
    arr_len -= half; // length -> ceil(length / 2)

    // We write 0 to explicitly refer to base
    const T mid_val = *(const T *)(arr + (0 + half) * arr_str);

    for (npy_intp i = 0; i < key_len; ++i) {
        const T key_val = *(const T *)(key + i * key_str);
        *(npy_intp *)(ret + i * ret_str) = cmp(mid_val, key_val) * half;
    }

    while (arr_len > 1) {
        npy_intp half = arr_len >> 1;
        arr_len -= half; // length -> ceil(length / 2)

        for (npy_intp i = 0; i < key_len; ++i) {
            npy_intp &base = *(npy_intp *)(ret + i * ret_str);
            const T mid_val = *(const T *)(arr + (base + half) * arr_str);
            const T key_val = *(const T *)(key + i * key_str);
            base += cmp(mid_val, key_val) * half;
        }

Review comment (suggestion):
```suggestion
    This initial block replaces the initialization loop that is used for the arr_len==0 case.
```

Code after:
    keys, we compute each step for all keys before continuing on the next
    step. All the computations at a given step are independent across
    different keys. Meaning that the CPU can execute multiple keys
    out-of-order in parallel.

    Invariant:
    - cmp(arr[i], key_val) == true  for all i < base
    - cmp(arr[i], key_val) == false for all i >= base + length

    The insertion index candidates are in range [base, base+length] and
    on each iteration we shrink the range into either
        [base, ceil(length / 2)]
    or
        [base + floor(length / 2), ceil(length / 2)]

    Optimization: we unroll the first iteration for the following reasons:
        1. ret is not initialized with the bases, so we save |keys| writes
        by not having to intialize it with 0s.
        2. By assuming the initial base for every key is 0, we also save
        |keys| reads.
        3. In the first iteration, all elements are compared against the
        median. So we can store it in a variable and use it for all keys.

    This initial block replaces the initialization loop that is used for the
    arr_len==0 case. Note that when arr_len = 1, then half is 0 so the
    following block initializes the array as with 0s.
    */
    npy_intp interval_length = arr_len;
    npy_intp half = interval_length >> 1;
    interval_length -= half; // length -> ceil(length / 2)

    npy_intp base = 0;
    const T mid_val = *(const T *)(arr + (base + half) * arr_str);

    for (npy_intp i = 0; i < key_len; ++i) {
        const T key_val = *(const T *)(key + i * key_str);
        *(npy_intp *)(ret + i * ret_str) = cmp(mid_val, key_val) * half;
    }

    while (interval_length > 1) {
        npy_intp half = interval_length >> 1;
        interval_length -= half; // length -> ceil(length / 2)

        for (npy_intp i = 0; i < key_len; ++i) {
            npy_intp &base = *(npy_intp *)(ret + i * ret_str);
            const T mid_val = *(const T *)(arr + (base + half) * arr_str);
            const T key_val = *(const T *)(key + i * key_str);
            base += cmp(mid_val, key_val) * half;
        }
    }