#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
 // looking for index i such that om[i]<x<om[i+1]
 // We know that the search is necessary only in the region om[ia]<x<om[N-1]!
 // Be carefull that 0<i<N! Otherwise coredump might occur.
  // ..........
  // return ...;
}

const int dN=10;

inline int find_ordered(double x, const double om[], int ia, int N)
{
  // The optimized search for correlated lookups
  // Index ia is the index from the previous lookup. Most probably,
  // the new point x is close to om[ia].
  //  ........................
}

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);
    if (ia<0 || ia>=N) {cerr<<"This code does not yet work!"<<endl; return 1;}
    else 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 (ia<0 || ia>=N) {cerr<<"This code does not yet work!"<<endl; return 1;}
      
      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;
}