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 &gt;= *cache &amp;&amp; x &lt; *(cache + 1)) {
    outside = false;
    position = cache;
} else if (x &lt; *begin || x &gt; *end) {
    outside = true;
} else {
    outside = false;

    if (n_ &lt; 64) {
[3%]    position = std::find_if(begin, end, std::bind2nd(std::greater_equal&lt;argument_type&gt;(), 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 &lt;iostream&gt;
#include &lt;ctime&gt;
#include &lt;algorithm&gt;
using namespace std;

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


double* binary_search16(double* arr, const double&amp; x)
{
    if (x&lt;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&lt;double&gt; test_vals;
  size_t num_elem = 10000000;
  for (size_t i=0; i&lt;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&lt;num_elem; i++) {   
    linavg += *(find_if(test, test+16, bind2nd(std::greater_equal&lt;double&gt;(), test_vals[i]))-1);
  }
  end = clock();
  elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
  cout &lt;&lt; "linear secs: " &lt;&lt;elapsed_secs &lt;&lt; ", testval = " &lt;&lt; linavg &lt;&lt; std::endl;


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

  // Binary search 2
  double binavg2 = 0;
  begin = clock();
  for (size_t i=0; i&lt;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 &lt;&lt; "binary2 secs: " &lt;&lt;elapsed_secs &lt;&lt; ", testval = " &lt;&lt; binavg2 &lt;&lt; std::endl;

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

if (x &gt;= *cache &amp;&amp; x &lt; *(cache + 1)) {


I'd try replacing this with:

if (static_cast&lt;unsigned&gt;(x-*cache) &lt; static_cast&lt;unsigned&gt;(*(cache+1) - *cache)) {


[Technically, unsigned is a place-holder here. It should really be something like std::make_unsigned&lt;input_type&gt;::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 &amp;&amp; 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 &lt; *begin || x &gt; *end) {
    outside = true;
} else {
    outside = false;


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

outside = x &lt; *begin || x &gt; *end;


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

outside = static_cast&lt;unsigned&gt;(x-*begin) &lt;= static_cast&lt;unsigned&gt;(*end-*begin);


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

if (n_ &lt; 64) {
[3%]    position = std::find_if(begin, end, std::bind2nd(std::greater_equal&lt;argument_type&gt;(), 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.