On 2011-03-23, at 9:32 PM, chevalma@iro.umontreal.ca wrote:
Ha! I hadn't thought of that! Indeed if fixnums are artificially limited to only have 52 significant bits as you are doing, then the radix must be <= 2^26.
That would probably be sufficient for the case I'm dealing with, but for the sake of future-proofing, would it really be that difficult to port your existing bignum division implementation?
Well, take a look. How long would it take to rewrite this algorithm in JavaScript and test it on all the corner cases? I just don't feel it is a top priority right now.
Can you explain where this is needed in the Tachyon compiler? Is it blocking us from progressing? Moreover, in JavaScript there is no integer division, so how can the bignum division help constant fold expressions like 101/4 or 2/3 ?
Marc
;;---------------------------------------------------------------------------
;; Bignum division
(define (bignum-div x y)
(define (single-digit-divisor-div x y lenx leny r)
;; simple algo for single digit divisor
(let ((d (bignum-digit-ref y 1))) (let loop1 ((i (- lenx 1)) (k 0)) (if (< 0 i) (let ((w (+ (* k (bignum-radix)) (bignum-digit-ref x i)))) (bignum-digit-set! r i (quotient w d)) (loop1 (- i 1) (remainder w d))) (cons (bignum-remove-leading-zeroes r) (fixnum->bignum (if (bignum-negative? x) (- 0 k) k)))))))
(define (multi-digit-divisor-div x y lenx leny r)
;; general algo from Knuth
;; STEP 1: normalize x and y
(let loop2 ((shift 1) (n (* (bignum-digit-ref y (- leny 1)) 2))) (if (< n (bignum-radix)) (loop2 (* shift 2) (* n 2))
(let ((nx (make-bignum (+ lenx 1))) (ny (make-bignum leny)))
(bignum-sign-set! nx (bignum-sign x))
(let loop3 ((i 1) (c 0)) (if (< i lenx) (let ((w (+ (* (bignum-digit-ref x i) shift) c))) (bignum-digit-set! nx i (modulo w (bignum-radix))) (loop3 (+ i 1) (quotient w (bignum-radix)))) (bignum-digit-set! nx i c)))
(let loop4 ((i 1) (c 0)) (if (< i leny) (let ((w (+ (* (bignum-digit-ref y i) shift) c))) (bignum-digit-set! ny i (modulo w (bignum-radix))) (loop4 (+ i 1) (quotient w (bignum-radix))))))
(let loop5 ((i lenx)) (if (not (< i leny))
;; STEP 2: calculate next digit in quotient
(let ((msd-of-ny (bignum-digit-ref ny (- leny 1))) (next-msd-of-ny (bignum-digit-ref ny (- leny 2))) (msd-of-nx (bignum-digit-ref nx i)) (next-msd-of-nx (bignum-digit-ref nx (- i 1))) (next-next-msd-of-nx (bignum-digit-ref nx (- i 2))))
(define (next-digit q u) (if (< u (bignum-radix)) (let* ((temp1 (* q next-msd-of-ny)) (temp2 (quotient temp1 (bignum-radix)))) (if (or (< u temp2) (and (= temp2 u) (< next-next-msd-of-nx (remainder temp1 (bignum-radix))))) (next-digit (- q 1) (+ u msd-of-ny)) q)) q))
(let ((q (if (= msd-of-nx msd-of-ny) (next-digit (bignum-radix-minus-1) (+ msd-of-ny next-msd-of-nx)) (let ((temp (+ (* msd-of-nx (bignum-radix)) next-msd-of-nx))) (next-digit (quotient temp msd-of-ny) (modulo temp msd-of-ny))))))
;; STEP 3: multiply and substract
(let loop7 ((j 1) (k (- i (- leny 1))) (b 0)) (if (< j leny)
(let ((w (- (+ (bignum-digit-ref nx k) b) (* (bignum-digit-ref ny j) q)))) (bignum-digit-set! nx k (modulo w (bignum-radix))) (loop7 (+ j 1) (+ k 1) (quotient (- w (bignum-radix-minus-1)) (bignum-radix))))
(let ((w (+ (bignum-digit-ref nx k) b))) (bignum-digit-set! nx k (modulo w (bignum-radix))) (if (< w 0) (begin (bignum-digit-set! r (- i (- leny 1)) (- q 1)) (let loop8 ((j 1) (k (- i (- leny 1))) (c 0)) (if (< j leny)
(let ((w (+ (+ (bignum-digit-ref nx k) (bignum-digit-ref ny j)) c))) (bignum-digit-set! nx k (modulo w (bignum-radix))) (loop8 (+ j 1) (+ k 1) (quotient w (bignum-radix)))) (bignum-digit-set! nx k (modulo (+ (bignum-digit-ref nx k) c) (bignum-radix)))))) (bignum-digit-set! r (- i (- leny 1)) q)) (loop5 (- i 1)))))))))
(let loop9 ((i (- leny 1)) (k 0)) (if (< 0 i) (let ((w (+ (* k (bignum-radix)) (bignum-digit-ref nx i)))) (bignum-digit-set! nx i (quotient w shift)) (loop9 (- i 1) (remainder w shift)))))
(cons (bignum-remove-leading-zeroes r) (bignum-remove-leading-zeroes nx))))))
(define (div x y lenx leny) (if (< lenx leny)
(cons bignum-zero x)
(let ((r (make-bignum (+ (- lenx leny) 2))))
(if (= (bignum-sign x) (bignum-sign y)) (bignum-set-nonneg! r) (bignum-set-neg! r))
(if (= leny 2) (single-digit-divisor-div x y lenx leny r) (multi-digit-divisor-div x y lenx leny r)))))
(if (bignum-zero? y) (snow-error "divide by zero") (div x y (bignum-length x) (bignum-length y))))
(define (bignum-quotient x y) (bignum-normalize (car (bignum-div x y))))
(define (bignum-remainder x y) (bignum-normalize (cdr (bignum-div x y))))
(define (bignum-modulo x y) (let ((r (cdr (bignum-div x y)))) (if (or (bignum-zero? r) (eqv? (bignum-negative? x) (bignum-negative? y))) (bignum-normalize r) (bignum+ r y))))
;;---------------------------------------------------------------------------