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

#ifdef __MINGW32__
#include <fenv.h>
#endif


// show the hex and decimal interpretations of the double
union U{double d;unsigned char c[8];};
void display(union U x)
{
  for(int b=7;b>=0;b--)printf("%02x",x.c[b]);
  printf(" %.21f",x.d);
}


int main(int argc, const char**argv)
{

#ifdef __MINGW32__
  if (argc >= 2)
    {
    int nrm, n;
        n = sscanf (argv[1], "%i", &nrm);
        if (n == 1) switch (nrm)
        {
            case 0:
            case 0x400:
            case 0x800:
            case 0xc00:
                fesetround (nrm);
        }
    }

  printf ("Rmode $%08x\n", fegetround ());
#endif

  printf("DECIMAL_DIG=%d\n",DECIMAL_DIG);

#define PASTE(x) {.s=#x,.d1=strtod(#x,0),.d2=x}
  struct S {const char*s;union U d1,d2;}q[]=
    {
      // test the rounding with the last digits...
      // first right on the boundary, to test round-to-even
      PASTE(1.0000000000000000000000000000000000000000000000000000000), // 1
      PASTE(1.0000000000000001110223024625156540423631668090820312500), // 1 + 2^-53
      PASTE(1.0000000000000002220446049250313080847263336181640625000), // 1 + 2^-52
      PASTE(1.0000000000000003330669073875469621270895004272460937500), // 1 + 2^-52 +2^-53
      {0},
      // now just above the boundary
      PASTE(1.00000000000000000001),
      PASTE(1.00000000000000011103),
      PASTE(1.00000000000000022205),
      PASTE(1.00000000000000033307),
      {0},
      // now just below the boundary
      PASTE(0.99999999999999999999),
      PASTE(1.00000000000000011101),
      PASTE(1.00000000000000022204),
      PASTE(1.00000000000000033301),
      {0},
      PASTE(0.105),
    };
 
  printf("                 1.23456789012345678901234567890\n"); // for counting sig figs
  for(int a=0;a<sizeof(q)/sizeof(struct S);a++)
    {
      if(!q[a].s){printf("\n");continue;}
      display(q[a].d1);
      printf(" ");
      display(q[a].d2);
      printf("\n");
    }
  
  printf("\n");

  union U w[]={0.105,0.115,0.125,0.135,0.145,0.155,0.375};
  for(int a=0;a<sizeof(w)/sizeof(union U);a++)
    {
      display(w[a]);
      printf("\n");
    }

  printf("%.40f\n",0.105);

  // here's some numbers with an exact representation in double
  printf("%.20f %.2f\n",.125,.125); 
  printf("%.20f %.2f\n",.375,.375); 
  printf("%.20f %.13f\n",0.95367431640625000000,0.95367431640625000000);
  printf("%.20f %.17f\n",3.5096855163574218750,3.5096855163574218750); // 1 + 657899/2^18
  printf("%.20f %.15f\n",-1.5035324096679687500,-1.5035324096679687500);

  // how many digits precision does gcc use when making double constants?
  display((union U){.d=1.000000000000000333066907387546962127089500427246});
}

// Useful Mathematica stuff...
/*
  Table[{x/1000.,
            BaseForm[FromDigits[RealDigits[#2, 2] // First // Rest, 2], 16],
            N[#2*2^(-53 + #[[2]]),
              41]} &[#,
        Ceiling[FromDigits[#[[1]], 2]/
            2]] &[RealDigits[x/1000, 2, 54]], {x, 105, 155, 10}]

*/