Ancient Algorithms

December 23, 2014

[ Today’s exercise is a Christmas stocking stuffed with six little exercises. I wish all my readers a Merry Christmas and a Happy New Year! I’ll be taking a few days with family, so the next exercise will appear on January 6, 2015. ]

We tend to think of algorithms as procedures used by computers, but algorithms actually have a long and distinguished history. Today’s exercise discusses six algorithms that pre-date the Christian era, arranged roughly in chronological order. Despite their antiquity, all these algorithms are still in use today.

Peasant Multiplication: This algorithm has been independently invented by most societies as they progress to numeracy, so it is called the Egyptian method, the Ethiopian method, the Russian method, the Ukranian method, and so on for many different peoples; we simply call it the Peasant method. It is known that the Egyptians used this algorithm during the construction of the pyramids.

The algorithm starts by writing the two numbers to be multiplied at the heads of two columns. Then, on successive rows, the number in the left column is halved, discarding any remainder, and the number in the right column is doubled, until the left column is reduced to one. Then the sum of the numbers in the right column that correspond to odd numbers in the left column is the desired product. The ancients use pairs of pebbles to halve numbers and determine parity; we find it simpler to sum the odds as we go along.

function product(left, right)
    prod := 0
    while left > 0
        if left is odd
            prod := prod + right
        left := halve(left)
        right := double(right)
    return prod

In the United States we torment eight-year-old kids with the times tables, but school children in the rest of the world learn this method, which is simpler and, if you’ve ever taught third-grade math, more likely to be done correctly. Additionally, all computers do arithmetic in binary, and use this algorithm deep in the microcode that is stamped into the silicon chip that runs them, shifting bits to perform the halving and doubling operations. The only improvement is that computers nowadays process multiple bits at the same time, instead of just a single bit, reducing the number of steps in the loop.

Babylonian Square Roots: Mathematical historians know that the Babylonians invented a method of calculating square roots because a clay tablet with the calculations etched into it (ancient “scrap paper”) exists, but the method was first described by Hero of Alexandria in the first century after Christ. Hero (his name is sometimes spelled Heron) invented the jet engine (a hollow metal ball filled with water, heated over a fire to convert the water to steam, which exited by two nozzles on opposite sides of the globe; the rotational motion could be captured by gears or belts), the windmill (he used it to power a variety of machines, including a grain grinder), and the vending machine (you dropped a coin into a slot, which caused a balance beam to rise, opening a valve and filling a small vessel with holy water; when the vessel was full, it tipped to deliver the holy water, the coin was tipped into the safe, the balance beam fell, the valve closed, and the machine reset for the next purchase).

Hero observed that if x is an approximation of the square root of n, a better approximation is given by the average of x and n /x; we say x ′ = (x + n /x) / 2. The initial value of x can be anything on the range 1 ≤ xn, though the process converges more quickly if x is somewhat closer to n. In the version of the function given below, we bound n on the range 1 ≤ n < 4 so we can fix the number of iterations required for double-precision accuracy; an alternative is to iterate until the difference between two successive approximations is less than some desired epsilon.

function sqrt(n)
    if n < 1 return sqrt(n * 4) / 2
    if 4 <= n return sqrt(n / 4) * 2
    x = (n + 1) / 2
    repeat 5 times
        x = (x + n / x) / 2
    return x

The only improvement since the Babylonians is the initialization of x, which nowadays is done by some kind of black magic operating on the internal representation of the number (for integers, it’s generally one less than the number of bits in the binary representation of the number, for floating point it’s generally half the mantissa of the number with the exponent reduced by one), usually reducing the number of iterative steps to two or three.

Pythagorean Triples: Pythagoras, who lived about five centures before Christ, was a mathematician, philosopher and mystic. When he discovered that the square root of 2 could not be expressed as the ratio of two integers (we now say that it is irrational), he thought he had discovered some impossible flaw in the perfection of the World, and swore his followers to secrecy; legend tells us that one of his disciples, Hippasus, was murdered when he divulged the secret (another legend has it that Hippasus, not Pythagoras, who discovered the irrationality of the square root of 2).

