337 lines
9.9 KiB
Scheme
337 lines
9.9 KiB
Scheme
|
; Copyright (c) 1993, 1994 Richard Kelsey and Jonathan Rees. See file COPYING.
|
|||
|
|
|
|||
|
|
; This is file bignum.scm.
|
|||
|
|
|
|||
|
|
; Integer arithmetic
|
|||
|
|
|
|||
|
|
(define-extended-number-type :bignum (:exact-integer)
|
|||
|
|
(make-bignum sign magnitude)
|
|||
|
|
bignum?
|
|||
|
|
(sign bignum-sign)
|
|||
|
|
(magnitude bignum-magnitude))
|
|||
|
|
|
|||
|
|
(define (integer->bignum m)
|
|||
|
|
(if (bignum? m)
|
|||
|
|
m
|
|||
|
|
(cond ((>= m 0)
|
|||
|
|
(make-bignum 1 (integer->magnitude m)))
|
|||
|
|
((= m least-non-bignum)
|
|||
|
|
(make-bignum -1 least-non-bignum-magnitude))
|
|||
|
|
(else
|
|||
|
|
(make-bignum -1 (integer->magnitude (- 0 m)))))))
|
|||
|
|
|
|||
|
|
;(define (bignum->integer n) ;For debugging
|
|||
|
|
; (* (bignum-sign n)
|
|||
|
|
; (reduce (lambda (d n) (+ d (* n radix)))
|
|||
|
|
; 0
|
|||
|
|
; (bignum-magnitude n))))
|
|||
|
|
|
|||
|
|
(define (make-integer sign mag)
|
|||
|
|
(if (> sign 0)
|
|||
|
|
(if (smaller-magnitude? greatest-non-bignum-magnitude mag)
|
|||
|
|
(make-bignum sign mag)
|
|||
|
|
(magnitude->integer mag))
|
|||
|
|
(if (smaller-magnitude? least-non-bignum-magnitude mag)
|
|||
|
|
(make-bignum sign mag)
|
|||
|
|
(if (same-magnitude? mag least-non-bignum-magnitude)
|
|||
|
|
least-non-bignum
|
|||
|
|
(- 0 (magnitude->integer mag))))))
|
|||
|
|
|
|||
|
|
; Arithmetic
|
|||
|
|
|
|||
|
|
(define (integer+ m n)
|
|||
|
|
(let ((m (integer->bignum m))
|
|||
|
|
(n (integer->bignum n)))
|
|||
|
|
(let ((m-sign (bignum-sign m))
|
|||
|
|
(m-mag (bignum-magnitude m))
|
|||
|
|
(n-sign (bignum-sign n))
|
|||
|
|
(n-mag (bignum-magnitude n)))
|
|||
|
|
(if (= m-sign n-sign)
|
|||
|
|
(make-integer m-sign (add-magnitudes m-mag n-mag))
|
|||
|
|
(if (smaller-magnitude? m-mag n-mag)
|
|||
|
|
(make-integer (- 0 m-sign) (subtract-magnitudes n-mag m-mag))
|
|||
|
|
(make-integer m-sign (subtract-magnitudes m-mag n-mag)))))))
|
|||
|
|
|
|||
|
|
(define (integer- m n)
|
|||
|
|
(integer+ m (integer-negate n)))
|
|||
|
|
|
|||
|
|
(define (integer-negate m)
|
|||
|
|
(cond ((bignum? m)
|
|||
|
|
(make-integer (- 0 (bignum-sign m))
|
|||
|
|
(bignum-magnitude m)))
|
|||
|
|
((= m least-non-bignum)
|
|||
|
|
(make-bignum 1 least-non-bignum-magnitude))
|
|||
|
|
(else (- 0 m))))
|
|||
|
|
|
|||
|
|
(define (integer* m n)
|
|||
|
|
(let ((m (integer->bignum m))
|
|||
|
|
(n (integer->bignum n)))
|
|||
|
|
(make-integer (* (bignum-sign m) (bignum-sign n))
|
|||
|
|
(multiply-magnitudes
|
|||
|
|
(bignum-magnitude m)
|
|||
|
|
(bignum-magnitude n)))))
|
|||
|
|
|
|||
|
|
(define (integer-divide m n cont)
|
|||
|
|
(let ((m (integer->bignum m))
|
|||
|
|
(n (integer->bignum n)))
|
|||
|
|
(divide-magnitudes
|
|||
|
|
(bignum-magnitude m)
|
|||
|
|
(bignum-magnitude n)
|
|||
|
|
(lambda (q r)
|
|||
|
|
(cont (make-integer (* (bignum-sign m) (bignum-sign n)) q)
|
|||
|
|
(make-integer (bignum-sign m) r))))))
|
|||
|
|
|
|||
|
|
(define (integer-quotient m n)
|
|||
|
|
(integer-divide m n (lambda (q r) q)))
|
|||
|
|
|
|||
|
|
(define (integer-remainder m n)
|
|||
|
|
(integer-divide m n (lambda (q r) r)))
|
|||
|
|
|
|||
|
|
(define integer=
|
|||
|
|
(lambda (m n)
|
|||
|
|
(let ((m (integer->bignum m))
|
|||
|
|
(n (integer->bignum n)))
|
|||
|
|
(and (= (bignum-sign m) (bignum-sign n))
|
|||
|
|
(same-magnitude? (bignum-magnitude m)
|
|||
|
|
(bignum-magnitude n))))))
|
|||
|
|
|
|||
|
|
(define integer<
|
|||
|
|
(lambda (m n)
|
|||
|
|
(let ((m (integer->bignum m))
|
|||
|
|
(n (integer->bignum n)))
|
|||
|
|
(let ((m-sign (bignum-sign m))
|
|||
|
|
(n-sign (bignum-sign n)))
|
|||
|
|
(or (< m-sign n-sign)
|
|||
|
|
(and (= m-sign n-sign)
|
|||
|
|
(if (< m-sign 0)
|
|||
|
|
(smaller-magnitude? (bignum-magnitude n)
|
|||
|
|
(bignum-magnitude m))
|
|||
|
|
(smaller-magnitude? (bignum-magnitude m)
|
|||
|
|
(bignum-magnitude n)))))))))
|
|||
|
|
|
|||
|
|
|
|||
|
|
; Magnitude (unsigned integer) arithmetic
|
|||
|
|
|
|||
|
|
(define log-radix 14) ;Cutting it close here...
|
|||
|
|
(define radix (expt 2 log-radix))
|
|||
|
|
|
|||
|
|
(define greatest-non-bignum (+ (expt 2 28) (- (expt 2 28) 1)))
|
|||
|
|
(define least-non-bignum (* (expt 2 28) -2))
|
|||
|
|
|
|||
|
|
(define zero-magnitude '())
|
|||
|
|
(define zero-magnitude? null?)
|
|||
|
|
|
|||
|
|
(define (low-digit m)
|
|||
|
|
(if (zero-magnitude? m)
|
|||
|
|
0
|
|||
|
|
(car m)))
|
|||
|
|
|
|||
|
|
(define (high-digits m)
|
|||
|
|
(if (zero-magnitude? m)
|
|||
|
|
m
|
|||
|
|
(cdr m)))
|
|||
|
|
|
|||
|
|
(define (adjoin-digit d m)
|
|||
|
|
(if (and (= d 0) (zero-magnitude? m))
|
|||
|
|
m
|
|||
|
|
(cons d m)))
|
|||
|
|
|
|||
|
|
(define (integer->magnitude n)
|
|||
|
|
(if (= n 0)
|
|||
|
|
zero-magnitude
|
|||
|
|
(let ((digit (remainder n radix)))
|
|||
|
|
(adjoin-digit digit
|
|||
|
|
(integer->magnitude (quotient n radix))))))
|
|||
|
|
|
|||
|
|
(define (magnitude->integer m)
|
|||
|
|
(if (zero-magnitude? m)
|
|||
|
|
0
|
|||
|
|
(+ (low-digit m)
|
|||
|
|
(* radix (magnitude->integer (high-digits m))))))
|
|||
|
|
|
|||
|
|
(define greatest-non-bignum-magnitude
|
|||
|
|
(integer->magnitude greatest-non-bignum))
|
|||
|
|
|
|||
|
|
(define least-non-bignum-magnitude
|
|||
|
|
(adjoin-digit (- 0 (remainder least-non-bignum radix))
|
|||
|
|
(integer->magnitude
|
|||
|
|
(- 0 (quotient least-non-bignum radix)))))
|
|||
|
|
|
|||
|
|
; Combine two numbers digitwise using op.
|
|||
|
|
|
|||
|
|
(define (combine-magnitudes m n op)
|
|||
|
|
(let recur ((m m) (n n) (carry 0))
|
|||
|
|
(if (and (zero-magnitude? m) (zero-magnitude? n))
|
|||
|
|
(integer->magnitude carry)
|
|||
|
|
(let ((result (+ carry (op (low-digit m) (low-digit n)))))
|
|||
|
|
(let ((q (quotient result radix))
|
|||
|
|
(r (remainder result radix)))
|
|||
|
|
(if (< r 0)
|
|||
|
|
(adjoin-digit (+ r radix)
|
|||
|
|
(recur (high-digits m)
|
|||
|
|
(high-digits n)
|
|||
|
|
(- q 1)))
|
|||
|
|
(adjoin-digit r
|
|||
|
|
(recur (high-digits m)
|
|||
|
|
(high-digits n)
|
|||
|
|
q))))))))
|
|||
|
|
|
|||
|
|
(define (add-magnitudes m n)
|
|||
|
|
(combine-magnitudes m n +))
|
|||
|
|
|
|||
|
|
(define (subtract-magnitudes m n)
|
|||
|
|
(combine-magnitudes m n -))
|
|||
|
|
|
|||
|
|
; Compare
|
|||
|
|
|
|||
|
|
(define same-magnitude? equal?)
|
|||
|
|
|
|||
|
|
(define (smaller-magnitude? m n)
|
|||
|
|
(let loop ((m m) (n n) (a #f))
|
|||
|
|
(cond ((zero-magnitude? m)
|
|||
|
|
(or (not (zero-magnitude? n)) a))
|
|||
|
|
((zero-magnitude? n) #f)
|
|||
|
|
(else
|
|||
|
|
(loop (high-digits m)
|
|||
|
|
(high-digits n)
|
|||
|
|
(or (< (low-digit m) (low-digit n))
|
|||
|
|
(and (= (low-digit m) (low-digit n)) a)))))))
|
|||
|
|
|
|||
|
|
; Multiply
|
|||
|
|
|
|||
|
|
(define (multiply-magnitudes m n)
|
|||
|
|
(let recur ((m m) (a zero-magnitude))
|
|||
|
|
(if (zero-magnitude? m)
|
|||
|
|
a
|
|||
|
|
(let ((a (combine-magnitudes a n (lambda (d e)
|
|||
|
|
(+ d (* e (low-digit m)))))))
|
|||
|
|
(adjoin-digit (low-digit a)
|
|||
|
|
(recur (high-digits m) (high-digits a)))))))
|
|||
|
|
|
|||
|
|
; Divide m/n: find q and r such that m = q*n + r, where 0 <= r < n.
|
|||
|
|
; Oh no... time to get out Knuth...
|
|||
|
|
; The main thing we don't do that Knuth does is to normalize the
|
|||
|
|
; divisor (n) by shifting it left.
|
|||
|
|
|
|||
|
|
(define (divide-magnitudes m n cont)
|
|||
|
|
(if (zero-magnitude? (high-digits n))
|
|||
|
|
(divide-by-digit m (low-digit n)
|
|||
|
|
(lambda (q r)
|
|||
|
|
(cont q (adjoin-digit r zero-magnitude))))
|
|||
|
|
(let recur ((m m) (cont cont))
|
|||
|
|
(if (smaller-magnitude? m n)
|
|||
|
|
(cont zero-magnitude m)
|
|||
|
|
(recur
|
|||
|
|
(high-digits m)
|
|||
|
|
(lambda (q r)
|
|||
|
|
;; 0 <= r < n and d < b
|
|||
|
|
;; so b*r + d < b*n.
|
|||
|
|
(divide-step (adjoin-digit (low-digit m) r)
|
|||
|
|
n
|
|||
|
|
(lambda (q1 r)
|
|||
|
|
(cont (adjoin-digit q1 q) r)))))))))
|
|||
|
|
|
|||
|
|
; Divide m by n, where n <= m < b*n, i.e. 1 <= quotient < b.
|
|||
|
|
; E.g. if n = 100 then 100 <= m <= 999
|
|||
|
|
; if n = 999 then 999 <= m <= 9989
|
|||
|
|
|
|||
|
|
(define (divide-step m n cont)
|
|||
|
|
(do ((m-high m (high-digits m-high))
|
|||
|
|
(n-high n (high-digits n-high)))
|
|||
|
|
((zero-magnitude? (high-digits (high-digits n-high)))
|
|||
|
|
;; Occasionally q^ is one larger than the actual first digit.
|
|||
|
|
;; This loop will never iterate more than once.
|
|||
|
|
(let loop ((q^ (min (guess-quotient-digit m-high n-high)
|
|||
|
|
(- radix 1))))
|
|||
|
|
(let ((r (combine-magnitudes m n (lambda (d e)
|
|||
|
|
(- d (* e q^))))))
|
|||
|
|
(if (improper-magnitude? r)
|
|||
|
|
;; (begin (write `(addback ,m ,n ,q^ ,r)) (newline) ...)
|
|||
|
|
(loop (- q^ 1))
|
|||
|
|
(cont q^ r)))))))
|
|||
|
|
|
|||
|
|
; Compute q such that [m1 m2 m3] = q*[n1 n2] + r with 0 <= r < [n1 n2]
|
|||
|
|
; Can assume b <= [0 n1 n2] <= [m1 m2 m3] <= [n1 n2 b-1]
|
|||
|
|
; Some examples:
|
|||
|
|
; m / n : 100[1] / 10[02], 099 / 10, 099[1] / 99[0], 999[8] / 99[99]
|
|||
|
|
; Various hacks are possible to improve performance. In particular, the
|
|||
|
|
; second division can be eliminated if the divisor is normalized.
|
|||
|
|
; See Knuth.
|
|||
|
|
; [m1 m2] = q0*[n1] + r0
|
|||
|
|
; [m1 m2 m3] = q0*[n1 n2] + r^
|
|||
|
|
; r^ = b*r0 + m3 - q0*n2
|
|||
|
|
|
|||
|
|
(define (guess-quotient-digit m n)
|
|||
|
|
(let ((n1 (low-digit (high-digits n)))
|
|||
|
|
(n2 (low-digit n))
|
|||
|
|
(m1 (low-digit (high-digits (high-digits m))))
|
|||
|
|
(m2 (low-digit (high-digits m)))
|
|||
|
|
(m3 (low-digit m)))
|
|||
|
|
(let ((m12 (+ (* m1 radix) m2)))
|
|||
|
|
(let ((q0 (quotient m12 n1))
|
|||
|
|
(r0 (remainder m12 n1)))
|
|||
|
|
(let ((r^ (- (+ (* radix r0) m3) (* q0 n2)))
|
|||
|
|
(n12 (+ (* n1 radix) n2)))
|
|||
|
|
(let ((q1 (quotient r^ n12))
|
|||
|
|
(r1 (remainder r^ n12)))
|
|||
|
|
(if (> q1 0)
|
|||
|
|
(begin (display "This should never happen: q1 = ")
|
|||
|
|
(write q1) (newline)))
|
|||
|
|
(let ((q (+ q0 q1)))
|
|||
|
|
(if (< r1 0) (- q 1) q))))))))
|
|||
|
|
|
|||
|
|
(define (improper-magnitude? m)
|
|||
|
|
(and (not (zero-magnitude? m))
|
|||
|
|
(or (< (low-digit m) 0)
|
|||
|
|
(improper-magnitude? (high-digits m)))))
|
|||
|
|
|
|||
|
|
; Special case of division algorithm for single-digit divisor.
|
|||
|
|
|
|||
|
|
(define (divide-by-digit m d cont)
|
|||
|
|
(if (= d 0)
|
|||
|
|
(error "integer division by zero" m d)
|
|||
|
|
(let recur ((m m) (cont cont))
|
|||
|
|
(if (and (zero-magnitude? (high-digits m))
|
|||
|
|
(< (low-digit m) d))
|
|||
|
|
(cont zero-magnitude (low-digit m))
|
|||
|
|
(recur (high-digits m)
|
|||
|
|
(lambda (q r)
|
|||
|
|
(let ((m1 (+ (low-digit m) (* radix r))))
|
|||
|
|
(cont (adjoin-digit (quotient m1 d) q)
|
|||
|
|
(remainder m1 d)))))))))
|
|||
|
|
|
|||
|
|
;(define (divide-test seed)
|
|||
|
|
; (let ((random (make-random seed)))
|
|||
|
|
; (let loop ()
|
|||
|
|
; (let* ((z1 (integer+ (random) (integer* (random) 10000000)))
|
|||
|
|
; (z2 (integer+ (random) (integer* (random) 10000000)))
|
|||
|
|
; (n (max z1 z2))
|
|||
|
|
; (r (min z1 z2))
|
|||
|
|
; (q (random))
|
|||
|
|
; (m (integer+ (integer* n q) r)))
|
|||
|
|
; (if (not (= n r))
|
|||
|
|
; (integer-divide m n
|
|||
|
|
; (lambda (q1 r1)
|
|||
|
|
; (if (and (= q q1) (= r r1))
|
|||
|
|
; (begin (display ".")
|
|||
|
|
; (force-output (current-output-port)))
|
|||
|
|
; (error "division error" m n q q1 r r1)))))
|
|||
|
|
; (loop)))))
|
|||
|
|
|
|||
|
|
|
|||
|
|
; Extend the generic arithmetic operators.
|
|||
|
|
|
|||
|
|
(define-method &integer? ((n :bignum)) #t)
|
|||
|
|
(define-method &exact? ((n :bignum)) #t)
|
|||
|
|
|
|||
|
|
(define-method &+ ((n1 :exact-integer) (n2 :exact-integer)) (integer+ n1 n2))
|
|||
|
|
(define-method &- ((n1 :exact-integer) (n2 :exact-integer)) (integer- n1 n2))
|
|||
|
|
(define-method &* ((n1 :exact-integer) (n2 :exact-integer)) (integer* n1 n2))
|
|||
|
|
(define-method &= ((n1 :exact-integer) (n2 :exact-integer)) (integer= n1 n2))
|
|||
|
|
(define-method &< ((n1 :exact-integer) (n2 :exact-integer)) (integer< n1 n2))
|
|||
|
|
|
|||
|
|
(define-method "ient ((n1 :exact-integer) (n2 :exact-integer))
|
|||
|
|
(integer-quotient n1 n2))
|
|||
|
|
(define-method &remainder ((n1 :exact-integer) (n2 :exact-integer))
|
|||
|
|
(integer-remainder n1 n2))
|