# Calculate median and quartile in C without sorting the array

Here's a piece of working code wrote in C to calculate the median without putting all values in an array and sorting them. It reads values from a file.

#include <stdio.h>
#include <stdlib.h>

double find_min(FILE *, double, double);
double find_max(FILE *, double, double);

int main(int argc, char *argv[]) {

int CLASSES = 2; // number of parts we create when splitting the set every single pass
double GRANULE = 0.001; // minimum precision desired, should be >=1 if numbers are integers

FILE * file_p; // data file
file_p = fopen(argv[1], "r");

int i, valueIndex1, valueIndex2, N;
// i is the index used in the for loops
// valueIndex1 is the 1st token we read from each line of the data file
// valueIndex2 is the 1st token we read from each line of the data file, in a different read
// N is the number of values in the initial set, we'll check if it's even or odd

double value;
// value is the 2nd token we read from each line of the data file, it's an observation

int count[CLASSES];
// count[n] is the number of elements in the nth part we create during a single pass

double min, max, delta;
// min is the minimum (inclusive) of the elements considered in the current pass
// max is the maximum (exclusive) of the elements considered in the current pass
// delta is the difference between the minimum and maximum values of every single part created during the current pass

double lowerBounds[CLASSES], upperBounds[CLASSES];
// lowerBounds[n] is the minimum value (inclusive) that can be in the nth part created during a pass
// upperBounds[n] is the maximum value (exclusive) that can be in the nth part created during a pass

int sum, tmpsum, limit;
sum = tmpsum = 0;
// sum is the number of elements in the parts before the ones we just create each pass
// tmpsum is the number of elements in the parts before the one we are considering during the pass
// limit is the number of elements the median should be greater than

int index_median;
// index_median is the number of the class containing the median element

short beyondPrecision; // beyondPrecision is raised to 1 if we are creating classes with a range < GRANULE

double median;
// median is where we store the value of the median before returning it

// find minimum and maximum values of the whole set

fscanf(file_p, "%d %lf", &valueIndex1, &min);

max = min;

while (fscanf(file_p, "%d %lf", &valueIndex1, &value) != EOF) {

if (value < min)
min = value;

if (value > max)
max = value;

}

N = ++valueIndex1;
// The number of values is calculated as the last of indexes found + 1
// It's equivalent to the number of lines in the file

max++; // make exclusive

// check if N is even or odd
if (N % 2)
limit = N / 2 + 1;
else
limit = N / 2;

delta = (max - min) / (CLASSES);

// define the borders of the initial frequency classes

for (i = 0; i < CLASSES; ++i) {

lowerBounds[i] = min + i*delta;
upperBounds[i] = min + (i + 1) * delta;
count[i] = 0;

}

// find the class (part) the median is in
// we keep dividing in smaller classes until we meet one of the 2 exit conditions

while (1) {

fclose(file_p);
file_p = fopen(argv[1], "r");

// fill the frequency classes
// during the first pass they are the initial classes, then they will be parts of the chosen part

while (fscanf(file_p, "%d %lf", &valueIndex2, &value) != EOF) {

for (i = 0; i < CLASSES; i++) {

if (value >= lowerBounds[i] && value < upperBounds[i])
count[i]++;

}

}

// select the class containing the median. It will be the first one with (#elements + #elements of classes on its left) > N/2

for (i = 0; i < CLASSES; ++i) {

tmpsum += count[i];

if (tmpsum > limit) {

index_median = i;
break;

}

sum += count[i];

}

min = lowerBounds[index_median];
max = upperBounds[index_median];

// if the number of elements on the left is exactly N/2, then we are done (exit condition satisfied)
// in that case, the median is calculated outside this cycle as the maximum element of the class (part) numbered index_median - 1

if (sum == limit)
break;

// if the median appears many times in the set, then we need a way to figure it out before we create classes that are too small,
// so we exit after reaching a certain fixed precision

if (max - min < GRANULE) {

beyondPrecision = 1;
index_median++;
break;

}

// set the borders of the new frequency classes (for next pass)
// they are created by dividing the class containing the median

delta = (max - min) / (CLASSES);

for (i = 0; i < CLASSES; ++i) {

lowerBounds[i] = min + i*delta;
upperBounds[i] = min + (i + 1) * delta;
count[i] = 0;

}

tmpsum = sum;

}

// now we now which class the median is in and we need to extract it
// we distinguish the 2 cases of having an ever or an odd number of observations (N)

fclose(file_p);
file_p = fopen(argv[1], "r");

if (N % 2 || beyondPrecision) { // N is odd or the median value appears more than one time
median = find_max(file_p, lowerBounds[index_median - 1], upperBounds[index_median - 1]);
} else { // N is even
double a, b;
// a is the maximum value in the left part, where the 1st value necessary to calculate the median
// b is the minimum value in the right part, where the 2nd value necessary to calculate the median

a = find_max(file_p, lowerBounds[index_median - 1], upperBounds[index_median - 1]);

fclose(file_p);
file_p = fopen(argv[1], "r");

b = find_min(file_p, lowerBounds[index_median], upperBounds[index_median]);

median = (a + b) / 2;

}

fclose(file_p);

printf("Our median is %f\n", median);

return (EXIT_SUCCESS);

}

