Below mathematical (not software) conventions are in effect!
The following expressions are exact reals and for positive y, x/y denotes infinite precision division ((x/y)y = x) and x÷y is the integer quotient of dividing x by y and discarding any remainder.
Theorems: y>0 → (0 ≤ x/y − x÷y < 1 and x−1 < (x÷y)y ≤ x).
“b” below is a design parameter, fixed at the factory and k depends on b and is likewise fixed at factory.
“ε” is fixed at the factory and is derived from design details not discussed on this page.
Increasing the hardware can make ε smaller.
“n” and “δ” are runtime arguments to a functional unit that performs 64 bit IEEE floating point division.
Conditional branch to same ideas as below but trying to clarify them:
abstract
Given at runtime
- n: 0 ≤ n < 253
- δ: 0 < δ < ε ≤ ¼.
We seek the integral part of n/(1−δ), i.e. n÷(1−δ).
We assume radix 2b for some positive integer b.
Let k = 54÷b.
0 ≤ n < 2(k+1)b.
We take t = n÷2kb as the trial divisor.
0 ≤ t < 2b
and 0 ≤ n − t2kb < 2kb.
Our ‘next’ remainder is n' = n − t(1−δ)2kb
= n − t2kb + tδ2kb.
We contribute t2kb to the quotient.
We need bounds on n'.
The least n' is when n − t2kb = 0
and δ = 0 whereupon n' = 0.
Since δ cannot reach 0, n' can be 0 only if n=0.
The largest n' is when n − t2kb = 2kb and δ = ε.
Then n' = n − t(1−δ)2kb
= n − t(1 − ε)2kb
= n − t2kb + tε2kb
= 2kb + tε2kb
= (1+tε)2kb.
t maxes out at 2b and so we have:
0 < n' < 2kb + tε2kb
= (1+ε2b)2kb.
The benefit of small ε is now clear.
Our upper bound on n' depends on t < 2b which we derived from n < 2(k+1)b.
The next super digit will be calculated similarly but with n' in the role of n and k−1 in the role of k.
Limiting t is critical to hardware design for it impacts the computation of t(1−δ)2kb.
I think that the bounds on successive super digits tj =
Bk = 2b
Bj = (1+ε2b)k−j.
t = tk and 0 ≤ tj < Bj.
See this.
It is time to corroborate these inequalities with code.
In routine div:
Parameter n there is 29n here.
Parameter d there is 263δ here.
Global ep there is 263ε here.