Programming Praxis

5 Responses to “Ancient Algorithms”

1. Globules said

   December 28, 2014 at 9:01 AM

   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

   December 29, 2014 at 4:49 AM

   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

   December 31, 2014 at 7:36 PM

   Just so everybody’s clear, 2^2 – 1^2 = 3
4. Algorithms | Erika M said

   January 6, 2015 at 3:13 PM

   […] Ancient Algorithms […]
5. My Mailbag | Programming Praxis said

   February 7, 2020 at 10:01 AM

   […] 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]). […]