[racket] Math library kudos
Hi Joe,
2013/2/20 Joe Gilray <jgilray at gmail.com>:
> 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
That's great to hear.
> 2) I have a couple of functions to donate if you want them:
Contributions to the library is very welcome. For example
it would be nice to have a better implementation of next-prime.
> 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))])))))))]))
As far as I can tell, this is the same algorithm as prime-strong-pseudo-single?
from number-theory.rkt. See line 134.
https://github.com/plt/racket/blob/master/collects/math/private/number-theory/number-theory.rkt
--
Jens Axel Søgaard