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

union u {
  unsigned long long u; /* assumed to be 64-bit */
  double d;
};

typedef unsigned long long boxed;

boxed smalldouble_to_boxed(double d)
{
  union u u;
  u.d = d;
  boxed rest = u.u & 0x7fffffffffffffff;
  boxed packed;
  if (rest > 0x7ff0000000000000)
    /* NaN */
    packed = u.u;
  else 
    {
      boxed sign = u.u & 0x8000000000000000;
      if (rest >= 0x3ff0000000000000)
    /* inf */
	rest = 0x3ff0000000000000;
      packed = sign | (rest << 1);
    }
  return packed | 1;
}

boxed double_to_boxed(double d)
{
  return smalldouble_to_boxed(ldexp(d, -512));
}

double smalldouble_from_boxed(boxed l)
{
  union u u;
  boxed sign = l & 0x8000000000000000;
  boxed rest = (l & 0x7fffffffffffffff) >> 1;
  if (rest >= 0x3ff0000000000000) rest <<= 1;
  u.u = sign | rest;
  return u.d;
}

double double_from_boxed(boxed l)
{
  double d = smalldouble_from_boxed(l);
  return ldexp(d, +512);
}

boxed smalldouble_add(boxed s1, boxed s2)
{
  double d1 = smalldouble_from_boxed(s1);
  double d2 = smalldouble_from_boxed(s2);
  return smalldouble_to_boxed(d1 + d2);
}

boxed smalldouble_sub(boxed s1, boxed s2)
{
  double d1 = smalldouble_from_boxed(s1);
  double d2 = smalldouble_from_boxed(s2);
  return smalldouble_to_boxed(d1 - d2);
}

boxed smalldouble_mul(boxed s1, boxed s2)
{
  double d1 = smalldouble_from_boxed(s1);
  double d2 = smalldouble_from_boxed(s2);
  return smalldouble_to_boxed(ldexp(d1, +512) * d2);
}

boxed smalldouble_div(boxed s1, boxed s2)
{
  double d1 = smalldouble_from_boxed(s1);
  double d2 = smalldouble_from_boxed(s2);
  return smalldouble_to_boxed(ldexp(d1/d2, -512));
}

int main()
{
  boxed l1, l2, l3;

  l1 = double_to_boxed(3.0);
  l2 = double_to_boxed(7.0);
  l3 = smalldouble_add(l1, l2);
  printf("3.0 + 7.0 = %f\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  l1 = double_to_boxed(nextafter(1.0, 2.0));
  l2 = double_to_boxed(nextafter(1.5, 0.0));
  l3 = smalldouble_add(l1, l2);
  printf("(1+) + (1.5-) = %f\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  l1 = double_to_boxed(1.0e-168);
  l2 = double_to_boxed(1.3e-168);
  l3 = smalldouble_add(l1, l2);
  printf("(1.0-168) + (1.3-168) = %e\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  l1 = double_to_boxed(1.00000000001e-155);
  l2 = double_to_boxed(1.0e-155);
  l3 = smalldouble_sub(l1, l2);
  printf("(1.0...01-155) + (1.0-155) = %.16e\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  l1 = double_to_boxed(1.0e-90);
  l2 = double_to_boxed(1.3e-78);
  l3 = smalldouble_mul(l1, l2);
  printf("(1.0-90) * (1.3-78) = %e\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  l1 = double_to_boxed(1.0e+70);
  l2 = double_to_boxed(1.3e+78);
  l3 = smalldouble_mul(l1, l2);
  printf("(1.0+70) * (1.3+78) = %e\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  l1 = double_to_boxed(1.0e+90);
  l2 = double_to_boxed(1.3e+78);
  l3 = smalldouble_mul(l1, l2);
  printf("(1.0+90) * (1.3+78) = %e\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  l1 = double_to_boxed(1.0e+90);
  l2 = double_to_boxed(1.3e+78);
  l3 = smalldouble_div(l1, l2);
  printf("(1.0+90) / (1.3+78) = %e\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  l1 = double_to_boxed(1.0e200);
  l2 = double_to_boxed(1.0e-200);
  l3 = smalldouble_div(l1, l2);
  printf("inf / 0 = %e\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  l1 = double_to_boxed(1.0e200);
  l2 = double_to_boxed(1.0e-200);
  l3 = smalldouble_mul(l1, l2);
  printf("inf * 0 = %e\n", double_from_boxed(l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  return 0;
}