Programming Praxis

The greatest common divisor is calculated by Euclid’s algorithm:

(define (big-gcd m n)
  (let loop ((m (big-abs m)) (n (big-abs n)))
    (cond ((big-lt? m n) (loop m (big-minus n m)))
          ((big-lt? n m) (loop (big-minus m n) n))
          (else m))))

Raising to a power is done by the square-and-multiply algorithm; so is modular exponentiation:

(define (big-expt b e)
  (if (big-zero? e) '(1 1)
    (let loop ((s b) (i e) (a '(1 1)))
      (let ((a (if (big-odd? i) (big-times a s) a))
            (i (big-quotient i '(1 2))))
        (if (big-zero? i) a (loop (big-times s s) i a))))))

(define (big-expm b e m)
  (define (m* x y) (big-modulo (big-times x y) m))
  (cond ((big-zero? e) '(1 1))
        ((big-even? e) (big-expm (m* b b) (big-quotient e '(1 2)) m))
        (else (m* b (big-expm (m* b b) (big-quotient (big-minus e '(1 1)) '(1 2)) m)))))

Square roots are calculated by Newton’s algorithm based on derivatives; the initial value is approximated by taking half the logarithm in the car of the big number:

(define (big-sqrt n)
  (if (not (big-positive? n))
      (error 'big-sqrt "must be positive")
      (let loop ((x (cons (quotient (+ (car n) 1) 2) (cdr n))))
        (let ((y (big-quotient (big-plus x (big-quotient n x)) '(1 2))))
          (if (big-lt? y x) (loop y) x)))))

Random numbers are generated by the Park/Miller algorithm. Since the seed must be between 0 and 2³², there are many big numbers that can never be returned by the random number generator.

(define big-rand ; Park/Miller algorithm
  (let ((a (string->big "16807"))
        (m (string->big "2147483647"))
        (s (string->big "1043618065")))
    (lambda args
      (if (eq? (car args) 'seed)
          (begin (set! s (cadr args)) s)
          (let ((first (if (pair? (cdr args)) (car args) '(0)))
                (past (if (pair? (cdr args)) (cadr args) (car args))))
            (set! s (big-modulo (big-times a s) m))
            (big-plus first (big-quotient (big-times s (big-minus past first)) m)))))))

Conversions between native big numbers and our big numbers are done a digit at a time; note the awkward method for determining the internal base of our big numbers:

(define integer->big
  (let ((base (string->number (big->string '(2 0 1)))))
    (lambda (n)
      (if (zero? n) '(0)
        (if (negative? n) (big-negate (integer->big (- n)))
          (let loop ((n n) (bs '()))
            (if (zero? n) (cons (length bs) (reverse bs))
              (loop (quotient n base) (cons (modulo n base) bs)))))))))

(define big->integer
  (let ((base (string->number (big->string '(2 0 1)))))
    (lambda (b)
      (if (big-zero? b) 0
        (if (big-negative? b) (* -1 (big->integer (big-negate b)))
          (let loop ((bs (reverse (cdr b))) (n 0))
            (if (null? bs) n
              (loop (cdr bs) (+ (car bs) (* n base))))))))))

The complete big number library is shown on the next page, and at http://programmingpraxis.codepad.org/hgy1QKi2.

Pages: 1 2 3