double find_min(FILE * file_p, double lowerBound, double upperBound) {

double value, min;
int valueIndex;

min = upperBound;

while (fscanf(file_p, "%d %lf", &valueIndex, &value) != EOF) {

if (value >= lowerBound && value < upperBound && value < min)
min = value;

}

return min;

}

double find_max(FILE * file_p, double lowerBound, double upperBound) {

double value, max;
int valueIndex;

max = lowerBound;

while (fscanf(file_p, "%d %lf", &valueIndex, &value) != EOF) {

if (value >= lowerBound && value < upperBound && value > max)
max = value;

}

return max;

}


And here is the code to calculate the quartile, based on the same principle.

#include <stdio.h>
#include <stdlib.h>

double find_min(FILE *, double, double);
double find_max(FILE *, double, double);

int main(int argc, char *argv[]) {

int CLASSES = 2; // number of parts we create when splitting the set every single pass
double GRANULE = 0.001; // minimum precision desired, should be >=1 if numbers are integers
double MEDIAN = 34.2275;

FILE * file_p; // data file
file_p = fopen(argv[1], "r");

int i, valueIndex1, valueIndex2, N = 0;
// i is the index used in the for loops
// valueIndex1 is the 1st token we read from each line of the data file
// valueIndex2 is the 1st token we read from each line of the data file, in a different read
// N is the number of values in the initial set, we'll check if it's even or odd

double value;
// value is the 2nd token we read from each line of the data file, it's an observation

int count[CLASSES];
// count[n] is the number of elements in the nth part we create during a single pass

double min, max, delta;
// min is the minimum (inclusive) of the elements considered in the current pass
// max is the maximum (exclusive) of the elements considered in the current pass
// delta is the difference between the minimum and maximum values of every single part created during the current pass

double lowerBounds[CLASSES], upperBounds[CLASSES];
// lowerBounds[n] is the minimum value (inclusive) that can be in the nth part created during a pass
// upperBounds[n] is the maximum value (exclusive) that can be in the nth part created during a pass

int sum, tmpsum, limit;
sum = tmpsum = 0;
// sum is the number of elements in the parts before the ones we just create each pass
// tmpsum is the number of elements in the parts before the one we are considering during the pass
// limit is the number of elements the quartile should be greater than

int index_quartile;
// index_quartile is the number of the class containing the quartile element

short beyondPrecision; // beyondPrecision is raised to 1 if we are creating classes with a range < GRANULE

double quartile;
// quartile is where we store the value of the quartile before returning it

// find minimum and maximum values of the whole set

while (N == 0) {
fscanf(file_p, "%d %lf", &valueIndex1, &value);

if (value < MEDIAN) {

min = value;
max = min;
N++; // N is now set as a counter

}
}

while (fscanf(file_p, "%d %lf", &valueIndex1, &value) != EOF) {

if (value < MEDIAN) {

if (value < min)
min = value;

if (value > max)
max = value;

N++;
}

}

max++; // make exclusive

// check if N is even or odd
if (N % 2)
limit = N / 2 + 1;
else
limit = N / 2;

delta = (max - min) / (CLASSES);

// define the borders of the initial frequency classes

for (i = 0; i < CLASSES; ++i) {

lowerBounds[i] = min + i*delta;
upperBounds[i] = min + (i + 1) * delta;
count[i] = 0;

}

// find the class (part) the quartile is in
// we keep dividing in smaller classes until we meet one of the 2 exit conditions

while (1) {

fclose(file_p);
file_p = fopen(argv[1], "r");

// fill the frequency classes
// during the first pass they are the initial classes, then they will be parts of the chosen part

while (fscanf(file_p, "%d %lf", &valueIndex2, &value) != EOF) {

if (value < MEDIAN) {

for (i = 0; i < CLASSES; i++) {

if (value >= lowerBounds[i] && value < upperBounds[i])
count[i]++;

}

}

}

// select the class containing the quartile. It will be the first one with (#elements + #elements of classes on its left) > N/2

for (i = 0; i < CLASSES; ++i) {

tmpsum += count[i];

if (tmpsum > limit) {

index_quartile = i;
break;

}

sum += count[i];

}

min = lowerBounds[index_quartile];
max = upperBounds[index_quartile];

// if the number of elements on the left is exactly N/2, then we are done (exit condition satisfied)
// in that case, the quartile is calculated outside this cycle as the maximum element of the class (part) numbered index_quartile - 1

if (sum == limit)
break;

// if the quartile appears many times in the set, then we need a way to figure it out before we create classes that are too small,
// so we exit after reaching a certain fixed precision

if (max - min < GRANULE) {

beyondPrecision = 1;
index_quartile++;
break;

}

// set the borders of the new frequency classes (for next pass)
// they are created by dividing the class containing the quartile

delta = (max - min) / (CLASSES);

for (i = 0; i < CLASSES; ++i) {

lowerBounds[i] = min + i*delta;
upperBounds[i] = min + (i + 1) * delta;
count[i] = 0;

}

tmpsum = sum;

}

// now we now which class the quartile is in and we need to extract it
// we distinguish the 2 cases of having an ever or an odd number of observations (N)

fclose(file_p);
file_p = fopen(argv[1], "r");

if (N % 2 || beyondPrecision) { // N is odd or the quartile value appears more than one time
quartile = find_max(file_p, lowerBounds[index_quartile - 1], upperBounds[index_quartile - 1]);
} else { // N is even
double a, b;
// a is the maximum value in the left part, where the 1st value necessary to calculate the quartile
// b is the minimum value in the right part, where the 2nd value necessary to calculate the quartile

a = find_max(file_p, lowerBounds[index_quartile - 1], upperBounds[index_quartile - 1]);

fclose(file_p);
file_p = fopen(argv[1], "r");

b = find_min(file_p, lowerBounds[index_quartile], upperBounds[index_quartile]);

quartile = (a + b) / 2;

}

fclose(file_p);

printf("Our quartile is %f\n", quartile);

return (EXIT_SUCCESS);

}

