/*
 * Binary search routines. 
 * Most of the ideas are taken from Bentleys book 
 * "Programming Pearls", Addison-Wesley, 1986.
 *
 * Copyright by Jyrki Katajainen and Tomi Pasanen
 * January - March 1998
 *
 * A driver program as well as the files included by it
 * are available from
 *
 * http://www.diku.dk/hjemmesider/ansatte/jyrki//Performance_Engineering/binary_driver.cc
 * http://www.diku.dk/hjemmesider/ansatte/jyrki//Performance_Engineering/heapsort.c
 * http://www.diku.dk/hjemmesider/ansatte/jyrki//Performance_Engineering/random_permutation.c
 */

const int null = 0;

/*
 * Source: Jon Bentley, Programming Pearls,
 * Addison-Wesley Publishing Company, 1986,
 * p. 38
 */

#define loop for( ; ; )

template<class T>
T* correct_search(T* first, T* last, const T& x) {

  T* l; T* m; T* p; T* u;

  /*
   * Preconditions: elements in [first,last) are 
   * in increasing order
   */
 
  l = first; u = last - 1;
  loop {

    /*
     * Invariant: if x is in [first,last), it is
     *            in [l,u]
     */

    if (l > u) {p = null; break;}
    m =  l + (u - l) / 2;
    if (*m < x) l = m + 1;
    else if (*m == x) {p = m; break;}
    else u = m - 1;
  }
  return p;
}

#undef loop

/*
 * Source: Gregory J. E. Rawlins, Compared to What?
 * An Introduction to the Analysis of Algorithms,
 * Computer Sceince Press, 1992, Algorithm 2.4
 */

template<class T>
T* recursive_search(T* first, T* last, const T& x) {

  T* mid;

  /*
   * Preconditions: elements in [first,last) are 
   * in increasing order
   */
 
  if (first >= last) return null;
  if (first + 1 == last) {
    if (*first == x) return first;
    else return null;
  }
  else {
    mid = first + (last - first) / 2;
    if (x < *mid) return recursive_search(first, mid, x);
    else return recursive_search(mid, last, x);
  }
}

/*
 * Source: Gregory J. E. Rawlins, Compared to What?
 * An Introduction to the Analysis of Algorithms,
 * Computer Sceince Press, 1992, Algorithm 2.5
 */

template<class T>
T* iterative_search(T* first, T* last, const T& x) {

  T* low; T* mid; T* high;
 
  /*
   * Preconditions: elements in [first,last) are 
   * in increasing order
   */

  low = first; high = last - 1;
  while (high >= low) {
    mid =  low + (high - low) / 2;
    if (*mid < x) low = mid + 1;
    else if (*mid == x) return mid;
    else  high = mid - 1;
  }
  return null;
}

/*
 * Source: Jon Bentley, Programming Pearls,
 * Addison-Wesley Publishing Company, 1986,
 * p. 85
 */

template<class T>
T* twoway_search(T* first, T* last, const T& x) {

  T* p; T* l; T* m; T* u;

  /*
   * Preconditions: elements in [first,last) are 
   *                in increasing order and
   *                *(first - 1) < x <= *last
   */

  l = first - 1; u = last;
  while (l + 1 != u) {

    /* assert(*l < x && *u >= x && l < u); */

    m = l + (u - l) / 2;
    if (*m < x) l = m;
    else u = m;
  }

  /* assert(l + 1 == u && *l < x && *u >= x); */

  p = u;
  if (p >= last || *p != x) p = null;
  return p;
}

/*
 * Source: Arne Andersson and Mikkel Thorup,
 * Implementing monotone priority queues,
 * Preliminary draft, 1997
 */

typedef struct {
  unsigned int signbit : 1;
  unsigned int exponentbias1023 : 11;
  /* unsigned int mantissa : 52; */
} doublemask;

typedef union {
  double as_double;
  doublemask as_mask;
} doubleandmask;

/* Takes a 32 bit word and finds the leftmost one bit */

