Using the assertions from
http://primes.utm.edu/prove/prove2_3.html
witnesses=: 3 : 0
r=. 0 1373653 25326001 3215031751 2152302898747 3474749600383 341550071728321x
p: i. 1+(r-1) I. y
)
dmr=: 3 : 0 " 0
if. 4>y do. y e. 2 3 return. end.
e=. ({: * 2 ^ [EMAIL PROTECTED]) huo y-1
for_a. witnesses y do. if. (+./c=y-1) +: 1={.c=. a y&|@^ e do. 0 return. end.
end.
1
)
(dmr -: 1&p:) i. 1e5
1
dmr is good for up to 3.4e14 (the last element of r in witnesses).
----- Original Message -----
From: John Randall <[EMAIL PROTECTED]>
Date: Thursday, October 30, 2008 16:48
Subject: Re: [Jprogramming] Miller-Rabin primality test
To: Programming forum <[email protected]>
> Roger Hui wrote:
> >> There are a lot of alternatives based on the range you want.
> >
> > Some popular ranges would be
> > y<2^31
> > y<2^32
> > y<2^63
> > y<2^64
>
> With no code conservation, here are a couple of examples.
>
> time=:6!:2
> huo=: <[EMAIL PROTECTED]:^:(0=2&|)^:a: NB. halve until odd
>
> DMR=: 3 : 0"0 NB. up to 2^32
> h=.huo y-1
> s=.#h
> d=.{: h
> for_a. (#~ y&>) 2 7 61 do.
> c=.a y&|@^d*2 y&|@^ i.s
> if. (1~: {. c) *. *./ (y-1)~: c do. 0 return. end.
> end.
> 1
> )
> DMR2=: 3 : 0"0 NB. up to 2^48
> h=.huo y-1
> s=.#h
> d=.{: h
> for_a. (#~ y&>) 2 3 5 7 11 13 do.
> c=.a y&|@^d*2 y&|@^ i.s
> if. (1~: {. c) *. *./ (y-1)~: c do. 0 return. end.
> end.
> 1
> )
>
> N=:(i.1e3)+10^12x
> time 'DMR2 N'
> 1.05322
> time '1&p: N'
> 1.49582
> (DMR2-:1&p:) N
> 1
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
