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
