Jump to content

Miller–Rabin primality test: Difference between revisions

Miller–Rabin primality test (view source)

Revision as of 05:14, 11 May 2014

250 bytes added ,  10 years ago
m
→‎{{header|REXX}}: added/changed comments, optimized the SUSPENDERS subroutine.
m (→‎{{header|REXX}}: added/changed comments, optimized the SUSPENDERS subroutine.)
Line 2,246:
 
This would be in the same vein as: 3 is prime, 5 is prime, 7 is prime, all odd numbers are prime.
<lang rexx>/*REXX program puts the Miller-Rabin primality test through its paces. */
parse arg limit accur . /*get some optional arguments (if any). */▼
 
if limit=='' | limit==',' then limit=1000 /*maybe assume LIMIT default.*/
▲arg limit accur . /*get some arguments (if any). */
if limitaccur=='' | limitaccur==',' then limitaccur=100010 /*maybe assume LIMIT" " ACCUR " default*/
numeric digits max(200, 2*limit) /*we're dealing with some biggies*/▼
if accur=='' | accur==',' then accur=10 /* " " ACCUR " */
tell= accur<0 /*show primes if K is negative.*/
▲numeric digits max(200,2*limit) /*we're dealing with some biggies*/
tellaccur=abs(accur<0 ) /*shownow, primesuse if absolute Kvalue of is negativeK.*/
accur=abs(accur) call suspenders limit /*suspenders now, make K positive. belt later ··· */
call suspenders /*suspenders now, belt later... */
primePi=# /*save the count of (real) primes*/
say "There are" primePi 'primes ≤' limit /*might as well crow a wee bit.*/
say /*nothing wrong with whitespace. */
do a=2 to accur /*(skipping 1) do range of K's.*/
say copies('─',79) /*show separator for the eyeballs*/
mrp=0 /*prime counter for this pass. */
do z=1 for limit /*now, let's get busy and crank. */▼
 
▲ do z=1 for limit /*now, let's get busy and crank. */
p=Miller_Rabin(z,a) /*invoke and pray... */
if p==0 then iterate /*Not prime? Then try another. */
mrp=mrp+1 /*well, found another one, by gum*/
if tell then say z, /*maybe should do a show & tell ?*/▼
 
▲ if tell then say z, /*maybe should do a show & tell ?*/
'is prime according to Miller-Rabin primality test with K='a
if !.z\==0 then iterate▼
 
▲ if !.z\==0 then iterate
say '[K='a"] " z "isn't prime !" /*oopsy-doopsy & whoopsy-daisy!*/
end /*z*/
 
say 'for 1──►'limit", K="a', Miller-Rabin primality test found' mrp,
'primes {out of' primePi"}"
Line 2,281 ⟶ 2,276:
/*─────────────────────────────────────Rabin─Miller (also known as)─────*/
Miller_Rabin: procedure; parse arg n,k
if n==2 then return 1 then return 1 /*special case of an even prime. */
if n<2 | n//2==0 then return 0 /*check for low, or even number.*/
d=n-1 /*just a temp variable for speed.*/
d=n-1
nL=n-1 /*saves a bit of time, down below*/
s=0 /*teeter-toter variable (see-saw)*/
s=0
do while d//2==0; d=d%2; s=s+1; end /*while d//2==0*/
 
do while d//2==0; d=d%2; do s=s+1k; end /*while d//a=random(2==0 */,nL)
x=(a**d) // n /*this number can get big fast. */
 
if x==1 | x==nL then iterate /*if so, try another power of A. */
do k▼
do r=1 for s-1; x=(x*x)//n /*compute new X ≡ X² modulus N.*/
a=random(2,nL)
if x=(a**d) // n =1 then return 0 /*thisit's numberdefinitely cannot prime. get big fast. */
if x==1nL | x==nL then iterateleave
▲ do kend /*r*/
 
if x\==nL do r=1then return 0 /*nope, it ain't prime nohows. for s-1*/
end x=(x*x) /*k*/ n
if x==1 then return 0 /*it's definitely not prime. */
if x==nL then leave
end /*r*/
 
if x\==nL then return 0 /*nope, it ain't prime nohows. */
end /*k*/
/*maybe it is, maybe it ain't ...*/
return 1 /*coulda/woulda/shoulda be prime.*/
/*──────────────────────────────────SUSPENDERS subroutine───────────────*/
suspenders: parse arg high; @.=0; !.=0 /*crank up the ole prime factory.*/
@.1_='2; @.2=3; @.3=5; 7 11 #=313 17 19 23 29 31 37 41 43 47 53 59 61 67 /*prime71 the73 pump79 with83 low89 primes.*/97'
do #=1 for words(_); p=word(_,#); @.#=p; s.#=p*p; !.p=1; end /*#*/
!.2=1; !.3=1; !.5=1 /*and don't forget the water jar.*/
#=#-1 /*adjust the # of primes so far. */
 
do j =@.#+2 by 2 to limithigh /*just process the odd integers. */
if j//3==0 do k=2 while @.k**2<=j then iterate /*let's dodivisible theby olethree? primalityYes, test¬prime*/
if right(j//@.k,1)==05 then iterate j /*the Greek way, in days" " five? " " of yore.*/
do end /*k*/ =4 while s.k<=j /*a useless comment,let's butdo hey!!the ole primality test*/
#=#+1 if j//@.k==0 then iterate j /*bump the prime counter. Greek way, in days of yore.*/
@.#=j end /*k*/ /*keepa priminguseless thecomment, primebut pump.hey!! */
#=#+1; @.#=j; s.#=j*j; !.j=1 /*andbump keepprime fillingctr, theprime, water jar.primeidx*/
end /*j*/ /*this comment not left blank. */
return /*whew! All done with the primes*/</lang>
{{out|Output when using the input of: <tt>10000 10</tt>}}
<pre style="height:30ex;overflow:scroll">
There are 1229 primes ≤ 10000
 
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/107869"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information