This Scheme program hews to the high level view modeling the values of n and d by Scheme integers by the same names bit multiplied by 253. Rounding logic goes into the expression “(* q (expt 2 c))” which now rounds towards 0.

To compute the fraction n/d to m bits beyond the binary point, we compute the integer quotient (2mn)/d and then move the binary point n places to the left.

The routine (dz n np) takes positive integer arguments and returns the integer └n∙2p/Nd┘.

We need to compute the bounds of the ‘trial divisors’ which are the values of td in the program. I think that the current program assures that td will never be negative. Each value of td depends on the previous value; we want to define an attractor basin. More accurately we want to find a contiguous set S of td’s which f maps into S.
Below we write “d” for trial divisor, “s” for log2 of table size, “z” for bits in table entry, “b” for bits of quotient per cycle.
The errors arise from: If b is too large for some pair of z and s, then the trial divisors will grow as the division proceeds. It would be good to find a formula for a safe b depending on z and s. For b smaller than that limit some table entries will not be used. Perhaps the table size need not be a power to two.
f(d) = (r + ¾
In grammar school long division we subtract multiples of the divisor from the numerator until that ‘remainder’ gets small enough to just call it the remainder and let it go at that. The multiples are by descending powers of 10. The cofactors of the divisor compose the quotient. If we guess the quotient digit wrong we may have to subtract again, or perhaps add and so adjust the quotient digit. It is the same here except that in place of quotient digit we adopt fixed numbers of bits. If we chose 4 bit it would be like dividing hexadecimal numbers. Amdahl’s trick works in decimal by making the numerator begin with leading 9's making guessing the next quotient digit especially easy.

The Iterated Part

We begin with a remainder and take some fixed bits from that remainder as an guess for the next batch of quotient bits. These batches will necessarily overlap and with care overlap only by one bit. Combining these batches into the ultimate quotient is a process of summation. There may be advantage to be taken of the fact that the overlaps are sparse. We examine next the bounds on guessing the quotient bits by this scheme.
There are three significant stages to this division:

This and this C code show signs of life but have become encrusted with debugging logic. I begin a new slate with built-in limit calculations that might be comprehensible. Here we crowd the significant bits to the right in the machine words. We seek 54 bits of quotient for the case that the numerator is greater than the denominator. The remainder is available to us, in internal form, to make a variety of rounding choices. Just now the code assumes that we will get 3 (=b) bits of quotient per iteration. The plan is to make that variable per cycle depending only on factory parameters. (The code uses this.)

Consider both arguments to dv and tsd to be fractions x such that ½ ≤ x < 1 scaled at −63. This means that the value ascribed by C semantics to the long int parameters, times 2−63 is the value that we reason about. The highest bit in each of these arguments is 1. Routine va verifies that values are in this range and that only 53 significant bits are present.

The arguments to rd are scaled to the right and allow for both numerator and denominator to be multiplied by the 9 (=rs) bit approximation to the reciprocal of the denominator. This requires 53+rs significant bits for the new numerator, nn, and the new denominator, nd. nn = num*tab[j] occupies 53+rs bits. nd = den*tab[j] occupies 53+rs bits except the high bit of this value is 0 considering how j was chosen. As integers, nn < 253+rs and 254+rs < nd < 254+rs+1. It is unnecessary to specify the scaling of the new arguments, except to say that they are scaled the same. The quotient q will be such that ½ < q < 2, scaled at −54.

Reading C Code

I read constructs such as ‘num>>(64-53)’ as “Put the left most 53 bits in the rightmost 53 places, inserting 0’s on the left.”, at least if num is of type uL.
1/(1+x) = 1 − x + x2 − x3 …