## Big Numbers: Functions

### June 28, 2011

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^{32}, 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.