TAOCP 1.2.2 Exercise 25
Let the value of $x$ at the beginning of an execution of step L3 be denoted by $x^{(t)}$, where $t = 0, 1, 2, \ldots$ counts the number of times step L4 has been performed.
Section 1.2.2: Numbers, Powers, and Logarithms
Exercise 25. [22] Suppose that we have a binary computer and a number $x$, $1 \le x < 2$. Show that the following algorithm, which uses only shifting, addition, and subtraction operations proportional to the number of places of accuracy desired, may be used to calculate an approximation to $y = \log_b x$:
L1. [Initialize.] Set $y \leftarrow 0$, $z \leftarrow x$ shifted right 1, $k \leftarrow 1$.
L2. [Test for end.] If $x = 1$, stop.
L3. [Compare.] If $x - z < 1$, set $z \leftarrow z$ shifted right 1, $k \leftarrow k + 1$, and repeat this step.
L4. [Reduce values.] Set $x \leftarrow x - z$, $z \leftarrow x$ shifted right $k$, $y \leftarrow y + \log_b (2^k/(2^k - 1))$, and go to L2. ▎
[Notes: This method is very similar to the method used for division in computer hardware. The idea goes back in essence to Henry Briggs, who used it (in decimal rather than binary form) to compute logarithm tables, published in 1624. We need an auxiliary table of the constants $\log_b 2$, $\log_b (4/3)$, $\log_b (8/7)$, etc., to as many values as the precision of the computer. The algorithm involves intentional computational errors, as numbers are shifted to the right, so that eventually $x$ will be reduced to 1 and the algorithm will terminate. The purpose of this exercise is to explain why it will terminate and why it computes an approximation to $\log_b x$.]
Verified: no
Solve time: 27s
Solution
Let the value of $x$ at the beginning of an execution of step L3 be denoted by $x^{(t)}$, where $t = 0, 1, 2, \ldots$ counts the number of times step L4 has been performed. Let the corresponding values of $y$ and $k$ be denoted by $y^{(t)}$ and $k_t$.
Since $1 \le x < 2$ initially, we have
$$ 0 \le x - (x \text{ shifted right } 1) < 1. $$
Hence the first execution of step L3 must eventually terminate, because repeated right shifts of $z$ produce values
$$ \frac{x}{2}, \frac{x}{4}, \frac{x}{8}, \ldots, $$
and therefore
$$ x - \frac{x}{2^k} < 1 $$
whenever
$$ \frac{x}{2^k} > x - 1. $$
Since $x - 1 < 1$, such a $k$ always exists.
Suppose step L3 terminates with value $k$. Since the previous value of $z$ failed the test, we have
$$ x - \frac{x}{2^{k-1}} \ge 1, $$
while the accepted value satisfies
$$ x - \frac{x}{2^k} < 1. $$
Therefore
$$ 1 \le x\left(1 - \frac1{2^k}\right) < 1, $$
after the subtraction in step L4. Hence the new value of $x$ is
$$ x' = x - \frac{x}{2^k} = x\left(1 - \frac1{2^k}\right) = x\frac{2^k-1}{2^k}, $$
and it satisfies
$$ 1 \le x' < 1 + \frac1{2^k}. $$
Now consider the quantity
$$ x^{(t)} b^{y^{(t)}}. $$
Initially,
$$ x^{(0)} b^{y^{(0)}} = x. $$
During one execution of step L4, the value of $x$ is multiplied by
$$ \frac{2^k-1}{2^k}, $$
while the value of $y$ is increased by
$$ \log_b \frac{2^k}{2^k-1}. $$
Hence
$$ b^{y'} = b^y \cdot \frac{2^k}{2^k-1}, $$
by equation (9). Therefore
$$ x' b^{y'} = x\frac{2^k-1}{2^k} \cdot b^y \frac{2^k}{2^k-1} = x b^y. $$
Thus the product $x^{(t)} b^{y^{(t)}}$ is invariant throughout the algorithm, and consequently
$$ x^{(t)} b^{y^{(t)}} = x $$
for all $t$.
When the algorithm terminates, $x = 1$. If the final value of $y$ is denoted by $Y$, the invariant gives
$$ 1 \cdot b^Y = x, $$
hence
$$ Y = \log_b x, $$
by equation (9).
This proves that the exact algorithm computes $\log_b x$.
The actual implementation uses shifted values truncated to the machine precision. Consequently the subtraction in step L4 is not exact, and eventually the computed value of $x$ becomes exactly $1$ after sufficiently many reductions. The algorithm therefore terminates after a finite number of steps.
To estimate the amount of work, suppose the machine retains $n$ binary places. Each right shift decreases the magnitude of $z$ by a factor of $2$, so after at most $n$ shifts the quantity $z$ becomes smaller than one unit in the last place. Therefore each execution of step L3 requires at most $n$ comparisons and shifts. Each execution of step L4 decreases $x - 1$, because
$$ x' - 1 = x\left(1-\frac1{2^k}\right)-1 < x-1. $$
Since the computation uses only $n$ binary places, at most $n$ significant reductions can occur before $x$ reaches the value $1$ exactly. Hence the total number of shifts, additions, and subtractions is proportional to $n$.
Therefore the algorithm computes an approximation to $\log_b x$ using only shifting, addition, and subtraction operations, and the total number of such operations is proportional to the number of binary places of accuracy desired.
This completes the proof.
∎