/* floating_point_precision.c
 *
 * This is a simple example how you can test the precision of a
 * floating point number on your system
 *
 * Copyright (C) 2009 by Heiko Voigt <hvoigt@hvoigt.net>
 *
 * Compile with:
 *
 *    gcc floating_point_precision.c -o floating_point_precision
 *
 * Description:
 *
 * A floating point number can only store its values with a certain
 * precision which depends on the size of its mantissa. You can simply
 * investigate this by performing the simple calculation:
 *
 *    x = 2^b + 1			(1)
 *
 * and then trying to substract this in the same order
 *
 *    0 = x - 2^b - 1			(2)
 *
 * where b is the number of bit precision which should be tested.
 *
 * This test relies on the property that when the storage
 * precision of double is exeeded the + 1 of calculation (1)
 * fails and because we substract 2^b in (2) first, the second
 * operation of ( - 1 ) can succeed.
 */

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

#define bitsof(x) sizeof(x) * 8

/* the macro below is a C type independent implementation of this function */
#if 0
double test_precision(int mantissa_bits)
{
	double d;
	double one = 1.0;
	double two = 2.0;
	double bits = mantissa_bits;

	/* calculate forward */
	d = one * pow(two, bits) + one;

	/* try to inverse it to zero */
	d = d - pow(two, bits) - one;

	return d;
}
#endif

#define test_precision(d, mantissa_bits) \
	typeof(d) one = 1.0; \
	typeof(d) two = 2.0; \
	typeof(d) bits = mantissa_bits; \
\
	/* calculate forward */ \
	d = one * pow(two, bits) + one; \
\
	/* try to inverse it to zero */ \
	d = d - pow(two, bits) - one;

double test_double_precision(int mantissa_bits)
{
	double d;
	test_precision(d, mantissa_bits);
	return d;
}

float test_float_precision(int mantissa_bits)
{
	float d;
	test_precision(d, mantissa_bits);
	return d;
}

int main(int argc, char *argv[])
{
	int i;
	double d;

	printf("Some system information:\n");
	printf("%-15s = %lu\n", "bitsof(int)", bitsof(int));
	printf("%-15s = %lu\n", "bitsof(long)", bitsof(long));
	printf("%-15s = %lu\n", "bitsof(float)", bitsof(float));
	printf("%-15s = %lu\n", "bitsof(double)", bitsof(double));
	printf("\n");

	printf("Testing float precision:\n");
	for (i=0; i < bitsof(float); i++) {
		d = test_float_precision(i);
		if (d < 0) {
			break;
		}
	}
	printf("Fails on Bit %d Result: %.1lf instead of 0.0\n", i, d);

	printf("Testing double precision:\n");
	for (i=0; i < bitsof(double); i++) {
		d = test_double_precision(i);

		if (d < 0) {
			break;
		}
	}
	printf("Fails on Bit %d Result: %.1lf instead of 0.0\n", i, d);

	return 0;
}