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
- I have included the switch after profiling since Linear search can faster for small arrays (and there are a lot of small arrays).
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;
}
