[gambit-list] Patch to speed bignum division in common case
Marc Feeley
feeley at iro.umontreal.ca
Wed Sep 25 21:12:02 EDT 2013
You seem to be on a 32 bit CPU. That is probably an important parameter. Can you try on another, 64 bit, machine?
Marc
On 2013-09-25, at 7:45 PM, Bradley Lucier <lucier at math.purdue.edu> wrote:
> 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