double find_min(FILE * file_p, double lowerBound, double upperBound) {

double value, min;
int valueIndex;

min = upperBound;

while (fscanf(file_p, "%d %lf", &valueIndex, &value) != EOF) {

if (value >= lowerBound && value < upperBound && value < min)
min = value;

}

return min;

}

double find_max(FILE * file_p, double lowerBound, double upperBound) {

double value, max;
int valueIndex;

max = lowerBound;

while (fscanf(file_p, "%d %lf", &valueIndex, &value) != EOF) {

if (value >= lowerBound && value < upperBound && value > max)
max = value;

}

return max;

}


I remember both of them working. I wrote them together with a friend less than a year ago. I forgot the logic since then (fortunately I commented it!), but now that I'm about to use it, I'd like to have it reviewed.

I'd like your opinion on this, correctness and all.

Maybe there's an elegant way to put them together?

• Interesting question, but I would recommend refamiliarizing yourself with the code and reviewing it yourself first. It may be better to clean it up before asking here. You're likely to get more meaningful answers and you'll be more prepared to answer questions reviewers may have. – RubberDuck Mar 27 '15 at 0:34
• I reviewed it back then. It took hours to figure out and comment. I also have a couple of pages describing the method that I could post here. I'd like to hear somebody else first, though. Hint, the trickiest part was figuring out all the corner cases. – Agostino Mar 27 '15 at 0:38
• Right. You reviewed it back then. I often put code aside for a few days and come back to it after a few days to give a self review. If you give yourself an earnest review now, you may be surprised at how much better you can make it all on your own. If you take it that first step yourself, reviewers can push you even farther beyond that. But suit yourself. It was just a friendly suggestion. – RubberDuck Mar 27 '15 at 0:44

