[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