#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
typedef long unsigned int uL;
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;
void pLL(char* n, res r){printf("%016lX.%016lX = %s:\n", r.q, r.r, n);}
res divq(res, uL);
res mulq(uL, uL);
res addq(res, res);
int min(int a, int b){return a<b?a:b;}
int max(int a, int b){return a<b?b:a;}
char p = 0; char const pz=0;
// llls long long left shift
res llls(res v, int n){return n<64?(res){(v.q<<n)|(n?(v.r>>64-n):0),v.r<<n}
     :(res){v.r<<(n-64), 0};}
int lleq(res a, res b){return a.q==b.q && a.r==b.r;}
res coer(uL a){return (res){0, a};}

int b; //calculated eventually ...
int const rs = 10; // number of bits in reciprocal approximation
int const lts = 9, ts = 1<<lts, bi = 8;

// reciprocal size and log table size.
short tab[ts];
// dv returns 2*(the quotient) using the proposed scheme.
// The quotient may not be normalized.

// the following variables would be neither global nor mutable
// if clang had nested fuctions!!
uL nd, nn, delta, Ng;
int epsilon, big;

uL rd(uL nnp, uL quo, int qb){
if(p) {pL("nnp", nnp); pL("Ng", Ng);}
if(1) assert(nnp%Ng==0);
if(p) {printf("qb=%d ", qb); pL("nnp", nnp); pL("quo", quo);}
   if(!qb) {if(pz){pL("enn", nnp); pL("equo", quo);
     printf("eqb=0 b=%d rs=%d\n", b, rs);}
   assert(nnp<(nd<<1));
   return nnp>=nd?quo+1:quo;}
   b = min(b, qb); // crude and safe estimate
   int splitpoint = 52+rs-b; if(p) printf("splitpoint=%d\n", splitpoint);
   int td = nnp>>splitpoint;
   uL nx = nnp - ((uL)td<<splitpoint); // split nnp.
   if(p) {printf("%d ", td); pL("znx", nx);}
   uL vv = (nx<<b)+td*delta;
   if(0) {pL("vv", vv);
     pL("(nnp<<b)", (nnp<<b)); pL("td*nd", td*nd);
     pL("(nnp<<b)-td*nd", (nnp<<b)-td*nd);}
   assert(vv==(nnp<<b)-td*nd);
   // below: if we were dividing nnp by nd would we be getting closer?
return rd(vv, (quo<<b)+td, qb-b);}

uL dv(uL const num, uL const den){
  {va(num); va(den);}
  big = num>=den; // not used soon; part of result
  if(num==den) return 1L<<52;
  b = bi;
  int N = tab[(den>>(64-(lts+1)))-ts]; Ng=N;
  if(p) printf("N=%d=0x%03x\n", N, N);
  assert(54+(rs-1) <= 64);
  nn = N*(num>>(64-53)); nd = N*(den>>(64-53));
  delta = (1L<<(53+rs-1))-nd;
  if(pz) {pL("nn", nn); pL("nd", nd); pL("delta", delta);}
  assert(nn < (1L<<(53+rs)));
  {uL tq = divq(llls((res){0, num}, 53-big), den).q,
      rq = divq(llls((res){0, nn}, 53-big), nd).q;
      if(pz) {pL("tq", tq); pL("rq", rq);} assert(tq==rq);
      uL xq = rd(nn, 0, 53-big); if(pz) pL("xq", xq);
      assert(xq==tq);
   return xq;}}

uL dvg(uL const num, uL const den){
  {va(num); va(den);}
  if(num<den) return divq(llls(coer(num>>11), 53), den>>11).q;
  if(den<num) return divq(llls(coer(num>>11), 52), den>>11).q;
  return 1L<<52;}
void tsd(uL n, uL d){if(pz) {pL("num", n); pL("den", d);}
  uL qg= dvg(n, d); if(pz) pL("qg", qg);
  uL qt = dv(n, d); if(pz) pL("qt", qt);
  uL qr = (divq((res){n>>1, 0}, d).q>>10)>>(n>=d);
  if(pz) pL("qr", qr);
  assert(qt == qr); assert(qt == qg);}

// This code is to compute the quotient to 53 bits without rounding
// but since position of leftmost 1 bit depends, by 1, onwhether num<den,
// we compute to 54 bits.
// The construct "&~(d<n)" above conditionally discards the extra low bit
// so as to always to deliver exatcly 53 significant bits.

uL r52(){return (random()^((((uL)random())<<21) | (1L<<52)))<<11;}
int main (){
{int const mg = lts+rs; int j=ts; int y = 1<<mg, z=-1<<mg;
  while(j--) {int r = (1<<mg)/(ts+1+j);
    y = min(y, (ts+j)*r); z = max(z, (ts+1+j)*r);
    tab[j] = r; if(0) printf("%4d %4d, %5x\n", j, r, r*(ts+j));}
  printf("epsilon = %d = 0x%x; mx = %d = 0x%x\n", y, y, z, z);}
  tsd(0x807FFFFFFFFFF800L, 0x8000000000000000L);
  if(1) {p=0; for(int i=0; i<ts; ++i) for(int j=0; j<ts; ++j) {
    uL n = (uL)i<<(64-(lts+1)), k=1L<<63,
       d = (uL)j<<(64-(lts+1)), m = (1L<<(64-(lts+1)))-(1<<11);
    if(0) {pL("d", d); pL("m", m); pL("k", k); pL("n", n);}
    tsd(k+n, k+d);   tsd(k+n+m, k+d);
    tsd(k+n, k+d+m); tsd(k+n+m, k+d+m);}
printf("Good patterns\n");}
{int j = 10000000; while(j--) tsd(r52(), r52());}
printf("Good random\n");
return 0;}