Posts Tagged ‘precision’

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