# How to speed up interval lookup for a piecewise defined function?

At the innermost loop of my software is a lookup of a function value from a piecewise defined function with interpolation between the sample points.

Illustration:

*                                                *
*                               *
*      *
|--|------|------|--------X---------|-------------|


There are several different interpolation methods in between the sample points.

My current implementation consists of two arrays, one for the interval endpoints and one for the function values. Searching is done via std::upper_bound if the array is large and find_if if the arrays are small. Based on the access patterns the speed improved tremendously by adding a cache value for the previous positions interval and checking if it is still valid for the new position

[1%] if (x >= *cache && x < *(cache + 1)) {
outside = false;
position = cache;
} else if (x < *begin || x > *end) {
outside = true;
} else {
outside = false;

if (n_ < 64) {
[3%]    position = std::find_if(begin, end, std::bind2nd(std::greater_equal<argument_type>(), x)) - 1;
} else {
position = std::upper_bound(begin, end, x) - 1;
}
cache = position ;
}


The lookup currently still takes about 5% of total execution time, which seems too much. I have indicated costs of the most expensive lines which is the comparison with the cached position [1%] and the linear search [3%].

Is there a faster way to lookup the interval which contains X?

Edit

Edit It seems the std::find_if implemenation is actually faster than binary search after all. I have made a test here: http://ideone.com/Fo5ZeW

#include <iostream>
#include <ctime>
#include <algorithm>
using namespace std;

double* binary_search8(double* arr, const double& x)
{
if (x<arr[4]) {
if(x<arr[2]) {
if(x<arr[1]) {
return arr;
} else {
return arr + 1;
}
}
else { //2..3
if(x<arr[3]){
return arr + 2;
} else {
return arr + 3;
}
}
} else { // 4..7
if(x<arr[6]) {
if(x<arr[5]) {
return arr + 4;
} else {
return arr + 5;
}
} else { //6..7
if(x<arr[7]) {
return arr + 6;
} else {
return arr + 7;
}
}
}
}

double* binary_search16(double* arr, const double& x)
{
if (x<arr[8]) return binary_search8(arr, x);
else return binary_search8(arr+8, x);
}

int main() {
double test[16] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150};

vector<double> test_vals;
size_t num_elem = 10000000;
for (size_t i=0; i<num_elem; i++)
test_vals.push_back(rand()%150);

clock_t begin, end;
double elapsed_secs;

// Linear search
double linavg = 0;
begin = clock();
for (size_t i=0; i<num_elem; i++) {
linavg += *(find_if(test, test+16, bind2nd(std::greater_equal<double>(), test_vals[i]))-1);
}
end = clock();
elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
cout << "linear secs: " <<elapsed_secs << ", testval = " << linavg << std::endl;

// Binary search
double binavg = 0;
begin = clock();
for (size_t i=0; i<num_elem; i++) {
binavg += *binary_search16(test, test_vals[i]);
}
end = clock();
elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
cout << "binary secs: " <<elapsed_secs << ", testval = " << binavg << std::endl;

// Binary search 2
double binavg2 = 0;
begin = clock();
for (size_t i=0; i<num_elem; i++) {
binavg2 += *(std::upper_bound(test, test+16, test_vals[i])-1);
}
end = clock();
elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
cout << "binary2 secs: " <<elapsed_secs << ", testval = " << binavg2 << std::endl;

return 0;
}

• Would be nice, if you post, if and how much you have improved the speed with witch attempt. maybe there are also some possible improvements in your interpolation between the sample points. Maybe make it part of another review question. – MrSmith42 Mar 10 '14 at 16:25
• Have you tested if linear search is still faster for your 64 elements? Binary search would for 64 cases reduce the average number of comparisons from 32 to 6. – MrSmith42 Mar 11 '14 at 20:13
• @MrSmith42 yes, it is. Branch prediction is phenomenally effective in this case it seems. – user35796 Mar 12 '14 at 5:45

It'll take some testing to be entirely certainly, but a few possibilities to consider:

if (x >= *cache && x < *(cache + 1)) {


I'd try replacing this with:

if (static_cast<unsigned>(x-*cache) < static_cast<unsigned>(*(cache+1) - *cache)) {


[Technically, unsigned is a place-holder here. It should really be something like std::make_unsigned<input_type>::type. That is, if the input is (say) long long, you want unsigned long long, not just unsigned int.]

This basically replaces two comparisons (that must be executed serially, thanks to && imposing a sequence point) with one plus a little extra math. The math can typically be done in parallel, and the branch becomes more predictable as well.

if (x < *begin || x > *end) {
outside = true;
} else {
outside = false;


I'd at least try just assigning the result of the comparison:

outside = x < *begin || x > *end;


Again, you can try the same "trick" as above, and get something like:

outside = static_cast<unsigned>(x-*begin) <= static_cast<unsigned>(*end-*begin);


As before, this replaces two relatively unpredictable comparisons with one relatively predictable one.

if (n_ < 64) {
[3%]    position = std::find_if(begin, end, std::bind2nd(std::greater_equal<argument_type>(), x)) - 1;
} else {
position = std::upper_bound(begin, end, x) - 1;
}


If you're spending significantly more time in the linear search than the binary search, I'd try reducing the cut-off between the two, perhaps to 32 instead of 64.

As with any micro-optimization, none of these is really guaranteed to be effective. It's even possible your compiler may be incorporating one or more of them into the code it generates, so they'll do no good at all, or even that its optimizer won't recognize the resulting patterns in the code as well, so they could actually hurt performance. The only real way to find out is some testing and profiling.

Why do you need a linear search [O(n)]?

You can determine in which interval the value belongs by a binary search [O(log n)] (do not implement is with a recursion! use a loop or hard-code it))).

Average number of comparisons for linear search is n/2. Binary search needs a constant number of comparisons of ceiling( ln(n) / ln(2) ).

           \ n|  1 |  2 |  4 |  8 | 16 | 32 | 64
--------------+----+----+----+----+----+----+---- ...
linear-search |  0 |  1 |  2 |  4 |  8 | 16 | 32
binary-search |  0 |  1 |  2 |  3 |  4 |  5 |  6


So binary search should be faster, if comparison is the main factor, for n > 4.

If n is constant (its seams 64 in your case) you can even hard-code the binary search to avoid the overhead of a loop. (that is also possible for linear search.)

As an example how to hard-code a binary search: [for a array with 8 entries] Code will get quite long for 64 entries, but it will be fast.

If your intervallBoundary array is intervals[0..7]:

  if(x<intervalls[4]){
if(x<intervalls[2]){
if(x<intervalls[1])
return intervall0;
else
return intervall1
}
else{ //2..3
if(x<intervalls[3])
return intervall2;
else
return intervall3
}
else{ // 4..7
if(x<intervalls[6]){
if(x<intervalls[4])
return intervall4;
else
return intervall5
}
else{ //6..7
if(x<intervalls[7])
return intervall6;
else
return intervall7
}
}


REMARK

If there is any rule for the distribution of the intervals, you might be able to calculate the correct interval (or at least calculate a good guess and search from there).

• The constant that gets dropped in Landau notation is the reason $c$ in $O(c*log(n))$ for all arrays smaller than 64 linear search is faster. – user35796 Mar 10 '14 at 11:58
• Added a hard-coded binary search example for 8 cases. – MrSmith42 Mar 10 '14 at 13:59