Rosetta Code

Miller–Rabin primality test: Difference between revisions

Page Discussion

Miller–Rabin primality test (view source)

Revision as of 00:16, 21 January 2019

110 bytes added ,  5 years ago
m
→‎{{header|REXX}}: changed word in the REXX section header, added whitespace and optimization, added option (random seed) to allow for repeatable Miller-Rabin primality tests.
m (Reordered properly)
m (→‎{{header|REXX}}: changed word in the REXX section header, added whitespace and optimization, added option (random seed) to allow for repeatable Miller-Rabin primality tests.)
Line 4,040:
The   '''K'''   (above) is signified by the REXX variable   '''times'''   in the REXX program below.
 
<br>To make the program small, the true &nbsp; ''true prime generator'' &nbsp; ('''GenP''') &nbsp; was coded to be small, but not particularly fast.
<lang rexx>/*REXX program puts the Miller─Rabin primality test through its paces. */
parse arg limit times seed . /*obtain optional arguments from the CL*/
if limit=='' | limit=="," then limit=1000 1000 /*Not specified? Then use the default.*/
if times=='' | times=="," then times= 10 10 /* " " " " " " */
if datatype(seed, 'W') then call random ,,seed /*If seed specified, use it for RANDOM.*/
numeric digits max(200, 2*limit) /*we're dealing with some ginormous #s.*/
tell= times<0 /*display primes only if times is neg.*/
times= abs(times); w= length(times) /*use absolute value of TIMES; get len.*/
call genP limit /*suspenders now, use a belt later ··· */
@MR= 'Miller─Rabin primality test' /*define a character literal for SAY. */
Line 4,053 ⟶ 4,054:
say /* [↓] (skipping unity); show sep line*/
do a=2 to times; say copies('─', 89) /*(skipping unity) do range of TIMEs.*/
p= 0 /*the counter of primes for this pass. */
do z=1 for limit /*now, let's get busy and crank primes.*/
if \M_Rt(z, a) then iterate /*Not prime? Then try another number.*/
p= p + 1 /*well, we found another one, by gum! */
if tell then say z 'is prime according to' @MR "with K="a
if !.z then iterate
say '[K='a"] " z "isn't prime !" /*oopsy─doopsy and/or whoopsy─daisy !*/
Line 4,065 ⟶ 4,066:
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
genP: parse arg high; @.=0; !.=0; @.1=2; @.2=3; !.=@.; !.2=1; !.3=1; #=2
do j=@.#+2 by 2 to high /*just examine odd integers from here. */
do k=2 while k*k<=j; if j//@.k==0 then iterate j; end /*k*/
#= # + 1; @.#= j; !.j=1 1 /*bump prime counter; add prime to the */
end /*j*/; return /*@. array; define a prime in !. array.*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
M_Rt: procedure; parse arg n,k; d= n-1; nL=d /*Miller─Rabin: A.K.A. Rabin─Miller.*/
if n==2 then return 1 /*special case of (the) even prime. */
if n<2 | n//2==0 then return 0 /*check for too low, or an even number.*/
s=-1; do while d//2==0; d=d%2; s=s+1 /*keep halving until a zero remainder.*/
end /*while*/▼
 
do s=-1 while do kd//2==0; ?d=random(2, nL) d % 2 /*keep [↓] perform the DO loophalving until Ka zero timesremainder.*/
▲ end /*while*/
x=?**d // n /*X can get real gihugeic really fast.*/▼
 
if x==1 | x==nL then iterate /*First or penultimate? Try another pow*/
do sk; x=x**2 // n ?= random(2, nL) /*compute new[↓] Xperform the DO ≡loop X²K modulus Ntimes. */
if x==1 ?**d // n then return 0 /*ifX unity,can get it'sreal definitelygihugeic notreally primefast.*/
if x==1 | x==nL then leave iterate /*if N-1, then it couldFirst beor prime.penultimate? Try another pow*/
end do /*r*/s; x= x**2 // n /*compute new X /* [↑]≡ // isX² modulus REXX's divisionN. remainder*/
if x\==nL1 then return 0 /*nopeif unity, it ainit'ts primedefinitely nohows,not nowayprime. */
end if /*k*/x==nL then leave /*if N-1, /*maybethen it's prime, maybecould it ain'tbe ···prime. */
▲ x=?**d // n end /*r*/ /*X [↑] can get// realis gihugeicREXX's reallydivision fast.remainder*/
if x\==nL then return 0 /*nope, it ain't prime nohows, noway. */
end /*k*/ /*maybe it's prime, maybe it ain't ··· */
return 1 /*coulda/woulda/shoulda be prime; yup.*/</lang>
{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 10000 &nbsp; 10 </tt>}}
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/107959"