[racket] Math library kudos

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 <jgilray at 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 <lvilnis at 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 <jgilray at 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
>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.racket-lang.org/users/archive/attachments/20130220/a091d32b/attachment-0001.html>