Big Numbers: Division

June 7, 2011

(define (big-divide ns ds) ; (values quotient remainder)
  (define (sign x) (if (negative? x) -1 (if (positive? x) 1 0)))
  (define (lt? xs ys)
    (let ((xlen (length xs)) (ylen (length ys)))
      (if (< xlen ylen) #t (if (< ylen xlen) #f
        (let loop ((xs (reverse xs)) (ys (reverse ys)))
          (cond ((null? xs) #f) ; equal
                ((< (car xs) (car ys)) #t)
                ((< (car ys) (car xs)) #f)
                (else (loop (cdr xs) (cdr ys)))))))))
  (define (reduce xs)
    (if (null? (cdr xs)) xs
      (if (positive? (car xs)) xs
        (reduce (cdr xs)))))
  (define (add b1 b2)
    (let loop ((b1 b1) (b2 b2) (c 0) (bs '()))
      (cond ((null? b1)
              (if (zero? c) (reverse bs) (reverse (cons 1 bs))))
            ((null? b2)
              (let* ((sum (+ (car b1) c))
                     (new-c (if (<= big-base sum) 1 0)))
                (loop (cdr b1) b2 new-c
                      (cons (modulo sum big-base) bs))))
            (else (let* ((sum (+ (car b1) (car b2) c))
                         (new-c (if (<= big-base sum) 1 0)))
                    (loop (cdr b1) (cdr b2) new-c
                          (cons (modulo sum big-base) bs)))))))
  (define (sub b1 b2)
    (let loop ((b1 b1) (b2 b2) (c 0) (bs '()))
      (cond ((null? b1) (reverse (reduce bs)))
            ((null? b2)
              (let* ((diff (- (car b1) c))
                     (new-c (if (< diff 0) 1 0)))
                (loop (cdr b1) b2 new-c
                      (cons (modulo diff big-base) bs))))
            (else (let* ((diff (- (car b1) (car b2) c))
                         (new-c (if (< diff 0) 1 0)))
                    (loop (cdr b1) (cdr b2) new-c
                          (cons (modulo diff big-base) bs)))))))
  (define (times big1 big2)
    (let loop ((b1 big1) (b2 big2) (zs '())
               (c 0) (ps '()) (bs '()))
      (cond ((null? b1) bs)
            ((null? b2) (let ((zs (cons 0 zs)))
              (loop (cdr b1) big2 zs 0 zs
                (add (reverse (if (zero? c) ps (cons c ps))) bs))))
            (else (let* ((x (+ (* (car b1) (car b2)) c))
                         (c (quotient x big-base))
                         (p (modulo x big-base)))
                    (loop b1 (cdr b2) zs c (cons p ps) bs))))))
  (define (mul1 ns d)
    (let loop ((ns ns) (c 0) (ps '()))
      (if (null? ns) (reverse (if (zero? c) ps (cons c ps)))
        (let* ((x (+ (* (car ns) d) c))
               (c (quotient x big-base))
               (p (modulo x big-base)))
          (loop (cdr ns) c (cons p ps))))))
  (define (div1 ns d)
    (let loop ((rev-ns (reverse ns)) (qs '()) (r 0))
      (if (null? rev-ns) (values (reverse (reduce (reverse qs))) (list r))
        (let* ((x (+ (* r big-base) (car rev-ns)))
               (q (quotient x d)) (r (modulo x d)))
          (loop (cdr rev-ns) (cons q qs) r)))))
  (define (nextq x d0 d1 ns)
    (let loop ((q (quotient x d0)))
      (if (< (* q d1) (+ (* (- x (* q d0)) big-base) (caddr ns))) q
        (loop (- q 1)))))
  (define (nextn j n rev-ns ds*q)
    (let ((zs (append (reverse (sub (reverse (take (+ n 1) rev-ns)) ds*q))
                      (drop (+ n 1) rev-ns))))
      (if (< (length zs) (+ n j)) (cons 0 zs) zs)))
  (define (div ns ds)
    (if (lt? ns ds) (values '(0) ns)
      (let* ((n (length ds)) (m (- (length ns) n)))
        (if (= n 1) (div1 ns (car ds))
          (let* ((d (quotient big-base (+ (car (reverse ds)) 1)))
                 (rev-ns (reverse (mul1 ns d)))
                 (rev-ns (if (= (length ns) (length rev-ns))
                           (cons 0 rev-ns) rev-ns))
                 (ds (mul1 ds d))
                 (d0 (car (reverse ds)))
                 (d1 (cadr (reverse ds))))
            (let loop ((j m) (rev-ns rev-ns) (qs '()))
              (if (negative? j)
                  (call-with-values
                    (lambda () (div1 (reverse rev-ns) d))
                    (lambda (q r) (values (reverse (reduce (reverse qs))) q)))
                  (let* ((x (+ (* (car rev-ns) big-base) (cadr rev-ns)))
                         (q (nextq x d0 d1 rev-ns))
                         (ds*q (mul1 ds q)))
                    (if (lt? (reverse (take (+ n 1) rev-ns)) ds*q)
                        (let* ((q (- q 1)) (ds*q (sub ds*q ds)))
                          (loop (- j 1) (nextn j n rev-ns ds*q) (cons q qs)))
                        (loop (- j 1) (nextn j n rev-ns ds*q) (cons q qs)))))))))))
  (if (big-zero? ds) (error 'big-divide "divide by zero")
    (let ((sn (sign (car ns))) (sd (sign (car ds))))
      (call-with-values
        (lambda () (div (cdr (big-abs ns)) (cdr (big-abs ds))))
        (lambda (qs rs)
          (values (cons (* (if (= sn sd) 1 -1) (length qs)) qs)
                  (cons (* sn (length rs)) rs)))))))

About these ads

Pages: 1 2 3

2 Responses to “Big Numbers: Division”

  1. […] today’s Programming Praxis exercise, our goal is to add division to our Big Number library. Let’s get […]

  2. 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]
    

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 637 other followers

%d bloggers like this: