Egyptian multiplication and some of its ramifications

TL;DR

This paper explores Egyptian multiplication and exponentiation, highlighting their historical significance and deriving algorithms for division, fractional powers, and logarithms based on these ancient methods.

Contribution

It presents a unified framework for understanding Egyptian multiplication and exponentiation, deriving classical algorithms and proposing new ones for fractional powers.

Findings

Egyptian multiplication can be used for efficient division algorithms.

The method yields a new approach to fractional power computation.

It provides historical context for algorithms like logarithms and division.

Abstract

Multiplication and exponentiation can be defined by equations in which one of the operands is written as the sum of powers of two. When these powers are non-negative integers, the operand is integer; without this restriction it is a fraction. The defining equation can be used in evaluation mode or in solving mode. In the former case we obtain "Egyptian" multiplication, dating from the 17th century BC. In solving mode we obtain an efficient algorithm for division by repeated subtraction, dating from the 20th century AD. In the exponentiation case we also distinguish between evaluation mode and solving mode. Evaluation mode yields a possibly new algorithm for raising to a fractional power. Solving mode yields the algorithm for logarithms invented by Briggs in the 17th century AD.

Equations

$$27\times 23=27+54+108+432$$

$$27\times 23=1\times 27+1\times 54+1\times 108+0\times 216+1\times 432.$$

$$c-r=a\times b=a\times \sum _{i=0}^{\infty }{d}_i{2}^i=\sum _{i=0}^{\infty }{d}_i{2}^{i}a$$

$$626-r=d_{0}27+d_{1}54+d_{2}108+d_{3}216+d_{4}432.$$

$$194-r=d_{0}27+d_{1}54+d_{2}108+d_{3}216$$

$$a^{2p/q}$$

$$a^{p/q}$$

$$x=\sum _{i=1}^{53}{d}_i{2}^{-i},$$

$$x=\sum _{i=1}^{k-1}{d}_i{2}^{-i}+d_k{2}^{-k}+\sum _{i=k+1}^{53}{d}_i{2}^{-i}.$$

$${\log }_{b}z=\sum _{i=1}^{k-1}{d}_i{2}^{-i}\text{mbox}and\text{mbox}frac=2^{-k}.$$

Full text

1 Egyptian multiplication

To begin at the beginning, let the starting point be the multiplication algorithm in the Rhind papyrus, a document found in Egypt. My source is an article by James R. Newman [4]. The document in question was a scroll of 18 feet long written in about 1700 B.C. Though not found intact, enough of the parts have been compiled to obtain a coherent whole. The scroll is a collection mathematical exercises and practical examples. Most of it is of no mathematical interest because of the unwieldy notation for fractions. However, one algorithm, for multiplying two integers, stands out for its clarity, simplicity, and relevance for computer programming.

Newman illustrates the algorithm by multiplying 23 and 27 as shown in Figure 1.

It is not clear from his article whether this example actually occurs in the papyrus or whether it is his way of presenting the idea gleaned from the papyrus. Let us translate the idea to modern mathematical notation.

The idea is that $a\times b$, with non-negative integer $b$, can be written as $a\times\sum_{i=0}^{\infty}d_{i}2^{i}$ where the $d_{i}$ are [math] or $1$. Of course, from a certain $i$ onward, all the $d_{i}$ are zero. We have $a\times b=a\times\sum_{i=0}^{\infty}d_{i}2^{i}=\sum_{i=0}^{\infty}d_{i}2^{i}a$. In the second equality, the multiplication is converted to additions (in the form of doublings).

In Newman’s example $a=27$ and $b=23$, so that we get

[TABLE]

corresponding to the binary expansion of $b$ being $10111$. When we determine this expansion in the usual way, the digits are generated in reverse order, so that we get

[TABLE]

This suggests evaluating $\sum_{i=0}^{\infty}d_{i}2^{i}a$ by combining the generation of the binary digits with the doublings of $a$ in a single loop, as in

r := 0;
while (b>0) {
  if (odd(b)) {r := r+a; b := b-1}
  a := a+a; b := b/2
}

