Studies in numerical drift

gistfile1.c

/*
Sometimes you get bored on flights between MDW and SFO and decide
to explore numerical drift and round off errors in doubles. 
Here's the difference between a linear fold over a list and 
using pyramidal/tree folds.

Theoretically the tree-like fold would minimize round-off error
as the operation propagates up the tree, but I'd need to double
check my numerical methods stuff to be sure :-).

Either way, floating point operations don't commute and this demonstrates
it reasonably well.

I've been building with this:
gcc -wAll -g -std=c89 -o drift drift.c

On my laptop it breaks at 184 with fairly significant results:

                   VVV
acc:  830425644985504.625000000000000 
acc2: 830425644985503.625000000000000 
									 ^^^
Total difference 1.0000000000000000000000000

*/

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

#define SEED 25

double pyramid_fold(double *numbers, int length){
	int rounds = log10(length) * 2;
	int i = 0;
	int prev_round = length;
	int round_size = 0;
	double result = 0.0;

	printf("Rounds : %d\n",rounds);
	double **mem_ref = calloc(rounds+1,sizeof(double*));
	int round_counter = 0;

	double *prev_holder;
	prev_holder = numbers;
	double *place_holder;
  do{
		round_size = floor(prev_round / 2.0);
		printf("Size: %d\n", round_size);

		place_holder = calloc(round_size,sizeof(double));
		mem_ref[round_counter++] = place_holder;

		int i = 0;
		int j = 0;
		for(i = 0; i < prev_round;){
			if((i+1) <= (prev_round - 1)){
				printf("i's %d %d\n", i, i+1);
				place_holder[j] = prev_holder[i] * prev_holder[i+1];
				printf("%5.15f = %5.15f %5.15f\n",place_holder[j],prev_holder[i],prev_holder[i+1]);
			}
			else{
				/* There's a trailing term... */
				printf("Odd sized round\n");
				printf("Merging in %5.15f \n",prev_holder[i]);
				if(j > 0){
					double debug_val = place_holder[j-1];
					place_holder[j-1] = place_holder[j-1] * prev_holder[i];
					printf("%5.15f = %5.15f %5.15f\n",place_holder[j-1],debug_val,prev_holder[i]);
				}
				else{
					/* Final round & a NOOP */
					printf("Final round + NOOP %5.15f\n", prev_holder[i]);
					place_holder[j] = prev_holder[i];
				}
			}

			i += 2;
			++j;
		}

		if(round_size == 0){
			printf("This was the final round.\n");
			printf("Round size %d\n", round_size);
			printf("prev_round %d\n", prev_round);
			printf("%5.15F\n",place_holder[0]);
			result = place_holder[0];
		}

		prev_holder = place_holder;
		prev_round = round_size;
	}while(round_size != 0);

	for(i = 0; i < rounds; i++){
		free(mem_ref[i]);
	}
	free(mem_ref);

	return result;
}

int main(){
	double acc = 0.0;
	double acc2 = 0.0;
	double threshold = 0.001;
	int i;
	int size = 1;

	while(fabs(acc - acc2) < threshold && size < 10000){
		printf("Size %d\n", size);
		srand(SEED);
		double *numbers = calloc(size,sizeof(double));
		double *numbers2 = calloc(size,sizeof(double));

		acc = 0.0;
		for(i = 0; i < size; ++i){
			numbers[i] = random();
			numbers2[i] = random();
			numbers[i] = numbers[i] / numbers2[i];
		}

		printf("Finished generating numbers...\n");
		acc = numbers[0];
		for(i = 1; i < size; ++i){
			acc = acc * numbers[i];
		}

		acc2 = pyramid_fold(numbers,size);

		printf("acc:  %5.15F \n", acc);
		printf("acc2: %5.15f \n", acc2);
		printf("Total difference %5.25f \n", fabs(acc - acc2));
		printf("Equality test? %d\n", acc == acc2);
		free(numbers);
		free(numbers2);
		size++;
	}
	printf("Size of array when test broke: %d\n", size);
	return 0;
}