blob: 995176afb7623bfa841eae8cb2245e124698eb7a [file] [log] [blame]
#include <stdint.h>
#include <stdio.h>
#include <math.h>
//the bits of 2/pi after the initial 1.
const char *pi = "45F306DC9C882A53F84EAFA3EA69BB81B6C52B3278872083FCA2C757BD778AC36E48DC74849BA5C00C925DD413A32439FC3BD63962534E7DD1046BEA5D768909D338E04D68BEFC827323AC7306A673E93908BF177BF250763FF12FFFBC0B301FDE5E2316B414DA3EDA6CFD9E4F96136E9E8C7ECD3CBFD45AEA4F758FD7CBE2F67A0E73EF14A525D4D7F6BF623F1ABA10AC06608DF8F6D757E19F784135E86C3B53C722C2BDCC3610CB330ABE2940D0811BFFB1009AE64E620C0C2AAD94E75192C1C4F78118D68F883386CF9BB9D0125506B388ED172C394DBB5E89A2AE320A7D4BFE0E0A7EFC67D06585BC9F3064FB77867A4DDED63CBDF13E74";
int get_bit( int offset )
{
int position = offset / 4;
int nibble = pi[ position ];
int value;
int shift = 3 - (offset % 4 );
if( nibble >= '0' && nibble <= '9' )
value = nibble - '0';
else
{
if( nibble >= 'A' && nibble <= 'F' )
value = nibble - 'A' + 10;
else
return -1;
}
return (value >> shift) & 1;
}
int main( void )
{
const int bits_per_float = 24;
const int mantissa_size = 52;
const int bias = 1023;
int currentp = -1; //-1 is 2/pi, -2 would be 1/pi
int nextp = currentp - 1;
int offset = 0;
uint64_t mantissa = 0;
int i, bit;
union
{
uint64_t u;
double d;
}value, exponent, startExp, expstep;
long double v = 0;
long double test = 2.0L/3.14159265358979323846264338327950288L;
startExp.u = mantissa_size + currentp - bits_per_float;
expstep.u = bias - bits_per_float;
startExp.u += bias;
startExp.u <<= mantissa_size;
expstep.u <<= mantissa_size;
// startExp.u &= 0x7ff0000000000000ULL;
startExp.u = 5; // stick a simple constant in there for gdb testing -- jsm
expstep.u &= 0x7ff0000000000000ULL; /* put a breakpoint here */
while( currentp + bias > 0 )
{
//set up exponent
value.u = 0;
for( i = 0; i < bits_per_float; i++ )
{
bit = get_bit( offset++ );
if( bit < 0 )
goto exit;
value.u = value.u << 1;
value.u |= bit;
nextp--;
}
value.u |= startExp.u;
value.d -= startExp.d;
v += value.d;
printf( "0x1.%13.13llxp%d, //%17.21g (%17.21g, %17.21g, %17.21g) %17.21g\n", value.u & 0x000FFFFFFFFFFFFFULL, ilogb(value.d), value.d, (double) test, (double) v, (double) (test - v), expstep.d );
startExp.d *= expstep.d;
}
exit:
return 0;
}