inline long unsigned int leftmost_one(long unsigned int x)
{
  doubleandmask fm;
  fm.as_double = (double) x;
  return (fm.as_mask.exponentbias1023 - 1023);
}

/*
 * Source: Jon Bentley, Programming Pearls,
 * Addison-Wesley Publishing Company, 1986,
 * p. 87
 */


template<class T>
T* alternative_search(T* first, T* last, const T& x) {

   long unsigned int i, j, n;
   T* l; T* m; T* p;

  /*
   * Preconditions: elements in [first,last) are 
   *                in increasing order and
   *                *(first - 1) < x <= *last
   */

   if (first >= last) return null;

   n = last - first;
   j = leftmost_one(n);
   i = 1 << j; /* i is the smallest power of 2 <= n */
   l = first - 1;
   if (*(last - i) < x) l = last - i;

   while (i != 1) {

     /* Invariant: *l < x && *(l + i) >= x && i = 2^j */

     i = i / 2;
     m = l + i;
     if (*m < x) l = m;
   }

   /* assert(i == 1 && *l < x && *(l + i) >= x) */
     
   p = l + 1;
   if (p >= last || *p != x) p = null;
   return p;
}

/*
 * Source: Jon Bentley, Programming Pearls,
 * Addison-Wesley Publishing Company, 1986,
 * p. 87
 */

template<class T>
T* unrolled_search(T* first, T* last, const T& x) {

  /*
   * Preconditions: elements in [first,last) are 
   *                in increasing order and
   *                last - first < 2^32
   */

  long unsigned int i, j, n;
  T* l; T* m; T* p;

  if (first >= last) return null;

  n = last - first;
  j = leftmost_one(n);
  i = 1 << j; /* i is the largest power of 2 <= n */
  l = first - 1;
  if (*(last - i) < x) l = last - i;

  switch(j) {
  case 31:
    m = l + 1073741824; if (*m < x) l = m;
  case 30:
    m = l + 536870912;  if (*m < x) l = m;
  case 29:
    m = l + 268435456;  if (*m < x) l = m;
  case 28:
    m = l + 134217728;  if (*m < x) l = m;
  case 27:
    m = l + 67108864;   if (*m < x) l = m;
  case 26:
    m = l + 33554432;   if (*m < x) l = m;
  case 25:
    m = l + 16777216;   if (*m < x) l = m;
  case 24:
    m = l + 8388608;    if (*m < x) l = m;
  case 23:
    m = l + 4194304;    if (*m < x) l = m;
  case 22:
    m = l + 2097152;    if (*m < x) l = m;
  case 21:
    m = l + 1048576;    if (*m < x) l = m;
  case 20:
    m = l + 524288;     if (*m < x) l = m;
  case 19:
    m = l + 262144;     if (*m < x) l = m;
  case 18:
    m = l + 131072;     if (*m < x) l = m;
  case 17:
    m = l + 65536;      if (*m < x) l = m;
  case 16:
    m = l + 32768;      if (*m < x) l = m;
  case 15:
    m = l + 16384;      if (*m < x) l = m;
  case 14:
    m = l + 8192;       if (*m < x) l = m;
  case 13:
    m = l + 4096;       if (*m < x) l = m;
  case 12:
    m = l + 2048;       if (*m < x) l = m;
  case 11:
    m = l + 1024;       if (*m < x) l = m;
  case 10:
    m = l + 512;        if (*m < x) l = m;
  case 9:
    m = l + 256;        if (*m < x) l = m;
  case 8:
    m = l + 128;        if (*m < x) l = m;
  case 7:
    m = l + 64;         if (*m < x) l = m;
  case 6:
    m = l + 32;         if (*m < x) l = m;
  case 5:
    m = l + 16;         if (*m < x) l = m;
  case 4:
    m = l + 8;          if (*m < x) l = m;
  case 3:
    m = l + 4;          if (*m < x) l = m;
  case 2:
    m = l + 2;          if (*m < x) l = m;
  case 1:
    m = l + 1;          if (*m < x) l = m;
  case 0:
    p = l + 1;
  }
  if ((p >= last) || *p != x) return null;
  return p;
}