## Ancient Algorithms

### December 23, 2014

```(define (product left right)
(let loop ((left left) (right right) (prod 0))
(if (zero? left) prod
(loop (quotient left 2) (* right 2)
(if (odd? left) (+ prod right) prod)))))
```
```(define (hero n)
(cond ((< n 1) (/ (hero (* n 4)) 2))
((<= 4 n) (* (hero (/ n 4)) 2))
(else (let* ((x (/ (+ n 1.0) 2))
(x (/ (+ x (/ n x)) 2))
(x (/ (+ x (/ n x)) 2))
(x (/ (+ x (/ n x)) 2))
(x (/ (+ x (/ n x)) 2))
(x (/ (+ x (/ n x)) 2)))
x))))
```
```(define (triple m n)
(values (- (* m m) (* n n))
(* 2 m n)
(+ (* m m) (* n n))))
```
```(define (ancient m n)
(cond ((< m n) (ancient m (- n m)))
((< n m) (ancient (- m n) n))
(else m)))
```
```(define (modern m n)
(if (zero? n) m
(modern n (modulo m n))))
```
```(define (pi n)
(let ((outer (* 3 (sqrt 3))))
(let loop ((n n) (outer outer) (inner (/ outer 2)))
(if (= n 1) (values inner outer)
(let ((outer (/ 2 (+ (/ outer) (/ inner)))))
(loop (- n 1) outer (sqrt (* outer inner))))))))
```
```(define (primes n)
(let ((sieve (make-vector (+ n 1) #t)))
(do ((p 2 (+ p 1))) ((< n p) (newline))
(when (vector-ref sieve p)
(display p) (display " ")
(do ((i (* p p) (+ i p))) ((< n i))
(vector-set! sieve i #f))))))
```

You can run the program at http://programmingpraxis.codepad.org/uMpTzHIU.

Pages: 1 2

### 5 Responses to “Ancient Algorithms”

1. Globules said
```import Control.Applicative ((<\$>))
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
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 ( ). […]