Posts Tagged ‘numbers’

Programming Praxis – Big Numbers: Getting Started

May 24, 2011

In today’s Programming Praxis exercise, our task is to implement the basics of a library for big numbers. Let’s get started, shall we?

A quick import:

import Data.List

We could represent a big number as a plain list like the Scheme version does. Using a custom data structure, however, has the advantage of being able to make it an instance of the standard numeric classes, which means shorter function names and easier literals. We also store the number of digits separately because I find it a bit cleaner.

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

We’re using 1000 as our base for now.

base :: Integer
base = 1000

The Num class gives us some of the required functions. Since addition and multiplication haven’t been implemented yet the compiler will throw some warnings.

instance Num BigNum where
    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)

I personally don’t see much use for three separate functions for the sign of a number, since you can either use signum or the appropriate comparsion, but we’ll stick to the assignment.

positive, negative, zero :: (Num a, Ord a) => a -> Bool
positive = (> 0)
negative = (< 0)
zero     = (== 0)

If we make BigNum an instance of the Integral class we could use the default even and odd functions, but this requires implementing modulo arithmetic. We’ll use these for now until we have to do so in a following exercise.

bigOdd, bigEven :: BigNum -> Bool
bigEven (B l ds) = l == 0 || even (head ds)
bigOdd           = not . bigEven

The Integral class also has a toInteger function, but for now we’ll use a separate function.

fromBig :: BigNum -> Integer
fromBig (B _ ds) = foldr (\x a -> fromIntegral x + base * a) 0 ds

Here’s another advantage of using the standard type classes: we only need to implement the compare function to get all the others for free.

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

Time for a whole bunch of tests.

main :: IO ()
main = do let a = 12345678
          let b = -87654321
          print $ 0               == B 0 []
          print $ 1               == B 1 [1]
          print $ (-1)            == B (-1) [1]
          print $ a               == B 3 [678,345,12]
          print $ b               == B (-3) [321,654,87]
          print $ fromBig a       == 12345678
          print $ fromBig (abs b) == 87654321
          print $ abs a           == a
          print $ positive a
          print $ negative b
          print $ zero (0 :: BigNum)
          print $ bigEven a
          print $ bigOdd b
          print $ b /= a
          print $ b < a

Looks like everything is working properly.

Programming Praxis – Partition Numbers

April 15, 2011

In today’s Programming Praxis exercise, our goal is to write a function to calculate partition numbers and to determine the 1000th partition number. Let’s get started, shall we?

Since the recursion in the algorithm is not as simple as in for example the Fibonacci numbers a simple zip will not suffice. To store previous results, we use an IntMap.

import Data.IntMap ((!), fromList, insert, findWithDefault)

Since we’re guaranteed to need every partition number below the one we’re looking for, we simply calculate all of them, starting from 1. This is done by folding over the numbers 1 through x with the IntMap as the accumulator.

partition :: Int -> Integer
partition x = if x < 0 then 0 else foldl p (fromList [(0,1)]) [1..x] ! x where
    p s n = insert n (sum [(-1)^(k+1) * (r pred + r succ) | k <- [1..n],
        let r f = findWithDefault 0 (n - div (k * f (3 * k)) 2) s]) s

Let’s see if we did everything correctly.

main :: IO ()
main = print $ partition 1000 == 24061467864032622473692149727991

Yup. After compiling it runs in about 0.4 seconds on my machine.

Programming Praxis – Sums Of Powers

February 11, 2011

In today’s Programming Praxis exercise, our goal is to implement an algorithm that calculates the bernoulli numbers and one that uses them to quickly calculate the sum of the mth power of numbers 1 through n. Let’s get started, shall we?

A quick import:

import Data.Ratio

To calculate the Bernouilli numbers I initially used the naive version, which simply uses the given mathematical formula. This is quick enough for the test case of 1000 numbers, but too slow for the test case that has a million, so we have to do some memoization. A closer look at the formula reveals that any row in the table depends only on the previous row. Since for the end result we are only interested in the last row, we can use iterate to produce the rows of the table. The value of a given column depends only on the number directly above it and the one to the upper right, so we can use a simple zip to calculate the new row.

a :: (Integral a, Integral b) => a -> a -> Ratio b
a i j = iterate (\xs -> zipWith (*) [1..] $ zipWith (-) xs (tail xs))
                (map (1 %) [1..]) !! fromIntegral i !! fromIntegral j

With this function calculating the Bernoulli numbers is trivial.

bernoullis :: (Integral a, Integral b) => a -> [Ratio b]
bernoullis upto = map (flip a 0) [0..upto]

For the algorithm we also need to calculate binomial coefficients, i.e. the amount of different ways you can choose k objects from a group of size n.

choose :: Integral a => a -> a -> Ratio a
choose n k = product [1..n] % (product [1..n-k] * product [1..k])

And some more executable math for the function that calculates the sum of powers.

s :: Integral a => a -> a -> Ratio a
s m n = 1 % (m+1) * sum [choose (m+1) k * a k 0 * (n%1)^(m+1-k) | k <- [0..m]]

We have one test case to test if the algorithm works correctly and one to judge the speed.

main :: IO ()
main = do print $ bernoullis 6 == [1, 1%2, 1%6, 0, -1%30, 0, 1%42]
          print $ s 10 1000 == 91409924241424243424241924242500
          print $ s 100 1000000

The program runs in about 150-170 ms, so we get the same speed as the Scheme version. Good enough for me.

Programming Praxis – Rational Numbers

January 25, 2011

