[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