#include <iostream>  // for printing
#include <cstdlib>   // for drand48
#include <algorithm> // for sort routine

inline int bisection(double x, const double om[], int ia, int N)
{// Bisection in the interval between ia and N
  int ib = N-1;
  while (ib-ia > 1){
    int k = (ib+ia) >> 1; // fast way of dividing by 2 is right shift
    if (om[k] > x) ib=k;  // x is smaller -> upper bound is ib
    else ia=k;            // x is bigger  -> upper bound becomes ia 
  }
  return ia;
}

const int dN=10;

inline int find_ordered(double x, const double om[], int ia, int N)
{
  if (x>om[ia]){
    if (ia>=N-1) return N-1;      // make sure we are not out of limits (avoid coredump)
    if (x<om[ia+1]) return ia;    // om[ia]<x<om[ia+1]! We are done.
    ia++;          	          // searching after om[ia]
    if (ia+dN>=N-1) return bisection(x, om, ia, N); // small interval is left, search in the remaining interval
    if (x<om[ia+dN]) return bisection(x, om, ia, ia+dN+1); // Need to search the first dN interval
    else return bisection(x, om, ia+dN, N);   // Need to search the rest of the table
  }else{
    if (ia<=0) return 0;          // check out of bounds
    if (x>om[ia-1]) return ia-1;  // point bracketed successfully
    ia--;
    if (ia-dN<=0) return bisection(x, om, 0, ia+1); // only small interval left, search it with bisection
    if (x>om[ia-dN]) return bisection(x, om, ia-dN, ia+1); //search small interval of dN points first
    else return bisection(x, om, 0, ia-dN+1);       // search the rest cause point is far
  }
}

int main()
{
  using namespace std;
  
  const int N=10000; // size of array is fixed at this value
  double om[N];      // real array of numbers is initialized

  srand48(time(0));  // random number generator from cstdlib is initialized to unique value
  
  for (int i=0; i<N; i++) om[i] = drand48(); // array is filled with random numbers
  
  sort(om, om+N); // array is sorted using STL call

  // Testing the first part of the assignment
  {
    double x0=0.5;
    int ia = bisection(x0, om, 0, N);
    cout<<"The number preceding "<<x0<<" is "<<om[ia]<<" and following "<<om[ia+1]<<endl;
  }

  // Testing the second part of the agorithm
  {
    int M=100*N; // testing sufficiently large number of casses to find possible bugs
    int ia=0;    // start searching table at the beginning in the absence of knowledge
    for (int i=0; i<M; i++){ 
      double x = drand48();            // random number to be searched in the table
      ia = find_ordered(x, om, ia, N); // actual searching routine
    
      if (x>om[0] && x<om[N-1]){        // interpolation possible only if number is inside the interval
	                                // otherwise extrapolation needs to be used
	if (!(om[ia]<x && x<om[ia+1]))  // if you see the following message, the code does not yet work!
	  cerr<<"Failed at point "<<om[ia]<<" "<<x<<" "<<om[ia+1]<<endl;
	
      }
    }
  }
  return 0;
}