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

4 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 […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: