[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