/*
 * pi.c
 * Copyright (C) 2008 Paul Kuliniewicz <paul@kuliniewicz.org>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02111-1301, USA.
 */

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

static const char *hex2binary[] = {
	"0000", "0001", "0010", "0011",
	"0100", "0101", "0110", "0111",
	"1000", "1001", "1010", "1011",
	"1100", "1101", "1110", "1111",
};

static char *
toBinary (const void *data, size_t len)
{
	char *result = malloc(8 * len + 1);
	const unsigned char *bytes = data;
	unsigned char byte;
	size_t i;

	if (result != NULL)
	{
		for (i = 0; i < len; i++)
		{
#if BYTE_ORDER == LITTLE_ENDIAN
			byte = bytes[len - i - 1];
#else
			byte = bytes[i];
#endif
			sprintf(result + 8 * i, "%s%s ", hex2binary[(byte >> 4) & 0x0f],
			                                 hex2binary[byte & 0x0f]);
		}
	}

	return result;
}

static void
dumpBinary (const void *data, size_t len)
{
	char *binary = toBinary(data, len);
	printf("%s", binary);
	free(binary);
}

static int
unbiasExponent (const char *binary, size_t exponentBits)
{
	int biased = 0;
	size_t i;

	for (i = 1; i < 1 + exponentBits; i++)
	{
		biased <<= 1;
		if (binary[i] == '1')
			biased++;
	}

	return biased - (1 << (exponentBits - 1)) + 1;
}

static void
dumpFloatingPoint (const void *data, size_t len, size_t exponentBits)
{
	char *binary = toBinary(data, len);
	int exponent = unbiasExponent(binary, exponentBits);
	int fractionBits = 8 * len - 1 - exponentBits;

	/* Note how we have to cheat for the long double's non-implied leading bit */

	int fudge = (len == sizeof(long double) - 2 ? 1 : 0);

	int errorExp2 = exponent - (fractionBits + 1 - fudge);
	double errorExp10 = errorExp2 * log10(2.0);
	double errorDecimal = exp10(errorExp10);

	printf("\tsign: %c\n", binary[0]);
	printf("\texponent: %.*s\n", exponentBits, binary + 1);
	printf("\tfraction: %s\n", binary + 1 + exponentBits);
	printf("\tvalue: %c1.%s * 2^(%d)\n",
		binary[0] == '1' ? '-' : ' ',
		binary + 1 + exponentBits + fudge,
		exponent);
	printf("\terror: +/- 2^(%d) = 10^(%.4f) = %.25f\n",
		errorExp2,
		errorExp10,
		errorDecimal);

	free(binary);
}

int
main (int argc, char *argv[])
{
	float piFloat = M_PI;
	double piDouble = M_PI;
	long double piLD = M_PIl;

	printf("float pi == %.20f\n", piFloat);
	puts("float pi:");
	dumpFloatingPoint(&piFloat, sizeof(piFloat), 8);

	printf("double pi == %.20f\n", piDouble);
	puts("double pi:");
	dumpFloatingPoint(&piDouble, sizeof(piDouble), 11);

	/* This probably only makes sense on x86. */

	printf("long double pi == %.20Lf\n", piLD);
	puts("long double pi:");
	/* sizeof(long double) == 12, but the 2 MSBytes are garbage */
#if BYTE_ORDER == LITTLE_ENDIAN
	dumpFloatingPoint(&piLD, sizeof(piLD) - 2, 15);
#else
	dumpFloatingPoint(((const unsigned char *) &piLD) + 2, sizeof(piLD) - 2, 15);
#endif

	puts("M_PI  == 3.14159265358979323846 exactly");
	puts("M_PIl == 3.1415926535897932384626433832795029L exactly");
	puts("Wikipedia pi == 1.10010010000111111011010101000100010000101101000110000100011010011 * 2^(1)");

	return EXIT_SUCCESS;
}