with the invariant $a_{0}\times b_{0}=r+a\times b$ where $a_{0}$ and $b_{0}$ are the initial values of $a$ and $b$.

Although Newman’s presentation suggests that a table of the successive doublings be kept, the above code shows that each result of doubling can overwrite the previous one.

This is an example of using the binary expansion to determine the value of $a\times b$ when both operands are known. Let us call this “evaluation mode”. We can also use the binary expansion in “solving mode”: to determine an unknown $b$ when $a$ and $c$ are given in $a\times b=c$. As this does not in general have a solution in integers, we solve instead $a\times b=c-r$ with $0\leq r<b$. In solving mode we compute the inverse, which in this case is division:

[TABLE]

.

For example

[TABLE]

Here it is apparent that $d_{4}=1$, so that we can expose $d_{3}$ by subtracting $432$ from both sides.

[TABLE]

from which we conclude the $d_{3}=0$ and $d_{2}=1$. Further subtractions give $86-r=d_{0}27+d_{1}54$ and $32-r=d_{0}27$ from which we conclude that $d_{0}=1$ and $r=5$, so that division of $626$ by $27$ gives quotient $11101$ (in binary, with most significant digit first) and leaves a remainder of 5.

The presentation by Newman suggests that a table of doublings of $a$ be kept. This is not necessary, as the next doubling can be written over the current one. Determining the binary expansion of $b$ tells one when to stop.

Most of the Rhind papyrus is concerned with division. No similarly simple and clean algorithm is found by Newman. This is probably also the case for the study [1] by Arnold Chace, on which Newman’s article is based.

The above numerical example for integer division shows that one needs to know in advance the largest doubling needed. Once this is known, the lesser doublings are efficiently obtained by halving; again no table is needed. See Figure 2.

This completes all that is said in this paper about multiplication and division. From now on the topic is exponentiation and its inverse. Just as multiplication can be done by repeated addition, exponentiation can be done by repeated multiplication. It can be speeded up in a similar way by using the binary expansion of one of the operands. We will be concerned with the case of a fractional operand, so that the binary expansion extends into negative powers of two. Again, algorithms can be obtained by using the defining equation in the two modes used earlier: evaluating mode and solving mode. We will need to take square roots because what halving is to doubling, taking the square root is to squaring. There is an algorithm for square root that is of interest in its own right.

2 The miracle of the square root

To obtain the square root $x$ of $a$, that is, to find $x$ such that $x^{2}=a$, observe that this equality implies $x=a/x$, hence $x=(x+a/x)/2$. This suggests considering the sequence $x_{0},x_{1},\ldots$ defined by $x_{n+1}=(x_{n}+a/x_{n})/2$ with some arbitrarily chosen $x_{0}$, say, 1. It is clear that if $x_{n}<\surd a$, then $a/x_{n}>\surd a$, and vice versa. If $x_{n}$ is regarded as a guess at the value of $\surd a$, then $x_{n+1}=(x_{n}+a/x_{n})/2$ is a plausible way of getting a better guess. It is not only plausible, but is guaranteed to converge to $\surd a$. Moreover, the number of correct figures doubles at every iteration. When I speak of the miracle of square root, I have in mind this combination of plausibility and algorithmic effectiveness.

The method is an instance of Newton’s method of finding roots of non-linear equations, and is often presented as such. This is a pity because it suggests that one needs calculus to understand a good method for the square root. This is not the case, as the algorithm is described by Heron of Alexandria who lived two thousand years ago111 There is “informed conjecture” that this algorithm was used by the Babylonians, and thus may be roughly as old as Egyptian multiplication. .

3 Fractional powers

The Heron algorithm might suggest looking for an equally elegant and effective algorithm for cube roots, fifth roots, $\ldots$. In the context of this paper a more attractive option is to exploit square roots for raising $a$ to any power of the form $p/q$, for non-negative integer $p$ and any positive integer $q$.

The mathematical basis of the algorithm is

[TABLE]

These equalities suggest mutual adjustments among $a$, $p$, and $q$ in such a way that after a sufficient number of steps the desired power is trivial to obtain. See Figure 4.

