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

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

void smalldouble_to_cell(void*p, double d)
{
  union u u;
  u.d = d;
  unsigned long long rest = u.u & 0x7fffffffffffffff;
  unsigned long long packed;
  if (rest > 0x7ff0000000000000)
    /* NaN */
    packed = u.u & 0xfffffffffffffffe;
  else 
    {
      unsigned long long sign = u.u & 0x8000000000000000;
      if (rest >= 0x3ff0000000000000)
    rest = 0x3ff0000000000000;
      packed = sign | (rest << 1);
    }
  memcpy(p, &packed, 8);
}

void double_to_cell(void *p, double d)
{
  smalldouble_to_cell(p, ldexp(d, -512));
}

double smalldouble_from_cell(void*p)
{
  union u u;
  unsigned long long l;
  memcpy(&l, p, 8);
  unsigned long long sign = l & 0x8000000000000000;
  unsigned long long rest = (l & 0x7fffffffffffffff) >> 1;
  if (rest >= 0x3ff0000000000000) rest <<= 1;
  u.u = sign | rest;
  return u.d;
}

double double_from_cell(void*p)
{
  double d = smalldouble_from_cell(p);
  return ldexp(d, +512);
}

void float_cell_add(void *d, void *s1, void *s2)
{
  double d1 = smalldouble_from_cell(s1);
  double d2 = smalldouble_from_cell(s2);
  smalldouble_to_cell(d, d1 + d2);
}

void float_cell_sub(void *d, void *s1, void *s2)
{
  double d1 = smalldouble_from_cell(s1);
  double d2 = smalldouble_from_cell(s2);
  smalldouble_to_cell(d, d1 - d2);
}

void float_cell_mul(void *d, void *s1, void *s2)
{
  double d1 = smalldouble_from_cell(s1);
  double d2 = smalldouble_from_cell(s2);
  smalldouble_to_cell(d, ldexp(d1, +512) * d2);
}

void float_cell_div(void *d, void *s1, void *s2)
{
  double d1 = smalldouble_from_cell(s1);
  double d2 = smalldouble_from_cell(s2);
  smalldouble_to_cell(d, ldexp(d1/d2, -512));
}

int main()
{
  unsigned long long l1, l2, l3;

  double_to_cell(&l1, 3.0);
  double_to_cell(&l2, 7.0);
  float_cell_add(&l3, &l1, &l2);
  printf("3.0 + 7.0 = %f\n", double_from_cell(&l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  double_to_cell(&l1, nextafter(1.0, 2.0));
  double_to_cell(&l2, nextafter(1.5, 0.0));
  float_cell_add(&l3, &l1, &l2);
  printf("(1+) + (1.5-) = %f\n", double_from_cell(&l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  double_to_cell(&l1, 1.0e-168);
  double_to_cell(&l2, 1.3e-168);
  float_cell_add(&l3, &l1, &l2);
  printf("(1.0-168) + (1.3-168) = %e\n", double_from_cell(&l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  double_to_cell(&l1, 1.0e-90);
  double_to_cell(&l2, 1.3e-78);
  float_cell_mul(&l3, &l1, &l2);
  printf("(1.0-90) * (1.3-78) = %e\n", double_from_cell(&l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  double_to_cell(&l1, 1.0e+70);
  double_to_cell(&l2, 1.3e+78);
  float_cell_mul(&l3, &l1, &l2);
  printf("(1.0+70) * (1.3+78) = %e\n", double_from_cell(&l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  double_to_cell(&l1, 1.0e+90);
  double_to_cell(&l2, 1.3e+78);
  float_cell_mul(&l3, &l1, &l2);
  printf("(1.0+90) * (1.3+78) = %e\n", double_from_cell(&l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  double_to_cell(&l1, 1.0e+90);
  double_to_cell(&l2, 1.3e+78);
  float_cell_div(&l3, &l1, &l2);
  printf("(1.0+90) / (1.3+78) = %e\n", double_from_cell(&l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  double_to_cell(&l1, 1.0e200);
  double_to_cell(&l2, 1.0e-200);
  float_cell_div(&l3, &l1, &l2);
  printf("inf / 0 = %e\n", double_from_cell(&l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  double_to_cell(&l1, 1.0e200);
  double_to_cell(&l2, 1.0e-200);
  float_cell_mul(&l3, &l1, &l2);
  printf("inf * 0 = %e\n", double_from_cell(&l3));
  printf("representations: %llx %llx %llx\n", l1, l2 , l3);

  return 0;
}