#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
typedef long unsigned int uL;
int const b = 6; // 1 2 3 4 5 6
int const epsilon = 0x1fd23;
void pL(char* n, uL x){printf("%016lX = %s\n", x, n);}
void va(uL a){assert(a>=(1L<<63) && !(a&((1<<11)-1)));}
char p = 1;

typedef struct{uL q; uL r;} res;
res divq(uL a, uL b, uL d);
res mulq(uL a, uL b);
short tab[256];
// dv returns 2*(the quotient) using the proposed scheme.
// The quotient may not be normalized.
uL dv(uL const num, uL const den){
  if(0<1){va(num); va(den);}
  int N = tab[255&(den>>(64-9))];
  if(p) printf("N=%d=0x%03x\n", N, N);
  uL nn = N*(num>>10), nd = N*(den>>8);
  // neither of above overflow!
  if(p) {pL("Inn", nn); pL("nd", nd);}
  assert ((nd>>(64-17)) >= epsilon);
  uL delta = ~nd+1; if(p) pL("delta", delta);
  uL vnum = num>>11; // A shadow value with which to argue and debug!!
{int const jc = 1+(54-1)/b; int j = jc;
 uL q = 0; while(j--) {int td = nn>>(63-b);
     if(p){if(td >= 1<<b) printf("Bang "); printf("%d ", td);}
     nn = ((nn<<b)&((1L<<63)-1)) + td*(delta>>1);
     vnum -= td*(den>>(10+b)); vnum <<= b;
     if(p) pL("vnum", vnum);
     assert(nn == N*(vnum<<1) && N*vnum<(1L<<63));
     // we need a more complex invariant.
     q = (q<<b) + td;}
   // do rounding here.
   // There is the hair to address that we don't know for sure
   // the position of the most significant quotient bit.
   pL("nn", nn); pL("q", q);
   if(0){
     printf("j=%d %d %d jc=%d sh=%d\n", j, (int)(nn>>63),
        (nn<<1)+delta==nn, jc, (64-b*jc));
     if((int)(nn>>62) && (nn<<2)+delta==(nn<<1)) ++q;}
   else {while (nn>=(nd>>1)) {nn-=(nd>>1); printf("+"); ++q;}}
// At this point q holds b*jc bits of quotient if num>den,
// and b*jc-1 otherwise.
   return (q<<(64-b*jc)) & -(1L<<(64-54));}}

void tsd(uL n, uL d){pL("num", n); pL("den", d);
  uL qt = dv(n, d); pL("qt", qt);
  uL qr = divq(n>>1, 0, d).q; pL("qr", qr);
  if(0) pL("-(1L<<(64-54))", -(1L<<(64-54)));
  if(qt!=(qr& -(1L<<(64-54)))) exit(printf("Foo\n"));}

int main (){assert(b<9);
{int j=256; int y = 1<<17;
  while(j--) {int t = (1<<17)/(257+j);
  if(0) if(8<=j && j<=10 || j==15 || j==21 ) ++t;
  {int z = (256+j)*t; if(z<y) y=z;}
   tab[j] = t; if(0) printf("%4d %4d, %5x\n", j, t, t*(256+j));}
  printf("epsilon = %x\n", y); assert (y==epsilon);}
  if(0) pL("quo", dv(5l<<61, 1l<<63));
  if(1) tsd(0x804L<<52, 0x840L<<52); // ++tab[8]
  {int j=54/b; double e=(1<<b)-1;
    while(j--) e *= ((1<<b)-e*(1.-((float)(1<<b)-epsilon); assert(e<2.);}
  tsd(0x800L<<52, 0x848L<<52); // ++tab[9]
  tsd(0x800L<<52, 0x850L<<52); // ++tab[10]
  tsd(0x804L<<52, 0x878L<<52); // ++tab[15]
  tsd(0x804L<<52, 0x8a8L<<52); // ++tab[21]
  tsd(0x804L<<52, 0x878L<<52); // ++tab[15]
  tsd(0x804L<<52, 0x928L<<52); // ++tab[21]
  tsd(5l<<61, 1l<<63); tsd(1l<<63, 5l<<61);
  if(1) {p=0; for(int i=0; i<256; ++i) for(int j=0; j<256; ++j) {
    uL n = (uL)i<<(63-8), d = (uL)j<<(63-8), m = 1L<<(63-8)-1, k=1L<<63;
    tsd(k+n, k+d); tsd(k+n, k+d+m); tsd(k+n+m, k+d);
    tsd(k+n+m, k+d+m);}}
return 0;}