There's a rare bug in binary search algorithms having to do with integer overflow in the calculation of the midpoint of the search range. See http://googleresearch.blogspot.com/2006/06/extra-extra-read-all-about-it-nearly.html
The offending code in PDL's vsearch is:
m = (jh+jl) >> 1;
where jh & jl are declared as having type PDL_Indx. At least on my 64bit linux box, PDL_Indx is declared as (signed) long. The fix (as per the above post) is to coerce them into unsigned longs, and then perform the shift, e.g.
m = ((unsigned PDL_Indx) jh + (unsigned PDL_Indx) jl) >> 1;
(which of course won't work).
It doesn't look like there's an unsigned version of PDL_Indx defined. To prevent future catastrophes, it'd be nice if that were done.
Log in to post a comment.