Marc:
Please find enclosed a patch that speeds bignum division in the common case where (quotient a b) is a relatively small integer.
Please apply this patch even though it's not a git fetch request. It's been going around in my head for a long time now, and I'd like to get rid of it before I learn about git.
In pidigit.scm, a spigot algorithm for calculating the digits of pi, all (or almost all, I haven't checked) of the quotients return a single decimal digit. So this new code improves this performance of this benchmark a bit: To calculate 10,000 decimal digits of pi:
Before the patch:
firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 5.560u 0.052s 0:05.63 99.6% 0+0k 0+0io 0pf+0w firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 5.564u 0.044s 0:05.62 99.6% 0+0k 0+0io 0pf+0w firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 5.588u 0.048s 0:05.65 99.4% 0+0k 0+0io 0pf+0w
After the patch:
firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 4.684u 0.024s 0:04.72 99.5% 0+0k 0+0io 0pf+0w firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 4.676u 0.040s 0:04.73 99.5% 0+0k 0+0io 0pf+0w firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 4.648u 0.032s 0:04.69 99.5% 0+0k 0+0io 0pf+0w
The code may also improve the performance of (very) bignum gcd, but I haven't checked that.
Thoroughly "white box" and "black box" tested.
The code includes better bignum shift code that I plan to use later in other places.
Brad
Afficher les réponses par date
On 2013-09-20, at 7:25 PM, Bradley Lucier lucier@math.purdue.edu wrote:
Marc:
Please find enclosed a patch that speeds bignum division in the common case where (quotient a b) is a relatively small integer.
Please apply this patch even though it's not a git fetch request. It's been going around in my head for a long time now, and I'd like to get rid of it before I learn about git.
In pidigit.scm, a spigot algorithm for calculating the digits of pi, all (or almost all, I haven't checked) of the quotients return a single decimal digit. So this new code improves this performance of this benchmark a bit: To calculate 10,000 decimal digits of pi:
Before the patch:
firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 5.560u 0.052s 0:05.63 99.6% 0+0k 0+0io 0pf+0w firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 5.564u 0.044s 0:05.62 99.6% 0+0k 0+0io 0pf+0w firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 5.588u 0.048s 0:05.65 99.4% 0+0k 0+0io 0pf+0w
After the patch:
firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 4.684u 0.024s 0:04.72 99.5% 0+0k 0+0io 0pf+0w firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 4.676u 0.040s 0:04.73 99.5% 0+0k 0+0io 0pf+0w firefly:~/programs/gambit/4.6.8/gambit> time ./pidigits2-gambit > /dev/null 4.648u 0.032s 0:04.69 99.5% 0+0k 0+0io 0pf+0w
The code may also improve the performance of (very) bignum gcd, but I haven't checked that.
Thoroughly "white box" and "black box" tested.
The code includes better bignum shift code that I plan to use later in other places.
Brad <pidigits2-gambit.scm><_num-scm.patch>
Unfortunately there's something wrong with your patch. Here's the output of "make check" on my Mac.
Marc
------------ TEST 3 (interpreter and library functions) ../gsi/gsi -:s,~~bin=../bin,~~lib=../lib,~~include=../include -f -e '(begin (load "r4rstest.scm") (test-cont) (test-sc4) (test-delay) (exit))' > test3.out diff test3.ok test3.out && rm -f test3.out tmp* 582,584c582,601 < ==> #t < ==> #t < ==> #t ---
Number readback failure for (+ 1. (* -97 2.220446049250313e-16)) .9999999999999785 Number readback failure for (+ 10. (* -94 1.7763568394002505e-15)) 9.999999999999833 Number readback failure for (+ 100. (* -98 1.4210854715202004e-14)) 99.99999999999861 ==> #f BUT EXPECTED #t Number readback failure for (+ 3. (* -85 4.440892098500626e-16)) 2.9999999999999623 Number readback failure for (+ 30. (* -83 3.552713678800501e-15)) 29.999999999999705 Number readback failure for (+ 300. (* -80 5.684341886080802e-14)) 299.99999999999545 ==> #f BUT EXPECTED #t Number readback failure for (+ 7. (* -94 8.881784197001252e-16)) 6.9999999999999165 ==> #f BUT EXPECTED #t
621c638,643 < Passed all tests ---
errors were: (SECTION (got expected (call))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f)))
631c653,658 < Passed all tests ---
errors were: (SECTION (got expected (call))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f)))
647c674,679 < Passed all tests ---
errors were: (SECTION (got expected (call))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f)))
658c690,695 < Passed all tests ---
errors were: (SECTION (got expected (call))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f)))
make[1]: *** [test3] Error 1 make: *** [check] Error 2
On 09/24/2013 05:00 PM, Marc Feeley wrote:
On 2013-09-20, at 7:25 PM, Bradley Lucier lucier@math.purdue.edu wrote:
Marc:
Please find enclosed a patch that speeds bignum division in the common case where (quotient a b) is a relatively small integer. gits2-gambit.scm><_num-scm.patch>
Unfortunately there's something wrong with your patch. Here's the output of "make check" on my Mac.
Hmmm. I don't see any problems when I run make check. I just regenerated the patch, but the file is the same.
I generated my patch on Ubuntu Linux; could this be a line ending mismatch problem?
I'm resending the patch in case it got corrupted the first time.
Brad
On 2013-09-25, at 11:13 AM, Bradley Lucier lucier@math.purdue.edu wrote:
On 09/24/2013 05:00 PM, Marc Feeley wrote:
On 2013-09-20, at 7:25 PM, Bradley Lucier lucier@math.purdue.edu wrote:
Marc:
Please find enclosed a patch that speeds bignum division in the common case where (quotient a b) is a relatively small integer. gits2-gambit.scm><_num-scm.patch>
Unfortunately there's something wrong with your patch. Here's the output of "make check" on my Mac.
Hmmm. I don't see any problems when I run make check. I just regenerated the patch, but the file is the same.
I generated my patch on Ubuntu Linux; could this be a line ending mismatch problem?
I'm resending the patch in case it got corrupted the first time.
Brad
<_num-scm.patch>
There's still a problem:
% uname -a Darwin neo.local 12.4.0 Darwin Kernel Version 12.4.0: Wed May 1 17:57:12 PDT 2013; root:xnu-2050.24.15~1/RELEASE_X86_64 x86_64 % patch -p1 < _num-scm.patch patching file lib/_num.scm Hunk #1 succeeded at 9043 (offset -2 lines). Hunk #2 succeeded at 9269 (offset -2 lines). % make check making all in include ... ------------ TEST 1 (debugging support) ../gsi/gsi -:~~bin=../bin,~~lib=../lib,~~include=../include -f debug.scm > test1.out diff test1.ok test1.out && rm -f test1.out ------------ TEST 2 (error handling) ../gsi/gsi -:m1,h4000,~~bin=../bin,~~lib=../lib,~~include=../include -f error.scm < error.scm > test2.out diff test2.ok test2.out && rm -f test2.out ------------ TEST 3 (interpreter and library functions) ../gsi/gsi -:s,~~bin=../bin,~~lib=../lib,~~include=../include -f -e '(begin (load "r4rstest.scm") (test-cont) (test-sc4) (test-delay) (exit))' > test3.out diff test3.ok test3.out && rm -f test3.out tmp* 582,584c582,601 < ==> #t < ==> #t < ==> #t ---
Number readback failure for (+ 1. (* -97 2.220446049250313e-16)) .9999999999999785 Number readback failure for (+ 10. (* -94 1.7763568394002505e-15)) 9.999999999999833 Number readback failure for (+ 100. (* -98 1.4210854715202004e-14)) 99.99999999999861 ==> #f BUT EXPECTED #t Number readback failure for (+ 3. (* -85 4.440892098500626e-16)) 2.9999999999999623 Number readback failure for (+ 30. (* -83 3.552713678800501e-15)) 29.999999999999705 Number readback failure for (+ 300. (* -80 5.684341886080802e-14)) 299.99999999999545 ==> #f BUT EXPECTED #t Number readback failure for (+ 7. (* -94 8.881784197001252e-16)) 6.9999999999999165 ==> #f BUT EXPECTED #t
621c638,643 < Passed all tests ---
errors were: (SECTION (got expected (call))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f)))
631c653,658 < Passed all tests ---
errors were: (SECTION (got expected (call))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f)))
647c674,679 < Passed all tests ---
errors were: (SECTION (got expected (call))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f)))
658c690,695 < Passed all tests ---
errors were: (SECTION (got expected (call))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f))) ((6 5 6) (#f #t (mult-float-print-test #f)))
make[1]: *** [test3] Error 1 make: *** [check] Error 2
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))))))))
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@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 bootstrap10 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))))))))
On 09/25/2013 09:58 PM, Bradley Lucier wrote:
The problem is with -m32, not -m64. I'll look at it, sorry.
Here is a fix. It passes "make check" on both 32-bit and 64-bit x86.
I've included a little quotient "torture test" in case you want to incorporate it into error.scm.
Brad
The times for 32-bit on my Core i7:
Original:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 13233 ms real time 13192 ms cpu time (13040 user, 152 system) 2840 collections accounting for 266 ms real time (224 user, 4 system) 12119389176 bytes allocated 192981 minor faults no major faults 13.040u 0.152s 0:13.23 99.6% 0+0k 0+0io 0pf+0w
New:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 11236 ms real time 11200 ms cpu time (11108 user, 92 system) 2300 collections accounting for 216 ms real time (232 user, 0 system) 9808248640 bytes allocated 127024 minor faults no major faults 11.120u 0.092s 0:11.24 99.7% 0+0k 0+0io 0pf+0w
For 64-bit:
Original:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 5632 ms real time 5612 ms cpu time (5552 user, 60 system) 2965 collections accounting for 203 ms real time (212 user, 4 system) 12173461160 bytes allocated 59591 minor faults no major faults 5.556u 0.060s 0:05.63 99.6% 0+0k 0+0io 0pf+0w
New:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 4645 ms real time 4632 ms cpu time (4592 user, 40 system) 2389 collections accounting for 159 ms real time (140 user, 0 system) 9831319240 bytes allocated 35307 minor faults no major faults 4.592u 0.040s 0:04.64 99.7% 0+0k 0+0io 0pf+0w
On 09/30/2013 08:52 PM, Bradley Lucier wrote:
The times for 32-bit on my Core i7:
<stuff deleted>
On my Q6600, which is the same, I believe, as as used here:
http://benchmarksgame.alioth.debian.org/u64/benchmark.php?test=pidigits&...
The timings are:
64 bit:
Old:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 16251 ms real time 16228 ms cpu time (15484 user, 744 system) 9987 collections accounting for 3023 ms real time (2888 user, 140 system) 12173870360 bytes allocated 263552 minor faults no major faults 15.488u 0.748s 0:16.26 99.7% 0+0k 0+0io 0pf+0w
New:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 13289 ms real time 13272 ms cpu time (12752 user, 520 system) 8065 collections accounting for 2421 ms real time (2268 user, 96 system) 9831650712 bytes allocated 205741 minor faults no major faults 12.760u 0.520s 0:13.30 99.8% 0+0k 0+0io 0pf+0w
32-bit:
Old:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 31351 ms real time 31308 ms cpu time (30512 user, 796 system) 11396 collections accounting for 2461 ms real time (2460 user, 68 system) 12119457704 bytes allocated 266083 minor faults no major faults 30.516u 0.800s 0:31.36 99.8% 0+0k 0+0io 0pf+0w
New:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 27339 ms real time 27300 ms cpu time (26824 user, 476 system) 9223 collections accounting for 2096 ms real time (2092 user, 56 system) 9808284504 bytes allocated 212081 minor faults no major faults 26.832u 0.476s 0:27.34 99.8% 0+0k 0+0io 0pf+0w
So current gmp kills us (all loops are in hand-coded assembler), but we don't do too badly compared to the other implementations.
Brad
Thanks for the patch. It is not committed and pushed to the repo.
What was the problem with your first patch?
Can you send me your pidigits2-gambit.scm file? It might make an interesting parallel benchmark.
Marc
On 2013-09-30, at 8:52 PM, Bradley Lucier lucier@math.purdue.edu wrote:
On 09/25/2013 09:58 PM, Bradley Lucier wrote:
The problem is with -m32, not -m64. I'll look at it, sorry.
Here is a fix. It passes "make check" on both 32-bit and 64-bit x86.
I've included a little quotient "torture test" in case you want to incorporate it into error.scm.
Brad
The times for 32-bit on my Core i7:
Original:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 13233 ms real time 13192 ms cpu time (13040 user, 152 system) 2840 collections accounting for 266 ms real time (224 user, 4 system) 12119389176 bytes allocated 192981 minor faults no major faults 13.040u 0.152s 0:13.23 99.6% 0+0k 0+0io 0pf+0w
New:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 11236 ms real time 11200 ms cpu time (11108 user, 92 system) 2300 collections accounting for 216 ms real time (232 user, 0 system) 9808248640 bytes allocated 127024 minor faults no major faults 11.120u 0.092s 0:11.24 99.7% 0+0k 0+0io 0pf+0w
For 64-bit:
Original:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 5632 ms real time 5612 ms cpu time (5552 user, 60 system) 2965 collections accounting for 203 ms real time (212 user, 4 system) 12173461160 bytes allocated 59591 minor faults no major faults 5.556u 0.060s 0:05.63 99.6% 0+0k 0+0io 0pf+0w
New:
time ./pidigits2-gambit > /dev/null (time (digits 10000)) 4645 ms real time 4632 ms cpu time (4592 user, 40 system) 2389 collections accounting for 159 ms real time (140 user, 0 system) 9831319240 bytes allocated 35307 minor faults no major faults 4.592u 0.040s 0:04.64 99.7% 0+0k 0+0io 0pf+0w
<_num-scm.patch><test-quotient.scm>
On 10/01/2013 01:38 PM, Marc Feeley wrote:
What was the problem with your first patch?
It's complicated. Basically, I didn't compute the last borrow correctly in certain situations. And if the shifted divisor had to be added back in because the guessed quotient digit was too big by 1, sometimes I didn't add the last carry correctly.
Can you send me your pidigits2-gambit.scm file? It might make an interesting parallel benchmark.
It was in the original email, but here it is again. What do you mean by "parallel"?
Brad
On 2013-10-01, at 1:38 PM, Marc Feeley feeley@iro.umontreal.ca wrote:
Thanks for the patch. It is not committed and pushed to the repo.
I meant "now" not "not"...
Marc
-
Bradley Lucier -
Marc Feeley