// Report running max and min of r.
#include <stdio.h>
#include <stdlib.h>
#define C const
typedef double R;
typedef unsigned long int ul;
typedef union{R r; ul i;} U;
static R rf(){return 0x1.p-31*random() + 0x1.p-62*random();}
static R l2R(ul x){return ((U)x).r;}
static ul R2l(R x){return ((U)x).i;}
static ul j(short h, ul a){return (ul)h<<52 | (a & ((1L<<52)-1));}
static ul rmin = 1L<<63, rmax = 0;
R fdiv(C R n, C R d){C ul a = R2l(n), b = R2l(d);
  C ul nf = j(1, a), df = j(1, b);
  C char ad = nf < df;
  C short N = (1<<20)/((df>>42)+1); // table lookup! M=8, z=2
  C ul Nd = N*df;
  ul r = (N*nf)<<ad, q = 0; short rb = 53;
  while(rb>=0) {C short nq = r>>54; r = (r<<8) - Nd*nq; q = (q<<8)+nq;
    if(rb ==  5){rmin = rmin<r?rmin:r; rmax = rmax<r?r:rmax;}
    rb -= 8;}
  while(r >= Nd){r -= Nd; ++q;}
  return l2R(j(1023 + (a>>52) - (b>>52) - ad, ((q>>-rb)+1)>>1));}
void t(R a, R b) {if (a/b != fdiv(a, b)) exit(printf("%20.17f %20.17f\n", a, b));}
int main(){for(int r=0; r<100000000; ++r) t(rf(), rf()); 
   printf("rb = -3, rmin = %016lx, rmax = %016lx\n", rmin, rmax); return 0;}
// rb = 53, rmin = 000b1d8367b0284d, rmax = 8c6c753785fbe5d4
// rb = 45, rmin = 0006aa311e7656f8, rmax = 91cd1936faef255c
// rb = 37, rmin = 0000ddeb1f590c16, rmax = 949b0cfa094ca740
// rb = 29, rmin = 0002e74199527000, rmax = 940f392f735de5ac
// rb = 21, rmin = 0006cbe16e9921cb, rmax = 961e3108f2731ad0
// rb = 13, rmin = 0002cf20ad1f9d3a, rmax = 96e6cf65610d7710
// rb =  5, rmin = 0003cddfabaecda0, rmax = 9425aac3c0bddf58