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.
Programming Praxis 
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 nThe output of this program is:
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)Just so everybody’s clear, 2^2 – 1^2 = 3
[…] Ancient Algorithms […]
[…] 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]). […]