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

Code before:

    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;
        }
    }

    /*
    At this point arr_len == 1, so the candidates are in [base, base + 1].

    We have two options:

Review comment (suggestion):
```suggestion
    npy_intp interval_length = arr_len;
    npy_intp half = interval_length >> 1;
    interval_length -= half; // length -> ceil(length / 2)
    ...
```?

The compiler is smart enough to optmize the `interval_length ` away. But the name is more readable.

Code after:
    - 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;
        }
    }

    /*
    At this point interval_length == 1, so the candidates are in the 
    interval [base, base + 1].

    We have two options:
        If cmp(arr[base], key_val) == true, insertion index is base + 1