Median value selection (Fixed)

In 2009 I presented a heap based median selection algorithm. It was original, and was apparently very fast when compiled with the Intel compiler (icc). Since I don't have the Intel compiler anymore, I can't test it's performance. It's slower than the nthElement code given below when compiled with g++-7 -O3.

Here is the fixed code.

float fixedHeapMedian (float *a) {
  const unsigned char HEAP_LEN = 13;
  float left[HEAP_LEN], right[HEAP_LEN], *p, median;
  unsigned char nLeft, nRight;

  // pick first value as median candidate
  p = a;
  median = *p++;
  nLeft = nRight = 0;

  for (;;) {
    //dumpState(left, nLeft, median, right, nRight, p, 27 - (p-a));
    //assert(stateIsValid(left, nLeft, median, right, nRight));

    // get next value
    float val = *p++;

    // if value is smaller than median, append to left heap
    if (val <= median) {
      // move biggest value to the top of left heap
      unsigned char child = nLeft++, parent = (child - 1) / 2;
      while (child && val > left[parent]) {
        left[child] = left[parent];
        child = parent;
        parent = (parent - 1) / 2;
      }
      left[child] = val;

      // if left heap is full
      if (nLeft == HEAP_LEN) {
        //cout << "---" << endl;
        // for each remaining value
        for (unsigned char nVal = 27-(p - a); nVal; --nVal) {
          //dumpState(left, nLeft, median, right, nRight, p, nVal);
          //assert(stateIsValid(left, nLeft, median, right, nRight));
          // get next value
          val = *p++;
          // discard values falling in other heap
          if (val >= median) {
            continue;
          }
          // if val is bigger than biggest in heap, val is new median
          if (val >= left[0]) {
            median = val;
            continue;
          }
          // biggest heap value becomes new median
          median = left[0];
          // insert val in heap
          parent = 0;
          child = 2;
          while (child < HEAP_LEN) {
            if (left[child-1] > left[child]) {
              child = child-1;
            }
            if (val >= left[child]) {
               break;
            }
            left[parent] = left[child];
            parent = child;
            child = (parent + 1) * 2;
          }
          left[parent] = val;
        }
        return median;
      }
    } else {
      // move smallest value to the top of right heap
      unsigned char child = nRight++, parent = (child - 1) / 2;
      while (child && val < right[parent]) {
        right[child] = right[parent];
        child = parent;
        parent = (parent - 1) / 2;
      }
      right[child] = val;

      // if right heap is full
      if (nRight == HEAP_LEN) {
        //cout << "---" << endl;
        // for each remaining value
        for (unsigned char nVal = 27-(p - a); nVal; --nVal) {
          //dumpState(left, nLeft, median, right, nRight, p, nVal);
          //assert(stateIsValid(left, nLeft, median, right, nRight));
          // get next value
          val = *p++;
          // discard values falling in other heap
          if (val <= median) {
            continue;
          }
          // if val is smaller than smallest in heap, val is new median
          if (val <= right[0]) {
            median = val;
            continue;
          }
          // heap top value becomes new median
          median = right[0];
          // insert val in heap
          parent = 0;
          child = 2;
          while (child < HEAP_LEN) {
            if (right[child-1] < right[child]) {
              child = child-1;
            }
            if (val <= right[child]) {
              break;
            }
            right[parent] = right[child];
            parent = child;
            child = (parent + 1) * 2;
          }
          right[parent] = val;
        }
        return median;
      }
    }
  }
}