Rosetta Code

Miller–Rabin primality test: Difference between revisions

Page Discussion

Miller–Rabin primality test (view source)

Revision as of 09:05, 12 March 2016

796 bytes added ,  8 years ago
m
→‎{{header|REXX}}: added/changed whitespace and comments, change indentations in program and output, simplified the genPrimes subroutine..
m (→‎{{header|Run BASIC}}: Fix opening comment)
m (→‎{{header|REXX}}: added/changed whitespace and comments, change indentations in program and output, simplified the genPrimes subroutine..)
Line 2,893:
 
=={{header|REXX}}==
With a   '''K'''   of   '''1''',   there seems to be a not uncommon number of failures, but
::: with a   '''K ≥ 2''',   the failures are rareseldom,
::: with a   '''K ≥ 3''',   the failures are rare as hen's teeth.
 
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.*/▼
<br>To make the program small, the true &nbsp; prime generator &nbsp; was coded to be small, but not particularly fast.
parse arg limit accur . /*get some optional arguments*/
▲<lang rexx>/*REXX program puts the Miller-RabinMiller─Rabin primality test through its paces. */
if limit=='' | limit==',' then limit=1000 /*maybe assume LIMIT default.*/
ifparse accur==''arg |limit accur==',' . then accur=10 /* " " ACCUR " /*obtain optional arguments from the CL*/
numericif digitslimit=='' max(200| limit==',' then 2*limit)=1000 /*we'reNot specified? dealingThen withuse somethe biggiesdefault.*/
tellif accur=='' | accur<0==',' then accur= 10 /* " " /*show primes if" " " " K is negative.*/
accur=abs(accur) numeric digits max(200, 2*limit) /*now, use we're absolutedealing valuewith ofsome ginormous K#s.*/
calltell= accur<0 suspenders limit /*suspendersshow primes only if now, beltK later ···is negative.*/
primePiaccur=#abs(accur) /*savenow, use the count absolute value of (real) primesK. */
call genPrimes limit /*suspenders now, use a belt later ··· */
say "There are" primePi 'primes ≤' limit /*might as well crow a wee bit*/▼
primePi=#; @MRpt='Miller─Rabin primality test' /*save count of actual primes; literal.*/
say /*nothing wrong with whitespace. */
▲say "There are" primePi 'primes ≤' limit limit /*might as well crowdisplay some stuff. a wee bit*/
do a=2 to accur /*(skipping 1) do range of K's.*/
say say copies('─',79) /*shownothing wrong with separatorsome forwhitespace. the eyeballs*/
mrpdo a=02 to accur /*prime(skipping counterunity) for this pass.do range of K's. */
say copies('─', 89) do z=1 for limit /*now,show let'sseparator getfor busythe andole crankeyeballs. */
mrp=0 p=Miller_Rabin(z,a) /*invoke and pray... /*counter of primes for this pass. */
if p= do z=01 for limit then iterate /*Notnow, prime?let's get busy Thenand trycrank anotherprimes. */
mrp=mrp+1 p=Miller_Rabin(z, a) /*well,invoke and hope for the best ··· found another one, by gum*/
if tellp==0 then sayiterate z, /*maybeNot shouldprime? do a showThen &try tellanother ?number.*/
end /*k*/ mrp=mrp+1 /*awell, found uselessanother commentone, butby hey!gum! */▼
'is prime according to Miller-Rabin primality test with K='a
if !.z\==0tell then iteratesay z 'is prime according to' @MRpt "with K="a
if !.z then iterate
say '[K='a"] " z "isn't prime !" /*oopsy-doopsyoopsy─doopsy and/or &whoopsy─daisy whoopsy-daisy!*/
end /*z*/
say 'for 1──►'limit", K="a', Miller-Rabin primality test found' mrp,
say ' for 1──►'limit", K="right(a,length(accur))',',
@MRpt "found" mrp 'primes {out of' primePi"}."
end /*a*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
/*─────────────────────────────────────Miller─Rabin primality test.─────*/
Miller_Rabin: procedure; parse arg n,k /*Miller─Rabin: A.K.A. Rabin─Miller.*/▼
/*─────────────────────────────────────Rabin─Miller (also known as)─────*/
return 1 if n==2 then return 1 /*coulda/woulda/shouldaspecial case of (the) even be prime. */▼
▲Miller_Rabin: procedure; parse arg n,k
if n==2 then return 1 if n<2 | n//2==0 then return 0 /*specialcheck casefor oftoo low, or an even primenumber. */
if n<2 | d=n//2=-1; nL=0d then return 0 /*checkjust fortwo lowtemp variables, orfor even numberspeed. */
d=n-1 s=0 /*just a tempteeter─toter variable for speed(see─saw). */
nL=n-1 do while d/*saves/2==0; a bitd=d%2; of time,s=s+1; down belowend /*while ···*/
s=0 /*teeter-toter variable (see-saw)*/
do while d//2==0; d=d%2; s=s+1; end /*while d//2==0*/
 
