No problem. They should be faster even for fairly small numbers since they usually require the evaluation of a polynomial (an approximation of (log)gamma) versus repeated multiplication/division. From memory the code should be something like:
(exp (fllog-gamma (+ 1.0 n)) - (fllog-gamma (+ 1.0 r)) - (fllog-gamma (+ 1.0 (- n r)))) fllog-gamma should also be faster than bflog-gamma or log-gamma if you don't need arbitrary precision. You're also right that this won't always give completely exact results - the Racket manual says that the only exact values are for log gamma of 1 and 2, but this usually is not a problem. PS. It looks like Racket's math collection has a built-in log-factorial function too, to avoid all the +1's, so you could try that. On Wed, Feb 20, 2013 at 1:44 AM, Joe Gilray <jgil...@gmail.com> wrote: > Hi Luke, > > Thanks for the knowledge. Do you have some code that I could try out. I > found gamma and bflog-gamma, but they work with floats and so I can't > imagine they are faster for exact answers... maybe for estimating nCr for > large numbers? > > -Joe > > > On Tue, Feb 19, 2013 at 8:26 PM, Luke Vilnis <lvil...@gmail.com> wrote: > >> FYI, log gamma is another fast way to calculate the number of >> combinations if you want to deal with really big numbers. >> >> On Tue, Feb 19, 2013 at 7:28 PM, Joe Gilray <jgil...@gmail.com> wrote: >> >>> Racketeers, >>> >>> Thanks for putting together the fantastic math library. It will be a >>> wonderful resource. Here are some quick impressions (after playing mostly >>> with math/number-theory) >>> >>> 1) The functions passed all my tests and were very fast. If you need >>> even more speed you can keep a list of primes around and write functions to >>> use that, but that should be rarely necessary >>> >>> 2) I have a couple of functions to donate if you want them: >>> >>> 2a) Probablistic primality test: >>> >>> ; function that performs a Miller-Rabin probabalistic primality test k >>> times, returns #t if n is probably prime >>> ; algorithm from http://rosettacode.org/wiki/Miller-Rabin_primality_test, >>> code adapted from Lisp example >>> ; (module+ test (check-equal? (is-mr-prime? 1000000000000037 8) #t)) >>> (define (is-mr-prime? n k) >>> ; function that returns two values r and e such that number = >>> divisor^e * r, and r is not divisible by divisor >>> (define (factor-out number divisor) >>> (do ([e 0 (add1 e)] [r number (/ r divisor)]) >>> ((not (zero? (remainder r divisor))) (values r e)))) >>> >>> ; function that performs fast modular exponentiation by repeated >>> squaring >>> (define (expt-mod base exponent modulus) >>> (let expt-mod-iter ([b base] [e exponent] [p 1]) >>> (cond >>> [(zero? e) p] >>> [(even? e) (expt-mod-iter (modulo (* b b) modulus) (/ e 2) p)] >>> [else (expt-mod-iter b (sub1 e) (modulo (* b p) modulus))]))) >>> >>> ; function to return a random, exact number in the passed range >>> (inclusive) >>> (define (shifted-rand lower upper) >>> (+ lower (random (add1 (- (modulo upper 4294967088) (modulo lower >>> 4294967088)))))) >>> >>> (cond >>> [(= n 1) #f] >>> [(< n 4) #t] >>> [(even? n) #f] >>> [else >>> (let-values ([(d s) (factor-out (- n 1) 2)]) ; represent n-1 as >>> 2^s-d >>> (let lp ([a (shifted-rand 2 (- n 2))] [cnt k]) >>> (if (zero? cnt) #t >>> (let ([x (expt-mod a d n)]) >>> (if (or (= x 1) (= x (sub1 n))) (lp (shifted-rand 2 (- n >>> 2)) (sub1 cnt)) >>> (let ctestlp ([r 1] [ctest (modulo (* x x) n)]) >>> (cond >>> [(>= r s) #f] >>> [(= ctest 1) #f] >>> [(= ctest (sub1 n)) (lp (shifted-rand 2 (- n 2)) >>> (sub1 cnt))] >>> [else (ctestlp (add1 r) (modulo (* ctest ctest) >>> n))])))))))])) >>> >>> 2b) combinations calculator >>> >>> ; function that returns the number of combinations, not the combinations >>> themselves >>> ; faster than using n! / (r! (n-r)!) >>> (define (combinations n r) >>> (cond >>> [(or (< n 0) (< r 0)) (error "combinations: illegal arguments, n and >>> r must be >= 0")] >>> [(> r n) 0] >>> [else >>> (let lp ([mord n] [total 1] [mult #t]) >>> (cond >>> [(or (= 0 mord) (= 1 mord)) total] >>> [(and mult (= mord (- n r))) (lp r total #f)] >>> [(and mult (= mord r)) (lp (- n r) total #f)] >>> [mult (lp (sub1 mord) (* total mord) #t)] >>> [else (lp (sub1 mord) (/ total mord) #f)]))])) >>> >>> Thanks again! >>> -Joe >>> >>> ____________________ >>> Racket Users list: >>> http://lists.racket-lang.org/users >>> >>> >> >
____________________ Racket Users list: http://lists.racket-lang.org/users