Pythagoras proved that for any right triangle the square of the length of the hypotenuse is equal to the sum of the squares of the other two sides. He also discovered that given any coprime m and n, a primitive right triangle (all sides coprime to each other) has sides m2n2, 2mn, m2 + n2; for instance, with m = 2 and n = 1, the triangle is 22 − 12 = 2, 2 × 2 × 1 = 4, and 22 + 12 = 5.

function triple(m, n)
    return m*m - n*n, 2*m*n, m*m + n*n

An alternative method is to take any primitive Pythagorean triple {a, b, c} and generate three new primitive Pythagorean triples as {a-2b+2c, 2ab+2c, 2a-2b+3c}, {a+2b+2c, 2a+b+2c, 2a+2b+3c}, {−a+2b+2c, −2a+2b+3c}; the process starts from {3, 4, 5}. Of course, non-primitive triples can be generated by multiplying all three elements of a primitive triple by some constant k.

Greatest Common Divisor: Euclid was either a Greek mathematician of the third century before Christ or, as some historians suggest, the conglomerate name of several mathematicians. Euclid wrote Elements, a textbook of geometry and number theory which has been in continuous publication since it was written, which is famous for its method of starting with five axioms and deriving all of geometry from them.

Euclid’s Elements, Book VII, Proposition 2, gives an algorithm for computing the greatest common divisor of two integers. Euclid’s algorithm works by repeatedly subtracting the smaller number from the larger, which reduces the magnitude of the numbers, thus guaranteeing termination of the algorithm, without affecting the greatest common divisor, since both the smaller number and the difference between the two numbers must be multiples of the greatest common divisor.

function gcd(m, n) # ancient
    if m < n return gcd(m, n-m)
    if n < m return gcd(m-n, n)
    return m
function gcd(m, n) # modern
    if n == 0 return m
    return gcd(n, m % n)

Euclid used repeated subtraction; nowadays we use modular arithmetic. Knuth calls this the “grand-daddy” of all algorithms, even though it isn’t the oldest, because it was the first to be written in rigorous form.

Approximation of Pi: Archimedes of Syracuse (on the island of Sicily) was a famous inventor of the third century before Christ: the screw pump, block-and-tackle pulley, odometer (it dropped a ball every time the wheel completed a turn), and numerous devices of war were among his inventions. He didn’t invent the lever, but was the first to explain how it worked, and he famously stated “Give me a lever and a place to stand, and I can move the world.” His most famous discovery was that objects float when they weigh less than the water they displace, which he used to prove that the king’s crown, thought to be pure gold, was actually an amalgam with silver; it is said that he discovered the principle while bathing, whereupon he ran naked through the street shouting “Eureka!”

The ratio of the circumference of a circle to its diameter is the mathematical constant known as π, which is approximately 3.141592654. Archimedes bounded the value of π in the range 223/71 < π < 22/7 by calculating the perimeters of two hexagons, one inscribed inside a circle and one circumscribing it, then successively doubling the number of sides of the inscribed and circumscribed regular polygons through the series 6, 12, 24, 48, and 96.

function pi(n) # archimedes used 6
    outer := 3 * sqrt(3)
    inner := outer / 2
    for i from 1 to n
        outer = 2 / (1/outer + 1/inner)
        inner = sqrt(outer * inner)
    return inner, outer

Archimedes’ approximation 22/7 is still frequently used today, but his method of computing the approximation has been replaced by a series expansion due to Srinivasa Ramanujan, or some similar series expansion, which converges in just a few steps; by contrast, Archimedes’ algorithm requires 27 steps for double-precision accuracy.

Prime Numbers: Eratosthenes, who measured the circumference of Earth, the distance from Earth to Sun, and the tilt of Earth’s axis, devised a system of latitude and longitude and was third chief librarian of the great library at Alexandria, succeeding Ptolemy and Apollonius, invented a method of enumerating prime numbers about 200BC; his writing has been lost, but we know of his invention through the writings of Nicomachus two centuries later.

