/* Some sample C code for the binmedian algorithm. 
 * Note: I haven't written the code for n even yet. */

#define SWAP(a,b) temp=a;a=b;b=temp;

float binmedian(int n, float *x) {
  // Compute the mean and standard deviation
  float sum = 0;
  int i;
  for (i = 0; i < n; i++) {
    sum += x[i];
  }
  float mu = sum/n;

  sum = 0;
  for (i = 0; i < n; i++) {
    sum += (x[i]-mu)*(x[i]-mu);
  }
  float sigma = sqrt(sum/n);

  // Bin x across the interval [mu-sigma, mu+sigma]
  int bottomcount = 0;
  int bincounts[1001];
  for (i = 0; i < 1001; i++) {
    bincounts[i] = 0;
  }
  float scalefactor = 1000/(2*sigma);
  float leftend =  mu-sigma;
  float rightend = mu+sigma;
  int bin;

  for (i = 0; i < n; i++) {
    if (x[i] < leftend) {
      bottomcount++;
    }
    else if (x[i] < rightend) {
      bin = (int)((x[i]-leftend) * scalefactor);
      bincounts[bin]++;
    }
  }

  // If n is odd
  if (n & 1) {
    // Recursive step
    int k, r, count, medbin;
    float oldscalefactor, oldleftend;
    int oldbin;
    float temp;

    k = (n+1)/2;
    r = 0;

    for (;;) {
      // Find the bin that contains the median, and the order
      // of the median within that bin
      count = bottomcount;
      for (i = 0; i < 1001; i++) {
        count += bincounts[i];

        if (count >= k) {
          medbin = i;
          k = k - (count-bincounts[i]);
          break;
        }
      }

      bottomcount = 0;
      for (i = 0; i < 1001; i++) {
        bincounts[i] = 0;
      }
      oldscalefactor = scalefactor;
      oldleftend = leftend;
      scalefactor = 1000*oldscalefactor;
      leftend = medbin/oldscalefactor + oldleftend;
      rightend = (medbin+1)/oldscalefactor + oldleftend;

      // Determine which points map to medbin, and put
      // them in spots r,...n-1
      i = r; r = n;
      while (i < r) {
        oldbin = (int)((x[i]-oldleftend) * oldscalefactor);

        if (oldbin == medbin) {
          r--;
          SWAP(x[i],x[r]);

          // Re-bin on a finer scale
          if (x[i] < leftend) {
            bottomcount++;
          }
          else if (x[i] < rightend) {
            bin = (int)((x[i]-leftend) * scalefactor);
            bincounts[bin]++;
          }
        }
        else {
          i++;
        }
      }

      // Stop if all points in medbin are the same
      int samepoints = 1;
      for (i = r+1; i < n; i++) {
        if (x[i] != x[r]) {
          samepoints = 0;
          break;
        }
      }
      if (samepoints) {
        return x[r];
      }

      // Stop if there's <= 20 points left
      if (n-r <= 20) {
        break;
      }
    }

    // Perform insertion sort on the remaining points,
    // and then pick the kth smallest
    float a;
    int j;
    for (i = r+1; i < n; i++) {
      a = x[i];
      for (j = i-1; j >= r; j--) {
        if (x[j] < a) {
          break;
        }
        x[j+1] = x[j];
      }
      x[j+1] = a;
    }

    return x[r-1+k];
  }

  // If n is even (not implemented yet)
  else {
    return 0;
  }
}