|
A method exists in this REXX program to also test Mersenne-type numbers <big>(2<sup>n</sup> - 1)</big>.
Since the majority of computing time is spent looking for primes, that part of the program was
<br>optimized somewhat (but could be extended if more optimization is wanted).
<syntaxhighlight lang="rexx"></*REXX pgm does prime decomposition of a range of positive integers (with a prime count)*/syntaxhighlight>
/*REXX pgm does prime decomposition of a range of positive integers (with a prime count)*/
numeric digits 1000 /*handle thousand digits for the powers*/
parseNumeric argDigits 1000 bot top step base add /*gethandle optionalthousand argumentsdigits fromFor the C.L. powers*/
ifParse Arg bot=='' top then do;step bot=1; top=100; end base add /*noget optional BOTarguments given? Then usefrom the defaultC.L. */
If bot='?' Then Do
if top=='' then top=bot /* " TOP? " " " " " */ ▼ Say 'rexx pfoo bot top step base add'
if step=='' then step= 1 /* " STEP? " " " " " */ ▼ Exit
if add =='' then add= -1 /* " ADD? " " " " " */ ▼ End
tell= top>0; top=abs(top) /*if TOP is negative, suppress displays*/
wIf bot==length(top)'' Then /* no arguments given /*get maximum width for aligned display*/
Parse Value 1 100 With bot top /* set default range .*/
if base\=='' then w=length(base**top) /*will be testing powers of two later? */
▲if addIf top=='' Then then addtop= -1 bot /* " ADD? " process one number " " " " */
@.=left('', 7); @.0="{unity}"; @.1='[prime]' /*some literals: pad; prime (or not).*/ ▼numericIf digitsstep=='' max(9,Then w+step=1) /* step=2 to process only odd numbers /*maybe increase the digits precision. */
#=0If add =='' Then add=-1 /* for Mersenne tests /*#: is the number of primes found. */
tell=top>0 do n=bot to top by step /*process a single numberIf TOP oris negative, asuppress range.displays*/
top=abs(top)
?=n; if base\=='' then ?=base**n + add /*should we perform a "Mercenne" test? */
w=length(top) pf=factr(?); f=words(pf) /*get primemaximum factors;width numberFor ofaligned factors.display*/
If base\=='' Then
if f==1 then #=#+1 /*Is N prime? Then bump prime counter.*/ ▼ w=length(base**top) if tell then say right(?,w) right('('f")",9) 'prime/*will factors:be ' @.ftesting powers of two later? pf*/
▲@tag.=left('', 7) ; @.0="{unity}"; @.1='[prime]' /*some literals: pad; prime (or not).*/
tag.0='{unity}'
say
tag.1='[prime]'
ps= 'primes'; if p==1 then ps= "prime" /*setup for proper English in sentence.*/ ▼sayNumeric rightDigits max(#9, w+9+1) ps 'found.' /*displaymaybe increase the numberdigits of primes foundprecision. */
exit np=0 /*sticknp: a fork in it,is the we'renumber allof doneprimes found. */
Do n=bot To top by step /*process a single number or a range. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
?=n
factr: procedure; parse arg x 1 d,$ /*set X, D to argument 1; $ to null.*/ ▼if x If base\==1 then return '' Then /*handleshould thewe specialperform casea of X = 1.'Mersenne' test? */
?=base**n+add
do while x//2==0; $=$ 2; x=x%2; end /*append all the 2 factors of new X.*/
pf=factr(?) do while x//3==0; $=$ 3; x=x%3; end /* " "/* get prime "factors 3 " " " " */
f=words(pf) do while x//5==0; $=$ 5; x=x%5; end /* " " " /* 5number of prime factors " " " " */
If f=1 Then do while x//7==0; $=$ 7; x=x%7; end /* " " /* " If the 7number is prime " " " " */
np=np+1 /* Then bump prime counter ___*/
If tell Then
q=1; do while q<=x; q=q*4; end /*these two lines compute integer √ X */
Say right(?,w) right('('f')',9) 'prime factors: ' tag.f pf
r=0; do while q>1; q=q%4; _=d-r-q; r=r%2; if _>=0 then do; d=_; r=r+q; end; end
End /*n*/
Say ''
ps='prime'
▲If f>1 Then ps= ps' primess' ; if p==1 then ps= "prime" /*setup forFor proper English in sentence.*/
Say right(np, w+9+1) ps 'found.' /*display the number of primes found. */
▲Exit if f==1 then #=#+1 /* Isstick Na prime?fork in Thenit, bump primewe're counterall done. */
/*---------------------------------------------------------------------------*/
factr: Procedure
▲factr: procedure; Parse parse argArg x 1 d, $pl /*set X, D to argument 1 ; $, pl to null . */
If x==1 Then Return '' /*handle the special case of X=1. */
primes=2 3 5 7
Do While primes>'' /* first check the small primes */
Parse Var primes prime primes
Do While x//prime==0
pl=pl prime
x=x%prime
End
End
r=isqrt(x)
Do if _ j= =311 by then6 iterateTo r /* Isinsure nextthat J Y isisn't divisible by 5?3. Skip.*/ ▼
Parse var j '' -1 _ /*obtain the last decimal digit of J. */
If _\==5 Then Do
pl=pl j
x=x%j
▲if top=='' then End top=bot /* " TOP? " /*maybe reduce by J. " " " " */
End
If _ ==3 Then Iterate /*If next Y is divisible by 5? Skip. */
y=j+2
Do While x//y==0
pl=pl y
x=x%y
▲if step=='' then step= End 1 /* " STEP? " /*maybe reduce by y. " " " " */
If do jx==11 by 6 to r 1 Then Return pl /*insureIs thatresidual=unity? Then J isndon't divisible by 3append.*/
parse var j '' -1 _ Return pl x /*obtainreturn the last decimalpl digit of with Jappended residual. */
if _\==5 then do while x//j==0; $=$ j; x=x%j; end /*maybe reduce by J. */
isqrt: Procedure
▲ if _ ==3 then iterate /*Is next Y is divisible by 5? Skip.*/
Parse Arg x
y=j+2; do while x//y==0; $=$ y; x=x%y; end /*maybe reduce by J. */
x=abs(x)
Parse Value 0 x with lo hi
/* [↓] The $ list has a leading blank.*/
Do While lo<=hi
if x==1 then return $ /*Is residual=unity? Then don't append.*/
t=(lo+hi)%2
return $ x /*return $ with appended residual. */</syntaxhighlight>
If t**2>x Then
hi=t-1
Else
lo=t+1
End
Return t
'''output''' when using the default input of: <tt> 1 100 </tt>
<pre style="font-size:75%;height:160ex">
1 (0) prime factors: {unity}
2 (1) prime factors: [prime] 2
<pre style="font-size:75%>
3 (1) prime factors: [prime] 3
7 (1) prime factors: [prime] 7
4095 (5) prime factors: 3 3 5 7 13
8191 (1) prime factors: [prime] 8191
16383 (23) prime factors: 3 546143 127
32767 (23) prime factors: 7 468131 151
65535 (4) prime factors: 3 5 17 257
131071 (1) prime factors: [prime] 131071
262143 (56) prime factors: 3 3 3 7 138719 73
524287 (1) prime factors: [prime] 524287
1048575 (6) prime factors: 3 5 5 11 41 31 41
2097151 (34) prime factors: 7 7 42799127 337
4194303 (4) prime factors: 3 23 89 683
8388607 (2) prime factors: 47 178481
16777215 (7) prime factors: 3 3 5 7 13 17 241
33554431 (13) prime factors: [prime] 33554431 31 601 1801
67108863 (23) prime factors: 3 223696212731 8191
134217727 (23) prime factors: 7 1917396173 262657
268435455 (56) prime factors: 3 5 29 43 113 5461127
536870911 (3) prime factors: 233 1103 2089
1073741823 (57) prime factors: 3 3 7 11 154941131 151 331
2147483647 (1) prime factors: [prime] 2147483647
4294967295 (5) prime factors: 3 5 17 257 65537
8589934591 (4) prime factors: 7 23 89 599479
17179869183 (3) prime factors: 3 43691 131071
34359738367 (34) prime factors: 31 71 122921127 3937122921
68719476735 (710) prime factors: 3 3 3 5 7 13 559377119 37 73 109
137438953471 (12) prime factors: [prime] 137438953471 223 616318177
274877906943 (23) prime factors: 3 91625968981174763 524287
549755813887 (24) prime factors: 7 7853654484179 8191 121369
1099511627775 (78) prime factors: 3 5 5 11 17 31 41 191211161681
2199023255551 (2) prime factors: 13367 164511353
4398046511103 (58) prime factors: 3 3 7 7 997289458343 127 337 5419
8796093022207 (3) prime factors: 431 9719 2099863
17592186044415 (67) prime factors: 3 5 23 89 397 683 8388612113
35184372088831 (26) prime factors: 7 502633886983331 73 151 631 23311
70368744177663 (4) prime factors: 3 47 178481 2796203
140737488355327 (23) prime factors: 2351 598628193774513 13264529
281474976710655 (810) prime factors: 3 3 5 7 13 17 97 241 257 15732721673
562949953421311 (12) prime factors: [prime] 562949953421311 127 4432676798593
1125899906842623 (47) prime factors: 3 11 31 251 135928999981601 1801 4051
11 8 primes found.
</pre>
'''output''' when using the input of: <tt> 1 50 1 2 +1 </tt>
<pre style="font-size:75%>
3 (1) prime factors: [prime] 3
57 (1) prime factors: [prime] 57
915 (2) prime factors: 3 35
1731 (1) prime factors: [prime] 1731
3363 (23) prime factors: 3 113 7
65127 (21) prime factors: [prime] 5 13127
129255 (23) prime factors: 3 435 17
257511 (12) prime factors: [prime] 257 7 73
5131023 (43) prime factors: 3 3 311 1931
10252047 (32) prime factors: 523 5 4189
20494095 (25) prime factors: 3 6833 5 7 13
40978191 (21) prime factors: [prime] 17 2418191
819316383 (23) prime factors: 3 273143 127
1638532767 (3) prime factors: 57 2931 113151
3276965535 (4) prime factors: 3 35 1117 331257
65537131071 (1) prime factors: [prime] 65537131071
131073262143 (26) prime factors: 3 436913 3 7 19 73
262145524287 (31) prime factors: [prime] 5 13 4033524287
5242891048575 (26) prime factors: 3 1747635 5 11 31 41
10485772097151 (24) prime factors: 177 616817 127 337
20971534194303 (34) prime factors: 3 323 89 233017683
41943058388607 (2) prime factors: 547 838861178481
838860916777215 (27) prime factors: 3 27962033 5 7 13 17 241
1677721733554431 (23) prime factors: 25731 65281601 1801
3355443367108863 (43) prime factors: 3 11 2512731 40518191
67108865134217727 (43) prime factors: 5 537 161373 157262657
134217729268435455 (56) prime factors: 3 35 329 343 113 1657009127
268435457536870911 (23) prime factors: 17233 157903211103 2089
5368709131073741823 (37) prime factors: 3 593 7 11 31 151 3033169331
10737418252147483647 (51) prime factors: [prime] 5 5 13 41 805812147483647
21474836494294967295 (25) prime factors: 3 7158278835 17 257 65537
42949672978589934591 (24) prime factors: 6417 670041723 89 599479
858993459317179869183 (43) prime factors: 3 3 68343691 1397419131071
1717986918534359738367 (4) prime factors: 531 13771 953127 26317122921
34359738369 68719476735 (510) prime factors: 3 113 3 5 7 13 19 28137 8617173 43109
68719476737137438953471 (2) prime factors: 17223 4042322161616318177
137438953473274877906943 (23) prime factors: 3 45812984491174763 524287
274877906945549755813887 (24) prime factors: 57 5497558138979 8191 121369
5497558138891099511627775 (38) prime factors: 3 35 5 11 17 31 41 6108397932161681
10995116277772199023255551 (2) prime factors: 25713367 4278255361164511353
21990232555534398046511103 (38) prime factors: 3 833 7 7 43 127 337 88314186975419
43980465111058796093022207 (53) prime factors: 5431 139719 29 113 206476212099863
879609302220917592186044415 (27) prime factors: 3 29320310074035 23 89 397 683 2113
1759218604441735184372088831 (36) prime factors: 177 35331 293154241773 151 631 23311
3518437208883370368744177663 (54) prime factors: 3 3 347 11178481 1184658992892796203
70368744177665 140737488355327 (43) prime factors: 5 10132351 302694513 45898913264529
140737488355329 281474976710655 (210) prime factors: 3 469124961184433 5 7 13 17 97 241 257 673
281474976710657562949953421311 (2) prime factors: 65537127 42949017614432676798593
5629499534213131125899906842623 (27) prime factors: 3 18764998447377111 31 251 601 1801 4051
1125899906842625 (6) prime factors: 5 5 5 41 101 2175126601
58 primes found.
</pre>
| 