Miller–Rabin primality test: Difference between revisions

Content added Content deleted
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: Line 2,893:


=={{header|REXX}}==
=={{header|REXX}}==
With a   K   of   1,   there seems to be a not uncommon number of failures, but
With a   '''K'''   of   '''1''',   there seems to be a not uncommon number of failures, but
::: with a   K ≥ 2,   the failures are rare,
::: with a   '''K ≥ 2''',   the failures are seldom,
::: with a   K ≥ 3,   rare as hen's teeth.
::: with a   '''K ≥ 3''',   the failures are rare as hen's teeth.


This would be in the same vein as:
This would be in the same vein as:
3 is prime, 5 is prime, 7 is prime, all odd numbers are prime.
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─Rabin primality test through its paces. */
if limit=='' | limit==',' then limit=1000 /*maybe assume LIMIT default.*/
if accur=='' | accur==',' then accur=10 /* " " ACCUR " */
parse arg limit accur . /*obtain optional arguments from the CL*/
numeric digits max(200, 2*limit) /*we're dealing with some biggies*/
if limit=='' | limit==',' then limit=1000 /*Not specified? Then use the default.*/
tell= accur<0 /*show primes if K is negative.*/
if accur=='' | accur==',' then accur= 10 /* " " " " " " */
accur=abs(accur) /*now, use absolute value of K.*/
numeric digits max(200, 2*limit) /*we're dealing with some ginormous #s.*/
call suspenders limit /*suspenders now, belt later ··· */
tell= accur<0 /*show primes only if K is negative.*/
primePi=# /*save the count of (real) primes*/
accur=abs(accur) /*now, use the absolute value of K. */
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 /*might as well display some stuff. */
do a=2 to accur /*(skipping 1) do range of K's.*/
say copies('─',79) /*show separator for the eyeballs*/
say /*nothing wrong with some whitespace. */
mrp=0 /*prime counter for this pass. */
do a=2 to accur /*(skipping unity) do range of K's. */
do z=1 for limit /*now, let's get busy and crank. */
say copies('─', 89) /*show separator for the ole eyeballs. */
p=Miller_Rabin(z,a) /*invoke and pray... */
mrp=0 /*counter of primes for this pass. */
if p==0 then iterate /*Not prime? Then try another. */
do z=1 for limit /*now, let's get busy and crank primes.*/
mrp=mrp+1 /*well, found another one, by gum*/
p=Miller_Rabin(z, a) /*invoke and hope for the best ··· */
if tell then say z, /*maybe should do a show & tell ?*/
if p==0 then iterate /*Not prime? Then try another number.*/
⚫
mrp=mrp+1 /*well, found another one, by gum! */
'is prime according to Miller-Rabin primality test with K='a
if !.z\==0 then iterate
if tell then say z 'is prime according to' @MRpt "with K="a
if !.z then iterate
say '[K='a"] " z "isn't prime !" /*oopsy-doopsy & whoopsy-daisy!*/
say '[K='a"] " z "isn't prime !" /*oopsy─doopsy and/or whoopsy─daisy !*/
end /*z*/
end /*z*/
say 'for 1──►'limit", K="a', Miller-Rabin primality test found' mrp,
say ' for 1──►'limit", K="right(a,length(accur))',',
'primes {out of' primePi"}"
@MRpt "found" mrp 'primes {out of' primePi"}."
end /*a*/
end /*a*/
exit /*stick a fork in it, we're done.*/
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)─────*/
⚫
if n==2 then return 1 /*special case of (the) even prime. */
⚫
Miller_Rabin: procedure; parse arg n,k
if n==2 then return 1 /*special case of an even prime. */
if n<2 | n//2==0 then return 0 /*check for too low, or an even number.*/
if n<2 | n//2==0 then return 0 /*check for low, or even number.*/
d=n-1; nL=d /*just two temp variables, for speed. */
d=n-1 /*just a temp variable for speed.*/
s=0 /*a teeter─toter variable (see─saw). */
nL=n-1 /*saves a bit of time, down below*/
do while d//2==0; d=d%2; s=s+1; end /*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)
do k; ?=random(2, nL) /* [↓] perform the DO loop K times.*/
x=(a**d) // n /*this number can get big fast. */
x=(?**d) // n /*X can get real gihugeic really fast.*/
if x==1 | x==nL then iterate /*if so, try another power of A. */
if x==1 | x==nL then iterate /*First or penultimate? Try another pow*/
do r=1 for s-1; x=(x*x)//n /*compute new X ≡ X² modulus N.*/
do r=1 for s-1; x=(x*x)//n /*compute new X ≡ X² modulus N. */
if x==1 then return 0 /*it's definitely not prime. */
if x==1 then return 0 /*if unity, it's definitely not prime.*/
if x==nL then leave
if x==nL then leave /*if N-1, then it could be prime. */
end /*r*/
end /*r*/ /* // is really REXX's ÷ remainder.*/
if x\==nL then return 0 /*nope, it ain't prime nohows. */
if x\==nL then return 0 /*nope, it ain't prime nohows, noway. */
end /*k*/
end /*k*/ /*maybe it's prime, maybe it ain't ··· */
/*maybe it is, maybe it ain't ...*/
return 1 /*coulda/woulda/shoulda be prime; 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───────────────*/
⚫
do j=@.#+2 by 2 to high /*just examine odd integers from 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; end /*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(_,#); @.#=p; s.#=p*p; !.p=1; end /*#*/
#=#+1; @.#=j; !.j=1 /*bump the prime counter; add prime to */
#=#-1 /*adjust the # of primes so far. */
end /*j*/ /* ··· array; define a prime in array·*/
do j=@.#+2 by 2 to high /*just process the odd integers. */
return /* [↓] generate 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>
'''output''' &nbsp; when using the input of: &nbsp; <tt> 10000 &nbsp; 10 </tt>
<pre>
<pre>
There are 1229 primes ≤ 10000
There are 1229 primes ≤ 10000


─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
[K=2] 2701 isn't prime !
[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 1230 primes {out of 1229}.
─────────────────────────────────────────────────────────────────────────────────────────
───────────────────────────────────────────────────────────────────────────────
for 1──►10000, K=2, Miller─Rabin primality test found 1229 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= 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= 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= 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= 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= 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= 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= 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}
for 1──►10000, K=10, Miller─Rabin primality test found 1229 primes {out of 1229}.
</pre>
</pre>


Retrieved from "https://rosettacode.org/wiki/Miller–Rabin_primality_test"