Fractional powers are defined by $a^{p/q}=b$. To compute such a power, $a$, $p$, and $q$ are given while $b$ is unknown. This defining equation is thus used in evaluation mode. Ideally, the same defining equation can be used in solving mode, with $a$ and $b$ given. But an unknown in the form $p/q$ is awkward. Therefore we use $a^{t}=b$ as defining equation, where $t$ is understood to be a fraction.

In Figure 5 we list a C function to compute $a^{t}$ even though this function is merely a predictable variant of the one listed in Figure 4.

4 In the footsteps of Henry Briggs

The Rhind papyrus uses a table of doublings of the multiplier. A modern version of this algorithm creates the table implicitly, each entry overwriting the previous one. In this section we will see that a similar table is useful in computing logarithms. This use occurred early in the history of computation; barely a quarter of a century after Simon Stevin taught the world to use decimal positional notation for fractions. I quote from The Feynman Lectures on Physics ([3], page 22-6):

This is how logarithms were originally computed by Mr Briggs of Halifax, in 1620. He said “I computed successively 54 square roots of 10.”

Briggs might easily have picked another number of square roots. On my computer (and on yours probably as well) the iteration

  double b = sqrt(10.0); while (b>1.0) b = sqrt(b);

stops after producing 53 distinct square roots. I don’t know how Briggs decided on 54, but this occurrence of “53” has an explanation. The double type in C is represented by the double-length format of the IEEE floating-point standard, which has a mantissa of 52 bits. This might suggest a relative precision of 52 bits, were it not for the fact that the first bit of the significand is suppressed in the format because it is always 1. Because that is taken into account, one gets effectively a precision of 53 bits. A curious near-coincidence with the choice of Briggs in 1620. For future reference we note that 53 bits is equivalent to about 16 decimal places.

We let the Feynman lectures ([3], page 22-7) continue with

$\ldots$ he calculated sixteen decimal places, and then reduced his answer to fourteen places when he published it, so that there were no rounding errors. He made tables of logarithms to fourteen decimal places by this method, which is quite tedious. But all logarithm tables for three hundred years were borrowed from Mr Briggs’s tables by reducing the number of decimal places. Only in modern times, with the WPA222 Works Progress Agency, a US government agency active in the 1930s to carry out public works, thus alleviating unemployment.

and computing machines, have new tables been independently computed.

The role of tables so far has been conceptual: no data structure was created to be filled with results of squarings or doublings, of square roots or halvings. Instead, the contents of the notional tables were constructed on the fly and overwritten promptly after use.

Rather than pursue an algorithm for base-10 logarithm, we note that Briggs’s method of computing a table of iterated square roots works equally well for any base $b$ (within a suitable range) as it does for base 10.

5 Logarithms to any base

So far, when we considered fractional powers, as in $b^{x}=a$ with fractional $x$, we only used this equality in evaluation mode. In solving mode, $b$ and $a$ are given and one is to find $x$. In other words, one is to find the logarithm of $a$ to base $b$. We determine the $x$ in $b^{x}=a$ when $a$ and $b$ are given, with $b>1$ and $1\leq a\leq b$. From these assumptions we conclude that $0\leq x\leq 1$. Accordingly, we have that

[TABLE]

with $d_{i}\in\{0,1\}$ and uniquely determined in the precision anticipated for $x$.

A plausible way to attain the goal of $b^{x}=a$ is to suppose we have a $k$ such that

[TABLE]

Initially this is easy to make true with $k=1$. To increase $k$ by one we maintain program variables $z$ and frac such that

[TABLE]

To discover whether $d_{k}=1$ or $d_{k}=0$ we test $z*b<a$. Truth indicates the former, falsity the latter. With $b=1$, $x$ cannot increase any more, and we know that its true value (to the anticipated precision) is representable as a sum of powers of 2, and we have included all of the powers that should be included.

Such considerations lead to the function listed in Figure 6.

Acknowledgements

Thanks to Paul McJones for corrections and discussions.