Posts Tagged ‘number’

Programming Praxis – Two Random Exercises

August 17, 2012

In today’s Programming Praxis exercise, our goal is to use convert two random number generators to a different range. Specifically, we have to convert an rng that produces numbers in the closed interval [1,3] to one that has the interval [1,9] and one with interval [1,5] to [1,7]. In both cases the resulting rngs must maintain an equal distribution across the interval. Let’s get started, shall we?

Some imports:

import Control.Applicative
import Control.Monad
import Data.List
import System.Random

These are the two functions we can use as imports.

rand3, rand5 :: IO Int
rand3 = randomRIO (1,3)
rand5 = randomRIO (1,5)

To test, we use a distribution function to see how many of each element occur in a sequence.

dist :: (Eq a, Ord a) => [a] -> [(a, Int)]
dist = map (\x -> (head x, length x)) . group . sort

Let’s start with the easiest one. Since 9 = 3 * 3, all we need to do scale up one input by a factor of 3 and “interpolate” using a second one. The reason for taking a source parameter is to facilitate testing: you can either feed it an IO Int for a random result, or you can give it a [Int] to see all possible results. We’ll get to why the testing is important in a bit.

rand9 :: (Applicative f, Num a, Enum a) => f a -> f a
rand9 src = (+) . (*3) . pred <$> src <*> src

The second one is more complicated, since 5*5 is not a multiple of 7. The correct solution is to do the same thing as above using a subrange that is a multiple of 7, i.e. 21 of the 25 numbers and rerolling when one of the four values outside the accepted range is chosen. The problem with this approach is that it isn’t guaranteed to terminate, since there is nothing preventing a true random number generator from generating the number 23 a few trillion times in a row. Granted, we’re only using pseudo-random number generators, so it probably will never come up. Here’s the implementation:

rand7 :: (Applicative m, Monad m, Integral a) => m a -> m a
rand7 src = (+) . (*5) . pred <$> src <*> src >>= \n -> 
    if n > 21 then rand7 src else return (mod n 7 + 1)

Since I didn’t like the chance of a non-terminating algorithm I came up with the following algorithm: sum seven random numbers in the range 1 to 5, and give the result modulo 7. Or in code:

rand7approx :: (Functor f, Monad f, Integral a) => f a -> f a
rand7approx src = succ . (`mod` 7) . sum <$> replicateM 7 src

Printing the distribution looked perfectly fine, with all of them getting a little over 14%. But something felt off to me. The code was shorter than the version with rerolling and didn’t suffer the inifinite loop problem, so there had to be reason it wasn’t the provided solution. I changed all the functions to take a source argument rather than calling rand3 and rand5 directly and printed the theoretical distributions as well (except for rand7, since the distribution for that one is non-terminating). This revealed that there is a 0.3% difference between the most and least common values. In retrospect, this should’ve been obvious, as 5^7 mod 7 isn’t equal to 0, which means it’s impossible for all 7 values to get the same distribution. Some more reading up revealed that there is indeed no solution that both provides an equal distribution and is guaranteed to terminate, since 7 is an infinitely repeating fraction in base 5.

Here are the tests I used:

main :: IO ()
main = do print $ dist (rand9 [1..3]) == zip [1..9] [1,1..]
          print . dist =<< replicateM 100000 (rand9 rand3)
          print . dist =<< replicateM 100000 (rand7 rand5)
          print $ dist (rand7approx [1..5])
          print . dist =<< replicateM 100000 (rand7approx rand5)

Programming Praxis – Big Numbers: Division

June 7, 2011

In today’s Programming Praxis exercise, our goal is to add division to our Big Number library. Let’s get started, shall we?

Because I found two bugs in the existing code and to avoid spreading the code around too much, this week I’ll post the full code for the module again.

I’ve added some imports:

import Control.Arrow
import Data.List
import Data.Ratio
import Test.QuickCheck

The data type, base and Ord instance are unchanged.

data BigNum = B Int [Int] deriving (Eq, Show)

base :: Integer
base = 1000

instance Ord BigNum where
    compare (B l1 ds1) (B l2 ds2) = case compare l1 l2 of
        EQ -> maybe EQ id . find (/= EQ) . reverse $ zipWith compare ds1 ds2
        c  -> c

While testing division I discovered that the result of a subtraction operation wasn’t trimmed of leading zeroes. The new definition for + solves this.

