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)))))))