In today’s Programming Praxis exercise, our goal is to implement common operations on rational numbers. Let’s get started, shall we?

We could of course just use the Data.Ratio module, which does everything we need, but that kind of defeats the purpose of the exercise.

Creating a rational number isn’t too difficult, but if you condense the various cases into a single expression like I did here you need to take care to get the signs right. Fortunately, that’s what unit testing is for. I’m using plain tuples here for the sake of brevity. In production it would be advisable to use a data type instead to prevent passing in invalid rationals such as (1,0).

ratio :: Integral a => a -> a -> (a, a)
ratio _ 0 = error "Division by zero"
ratio n d = (signum d * div n g, abs $ div d g) where g = gcd n d

For adding we use the given formula. Subtracting is the same as adding the negative of the second number.

plus :: Integral a => (a, a) -> (a, a) -> (a, a)
plus (nx, dx) (ny, dy) = ratio (nx * dy + dx * ny) (dx * dy)

minus :: Integral a => (a, a) -> (a, a) -> (a, a)
minus x (ny, dy) = plus x (-ny, dy)

Multiplication didn’t pass the unit tests at first. Turns out the formula in the Scheme solution is wrong, so I replaced it with the correct one. Division is just multiplying by the inverse of the second number.

times :: Integral a => (a, a) -> (a, a) -> (a, a)
times (nx, dx) (ny, dy) = ratio (nx * ny) (dx * dy)

divide :: Integral a => (a, a) -> (a, a) -> (a, a)
divide x (ny, dy) = times x (dy, ny)

Finally, there’s the comparison operator.

lessThan :: Integral a => (a, a) -> (a, a) -> Bool
lessThan (nx, dx) (ny, dy) = nx * dy < dx * ny

I used a decent number of unit tests to cover all of the potentially problematic cases.

main :: IO ()
main = do print $ ratio 1 2           == (1,2)
          print $ ratio 1 (-2)        == (-1,2)
          print $ ratio (-1) 2        == (-1,2)
          print $ ratio (-1) (-2)     == (1,2)
          print $ ratio 2 4           == (1,2)
          print $ plus (1,2) (-1,6)   == (1,3)
          print $ minus (1,2) (1,6)   == (1,3)
          print $ times (3,5) (5,3)   == (1,1)
          print $ times (2,5) (3,7)   == (6,35)
          print $ times (2,5) (-3,7)  == (-6,35)
          print $ divide (3,4) (3,2)  == (1,2)
          print $ divide (1,3) (-2,3) == (-1,2)
          print $ lessThan (1,3) (1,2)
          print $ lessThan (-1,2) (1,6)
          print $ lessThan (-1,2) (-1,6)

Everything passes, so it looks like things are working correctly.

Programming Praxis – Polite Numbers

December 17, 2010

In today’s Programming Praxis exercise, our task is to list all the ways that a number can be written as the sum of a consecutive series of integers. Let’s get started, shall we?

Some imports:

import Data.List
import Data.Numbers.Primes

We need a function to calculate the divisors of a number, which we recycle from a previous exercise:

divisors :: Integral a => a -> [a]
divisors = nub . sort . map product . subsequences . primeFactors

We use a somewhat more compact version of the algorithm than the Scheme solution. The main savings are replacing the if statement with simply taking the maximum of the two numbers and not including an explicit check for powers of two, since the basic algorithm already produces the correct answer and powers of 2 are rare enough that in my opinion the extra algorithm size is not warranted.

polites :: Integral a => a -> [[a]]
polites n = [ [max (d//2 - n//d + 1) (n//d - d//2)..n//d + d//2]
            | d <- tail $ divisors n, odd d, let (//) = div]

Thanks to laziness we don’t have to duplicate code when calculating the politeness of a number (the length of the result set). The result sets themselves are never evaluated, so all this does is count the odd divisors.

politeness :: Integral a => a -> Int
politeness = length . polites

Some tests to see if everything is working properly:

main :: IO ()
main = do print $ polites 15 == [[4..6],[1..5],[7,8]]
          print $ politeness 15 == 3
          print $ polites 28 == [[1..7]]
          print $ politeness 28 == 1
          print $ polites 33 == [[10..12],[3..8],[16,17]]
          print $ politeness 33 == 3
          print . all (null . polites) . take 10 $ iterate (* 2) 1

Looks like it is.

Programming Praxis – Happy Numbers

July 23, 2010

Today’s Programming Praxis problem is a pretty simple one, but it comes with a time limit. We have 15 minutes to write a function that finds all the happy numbers up to a given limit. The version below took me around 4 minutes. Let’s go into the explanation.

While we could write the function directly, we’ll split it up in two: one to determine if a single number is happy and one to get all the required happy numbers. We model the happy number algorithm with basic recursion. We keep a list of all previously “visited” numbers to detect loops.

isHappy :: (Read a, Integral a) => a -> Bool
isHappy = f [] where
    f _  1 = True
    f xs n = notElem n xs && f (n:xs) (sum . map (^ 2) $ digits n)
    digits = map (read . return) . show

Once we know if a single number is happy, determining all of them is just a simple matter of filtering all the numbers less than the limit.

happyUpto :: (Read a, Integral a) => a -> [a]
happyUpto n = filter isHappy [1..n - 1]

All that’s left is actually running the algorithm.

main :: IO ()
main = print $ happyUpto 50

Looks like everything’s working correctly. Of course, it would be somewhat embarrassing if it didn’t, since this is only a fraction more complicated than FizzBuzz.