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