/* quick hack to deal with IEEE-754 Doubles
 * September 16, 2008
 */
#include <stdio.h>
#include <math.h>

/* WRite Big-Endian 64-bit value */
#define WR_BE64(dest, offset, value) do { \
		unsigned long long _tmp=value; \
		(dest)[offset]=((_tmp)>>56)&255; \
		(dest)[(offset)+1]=((_tmp)>>48)&255; \
		(dest)[(offset)+2]=((_tmp)>>40)&255; \
		(dest)[(offset)+3]=((_tmp)>>32)&255; \
		(dest)[(offset)+4]=((_tmp)>>24)&255; \
		(dest)[(offset)+5]=((_tmp)>>16)&255; \
		(dest)[(offset)+6]=((_tmp)>>8)&255; \
		(dest)[(offset)+7]=(_tmp)&255; \
	} while(0)


/* ReaD Big-Endian 64-bit value */
#define RD_BE64(src, offset) ( \
	((src)[(offset)]*0x100000000000000ULL)| \
	((src)[(offset)+1]*0x1000000000000ULL)| \
	((src)[(offset)+2]*0x10000000000ULL)| \
	((src)[(offset)+3]*0x100000000ULL)| \
	((src)[(offset)+4]*0x1000000U)| \
	((src)[(offset)+5]*0x10000U)| \
	((src)[(offset)+6]*0x100U)| \
	((src)[(offset)+7]))

void ieee754_encode64(void *dest, double f) {
	const unsigned
		bits=64,
		expbits=11,
		sigbits=bits-expbits-1;
	const int bias=(1<<(expbits-1))-1;
	double norm;
	int shift;
	unsigned sign, exponent;
	unsigned long long sig;

	sign=f<0;
	norm=fabs(f);
	for(shift=0;norm>=2.;norm/=2.,shift++) ;
	for(;norm<1.;norm*=2.,shift--) ;
	norm-=1.;
	sig=norm*((1ll<<sigbits)+.5);
	exponent=shift+bias;

	sig|=((unsigned long long)sign<<(bits-1))|((unsigned long long)exponent<<sigbits);
	WR_BE64((unsigned char*)dest, 0, sig);
}

double ieee754_decode64(void *src) {
	const unsigned
		bits=64,
		expbits=11,
		sigbits=bits-expbits-1;
	const int bias=(1<<(expbits-1))-1;
	unsigned long long i;
	int shift;
	double res;
	i=RD_BE64((unsigned char*)src, 0);
	if(i==0) return 0.;
	printf("hex=0x%llX i=0x%llX\n", *(long long*)src, i);
	res=i&((1LL<<sigbits)-1); /* mask off significant */	
	printf("res=%f hex=0x%llX\n", res, *(long long*)&res);
	res/=(1LL<<sigbits);
	res+=1.;
	shift=((i>>sigbits)&((1LL<<expbits)-1))-bias;
	printf("shift=%d res=%f hex=0x%llX\n", shift, res, *(long long*)&res);
	while(shift>0) { res*=2.0; shift--; }
	while(shift<0) { res/=2.0; shift++; }
	if(i>>(bits-1)) {
		res=-res;
	}
	return res;
}

int main(void) {
	char buf[8];
	double tmp, pi;

	ieee754_encode64(buf, pi=M_PI);
	tmp=ieee754_decode64(buf);
	printf("in=%f hex=0x%llX out=%f hex=0x%llX buf=0x%llX ok? %d\n", pi, *(long long*)&pi, tmp, *(long long*)&tmp, *(long long*)buf, tmp==pi);
	return 0;
}