<div dir="ltr">Hi Jens,<div><br></div><div style>You are probably right (no pun intended!), but I was only looking at the docs and since prime-strong-pseudo-single? is unexported and not in the docs, I didn't see it.</div> <div style><br></div><div style>I may not fully understand your code, but where do you put in the number of trials (it is not in the signature of prime-strong-pseudo-single?)?</div><div style><br></div><div style>Thanks again for a very useful set of functions!</div> <div style>-joe</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Feb 20, 2013 at 11:24 AM, Jens Axel Søgaard <span dir="ltr"><<a href="mailto:jensaxel@soegaard.net" target="_blank">jensaxel@soegaard.net</a>></span> wrote:<br> <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Joe,<br> <br> 2013/2/20 Joe Gilray <<a href="mailto:jgilray@gmail.com">jgilray@gmail.com</a>>:<br> <div class="im">> Racketeers,<br> ><br> > Thanks for putting together the fantastic math library. It will be a<br> > wonderful resource. Here are some quick impressions (after playing mostly<br> > with math/number-theory)<br> ><br> > 1) The functions passed all my tests and were very fast. If you need even<br> > more speed you can keep a list of primes around and write functions to use<br> > that, but that should be rarely necessary<br> <br> </div>That's great to hear.<br> <div class="im"><br> > 2) I have a couple of functions to donate if you want them:<br> <br> </div>Contributions to the library is very welcome. For example<br> it would be nice to have a better implementation of next-prime.<br> <div><div class="h5"><br> > 2a) Probablistic primality test:<br> ><br> > ; function that performs a Miller-Rabin probabalistic primality test k<br> > times, returns #t if n is probably prime<br> > ; algorithm from <a href="http://rosettacode.org/wiki/Miller-Rabin_primality_test" target="_blank">http://rosettacode.org/wiki/Miller-Rabin_primality_test</a>,<br> > code adapted from Lisp example<br> > ; (module+ test (check-equal? (is-mr-prime? 1000000000000037 8) #t))<br> > (define (is-mr-prime? n k)<br> > ; function that returns two values r and e such that number = divisor^e *<br> > r, and r is not divisible by divisor<br> > (define (factor-out number divisor)<br> > (do ([e 0 (add1 e)] [r number (/ r divisor)])<br> > ((not (zero? (remainder r divisor))) (values r e))))<br> ><br> > ; function that performs fast modular exponentiation by repeated squaring<br> > (define (expt-mod base exponent modulus)<br> > (let expt-mod-iter ([b base] [e exponent] [p 1])<br> > (cond<br> > [(zero? e) p]<br> > [(even? e) (expt-mod-iter (modulo (* b b) modulus) (/ e 2) p)]<br> > [else (expt-mod-iter b (sub1 e) (modulo (* b p) modulus))])))<br> ><br> > ; function to return a random, exact number in the passed range<br> > (inclusive)<br> > (define (shifted-rand lower upper)<br> > (+ lower (random (add1 (- (modulo upper 4294967088) (modulo lower<br> > 4294967088))))))<br> ><br> > (cond<br> > [(= n 1) #f]<br> > [(< n 4) #t]<br> > [(even? n) #f]<br> > [else<br> > (let-values ([(d s) (factor-out (- n 1) 2)]) ; represent n-1 as 2^s-d<br> > (let lp ([a (shifted-rand 2 (- n 2))] [cnt k])<br> > (if (zero? cnt) #t<br> > (let ([x (expt-mod a d n)])<br> > (if (or (= x 1) (= x (sub1 n))) (lp (shifted-rand 2 (- n 2))<br> > (sub1 cnt))<br> > (let ctestlp ([r 1] [ctest (modulo (* x x) n)])<br> > (cond<br> > [(>= r s) #f]<br> > [(= ctest 1) #f]<br> > [(= ctest (sub1 n)) (lp (shifted-rand 2 (- n 2))<br> > (sub1 cnt))]<br> > [else (ctestlp (add1 r) (modulo (* ctest ctest)<br> > n))])))))))]))<br> <br> </div></div>As far as I can tell, this is the same algorithm as prime-strong-pseudo-single?<br> from number-theory.rkt. See line 134.<br> <br> <a href="https://github.com/plt/racket/blob/master/collects/math/private/number-theory/number-theory.rkt" target="_blank">https://github.com/plt/racket/blob/master/collects/math/private/number-theory/number-theory.rkt</a><br> <br> --<br> Jens Axel Søgaard<br> </blockquote></div><br></div>