Since ({: * 2 ^ [EMAIL PROTECTED]) huo y-1 is the same as |.huo y-1 ,
dmr can be simplified as follows:
dmr=: 3 : 0 " 0
if. 4>y do. y e. 2 3 return. end.
e=. huo y-1
for_a. witnesses y do. if. (+./c=y-1) +: 1={:c=. a y&|@^ e do. 0 return. end.
end.
1
)
----- Original Message -----
From: Roger Hui <[EMAIL PROTECTED]>
Date: Sunday, November 2, 2008 8:24
Subject: Re: [Jprogramming] Miller-Rabin primality test
To: Programming forum <[email protected]>
> 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
