The following is a concise implementation of
the Miller-Rabin probabilistic primality test.
m MillerRabin n is 0 or 1. If 0, then n is
not prime; if 1, then n is prime with an error
probability of 0.25^m .
huo=: <[EMAIL PROTECTED]:^:(0=2&|)^:a: NB. halve until odd
MillerRabin1=: 100&$: : (4 : 0) " 0
*./1=((1+x [EMAIL PROTECTED] y-1) y&|@^ {:s) y&|@^ 2x^#s=. huo y-1
)
A more verbose but more efficient implementation
is as follows:
MillerRabin=: 100&$: : (4 : 0) " 0
d=. {:s=. huo y-1
e=. 2x^#s
for_a. 1+x [EMAIL PROTECTED] y-1 do.
if. 1~:(a y&|@^ d) y&|@^ e do. 0 return. end.
end.
1
)
(MillerRabin -: 1&p:) 1000 [EMAIL PROTECTED] 2e9
1
MillerRabin n=: */ p: 1e8+0 100x
0
n
4153752848336204041
MillerRabin 4 p: n
1
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
