http://www.jsoftware.com/jwiki/PascalJasmin/Miller%20Rabin%20related%20code
with some updates. There were bugs in the code I posted in message. ----- Original Message ----- From: Raul Miller <[email protected]> To: Programming forum <[email protected]> Cc: Sent: Tuesday, February 18, 2014 11:39:52 PM Subject: Re: [Jprogramming] some Miller Rabin code This looks fun, and maybe even impressive. I want to be able to come back to it later, when I feel comfortable focusing on primality issues - can I convince you to put it on a J wiki page? (This does not have to be today - if you are still working through your ideas maybe next week would have more perspectives - but this looks promising.) Thanks, -- Raul On Tue, Feb 18, 2014 at 11:30 PM, Pascal Jasmin <[email protected]>wrote: > based on Roger's code (MillerRabin1) here are a couple of variations: > > huo=: <. @ -:^:(0=2&|)^:a: NB. halve until odd > NB. magic numbers from http://primes.utm.edu/prove/prove2_3.html > witnesses=: 4 : 0 > r=. 9080191 4759123141 2152302898747 3474749600383 341550071728321x > if. y>:{:r do. > 1+x ?@$ y-1 > else. > ((r-1) I. y){:: 31 73 ; 2 7 61 ; 2 3 5 7 11 ; 2 3 5 7 11 13 ; 2 3 5 7 11 > 13 17 > end. > ) > > NB. Miller-Rabin primality test > NB. deterministic for y< {:r (3.4e14) in witnesses; probabilistic for > y>:{:r (0 answer is always correct; > NB. 1 answer is wrong with error probability at most 0.25^x) > > MillerRabin1=: 100&$: : (4 : 0) " 0 > if. 0=2|y do. 2=y return. end. > if. 74>y do. y e. i.&.(p:^:_1) 74 return. end. > e=. huo y-1 > for_a. x witnesses y do. if. (+./c=y-1) +: 1={:c=. a y&|@^ e do. 0 > return. end. end. > 1 > ) > MillerRabin=: 100&$: : (4 : 0) " 0 > e=. huo y-1 > for_a. x witnesses y do. if. (+./c=y-1) +: 1={:c=. a y&|@^ e do. 0 > return. end. end. > 1 > ) > MillerRabinX=: 4&$: : (4 : 0) " 0 > e=. huo y-1 > for_a. x witnesses y do. if. (c=y-1) +: 1=c=. a y&|@^ {:e do. > if. (+./c=y-1) +: +./ 1=c y&|@^ 2 ^ >:i.<:#e do. 0 return. end. end. > end. > 1 > ) > > > MillerRabinX uses NIST's fips 186-3 algorithm, and appears faster at > finding composites. NIST also recommends much fewer rounds (x parameter) > than the original algorithm, though I seem to have a hard time finding > pseudoprimes that fail few rounds, and pass many rounds. And when I do > find them, I forget to write them down :P > > Anyways, this code is sometimes faster than the builtin p: command (for > medium to large y). Both 1 p: and 4 p: replacements are provided. > > maybePrimeMod210 =: 211 ,~ 2 3 5 7 (] #~ [: -. [: +./ 0 = |"0 1) 1+i.210 > quicknextcandidate =: ((- 210&|@:]) + maybePrimeMod210&([ {~ [ > I.`(>:@:I.)@.(e.~) 210&|@:])"1 ) > nextprime =: quicknextcandidate^:(-.@(1&p:))^:_@:quicknextcandidate > nextprime2 =: > quicknextcandidate^:(-.@(1&MillerRabin))^:_@:quicknextcandidate > nextprime3 =: > quicknextcandidate^:(-.@(3&MillerRabinX))^:_@:quicknextcandidate > nextprime4 =: > quicknextcandidate^:(-.@(1&MillerRabinX))^:_@:quicknextcandidate > nextprime40 =: > quicknextcandidate^:(-.@(40&MillerRabin))^:_@:quicknextcandidate > > > ts 'nextprime^:25 ] 35123443234232312332154433466451x' > 0.311963/sec 0.234368MB > ts 'nextprime4^:25 ] 35123443234232312332154433466451x' > 1.3341/sec 0.204416MB > ts 'nextprime2^:25 ] 35123443234232312332154433466451x' > 0.460181/sec 0.206208MB > > above tests finding 25 next primes from the same starting "random" number > nextprime4 is both faster than p: and supposedly more accurate. > > ts'nextprime^:50 ] 35123443234232312332451x' > 0.313615/sec 0.143744MB > ts'4 p:^:50 ] 35123443234232312332451x' > 0.334359/sec 0.147456MB > > ts'nextprime4^:50 ] 35123443234232312332451x' > 2.53272/sec 0.122752MB > ts'nextprime40^:50 ] 35123443234232312332451x' > 0.200787/sec 0.131456MB > ts'nextprime3^:50 ] 35123443234232312332451x' > 2.17427/sec 0.12352MB > > nextprime4 and 3 (both using MillerRabinX) get even faster than p: as y > is lowered > > Perhaps part of the speed increase is due to maybePrimeMod210, which is a > list of 48 possible results that a 210 | candidate could have and be prime > (simultaneously tests for divisibility by 2 3 5 7) > > quicknextcandidate is super fast if you only care about finding the next > number that is not divisible by 2 3 5 7. > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
