4 ms·

Beating Bisect with Branchless Binary Search

62 points by juliusgeo 3y ago

juliusgeo 3y ago
I was inspired by this article[0] to see if I could beat the Python standard library implementation of "bisect_left" using C-extensions. [0] https://probablydance.com/2023/04/27/beautiful-branchless-binary-search/ https://probablydance.com/2023/04/27/beautiful-branchless-bi...
- cycomanic 3y ago
  > This is the second time that I came across an algorithm in Knuth’s books that is brilliant and should be used more widely but somehow was forgotten. Maybe I should actually read the book… It’s just really hard to see which ideas are good and which ones aren’t. For example immediately after sketching out Shar’s algorithm, Knuth spends far more time going over a binary search based on the Fibonacci sequence. It’s faster if you can’t quickly divide integers by 2, and instead only have addition and subtraction. So it’s probably useless, but who knows? When reading Knuth’s book, you have to assume that most algorithms are useless, and that the good things have been highlighted by someone already. Luckily for people like me, there seem to still be a few hidden gems. Oh the irony, considering that fibonacci hashing is another of those forgotten algorithms, see yesterday's discussion: https://news.ycombinator.com/item?id=35739629 https://news.ycombinator.com/item?id=35739629
- dang 3y ago
  Discussed here a couple days ago: Beautiful branchless binary search - https://news.ycombinator.com/item?id=35737862 https://news.ycombinator.com/item?id=35737862 - April 2023 (143 comments)
xeromal 3y ago
I first read this as "Eating Bistec with..." and got hungry.
- juliusgeo 3y ago
  Now I’m hungry too! Darn.
  - xeromal 3y ago
    Hope you find something good to eat! I made some taters and chicken.
alexmolas 3y ago
I don't know if it's a problem on my side, but images are not rendering. Besides from that the implementation is very cool!
- juliusgeo 3y ago
  I just checked on my phone and images are rendering fine for me—it might be because they’re high res or otherwise delayed in loading.
  - cycomanic 3y ago
    I also can't load the images. I get a "image can't be loaded because it contains errors" error on Firefox mobile nightly and and just a failure on brave on mobile.
    - juliusgeo 3y ago
      Weirdly enough, when I go to this url on Safari it loads fine: https://raw.githubusercontent.com/juliusgeo/branchless_bisect/master/images/1.png https://raw.githubusercontent.com/juliusgeo/branchless_bisec... When I go to it on Chrome, I get an error.
    - raimue 3y ago
      This is most probably because they are actually TIFF, but with the .png file extension they are served as image/png.
- boywitharupee 3y ago
  The images aren't loading for me either
  - juliusgeo 3y ago
    That is concerning. I just tried loading it on an incognito tab and the images loaded after a short delay.
exebook 3y ago
I guess while Python interprets this branchless code it still can do some branching? Or am I missing something here?
- juliusgeo 3y ago
  There are some provisos to the “Branchless” tag that are covered in the article I drew inspiration from. One of which is that a “CMOVE” is not technically a branch :)
  - yarg 3y ago
    You have to do some really weird shit if you want generalised branchless ternary statements. begin += (arr[step+begin] < value)?step:0; Something like: int mask = ((arr[step + begin] - value) >> 31); //Depends on signed shift begin += (step & mask) | (0 & ~mask); (Obviously in this case, it's simplifiable.)
    - juliusgeo 3y ago
      I was considering doing it all in a ternary statement, but I feel that the current form is also branchless because it is simply a multiply and add. The extra bounds-checking condition can probably be omitted, but I haven't tested that. for (step >>= 1; step != 0; step >>=1) { if ((next = begin + step) < size) { begin += PyObject_RichCompareBool(PyList_GetItem(list_obj, next), value, Py_LT) * step; } }
      - yarg 3y ago
        My point was, anywhere there's a hidden 'if' can be branching. If there's no calculation being done, it'll simplify. value = (test) ? const0 : const1; But if calculations are being done, it won't. value = (test) ? calc0() : calc1(); If you want non-branching where the ternary options are calculated, you need to calculate both. This matters most with SIMD operations. Look at section 2.5.1 (Branch-Equivalent SIMD Processing) http://ftp.cvut.cz/kernel/people/geoff/cell/ps3-linux-docs/CellProgrammingTutorial/BasicsOfSIMDProgramming.html http://ftp.cvut.cz/kernel/people/geoff/cell/ps3-linux-docs/C...