In the beginning of your code where you delcare your find_min and find_max functions, I recommend that you put the parameter names along with their type.

If your code were to be longer, it would make quickly glancing at your code easier as the code already shows at the top which argument means what in a function call.

To me, a function called find_min or find_max should return which of two number arguments passed is either the greatest or the least.

To better describe the functions, I recommend calling them something like ffind_max or f_find_min.

To be safe, you should check that the return of fopen is not NULL.

Also to be safe, you should make sure that there is even a argv[1] to begin with.

(I think) A faster and more efficient way to find out if a number is odd or even than doing % 2 would be to & 1.

This will result in a 1 for an odd number and a 0 for an even number.

Here is why:

0111    = 7
&
0001    = 1


ANDing with 0 will make a bit 0, but ANDing with 1 will return the bit being ANDed.

And, to enhance understanding, I would create a macro for finding out if a number is even or odd:

#define IS_ODD(n) n & 1


That way, every time you are checking if an number is odd, you don't have to write a comment about it.

These are just small points about your code. I'll look more into the content of your code tomorrow.

• Thanks. Names for the find_ functions could be better, yes, although the parameter names do help understanding. Checking on fopen is a good idea. The trick for finding if a number is even is nice, maybe a just a little overkill. – Agostino Jun 16 '15 at 3:09
• I wouldn't consider the trick to be overkill at all; it's going to be easier on your processor. Think about it: by doing %, on a x86, the processor is going to run the div on either a 32 bit value and a 16 bit value, or a 64 bit value and a 32 bit value (probably the latter). The div instruction has a latency of 40. However, if the processor uses the or instruction, it will have an easier time since the or instruction only has a latency of 0.5 or 1. – SirPython Jun 20 '15 at 21:58
• OK, maybe this is a scenario that can benefit from the use of bitwise operators. I can use the extra speed. Other suggestions? – Agostino Jun 21 '15 at 1:29

## Segmentation violation

I got the program to crash with following input:

0 0
1 1
2 2
3 3
4 -1.0e+17


The problem was that the max-min calculation overflowed the precision of a double, and then the two buckets no longer covered the entire range. So in the loop that counted how many numbers were in each bucket, some numbers didn't belong to any buckets. After that, index_median didn't get set (it had an uninitialized value), and then index_median was used to index into the lowerBounds array. If you want to reproduce it every time, initialize index_median to some large value such as 0x10000000 and then run the program on that input file.

I think one of the flaws of this algorithm is using floating point arithmetic without checking for overflows and precision loss. There may be other unhandled corner cases other than the one I found.

• Thanks. Is this more of matter of input sanitation or an actual loss of precision during calculations? What would you suggest to solve/mitigate this? – Agostino Jun 16 '15 at 3:05
• @Agostino Well to begin with, I'm not really sure why you are choosing to find the median this way instead of just using an array. This way is going to be a lot slower. But if you insist on not creating an array, I could think of a way that can solve the problem without requiring so many floating point operations. Basically it would be something like a binary search for the median, similar to what you are doing but with mostly compare operations and very few add/subtract/divide operations. – JS1 Jun 16 '15 at 3:28
• Too much data to store it in memory. That's the reason. These files are huge data logs, and this is to cross-check some Apache Hadoop results. There's aren't many division operations in my code. Again, is it an input problem, or an actual loss of precision during run-time? – Agostino Jun 16 '15 at 17:03
• @Agostino It's a loss of precision during runtime. I think to fix you should set upperBounds[CLASSES-1] = max. That will fix the problem I found. Although as I said before, I'm not sure if that is the only problem. – JS1 Jun 16 '15 at 17:06
• Thanks. If you want to add some more comments, you're welcome to do it. – Agostino Jun 16 '15 at 19:06