Carmichael 3 strong pseudoprimes: Difference between revisions
Content added Content deleted
m (→{{header|REXX}}: optimized the isPrime function a wee bit.) |
(→{{header|REXX}}: added a 2nd REXX version, added/changed comments and whitespace, aligned some statements, 2nd version has many more primes.) |
||
Line 1,061: | Line 1,061: | ||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
==optimized== |
|||
Note that REXX's version of '''modulus''' (<big><code>'''//'''</code></big>) is really a ''remainder'' function, |
|||
Note that REXX's version of '''modulus''' (<big><code>'''//'''</code></big>) is really a ''remainder'' function. |
|||
<br>(It was necessary to use '''modulus''' instead of '''remainder''' when using a negative value.) |
|||
<br>The Carmichael numbers are shown in numerical order. |
<br>The Carmichael numbers are shown in numerical order. |
||
<br>Some code optimization was done, while not necessary for the small default number ('''61'''), it was significant for larger numbers. |
<br>Some code optimization was done, while not necessary for the small default number ('''61'''), it was significant for larger numbers. |
||
<lang rexx>/*REXX program calculates Carmichael |
<lang rexx>/*REXX program calculates Carmichael 3─strong pseudoprimes (up to and including N). */ |
||
numeric digits 30 /*handle big dig #s (9 is the default).*/ |
numeric digits 30 /*handle big dig #s (9 is the default).*/ |
||
parse arg N .; if N=='' then N=61 /*allow user to specify for the search.*/ |
parse arg N .; if N=='' then N=61 /*allow user to specify for the search.*/ |
||
Line 1,079: | Line 1,078: | ||
do h3=2 to pm; g=h3+p /*find Carmichael #s for this prime. */ |
do h3=2 to pm; g=h3+p /*find Carmichael #s for this prime. */ |
||
gPM=g*pm; npsH3=((nps//h3)+h3)//h3 /*define a couple of shortcuts for pgm.*/ |
gPM=g*pm; npsH3=((nps//h3)+h3)//h3 /*define a couple of shortcuts for pgm.*/ |
||
/*perform some weeding of |
/* [↓] perform some weeding of D values*/ |
||
do d=1 for g-1; |
do d=1 for g-1; if gPM//d \== 0 then iterate |
||
if npsH3 \== d//h3 then iterate |
|||
q=1+gPM%d; if \isPrime(q) then iterate |
|||
r=1+p*q%h3; if q*r//pm\==1 then iterate |
|||
if \isPrime(r) then iterate |
if \isPrime(r) then iterate |
||
carms=carms+1; @.q=r |
carms=carms+1; @.q=r /*bump Carmichael counter; add to array*/ |
||
if bot==0 then bot=q; bot=min(bot,q); top=max(top,q) |
if bot==0 then bot=q; bot=min(bot,q); top=max(top,q) |
||
end /*d*/ |
end /*d*/ /* [↑] find the minimum & the maximum.*/ |
||
end /*h3*/ |
end /*h3*/ |
||
$=0 /*display a list of some Carmichael #s.*/ |
$=0 /*display a list of some Carmichael #s.*/ |
||
do j=bot to top by 2 while tell; if @.j==0 then iterate; $=1 |
do j=bot to top by 2 while tell; if @.j==0 then iterate; $=1 |
||
say '──────── a Carmichael number: ' p "∙" j |
say '──────── a Carmichael number: ' p "∙" j '∙' @.j |
||
end /*j*/ |
end /*j*/ |
||
if $ then say /*show a blank line for beautification.*/ |
if $ then say /*show a blank line for beautification.*/ |
||
end /*p*/ |
end /*p*/ |
||
say; say carms ' Carmichael numbers found.' |
say; say carms ' Carmichael numbers found.' |
||
exit /*stick a fork in it, we're all done. */ |
exit /*stick a fork in it, we're all done. */ |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
isPrime: |
isPrime: parse arg x; if !.x then return 1 /*X a known prime?*/ |
||
if x<37 then return 0; if x//2==0 then return 0; if x// 3==0 then return 0 |
if x<37 then return 0; if x//2==0 then return 0; if x// 3==0 then return 0 |
||
parse var x '' -1 _; if _==5 then return 0; if x// 7==0 then return 0 |
parse var x '' -1 _; if _==5 then return 0; if x// 7==0 then return 0 |
||
Line 1,187: | Line 1,186: | ||
1038 Carmichael numbers found. |
1038 Carmichael numbers found. |
||
</pre> |
</pre> |
||
==more robust== |
|||
This REXX version (pre-)generates around 5,000 primes to assist the '''isPrime''' function. |
|||
<lang rexx>/*REXX program calculates Carmichael 3─strong pseudoprimes (up to and including N). */ |
|||
numeric digits 30 /*handle big dig #s (9 is the default).*/ |
|||
parse arg N .; if N=='' then N=61 /*allow user to specify for the search.*/ |
|||
tell= N>0; N=abs(N) /*N>0? Then display Carmichael numbers*/ |
|||
carms=0 /*number of Carmichael numbers so far. */ |
|||
!.=0; !.2=1; !.3=1; !.5=1; !.7=1; !.11=1; !.13=1; !.17=1; !.19=1; !.23=1; !.29=1; !.31=1 |
|||
do i=23 by 2 for 5000; if isPrime(i) then do; !.i=1; HI=i; end; end /*i*/ |
|||
HI=HI+2 |
|||
/*[↑] prime number memoization array. */ |
|||
do p=3 to N by 2; pm=p-1; bot=0; top=0 /*step through some (odd) prime numbers*/ |
|||
nps=-p*p; if \isPrime(p) then iterate /*is P a prime? No, then skip it.*/ |
|||
@.=0 |
|||
do h3=2 to pm; g=h3+p /*find Carmichael #s for this prime. */ |
|||
gPM=g*pm; npsH3=((nps//h3)+h3)//h3 /*define a couple of shortcuts for pgm.*/ |
|||
/* [↓] perform some weeding of D values*/ |
|||
do d=1 for g-1; if gPM//d \== 0 then iterate |
|||
if npsH3 \== d//h3 then iterate |
|||
q=1+gPM%d; if \isPrime(q) then iterate |
|||
r=1+p*q%h3; if q*r//pm\==1 then iterate |
|||
if \isPrime(r) then iterate |
|||
carms=carms+1; @.q=r /*bump Carmichael counter; add to array*/ |
|||
if bot==0 then bot=q; bot=min(bot,q); top=max(top,q) |
|||
end /*d*/ /* [↑] find the minimum & the maximum.*/ |
|||
end /*h3*/ |
|||
$=0 /*display a list of some Carmichael #s.*/ |
|||
do j=bot to top by 2 while tell; if @.j==0 then iterate; $=1 |
|||
say '──────── a Carmichael number: ' p "∙" j '∙' @.j |
|||
end /*j*/ |
|||
if $ then say /*show a blank line for beautification.*/ |
|||
end /*p*/ |
|||
say; say carms ' Carmichael numbers found.' |
|||
exit /*stick a fork in it, we're all done. */ |
|||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
|||
isPrime: parse arg x; if !.x then return 1 /*X a known prime?*/ |
|||
if x<HI then return 0; if x//2==0 then return 0; if x// 3==0 then return 0 |
|||
parse var x '' -1 _; if _==5 then return 0; if x// 7==0 then return 0 |
|||
if x//11==0 then return 0; if x//13==0 then return 0 |
|||
if x//17==0 then return 0; if x//19==0 then return 0 |
|||
do k=23 by 6 until k*k>x; if x// k ==0 then return 0 |
|||
if x//(k+2)==0 then return 0 |
|||
end /*k*/ /*K will never be divisible by three.*/ |
|||
!.x=1; return 1 /*Define a new prime (X). Indicate so.*/</lang> |
|||
'''output''' is identical to the 1<sup>st</sup> REXX version. <br><br> |
|||
=={{header|Ruby}}== |
=={{header|Ruby}}== |