The Sieve of Eratosthenes identifies primes by striking out multiples and keeping what’s left; the image is that composites fall through the sieve but primes remain. Initialize a list of numbers, from 2 to the desired limit n all marked as prime, then consider each number in turn; if it has not been marked as composite, report it as prime, then mark all its multiples as composite, starting from the square of the prime (since all smaller multiples will have already been marked composite by smaller primes).

function primes(n)
    sieve := makeArray(2..n, True)
    for p from 2 to n
        if sieve[p]
            output p
            for i from p*p to n step p
                sieve[i] = False

Modern mathematicians have optimized the basic Sieve of Eratosthenes shown above by eliminating the small primes that, because they have so many multiples, take most of the time; eliminating multiples of 2 is easy, and “wheels” can be used to eliminate multiples of 3, 5, 7, and so on, though increasingly large wheels quickly become cumbersome to operate. Still, an optimized Sieve of Eratosthenes is faster than any of its many modern alternatives.

Your task it to translate these seven programs into your favorite programming language. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

5 Responses to “Ancient Algorithms”

  1. Globules said
    import Control.Applicative ((<$>))
    import Control.Monad (when)
    import Control.Monad.ST.Safe
    import Data.List (intercalate, unfoldr)
    import Data.Ratio
    import qualified Data.Vector.Generic as V
    import qualified Data.Vector.Unboxed as UV
    import qualified Data.Vector.Unboxed.Mutable as MV
    import Text.Printf (printf)
    
    -- The product of x and y.  Assumes x and y are non-negative.
    peasant :: Integral a => a -> a -> a
    peasant x y = sum $ unfoldr step (x, y)
      where step (r, s) | r == 0 = Nothing
                        | otherwise = let next = (r `quot` 2, s*2)
                                      in Just (if even r then 0 else s, next)
    
    -- An approximation to the square root of x, with tolerance eps.
    hero :: (Fractional a, Ord a) => a -> a -> a
    hero eps x = let rs = iterate (\y -> (y + x/y)/2) x
                 in snd . head . dropWhile tooFar $ zip rs (tail rs)
      where tooFar (r, s) = abs (r - s) >= eps
    
    -- A Pythagorean triple corresponding to m and n.  Assumes m and n are positive
    -- and coprime, with m > n.
    pyth :: Integral a => a -> a -> (a, a, a)
    pyth m n = let (mm, nn) = (m*m, n*n) in (mm - nn, 2*m*n, mm + nn)
    
    -- An infinite list of primitive Pythagorean triples, starting with the seed
    -- (3, 4, 5).  (The list is the "Tree of primitive Pythagorean triples" in depth
    -- first order.  See Wikipedia.)
    pyths :: [(Integer, Integer, Integer)]
    pyths = concat $ iterate (concatMap step) [(3, 4, 5)]
      where step :: (Integer, Integer, Integer) -> [(Integer, Integer, Integer)]
            step (a, b, c) =  [( a-2*b+2*c,  2*a-b+2*c,  2*a-2*b+3*c),
                               ( a+2*b+2*c,  2*a+b+2*c,  2*a+2*b+3*c),
                               (-a+2*b+2*c, -2*a+b+2*c, -2*a+2*b+3*c)]
    
    -- Greatest common divisor.  (Ancient version.)
    gcdA :: Integral a => a -> a -> a
    gcdA m n | m < n     = gcd  m    (n-m)
             | n < m     = gcd (m-n)  n
             | otherwise = m
    
    -- Greatest common divisor.  (Modern version.)
    gcdM :: Integral a => a -> a -> a
    gcdM m n | n == 0    = m
             | otherwise = gcd n (m `rem` n)
    
    -- Greatest common divisor.  (Modern version without explicit recursion.)
    gcdM' :: Integral a => a -> a -> a
    gcdM' m n = fst . head . dropWhile ((/=0) . snd) $ iterate next (m, n)
      where next (a, b) = (b, a `rem` b)
    
    -- Archimedes approximation to π using n-gons.
    pi' :: Floating a => Int -> a
    pi' n = (\(o, i) -> (o + i)/2) . head . drop n $ iterate step oi
      where step (o, i) = let oi'@(o', _) = (2/(1/o + 1/i), sqrt (o'*i))
                          in oi'
            oi@(o, _) = (3*sqrt 3, o/2)
    
    -- A vector of primes less then or equal to n.  For fun we follow the imperative
    -- Python version more closely than usual, using a mutable vector of booleans
    -- (in the ST monad) to flag the primes.
    primes :: Int -> UV.Vector Int
    primes n = runST $ do
      pr <- MV.replicate (n + 1) True
      mapM_ (markMultiples pr) [2..MV.length pr - 1]
      (V.map fst . V.drop 2 . V.filter snd . V.indexed) <$> V.unsafeFreeze pr
        where markMultiples v m = do
                b <- MV.read v m
                when b $ mapM_ setFalse [mm, mm+m..MV.length v - 1]
                  where mm = m*m
                        setFalse i = MV.write v i False
    
    main :: IO ()
    main = do
      let (a, b) = (7, 11 :: Int)
      printf "Peasant multiplication: %d * %d = %d\n" a b $ peasant a b
      let x = 2 :: Int
          y = hero (1 % 100000) 2 :: Rational
      printf "Hero's square root using rationals: sqrt %d = %s\n" x $ show y
      let x = 2 :: Double
          y = hero (1/100000) x
      printf "Hero's square root using doubles: sqrt %f = %f\n" x y
      let (a, b) = (2, 1 :: Int)
      printf "Pythagorean triple for coprimes %d and %d = %s\n"
        a b $ show $ pyth a b
      printf "Pythagorean triples generated from (3, 4, 5):\n  %s\n"
        $ intercalate "\n  " $ take 10 $ map show pyths
      let (a, b) = (12345, 75 :: Int)
      printf "GCD of %d and %d: Euclid = %d, modern = %d, modern' = %d\n"
        a b (gcdA a b) (gcdM a b) (gcdM' a b)
      let n = 6
      printf "Archimedes approximation to π using %d-gons = %f\n"
        n (pi' n :: Double)
      let n = 80
      printf "Primes ≤ %d = %s\n" n $ show $ V.toList $ primes n
    

    The output of this program is:

    Peasant multiplication: 7 * 11 = 77
    Hero's square root using rationals: sqrt 2 = 665857 % 470832
    Hero's square root using doubles: sqrt 2.0 = 1.4142135623746899
    Pythagorean triple for coprimes 2 and 1 = (3,4,5)
    Pythagorean triples generated from (3, 4, 5):
      (3,4,5)
      (5,12,13)
      (21,20,29)
      (15,8,17)
      (7,24,25)
      (55,48,73)
      (45,28,53)
      (39,80,89)
      (119,120,169)
      (77,36,85)
    GCD of 12345 and 75: Euclid = 15, modern = 15, modern' = 15
    Archimedes approximation to π using 6-gons = 3.1416627611326433
    Primes ≤ 80 = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79]
    
  2. Globules said

    Just noticed that gcdA and gcdM are mistakenly calling the system gcd function. They should be:

    -- Greatest common divisor.  (Ancient version.)
    gcdA :: Integral a => a -> a -> a
    gcdA m n | m < n     = gcdA  m    (n-m)
             | n < m     = gcdA (m-n)  n
             | otherwise = m
    
    -- Greatest common divisor.  (Modern version.)
    gcdM :: Integral a => a -> a -> a
    gcdM m n | n == 0    = m
             | otherwise = gcdM n (m `rem` n)
    
  3. Josh said

    Just so everybody’s clear, 2^2 – 1^2 = 3

  4. […] Ancient Algorithms […]

  5. […] One writer saw Johnny Ball’s video about Russian Multiplication on Numberphile and suggested it would make a good exercise. Indeed it would; in fact, we have already done it twice ([1] [2]). […]

Leave a comment