tomsmeding 3y ago
I would love some logarithmic y axes on these plots. As is, the faster versions are basically horizontal lines, making one think whether perhaps something might be wrong with the benchmark and the compiler just optimised everything away to an empty function. Sinilarly, what's going on with the performance of bisect_left_c? (Second graph.) Why is the graph completely flat at first, only to then ramp up with a perfectly straight line? If you add some measurement points in between, will it turn out the high measurement on the right of the graph is actually a measurement error? If this works it would make a nice python lib.
- juliusgeo 3y ago
  Thanks for the feedback! I experimented with using the matplotlib `plt.xscale("log")` function, but because the x axis already distributed by powers of 2, it to some extent is already a log axis (unless I'm making a very wrong assumption). I was curious about the behavior of the lines, and how they seem to ramp up in weird ways. I tried to mitigate by increasing the number of samples in the benchmark, but I haven't been able to get graphs that look as nice as the original article. I think to some extent though, that is because the `sortedcontainers` implementation is so bad. If the `y` axis was scaled better, you would be able to see a more linear increase in execution time. Also, if you want to, the repo provides a build script for the C-extension so you can mess around with the benchmarks. I would love for you to find some bugs in my implementation :).
  - tomsmeding 3y ago
    > > I would love some logarithmic y axes on these plots. > I experimented with using the matplotlib `plt.xscale("log")` function, but because the x axis already distributed by powers of 2, [...] I said y, not x ;) Also, the x axis is not distributed by powers of 2, right? I see 0 through roughly 2.7*10^8 — either that or the axis is confusingly labeled. :) I would want two graphs here: 1. High-res linear-linear to show the expected logarithmic behaviour (surely binary search is logarithmic, right?), as well as possibly some performance cliff as you run out of cache and start hitting RAM for a decent part of the search. 2. Linear x but logarithmic y, on which you compare the various implementations. The point is that because all implementations presumably have the same complexity, and furthermore probably have caching issues starting from similar (though not necessarily equal) points, the most interesting point of info is their _ratio_. And since log(2*f(x)) = log(2) + log(f(x)), the ratio between two functions is simply their _distance_ on an yscale("log") plot. Much easier to see, especially if the ratio is like 100x. Perhaps you'll even see that for sufficiently large inputs, the ratio shrinks because the search becomes memory-bound. :)
kevingadd 3y ago
If "execution time (s)" reads 0 for one of your measured cases, I become extremely skeptical of your measurements. I get that the measurement probably just rounded down but that alongside a perfectly straight horizontal line just makes me go "what's going on here? What conclusion can I draw from this other than that the algorithm is somehow O(1)?" Is the value being plotted the average time per execution? You ran each test scenario for a while, not once, I assume. It's also worth considering whether you should measure against the same dataset and value every time, or have a bunch of different ones. If it's the same dataset and value it's possible you're priming the branch predictor tables, or worse, favoring one algorithm by accident due to the layout of the data.
- juliusgeo 3y ago
  Yeah it’s difficult because all the times are by default in seconds from ‘timeit’. I think the sortedcontainers impl is just really, really bad, and that’s why that graph looks so flat in comparison. The datasets and search values are generated on each loop in the benchmark, so hopefully no priming issues.
shoo 3y ago
If anyone has the time & energy to put into figuring out a production quality open source version of this, one place that could benefit is numpy's searchsorted function https://github.com/numpy/numpy/blame/main/numpy/core/src/npysort/binsearch.cpp https://github.com/numpy/numpy/blame/main/numpy/core/src/npy...