instance Num BigNum where
    a@(B l1 ds1) + b@(B l2 ds2) = B (length t * signum l) (reverse t) where
        B l _ = if abs b > abs a then b else a
        (_,t) = span (== 0) . reverse . f 0 $ (if abs b > abs a then flip else id)
             (prep $ if signum l1 == -signum l2 then (-) else (+)) ds1 ds2
        prep op (x:xs) (y:ys) = op (toInteger x) (toInteger y) : prep op xs ys
        prep _  xs     ys     = map toInteger $ xs ++ ys
        f r (x:xs) = let (d,m) = divMod (r + x) base in fromIntegral m : f d xs
        f r []     = if r == 0 then [] else [fromIntegral r]
    (B l1 ds1) * (B l2 ds2) = B (signum l1 * signum l2 * sl) sds where
        B sl sds = sum $ mult ds1 ds2
        mult (x:xs) (y:ys) = fromIntegral (toInteger x * toInteger y) :
                             map shift (mult xs (y:ys)) ++
                             map shift (mult [x] ys)
        mult _     _  = []
        shift (B l ds) = B (l + 1) (0 : ds)
    negate (B l ds) = B (-l) ds
    abs (B l ds)    = B (abs l) ds
    signum (B l _)  = fromIntegral $ signum l
    fromInteger n | n < 0     = negate $ fromInteger (-n)
                  | otherwise = B (length ds) (map fromIntegral ds)
                  where ds = tail $ f (n,0)
                        f (0,m) = [m]
                        f (d,m) = m : f (divMod d base)

These two instances aren’t strictly necessary, but making a data type an instance of Integral requires it instances Enum and Real as well. The definitions are trivial enough, so I added them.

instance Enum BigNum where
    fromEnum = fromIntegral . toInteger
    toEnum   = fromInteger . fromIntegral

instance Real BigNum where
    toRational = (% 1) . toInteger

The second bug I found was that converting a Big Number to an Integer didn’t keep the sign of the number. This is now fixed. Also, since we’re now instancing Integral I could give it the proper name.

instance Integral BigNum where
    toInteger (B l ds) = foldr (\x a -> fromIntegral (signum l * x) + base * a) 0 ds

For the actual division, I developed my own algorithm since I couldn’t find a good explanation of Knuth’s algorithm. It fairly simple, though probably a bit less efficient than Knuth’s version. Determine how often the denominator can go in the numerator based on the first digit group of the denominator (d1) and the first digit group of the numerator (n1). If n1 is less than d1, the second digit group of the numerator (n2) is used instead, and a value equal to how often d1 goes into the base is added. The resulting value is multiplied by the denominator, subtracted from the numerator and the algorithm is called recursively. Thanks to the Integral typeclass we get the other functions (div, mod, etc.) for free.

    quotRem n@(B l1 ds1) d@(B l2 ds2)
        | d == 0         = error "Division by zero"
        | n < 0 || d < 0 = let (q',r') = quotRem (abs n) (abs d)
                           in (signum n * signum d * q', signum n * r')
        | n < d          = (0, n)
        | otherwise      = first (+ q) $ quotRem (n - d*q) d where
        (n1,n2,d1) = (last ds1, last $ tail ds1, last ds2)
        q = if n1 < d1 then shift (l1 - l2 - 1) $ div n2 d1 +
                            fromIntegral (div base (fromIntegral d1))
                       else shift (l1 - l2) $ div n1 d1
        shift s i = B (s + 1) $ (replicate s 0) ++ [i]

Since there are plenty of possible corner cases I automated the testing using QuickCheck to make sure there is no difference between Integer and BigNum division.

divTest :: Integer -> Integer -> Property
divTest a b = b /= 0 ==> quotRem a b == (toInteger q, toInteger r)
    where (q, r) = quotRem (fromInteger a :: BigNum) (fromInteger b)

main :: IO ()
main = do print $ div 12345678 (3456 :: BigNum) == 3572
          print $ mod 12345678 (3456 :: BigNum) == 846
          quickCheck divTest

Programming Praxis – ISBN Validation

May 20, 2011

In today’s Programming Praxis exercise, our goal is to write a number of functions related to ISBN numbers. Let’s get started, shall we?

Some imports:

import Control.Applicative hiding ((<|>), optional)
import Data.Char
import Data.List
import Data.Map (elems)
import Network.HTTP
import Text.HJson
import Text.HJson.Query
import Text.Parsec

First, we need some parsers for ISBN and EAN numbers.

isbn = (++) <$> (concat <$> sepEndBy1 (many1 d) (oneOf " -"))
            <*> option [] ([10] <$ char 'X') where
    d = read . return <$> digit
ean = string "978" *> optional (oneOf " -") *> isbn

Since we need the check digits both for validation and conversion we make separate functions for them.

isbnCheck, eanCheck :: Integral a => [a] -> a
isbnCheck n = 11 - mod (sum $ zipWith (*) [10,9..] (take 9 n)) 11
eanCheck n = mod (sum $ zipWith (*) (cycle [1,3]) (take 9 n)) 10

Checking whether a number is valid is a matter of checking if the length and last digit are correct.

validISBN, validEAN :: String -> Bool
validISBN = valid isbn isbnCheck
validEAN = valid ean eanCheck

valid p c = either (const False) v . parse p "" where
    v ds = length ds == 10 && c ds == last ds

