[gambit-list] New fft bignum multiplication

Bradley Lucier lucier at math.purdue.edu
Thu May 17 19:48:21 EDT 2007


Marc:

After many years, I finally understand what Percival, Crandall,  
Brandt, et al. are trying to tell me about how to do fft bignum  
multiplication, and I implemented it.  Here are the timing results  
(in ms) on my 2GHz G5 comparing the code I sent you recently with the  
new code and mzscheme 360.  Results that had an inordinate gc time  
are noted.

                                                   old                
new         mzscheme

(expt 3 1000000)  ; a                             175                
101            273
(expt 3 1000001)  ; b                             176                
103            277
(* a a)           ; c                             204                
116            478
(* a b)           ;                               302                
165            667
(quotient c a)    ;                              1433                
825           1670
(sqrt c)          ;                              1430                
865           1176
(gcd a b)         ; note 1                       5385               
3979           6759
(* a b)                                           145                 
84            191
(expt1 3 1000000) ; note 2                        179                
106            277
(expt2 3 1000000) ; note 3                        763                
433            957
(* a a)           ; a=3^1000000                   202                
160 (48 gc)    473
(expt 10 10000000); a                             450                
227            492
(fib 10000)        ; note 4                        25                 
22            130 (89 gc)
(factorial 10000) ; note 5                        311                
310            666 (486 gc)
(partial-factorial 0 10000) ; note 6               29                 
19             30
(binary-splitting-compute-e 10000) ; note 7      1692               
1469           1857
(binary-splitting-compute-pi 10000) ; note 9     2093               
1820           2334
(pi-brent-salamin) ; n. 10, beta^k=10^100000    10325 (764 gc)      
5918 (473 gc)  8051
(pi-brent-salamin) ; beta^k=2^332193             7703 (683 gc)      
4497 (462 gc)  3985

MzScheme uses its own hybrid version (part of gmp 3, part gmp 4) of  
gmp without fft multiplication and some of the more-recently-added  
algorithms.

The new code uses less memory, so the number of gc's were reduced,  
too.  We really need to improve our rational arithmetic if such a big  
improvement in basic large operations results in such a small  
improvement in code using rationals (binary-splitting-compute-e and  
binary-splitting-compute-pi).  These times scale to arguments up to  
half a billion bits in size (after that, we slow down because we  
switch to Karatsuba multiply again because the floating-point error  
bounds get too bad).

Some timings for similar things from December 2003 can be found at

<http://groups.google.com/group/comp.lang.scheme/msg/4306a6bd50004f9a>

Note the different parameters for many of the timings.

I'm going to look at its packaging a bit more before sending it to you.

So now, instead of gambit taking about 5 times as long as gmp-4.2.1  
for multiplication/division/sqrt and being a bit faster for gcd,  
gambit takes "only" a bit less than 3 times as long as gmp-4.2.1 for  
multiplication/division/sqrt and is noticeably faster for gcd.  (That  
last bit will change with version 5 of gmp, which will use a faster  
version of gambit's algorithm that was written by Niels Möller of  
Sweden.)  The basic fft code in _num.scm is about half as fast as,  
e.g., FFTW so this isn't too bad.

Brad

note 1:

a=(fib 1000000) b=(fib 1000001)

note 2:

(define (expt1 a b)
   (define (square x) (* x x))
   (cond ((= b 0) 1)
         ((even? b)
          (square (expt1 a (quotient b 2))))
         (else
          (* a (square (expt1 a (quotient b 2)))))))

note 3:

(define (expt2 a b)
   (define (square x) (* x x))
   (cond ((= b 0) 1)
         ((even? b)
          (expt2 (square a) (quotient b 2)))
         (else
          (* a (expt2 (square a) (quotient b 2))))))

note 4:

(define (fib-ratio n)
   (if (= n 1)
       1
       (+ 1 (/ (fib-ratio (- n 1))))))
