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:<div>
<br></div><div><font face="courier new, monospace">(exp (fllog-gamma (+ 1.0 n)) - (fllog-gamma (+ 1.0 r)) - (fllog-gamma (+ 1.0 (- n r))))</font></div><div><br></div><div><font face="courier new, monospace">fllog-gamma</font> should also be faster than <font face="courier new, monospace">bflog-gamma</font> or <font face="courier new, monospace">log-gamma</font> 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.<br>
<br>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.<br><br><div class="gmail_quote">On Wed, Feb 20, 2013 at 1:44 AM, Joe Gilray <span dir="ltr"><<a href="mailto:jgilray@gmail.com" target="_blank">jgilray@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi Luke,<div><br></div><div>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?</div>
<span class="HOEnZb"><font color="#888888">
<div><br></div><div>-Joe</div></font></span></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, Feb 19, 2013 at 8:26 PM, Luke Vilnis <span dir="ltr"><<a href="mailto:lvilnis@gmail.com" target="_blank">lvilnis@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">FYI, log gamma is another fast way to calculate the number of combinations if you want to deal with really big numbers.<br>
<br><div class="gmail_quote"><div><div>On Tue, Feb 19, 2013 at 7:28 PM, Joe Gilray <span dir="ltr"><<a href="mailto:jgilray@gmail.com" target="_blank">jgilray@gmail.com</a>></span> wrote:<br>
</div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div><div dir="ltr">Racketeers,<div><br></div><div>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)<div>
<br>
</div><div>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</div><div>
<br></div><div>2) I have a couple of functions to donate if you want them:</div><div><br></div><div>2a) Probablistic primality test:</div><div><br></div><div><div><font face="courier new, monospace">; function that performs a Miller-Rabin probabalistic primality test k times, returns #t if n is probably prime</font></div>
<div><font face="courier new, monospace">; algorithm from <a href="http://rosettacode.org/wiki/Miller-Rabin_primality_test" target="_blank">http://rosettacode.org/wiki/Miller-Rabin_primality_test</a>, code adapted from Lisp example</font></div>
<div><font face="courier new, monospace">; (module+ test (check-equal? (is-mr-prime? 1000000000000037 8) #t))</font></div><div><font face="courier new, monospace">(define (is-mr-prime? n k)</font></div><div><font face="courier new, monospace"> ; function that returns two values r and e such that number = divisor^e * r, and r is not divisible by divisor</font></div>
<div><font face="courier new, monospace"> (define (factor-out number divisor)</font></div><div><font face="courier new, monospace"> (do ([e 0 (add1 e)] [r number (/ r divisor)])</font></div><div><font face="courier new, monospace"> ((not (zero? (remainder r divisor))) (values r e))))</font></div>
<div><font face="courier new, monospace"> </font></div><div><font face="courier new, monospace"> ; function that performs fast modular exponentiation by repeated squaring</font></div><div><font face="courier new, monospace"> (define (expt-mod base exponent modulus)</font></div>
<div><font face="courier new, monospace"> (let expt-mod-iter ([b base] [e exponent] [p 1])</font></div><div><font face="courier new, monospace"> (cond </font></div><div><font face="courier new, monospace"> [(zero? e) p]</font></div>
<div><font face="courier new, monospace"> [(even? e) (expt-mod-iter (modulo (* b b) modulus) (/ e 2) p)]</font></div><div><font face="courier new, monospace"> [else (expt-mod-iter b (sub1 e) (modulo (* b p) modulus))])))</font></div>
<div><font face="courier new, monospace"> </font></div><div><font face="courier new, monospace"> ; function to return a random, exact number in the passed range (inclusive)</font></div><div><font face="courier new, monospace"> (define (shifted-rand lower upper)</font></div>
<div><font face="courier new, monospace"> (+ lower (random (add1 (- (modulo upper 4294967088) (modulo lower 4294967088))))))</font></div><div><font face="courier new, monospace"> </font></div><div><font face="courier new, monospace"> (cond </font></div>
<div><font face="courier new, monospace"> [(= n 1) #f]</font></div><div><font face="courier new, monospace"> [(< n 4) #t]</font></div><div><font face="courier new, monospace"> [(even? n) #f]</font></div><div>
<font face="courier new, monospace"> [else </font></div>
<div><font face="courier new, monospace"> (let-values ([(d s) (factor-out (- n 1) 2)]) ; represent n-1 as 2^s-d</font></div><div><font face="courier new, monospace"> (let lp ([a (shifted-rand 2 (- n 2))] [cnt k])</font></div>
<div><font face="courier new, monospace"> (if (zero? cnt) #t</font></div><div><font face="courier new, monospace"> (let ([x (expt-mod a d n)])</font></div><div><font face="courier new, monospace"> (if (or (= x 1) (= x (sub1 n))) (lp (shifted-rand 2 (- n 2)) (sub1 cnt))</font></div>
<div><font face="courier new, monospace"> (let ctestlp ([r 1] [ctest (modulo (* x x) n)])</font></div><div><font face="courier new, monospace"> (cond </font></div><div><font face="courier new, monospace"> [(>= r s) #f]</font></div>
<div><font face="courier new, monospace"> [(= ctest 1) #f]</font></div><div><font face="courier new, monospace"> [(= ctest (sub1 n)) (lp (shifted-rand 2 (- n 2)) (sub1 cnt))]</font></div>
<div><font face="courier new, monospace"> [else (ctestlp (add1 r) (modulo (* ctest ctest) n))])))))))]))</font></div><div><br></div><div>2b) combinations calculator</div><div><br></div><div>
<div><font face="courier new, monospace">; function that returns the number of combinations, not the combinations themselves</font></div><div><font face="courier new, monospace">; faster than using n! / (r! (n-r)!)</font></div>
<div><font face="courier new, monospace">(define (combinations n r)</font></div><div><font face="courier new, monospace"> (cond </font></div><div><font face="courier new, monospace"> [(or (< n 0) (< r 0)) (error "combinations: illegal arguments, n and r must be >= 0")]</font></div>
<div><font face="courier new, monospace"> [(> r n) 0]</font></div><div><font face="courier new, monospace"> [else </font></div><div><font face="courier new, monospace"> (let lp ([mord n] [total 1] [mult #t])</font></div>
<div><font face="courier new, monospace"> (cond </font></div><div><font face="courier new, monospace"> [(or (= 0 mord) (= 1 mord)) total]</font></div><div><font face="courier new, monospace"> [(and mult (= mord (- n r))) (lp r total #f)]</font></div>
<div><font face="courier new, monospace"> [(and mult (= mord r)) (lp (- n r) total #f)]</font></div><div><font face="courier new, monospace"> [mult (lp (sub1 mord) (* total mord) #t)]</font></div><div><font face="courier new, monospace"> [else (lp (sub1 mord) (/ total mord) #f)]))]))</font></div>
<div><br></div><div>Thanks again!</div><span><font color="#888888"><div>-Joe</div></font></span></div></div></div></div>
<br></div></div>____________________<br>
Racket Users list:<br>
<a href="http://lists.racket-lang.org/users" target="_blank">http://lists.racket-lang.org/users</a><br>
<br></blockquote></div><br>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div>