Jump to content

Miller–Rabin primality test: Difference between revisions

→‎{{header|REXX}}: added/changed comments and whitespace, changed indentations, optimized the inner DO loop in the M_Rt function.
(→‎{{header|REXX}}: added/changed comments and whitespace, changed indentations, optimized the inner DO loop in the M_Rt function.)
Line 3,816:
if times=='' | times=="," then times= 10 /* " " " " " " */
numeric digits max(200, 2*limit) /*we're dealing with some ginormous #s.*/
tell= times<0; pad=left('', 7) /*showdisplay primes only if times is negativeneg.*/
times=abs(times); w=length(times) w=length(times) /*use absolute value of TIMES; get len.*/
call genPrimesgenP limit /*suspenders now, use a belt later ··· */
@MRptMR= 'Miller─Rabin primality test' /*define a character literal for SAY. */
say "There are" # 'primes ≤' limit limit /*might as well display some stuff. */
say /* [↓] (skipping unity); show sep line*/
do a=2 to times; say copies('─', 89) /*(skipping unity) do range of TIMEs.*/
mrpp=0 /*counter of primes for this pass. */
do z=1 for limit /*now, let's get busy and crank primes.*/
if p=Miller_Rabin\M_Rt(z, a) then iterate /*invoke and hope for theNot bestprime? ··· Then try another number.*/
p=p + 1 if p==0 then iterate /*Notwell, prime?we found another Thenone, tryby gum! another number.*/
if tell then say z mrp=mrp+1 'is prime according to' @MR "with /*well, found another one, by gum! */K="a
if !.z then end /*z*/iterate
if tell then say z 'is prime according to' @MRpt "with K="a
say '[K='a"] " ifz "isn't prime !.z" /*oopsy─doopsy then iterateand/or whoopsy─daisy !*/
end end /*az*/
say '[K='a"] " z "isn't prime !" /*oopsy─doopsy and/or whoopsy─daisy !*/
say pad ' for 1──►'limit", K="right(a,w)',' @MRptMR "found" mrp p 'primes {out of' #"}."
end /*z*/
end /*a*/
 
say pad 'for 1──►'limit", K="right(a,w)',' @MRpt "found" mrp 'primes {out of' #"}."
end /*a*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
genPrimesgenP: parse arg high; @.=0; !.=0; @.1=2; @.2=3; !.2=1; !.3=1; #=2
Miller_Rabin: procedure; parse arg n,k /*Miller─Rabin: A.K.A. Rabin─Miller.*/
do j=@.#+2 by 2 ifto n==2high then return 1 /*specialjust caseexamine ofodd (the)integers from even primehere. */
do k=2 while if nk*k<2=j; | nif j//2@.k==0 then returniterate j; 0 end /*check for too low, or an even number.k*/
d#=n-#+1; nL @.#=d j; !.j=1 /*justbump twoprime tempcounter; variables,add prime forto speed.the */
end /*j*/; s=0 return /*@. array; define /*a teeter─toterprime variablein !. (see─saw)array. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
do while d//2==0; d=d%2; s=s+1; end /*while ···*/
Miller_RabinM_Rt: procedure; parse arg n,k; d=n-1; nL=d /*Miller─Rabin: A.K.A. Rabin─Miller.*/
if n==2 #=#+1;then return 1 @.#=j; !.j=1 /*bumpspecial case of (the) prime counter;even add prime. to */
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 endhalving /*whileuntil a zero ···remainder.*/
end /*while*/
 
do k; ?=random(2, nL) /* [↓] perform the DO loop K times.*/
x=(?**d) // n /*X can get real gihugeic really fast.*/
if x==1 | x==nL then iterate /*First or penultimate? Try another pow*/
do r=1s; for s-1; x=x**2 // 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 ÷ division 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 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>
/*──────────────────────────────────────────────────────────────────────────────────────*/
genPrimes: 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 /*bump the prime counter; add prime to */
end /*j*/ /* ··· array; define a prime in array·*/
return /* [↓] generate primes (not optimal).*//</lang>
'''output''' &nbsp; when using the input of: &nbsp; <tt> 10000 &nbsp; 10 </tt>
<pre>
There are 1229 primes ≤ 10000
Cookies help us deliver our services. By using our services, you agree to our use of cookies.