62 points by juliusgeo 3y ago

I was inspired by this article[0] to see if I could beat the Python standard library implementation of "bisect_left" using C-extensions. [0] https://probablydance.com/2023/04/27/beautiful-branchless-binary-search/ https://probablydance.com/2023/04/27/beautiful-branchless-bi...

cycomanic 3y ago
> This is the second time that I came across an algorithm in Knuth’s books that is brilliant and should be used more widely but somehow was forgotten. Maybe I should actually read the book… It’s just really hard to see which ideas are good and which ones aren’t. For example immediately after sketching out Shar’s algorithm, Knuth spends far more time going over a binary search based on the Fibonacci sequence. It’s faster if you can’t quickly divide integers by 2, and instead only have addition and subtraction. So it’s probably useless, but who knows? When reading Knuth’s book, you have to assume that most algorithms are useless, and that the good things have been highlighted by someone already. Luckily for people like me, there seem to still be a few hidden gems. Oh the irony, considering that fibonacci hashing is another of those forgotten algorithms, see yesterday's discussion: https://news.ycombinator.com/item?id=35739629 https://news.ycombinator.com/item?id=35739629
dang 3y ago
Discussed here a couple days ago: Beautiful branchless binary search - https://news.ycombinator.com/item?id=35737862 https://news.ycombinator.com/item?id=35737862 - April 2023 (143 comments)

xeromal 3y ago

I first read this as "Eating Bistec with..." and got hungry.

juliusgeo 3y ago
Now I’m hungry too! Darn.
- xeromal 3y ago
  Hope you find something good to eat! I made some taters and chicken.

alexmolas 3y ago

I don't know if it's a problem on my side, but images are not rendering. Besides from that the implementation is very cool!

juliusgeo 3y ago
I just checked on my phone and images are rendering fine for me—it might be because they’re high res or otherwise delayed in loading.
- cycomanic 3y ago
  I also can't load the images. I get a "image can't be loaded because it contains errors" error on Firefox mobile nightly and and just a failure on brave on mobile.
  - juliusgeo 3y ago
    Weirdly enough, when I go to this url on Safari it loads fine: https://raw.githubusercontent.com/juliusgeo/branchless_bisect/master/images/1.png https://raw.githubusercontent.com/juliusgeo/branchless_bisec... When I go to it on Chrome, I get an error.
  - raimue 3y ago
    This is most probably because they are actually TIFF, but with the .png file extension they are served as image/png.
boywitharupee 3y ago
The images aren't loading for me either
- juliusgeo 3y ago
  That is concerning. I just tried loading it on an incognito tab and the images loaded after a short delay.

exebook 3y ago

I guess while Python interprets this branchless code it still can do some branching? Or am I missing something here?

juliusgeo 3y ago
There are some provisos to the “Branchless” tag that are covered in the article I drew inspiration from. One of which is that a “CMOVE” is not technically a branch :)
- yarg 3y ago
  You have to do some really weird shit if you want generalised branchless ternary statements. begin += (arr[step+begin] < value)?step:0; Something like: int mask = ((arr[step + begin] - value) >> 31); //Depends on signed shift begin += (step & mask) | (0 & ~mask); (Obviously in this case, it's simplifiable.)
  - juliusgeo 3y ago
    I was considering doing it all in a ternary statement, but I feel that the current form is also branchless because it is simply a multiply and add. The extra bounds-checking condition can probably be omitted, but I haven't tested that. for (step >>= 1; step != 0; step >>=1) { if ((next = begin + step) < size) { begin += PyObject_RichCompareBool(PyList_GetItem(list_obj, next), value, Py_LT) * step; } }
    - yarg 3y ago
      My point was, anywhere there's a hidden 'if' can be branching. If there's no calculation being done, it'll simplify. value = (test) ? const0 : const1; But if calculations are being done, it won't. value = (test) ? calc0() : calc1(); If you want non-branching where the ternary options are calculated, you need to calculate both. This matters most with SIMD operations. Look at section 2.5.1 (Branch-Equivalent SIMD Processing) http://ftp.cvut.cz/kernel/people/geoff/cell/ps3-linux-docs/CellProgrammingTutorial/BasicsOfSIMDProgramming.html http://ftp.cvut.cz/kernel/people/geoff/cell/ps3-linux-docs/C...