Conversion just requires changing the last digit.

toISBN, toEAN :: String -> Maybe String
toISBN = convert ean isbnCheck
toEAN = fmap ("978-" ++) . convert isbn eanCheck

convert p c = either (const Nothing) (Just . fixCheck) . parse p ""
    where fixCheck n = map intToDigit (init n) ++ [check $ c n]
          check n = if n == 10 then 'X' else intToDigit n

Since I don’t like APIs that require an access key, we’ll be using openlibrary instead of isbndb.

lookupISBN :: String -> IO [(String, [String])]
lookupISBN = get . ("http://openlibrary.org/api/books?format=json&\
                    \jscmd=data&bibkeys=ISBN:" ++) where
    f ~(JObject j) = map (\b -> (unjs $ key "title" b,
        map (unjs . key "name") . getFromArr $ key "authors" b)) $ elems j
    key k = head . getFromKey k
    unjs ~(JString s) = s
    get url = fmap (either (const undefined) f . fromString) .
              getResponseBody =<< simpleHTTP (getRequest url)

Some tests to see if everything is working properly:

main :: IO ()
main = do print $ validISBN "99921-58-10-7"
          print $ validISBN "80-902734-1-6"
          print $ validISBN "0-943396-04-2"
          print $ validISBN "0-9752298-0-X"
          print $ validISBN "0943396042"
          print $ not $ validISBN "99921-58-10-8"
          print $ not $ validISBN "99921-58-10-"
          print $ not $ validISBN "9"
          print $ validEAN "978-0-0700048-4-9"
          print $ validEAN "9780070004849"
          print $ not $ validEAN "9780070004848"
          print $ toISBN "9780070004849" == Just "0070004846"
          print $ toEAN "0070004846" == Just "978-0070004849"

          mapM_ (\(t,a) -> putStrLn ("Title: " ++ t) >>
                           putStrLn ("Authors: " ++ intercalate ", " a)) =<<
              lookupISBN "0070004846"

Programming Praxis – Rule 30 RNG

April 29, 2011

In today’s Programming Praxis exercise, our goal is to write a random number generator based on the Rule 30 cellular automaton. Let’s get started, shall we?

Some imports:

import Control.Monad.State
import Data.Bits
import Data.Word

First we define the amount of bits in our RNG state.

size :: Int
size = 7 * 43

We’re going to be needing a function to convert a list of bits to a number.

fromBits :: (Bits a, Integral a) => [Bool] -> a
fromBits = foldl (\a x -> shift a 1 .|. fromIntegral (fromEnum x)) 0

The step function replaces the state with the next generation of cells for a given rule and returns the middle bit.

step :: Monad m => Int -> StateT Integer m Bool
step r = modify (\s -> let b i = testBit s (mod i size) in
    fromBits [ testBit r $ fromBits [b (n+1), b n, b (n - 1)]
             | n <- [size - 1, size - 2..0]]) >> gets (`testBit` div size 2)

To produce a random number, we take the middle bit of the desired number of generations and convert that to an unsigned integer.

randomBits :: Monad m => Int -> Int -> StateT Integer m Word
randomBits n r = return . fromBits =<< replicateM n (step r)

To test the algorithm, we see if we get the correct result when we only set the bit in the middle. Afterwards, we produce a bunch of 8-bit numbers.

main :: IO ()
main = evalStateT (do r <- randomBits 32 30
                      liftIO $ print (r == 3112904540)
                      liftIO . print =<< replicateM 100 (randomBits 8 30)
                  ) $ bit (div size 2)

Looks like everything is working properly.

Programming Praxis – Two Kaprekar Exercises

March 22, 2011

In today’s Programming Praxis exercise,our goal is to determine the longest possible Kaprekar chain and the Kaprekar numbers below 1000. Let’s get started, shall we?

Some imports:

import Data.List
import qualified Data.List.Key as K
import Text.Printf

The Kaprekar chain algorithm is fairly simple. We use printf for easy padding of the number.

chain :: Int -> [Int]
chain 0 = []
chain 6174 = [6174]
chain n = n : chain (read (reverse p) - read p) where
    p = sort $ printf "%04d" n

Determining whether a number is a Kaprekar number is made easier by realizing that the whole ‘take the first n-1 or n and the last n digits’ bit can be replaced with the divMod function.

isKaprekar :: Integral a => a -> Bool
isKaprekar n = n == uncurry (+) (divMod (n^2) $ 10 ^ length (show n))

Some tests to see if everything is working properly:

main :: IO ()
main = do print $ chain 2011 == [2011, 1998, 8082, 8532, 6174]
          print . K.maximum length $ map chain [0..9999]
          print $ filter isKaprekar [1..999] ==
                  [1, 9, 45, 55, 99, 297, 703, 999]

The chain shown is just one of 2184 possible chains of length 8. The first chain of length 8 starts with 14.