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 <[email protected]> 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 <[email protected]> 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