tomsmeding 3y ago

I would love some logarithmic y axes on these plots. As is, the faster versions are basically horizontal lines, making one think whether perhaps something might be wrong with the benchmark and the compiler just optimised everything away to an empty function. Sinilarly, what's going on with the performance of bisect_left_c? (Second graph.) Why is the graph completely flat at first, only to then ramp up with a perfectly straight line? If you add some measurement points in between, will it turn out the high measurement on the right of the graph is actually a measurement error? If this works it would make a nice python lib.

juliusgeo 3y ago
Thanks for the feedback! I experimented with using the matplotlib `plt.xscale("log")` function, but because the x axis already distributed by powers of 2, it to some extent is already a log axis (unless I'm making a very wrong assumption). I was curious about the behavior of the lines, and how they seem to ramp up in weird ways. I tried to mitigate by increasing the number of samples in the benchmark, but I haven't been able to get graphs that look as nice as the original article. I think to some extent though, that is because the `sortedcontainers` implementation is so bad. If the `y` axis was scaled better, you would be able to see a more linear increase in execution time. Also, if you want to, the repo provides a build script for the C-extension so you can mess around with the benchmarks. I would love for you to find some bugs in my implementation :).
- tomsmeding 3y ago
  > > I would love some logarithmic y axes on these plots. > I experimented with using the matplotlib `plt.xscale("log")` function, but because the x axis already distributed by powers of 2, [...] I said y, not x ;) Also, the x axis is not distributed by powers of 2, right? I see 0 through roughly 2.7*10^8 — either that or the axis is confusingly labeled. :) I would want two graphs here: 1. High-res linear-linear to show the expected logarithmic behaviour (surely binary search is logarithmic, right?), as well as possibly some performance cliff as you run out of cache and start hitting RAM for a decent part of the search. 2. Linear x but logarithmic y, on which you compare the various implementations. The point is that because all implementations presumably have the same complexity, and furthermore probably have caching issues starting from similar (though not necessarily equal) points, the most interesting point of info is their _ratio_. And since log(2*f(x)) = log(2) + log(f(x)), the ratio between two functions is simply their _distance_ on an yscale("log") plot. Much easier to see, especially if the ratio is like 100x. Perhaps you'll even see that for sufficiently large inputs, the ratio shrinks because the search becomes memory-bound. :)

kevingadd 3y ago

If "execution time (s)" reads 0 for one of your measured cases, I become extremely skeptical of your measurements. I get that the measurement probably just rounded down but that alongside a perfectly straight horizontal line just makes me go "what's going on here? What conclusion can I draw from this other than that the algorithm is somehow O(1)?" Is the value being plotted the average time per execution? You ran each test scenario for a while, not once, I assume. It's also worth considering whether you should measure against the same dataset and value every time, or have a bunch of different ones. If it's the same dataset and value it's possible you're priming the branch predictor tables, or worse, favoring one algorithm by accident due to the layout of the data.

juliusgeo 3y ago
Yeah it’s difficult because all the times are by default in seconds from ‘timeit’. I think the sortedcontainers impl is just really, really bad, and that’s why that graph looks so flat in comparison. The datasets and search values are generated on each loop in the benchmark, so hopefully no priming issues.

shoo 3y ago

If anyone has the time & energy to put into figuring out a production quality open source version of this, one place that could benefit is numpy's searchsorted function https://github.com/numpy/numpy/blame/main/numpy/core/src/npysort/binsearch.cpp https://github.com/numpy/numpy/blame/main/numpy/core/src/npy...