[racket] Math library kudos
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/20130219/dd70d124/attachment.html>