[gambit-list] Patch to speed bignum division in common case
Bradley Lucier
lucier at math.purdue.edu
Wed Sep 25 19:45:35 EDT 2013
Here's what I did on my Mac portable, and didn't find any problems:
1 19:16 pu programs/gambit-devel/
3 19:16 git clone https://github.com/feeley/gambit.git
4 19:19 cd gambit
6 19:21 ./configure 'CC=/pkgs/gcc-4.8.1/bin/gcc -march=native' '--enable-multiple-versions' '--enable-single-host'
7 19:21 make -j4 bootstrap
8 19:23 make -j4 bootclean
9 19:23 make -j4 bootstrap
10 19:25 make -j4 bootclean
11 19:25 make -j4
12 19:28 make check
13 19:33 git apply --stat _num-scm.patch
14 19:33 git apply --check _num-scm.patch
15 19:34 git apply _num-scm.patch
16 19:34 vi lib/_num.scm
17 19:35 make
18 19:36 make check
The 'git apply _num-scm.patch' command gave the following output
[Bradley-Luciers-MacBook-Pro:~/programs/gambit-devel/gambit] lucier% git apply _num-scm.patch
_num-scm.patch:14: trailing whitespace.
_num-scm.patch:18: trailing whitespace.
_num-scm.patch:229: trailing whitespace.
_num-scm.patch:230: trailing whitespace.
_num-scm.patch:259: trailing whitespace.
warning: squelched 3 whitespace errors
warning: 8 lines add whitespace errors.
Here's my environment
[Bradley-Luciers-MacBook-Pro:~/programs/gambit-devel/gambit] lucier% gsi/gsi -v
v4.7.0 20130924213249 i386-apple-darwin10.8.0 "./configure 'CC=/pkgs/gcc-4.8.1/bin/gcc -march=native' '--enable-multiple-versions' '--enable-single-host'"
[Bradley-Luciers-MacBook-Pro:~/programs/gambit-devel/gambit] lucier% uname -a
Darwin Bradley-Luciers-MacBook-Pro.local 10.8.0 Darwin Kernel Version 10.8.0: Tue Jun 7 16:33:36 PDT 2011; root:xnu-1504.15.3~1/RELEASE_I386 i386
So I don't know what's happening. You seem to be running 10.8; I'm running 10.6.8, but I don't see how that could matter.
Here are the two routines inlined to this message in case you want to apply them by hand. The first routine is a new routine at top level, the second replaces the existing naive-div.
Brad
(define-prim (##bignum.arithmetic-shift-into! x shift result)
#|
Shifts x by shift bits into result.
Will eventually replace other "shift"ing code.
Left pads by sign bit as necessary, right pads by zeros as necessary.
Makes *no* error checks.
|#
;; allocates nothing
(declare (not interrupts-enabled))
(let* ((bit-shift
(##fxmodulo shift ##bignum.adigit-width))
(digit-shift
(##fxquotient (##fx- shift bit-shift)
##bignum.adigit-width))
(x-length
(##bignum.adigit-length x))
(result-length
(##bignum.adigit-length result))
(zeros
##bignum.adigit-zeros)
(left-fill
(if (##bignum.negative? x)
##bignum.adigit-ones
##bignum.adigit-zeros)))
(if (##fxzero? bit-shift)
;; Copy left-fill into leftmost digits of result as needed.
(let loop1 ((i (##fx- result-length 1)) ; index for adigit in result
(j (##fx- result-length 1 digit-shift))) ; index for adigit in x
(if (and (##fx>= i 0) (##fx>= j x-length))
(begin (##bignum.adigit-copy! result i left-fill 0)
(loop1 (##fx- i 1) (##fx- j 1)))
;; Copy the digits from x into result as needed.
(let loop2 ((i i)
(j j))
(if (and (##fx>= i 0) (##fx>= j 0))
(begin (##bignum.adigit-copy! result i x j)
(loop2 (##fx- i 1) (##fx- j 1)))
;; copy zero into digits of result as needed.
(let loop3 ((i i))
(if (##fx>= i 0)
(begin (##bignum.adigit-copy! result i zeros 0)
(loop3 (##fx- i 1)))))))))
(let ()
;; copy left-fill into leftmost digits of result as needed,
;; then concatenate left-fill with leftmost digit of x if needed.
(define (loop4 i j)
(if (and (##fx>= i 0) (##fx>= j x-length))
(begin (##bignum.adigit-copy! result i left-fill 0)
(loop4 (##fx- i 1) (##fx- j 1)))
(if (##fx>= i 0)
(if (##fx= (##fx+ j 1) x-length)
(begin (##bignum.adigit-cat! result i left-fill 0 x j bit-shift)
(loop5 (##fx- i 1) (##fx- j 1)))
(loop5 i j)))))
;; concatenate adjacent digits of x into result as needed,
;; then concatenate rightmost digit of x with 0 if needed.
(define (loop5 i j)
(if (and (##fx>= i 0) (##fx>= j 0))
(begin (##bignum.adigit-cat! result i x (##fx+ j 1) x j bit-shift)
(loop5 (##fx- i 1) (##fx- j 1)))
(if (##fx>= i 0)
(if (##fx= (##fx+ j 1) 0)
(begin (##bignum.adigit-cat! result i x 0 zeros 0 bit-shift)
(loop6 (##fx- i 1)))
(loop6 i)))))
;; copy 0 into rightmost digits of x as needed.
(define (loop6 i)
(if (##fx>= i 0)
(begin (##bignum.adigit-copy! result i zeros 0)
(loop6 (##fx- i 1)))))
(loop4 (##fx- result-length 1) ; index for adigit in result
(##fx- result-length digit-shift 2)))) ; index for adigit in x
;; return something useful
result))
(define (naive-div u v)
;; u is a normalized bignum, v is a possibly unnormalized bignum
;; u >= v >= ##bignum.mdigit-base
;; on a 64-bit machine, allocates three words (24 bytes) for temp
(let ((need-quotient? #t) ;; will be made an optional argument later
(keep-dividend? #t) ;; will be made an optional argument later
(n
(let loop1 ((i (##fx- (##bignum.mdigit-length v) 1)))
(if (##fx< 0 (##bignum.mdigit-ref v i))
(##fx+ i 1)
(loop1 (##fx- i 1)))))
(u-bits
(##integer-length u))
(v-bits
(##integer-length v)))
(let* ((temp
(##bignum.make (if (and (##fx= ##bignum.adigit-width 64)
(##fx= ##bignum.mdigit-width 16))
;; need three mdigits for top-bits-of-u, which will fit into one adigit
1
;; in the other cases, we need two adigits
2)
#f #f))
(top-2*mdigit-width-bits-of-v
(##bignum.arithmetic-shift-into! v
(##fx- (##fx* ##bignum.mdigit-width 2)
v-bits)
temp))
(v_n-1
(##bignum.mdigit-ref top-2*mdigit-width-bits-of-v 1))
(v_n-2
(##bignum.mdigit-ref top-2*mdigit-width-bits-of-v 0)))
;; Knuth says to simplify things by shifting u and v so that
;; the top nonzero mdigit of v is >= mdigit-base/2
;; We're not going to do the shift, but we're going to use that
;; idea to do the next calculation.
;; this strategy does more work, but generates less garbage.
(let* ((conceptual-shift
(##fx- ##bignum.mdigit-width
(##fxlength (##bignum.mdigit-ref v (##fx- n 1)))))
(shifted-v-adigits
(##fxquotient (##fx+ v-bits
conceptual-shift
64)
##bignum.adigit-width))
(shifted-u-adigits
(##fxquotient (##fx+ u-bits
conceptual-shift
64)
##bignum.adigit-width))
(q-adigits
(##fx+ (##fx- shifted-u-adigits
shifted-v-adigits)
2)) ; 1 is not always sufficient...
(q-mdigits
(##fxquotient (##fx* q-adigits ##bignum.adigit-width)
##bignum.mdigit-width))
(q
(and need-quotient? (##bignum.make q-adigits #f #f)))
(u
(if keep-dividend?
;; copy u
(##bignum.make (##bignum.adigit-length u) u #f)
;; overwrite u with remainder
u)))
(let loop3 ((j (##fx- q-mdigits 1)))
(if (##not (##fx< j 0))
(let* ((top-bits-of-u
(##bignum.arithmetic-shift-into!
u
(##fx- (##fx* (##fx- 2 j) ##bignum.mdigit-width)
v-bits)
temp))
(q-hat
(let ((q-hat
(##bignum.mdigit-quotient
top-bits-of-u
2
v_n-1))
(u_n+j-2
(##bignum.mdigit-ref
top-bits-of-u
0
)))
(let ((r-hat
(##bignum.mdigit-remainder
top-bits-of-u
2
v_n-1
q-hat)))
(if (or (##fx= q-hat ##bignum.mdigit-base)
(##bignum.mdigit-test?
q-hat
v_n-2
r-hat
u_n+j-2))
(let ((q-hat
(##fx- q-hat 1))
(r-hat
(##fx+ r-hat v_n-1)))
(if (and (##fx< r-hat ##bignum.mdigit-base)
(or (##fx= q-hat ##bignum.mdigit-base)
(##bignum.mdigit-test?
q-hat
v_n-2
r-hat
u_n+j-2)))
(##fx- q-hat 1)
q-hat))
q-hat)))))
(##declare (not interrupts-enabled))
(if (##fx= q-hat 0)
(begin
(and need-quotient?
(##bignum.mdigit-set! q j q-hat))
(loop3 (##fx- j 1)))
(let loop4 ((i j)
(k 0)
(borrow 0))
(if (##fx< k n)
(loop4 (##fx+ i 2)
(##fx+ k 2)
(##bignum.mdigit-div!
u
(##fx+ i 1)
v
(##fx+ k 1)
q-hat
(##bignum.mdigit-div!
u
i
v
k
q-hat
borrow)))
(let ((borrow
(if (or (##fx< n k)
(##fx= i (##bignum.mdigit-length u))
(##fx= (##bignum.mdigit-ref u i) 0))
borrow
(##bignum.mdigit-div!
u
i
##bignum.adigit-zeros
0
q-hat
borrow))))
(if (##fx< borrow 0)
(let loop5 ((i j)
(l 0)
(carry 0))
(if (##fx>= l n)
(begin
(if (##not (or (##fx< n k)
(##fx= i (##bignum.mdigit-length u))
(##fx= (##bignum.mdigit-ref u i) 0)))
(##bignum.mdigit-mul!
u
i
##bignum.adigit-zeros
0
1
carry))
(and need-quotient?
(##bignum.mdigit-set! q j (##fx- q-hat 1)))
(loop3 (##fx- j 1)))
(loop5 (##fx+ i 2)
(##fx+ l 2)
(##bignum.mdigit-mul!
u
(##fx+ i 1)
v
(##fx+ l 1)
1
(##bignum.mdigit-mul!
u
i
v
l
1
carry)))))
(begin
(and need-quotient?
(##bignum.mdigit-set! q j q-hat))
(loop3 (##fx- j 1)))))))))
(##cons (and need-quotient?
(##bignum.normalize! q))
(##bignum.normalize! u))))))))
More information about the Gambit-list
mailing list