Big Numbers: Division
June 7, 2011
We continue our examination of big-integer arithmetic today with a look at division. Long division was hard in grade school, and it’s hard for computers, too, with tricky algorithms and lots of special cases. Fortunately, Donald Knuth has made things easier for us, and we will be following his Algorithm 4.3.1 D.
Division takes as input two numbers, n (the dividend, also called the numerator) and d (the divisor, also called the denominator) and returns two numbers q (the quotient) and r (the remainder) such that q · d + r = n, with 0 ≤ r < d. The basic idea of Algorithm D is to take successive partial divisions, where in each case the partial dividend has one more digit than the divisor, each successive partial division revealing one more digit of the quotient. As with the other arithmetic operators, there is a notion of carry from one partial division to the next.
We’re not going to give a detailed explanation here, because Knuth says it far better; you can run to your nearest copy of Knuth, or peek at the solution. Beware that the code is lengthy, and there are lots of tricky bits, and plan from the start to do lots of debugging.
Your task is to write a function that performs division on big integers. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
[…] today’s Programming Praxis exercise, our goal is to add division to our Big Number library. Let’s get […]
My Haskell solution (see http://bonsaicode.wordpress.com/2011/06/07/programming-praxis-big-numbers-division/ for a version with comments):
import Control.Arrow instance Integral BigNum where 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]