do k; a ?=random(2, nL) /* [↓] perform the DO loop K times.*/
x=(a?**d) // n /*thisX number can get bigreal fast.gihugeic really fast.*/
if x==1 | x==nL then iterate /*ifFirst so,or trypenultimate? anotherTry power of A.another pow*/
do r=1 for s-1; x=(x*x)//n /*compute new X ≡ X² modulus N. */
if x==1 then return 0 /*if unity, it's definitely not prime. */
if x==nL then leave /*if N-1, then it could be prime. */
end /*r*/ /* // is really REXX's ÷ remainder.*/
if x\==nL then return 0 /*nope, it ain't prime nohows., noway. */
end end /*k*/ /*maybe it's prime, maybe it ain't ··· */
return 1 /*maybe it is, maybe/*coulda/woulda/shoulda be itprime; ain't ..yup.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
▲return 1 /*coulda/woulda/shoulda be prime.*/
genPrimes: parse arg high; @.=0; !.=0; @.1=2; @.2=3; !.2=1; !.3=1; #=2
/*──────────────────────────────────SUSPENDERS subroutine───────────────*/
end /*j*/ do j=@.#+2 by 2 to high /*thisjust commentexamine notodd leftintegers blankfrom here. */▼
suspenders: parse arg high; @.=0; !.=0 /*crank up the ole prime factory.*/
do k=2 while k*k<=j; if j//@.k==0 then iterate j; /*the Greek way, in daysend of yore./*k*/▼
_='2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97'
do #=1 for words(_); p #=word(_,#)+1; @.#=pj; s.#=p*p; !.pj=1; end /*#bump the prime counter; add prime to */
#=#-1 end /*j*/ /*adjust the··· array; # ofdefine primesa soprime far.in array·*/
do j=@.#+2 by 2return to high /*just process[↓] the oddgenerate integers.primes (not optimal).*//</lang>
if j//3==0 then iterate /*divisible by three? Yes, ¬prime*/
if right(j,1)==5 then iterate /* " " five? " " */
do k=4 while s.k<=j /*let's do the ole primality test*/
▲ if j//@.k==0 then iterate j /*the Greek way, in days of yore.*/
▲ end /*k*/ /*a useless comment, but hey!! */
#=#+1; @.#=j; s.#=j*j; !.j=1 /*bump prime ctr, prime, primeidx*/
▲ end /*j*/ /*this comment not left blank. */
return /*whew! All done with the primes*/</lang>
'''output''' &nbsp; when using the input of: &nbsp; <tt> 10000 &nbsp; 10 </tt>
<pre>
There are 1229 primes ≤ 10000
 
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
[K=2] 2701 isn't prime !
for 1──►10000, K= 2, Miller─Rabin primality test found 1230 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K= 2, Miller─Rabin primality test found 1229 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K= 3, Miller─Rabin primality test found 1229 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K= 4, Miller─Rabin primality test found 1229 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K= 5, Miller─Rabin primality test found 1229 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K= 6, Miller─Rabin primality test found 1229 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K= 7, Miller─Rabin primality test found 1229 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K= 8, Miller─Rabin primality test found 1229 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K= 9, Miller─Rabin primality test found 1229 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=10, Miller─Rabin primality test found 1229 primes {out of 1229}.
</pre>
 
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/107909"