(define (fib n)
   (numerator (fib-ratio n)))

note 5:

(define (factorial n)
   (let loop ((i 1)
              (result 1))
     (if (> i n)
         result
         (loop (+ i 1)
               (* i result)))))
note 6:

(define (partial-factorial m n)
   ;; computes the product (m+1) * ... * (n-1) * n
   (if (< (- n m) 10)
       (do ((i (+ m 1) (+ i 1))
            (result 1 (* result i)))
           ((> i n) result))
       (* (partial-factorial m (quotient (+ m n) 2))
          (partial-factorial (quotient (+ m n) 2) n))))

note 6:

(define (partial-factorial m n)
   ;; computes the product (m+1) * ... * (n-1) * n
   (if (< (- n m) 10)
       (do ((i (+ m 1) (+ i 1))
            (result 1 (* result i)))
           ((> i n) result))
       (* (partial-factorial m (quotient (+ m n) 2))
          (partial-factorial (quotient (+ m n) 2) n))))

note 7:

(define (binary-splitting-partial-sum m n
                                       partial-term
                                       common-factor-ratio)
   ;; sums (partial) terms from m to n-1
   ;; (partial-term n m) is the term at n with the common factors of  
terms >= m removed
   ;; (common-factor-ratio m n) is the ratio of the common factor of  
terms >= n divided by
   ;; the common factors of terms >= m
   (if (< (- n m) 10)
       (do ((i m (+ i 1))
            (result 0 (+ result (partial-term m i))))
           ((= i n) result))
       (+ (binary-splitting-partial-sum m (quotient (+ m n) 2)  
partial-term common-factor-ratio)
          (* (common-factor-ratio m (quotient (+ m n) 2))
             (binary-splitting-partial-sum (quotient (+ m n) 2) n  
partial-term common-factor-ratio)))))

(define (binary-splitting-sum n partial-term common-factor)
   (binary-splitting-partial-sum 0 n partial-term common-factor))

(define (binary-splitting-compute-e n)
   (binary-splitting-sum n
                        (lambda (m n) (/ (partial-factorial m n)))
                        (lambda (m n) (/ (partial-factorial m n)))))

note 9:

(define (binary-splitting-compute-atan n x)
   ;; here we just consider the common factor to be x^(2n+1)
   (* x    ;; common factor for all terms
      (binary-splitting-sum n
                            (lambda (m n) (/ (expt x (* 2 (- n m)))
                                          (* (if (odd? n) -1 1) ( +  
(* 2 n) 1))))
                            (lambda (m n) (expt x (* 2 (- n m)))))))

(define (binary-splitting-compute-pi n)
   (* 4 (- (* 4 (binary-splitting-compute-atan n 1/5))
           (binary-splitting-compute-atan (quotient (* n 10) 34)  
1/239))))

note 10:

(define (fixed.+ x y)
   (+ x y))
(define (fixed.- x y)
   (- x y))
(define (fixed.* x y)
   (quotient (* x y) beta^k))
(define (fixed.square x)
   (fixed.* x x))
(define (fixed./ x y)
   (quotient (* x beta^k) y))
(define (fixed.sqrt x)
   (##exact-int.sqrt (* x beta^k)))
(define (number->fixed x)
   (round (* x beta^k)))
(define (fixed->number x)
   (/ x beta^k))

(define (pi-brent-salamin)
   (let ((one (number->fixed 1)))
     (let loop ((a one)
                (b (fixed.sqrt (quotient one 2)))
                (t (quotient one 4))
                (x 1))
       (if (= a b)
           (fixed./ (fixed.square a) t)
           (let ((new-a (quotient (fixed.+ a b) 2)))
             (loop new-a
                   (fixed.sqrt (fixed.* a b))
                   (fixed.- t (* x (fixed.square (fixed.- new-a a))))
                   (* 2 x)))))))




More information about the Gambit-list mailing list