#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
typedef long unsigned int uL;
int const b = 2; // 1 2 3 4 5 6 7 8
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)));}

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))]; printf("N=%d=0x%03x\n", N, N);
  uL nn = N*(num>>10), nd = N*(den>>8);
  // neither of above overflow!
  pL("Inn", nn); pL("nd", nd); assert ((nd>>(64-17)) >= epsilon);
  uL delta = ~nd+1; pL("delta", delta);
  uL vnum = num/2; // 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(td >= 1<<b) printf("Bang "); printf("%d ", td);
    // {res r = mulq(nd, q); ...}
     nn = ((nn<<b)&((1L<<63)-1)) + td*(delta>>1);
     pL("nn", nn);
     vnum -= td*(den>>1); vnum <<= b;
     if(1) pL("vnum", vnum);
     if(1){pL("nn", nn); pL("N*(vnum>>9)", N*(vnum>>9));}
     assert(nn == N*(vnum>>9));
     // 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.
   if(p) {pL("nnf", nn); pL("q", q); pL("vnumf", vnum);}
   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;
// 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); {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(0) tsd(0x800L<<52, 0x804L<<52);
  if(1) 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);
    uL m = (1L<<(63-8)) - (1L<<(64-53)), 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;}