/* SIEVE prime numbers or determine smallest PRIME divisors */ /* The following lines upto procedure PRIME demonstrate the */ /* usage of FACTOR() and PRIME(). The lines below procedure */ /* FACTOR could be copied to other scripts. */ signal on novalue ; signal on halt ; numeric digits 40 parse upper source . T F ; arg N ; LINE = '' F = substr( F, 1+ lastpos( '\', F )) F = substr( F, 1+ lastpos( '/', F )) if \ verify( '.', F ) then F = left( F, pos( '.', F ) - 1 ) if verify( '(', N ) then N = strip( N ) else interpret 'N =' N select /* ----------------- */ when N = '+' then do /* undocumented '+': */ M = 1 ; numeric digits 50 /* efficiency limit */ do N = 1 /* for prime sieves */ PRIME = PRIME( N ) M = M * ( PRIME - 1 ) / PRIME say right( N, 10 ) right( PRIME, 12 ) || ':' M end N end /* ----------------- */ when N = '-' then do /* undocumented '-': */ S = 0 do N = 1 to 10000 P = right( PRIME( N ), 9 ) S = right( S + P, 9 ) T = right( P * 13 + 28, 9 ) if S <= T then say right( N, 6 ) P S ' ' T else say right( N, 6 ) P S '>' T end N end /* ----------------- */ when N <> '.' & \ datatype( N, 'N' ) then do if T <> 'FUNCTION' then do say F '0 table of primes : if not used' say F '1 table of divisors : as function' say F '. interactive usage : ' F || '( . )' say F 'n all resp. 1st divisor(s) of n > 1' say F "+r Euler's phi( r ) " say F '-p p-th prime number' say F '.c count primes <= c' say F '(?) interprets (?) as its argument' end else say F || '(' N ') ? Use [+|-|.] N or CALL' F end /* ----------------- */ when N > '.' & abbrev( N, '.' ) then do N = N * ( 10 ** length( N )) % 10 if T = 'FUNCTION' then return PI( N ) else say PI( N ) 'primes <=' N end /* ----------------- */ when \ datatype( N, 'W' ) then do if T <> 'FUNCTION' then do N = 'EXIT : exit' F say 'SAY ... : shows result of ... ,' N N = 'MU( N ) : Moebius function' say '... : REXX interprets ... ,' N N = 'PI( N ) : counts primes <= N' say 'FACTOR(N) : prime divisor of N ,' N N = 'PHI( N ) : totient function' say 'PRIME( N ): N-th prime number ,' N say ; signal off novalue do forever ; interpret linein() ; end end else say F || '(' N ') not supported, use CALL' F N end /* ----------------- */ when abbrev( N, '+' ) then do if T = 'FUNCTION' then return PHI( N ) else say 'phi(' N ') =' PHI( N ) end /* ----------------- */ when N = 0 & T <> 'FUNCTION' then do say 'sieve first prime numbers (ESC to stop):' do N = 0 by 10 until BREAK() LINE = right( N, 7 ) || ':' do M = N to N + 9 LINE = LINE || right( PRIME( M ) , 7 ) end M /* BREAK look ahead: */ say LINE ; LINE = '' ; call PRIME M + 8 end N end /* ----------------- */ when N = 1 & T <> 'FUNCTION' then do say 'find smallest divisors > 1 (ESC to stop):' do M = 0 to 5 LINE = LINE || right( M, 6 ) LINE = LINE || right( FACTOR( M ), 6 ) || ',' end M say LINE ; LINE = '' do N = 7 by 12 until BREAK() do M = N to N + 10 by 2 /* 7..17. 19..29, .. */ LINE = LINE || right( M, 6 ) LINE = LINE || right( FACTOR( M ), 6 ) || ',' end M say LINE ; LINE = '' end N end /* ----------------- */ when N > 1 & T <> 'FUNCTION' then do call charout /**/, N '=' sign( N ) do while abs( N ) > 1 /* show ALL divisors */ M = FACTOR( N ) ; N = N % M call charout /**/, ' *' M end end /* ----------------- */ when N > 0 then return FACTOR( N ) when T = 'FUNCTION' then return PRIME( -N ) when N = -1 then say '1st prime is' PRIME( -N ) when N = -2 then say '2nd prime is' PRIME( -N ) when N = -3 then say '3rd prime is' PRIME( -N ) when N > -199 then say -N || '-th prime is' PRIME( -N ) otherwise /* ----------------- */ call charout /**/, left( -N, 16 ) do M = 1 to 10 /* show progress bar */ call charout /**/, '0....' || M end M call charout /**/, '0%' || d2c( 13 ) call charout /**/, left( -N, 11 ) do M = 1 to 64 /* 64 = 10*6 + 15-11 */ call PRIME -N * M % 64 call charout /**/, '#' /* show progress dot */ end M say ; say -N || '-th prime is' PRIME( -N ) pull end /* ----------------- */ exit 0 HALT: say ' *...' ; exit 1 /* line end divisors */ BREAK: procedure /* await pressed key */ parse version X . . /* (depending on OS) */ select when X = 'REXXSAA' then do parse source X . . if X = 'DOS' then do /* external PC DOS 7 */ X = RxGetKey( 'NoEcho' ) /* RXGETKEY.RX */ if X = d2c( 0 ) | X = d2c( 224 ) then X = X || RxGetKey( 'NoEcho' ) end else do /* assume loaded IBM */ X = SysGetKey( 'NoEcho' ) /* OS/2 REXXUTIL.DLL */ if X = d2c( 0 ) | X = d2c( 224 ) then X = X || SysGetKey( 'NoEcho' ) end end when X = 'REXX/Personal' then X = InKey() otherwise X = charin() end X = c2d( X ) + 256 * ( length( X ) - 1 ) return X = 27 | X = 317 | X = 363 /* ESC | F3 | Alt-F4 */ /* -------------------------------------------------------- */ /* PI( 10000000000) = 455052512, PI( 316241) = 27294 */ /* PI( 1000000000) = 50847534, PRIME( 27293) = 316223 */ /* PI( 100000000) = 5761455, 316241 * 316241 > 10 ** 11 */ /* PI( 10000000) = 664579, 316223 * 316223 < 10 ** 11 */ /* PI( 1000000) = 78498, PRIME( 1000000) = 15485863 */ /* PI( 100000) = 9592, PRIME( 100000) = 1299709 */ /* PI( 10000) = 1229, PRIME( 10000) = 104729 */ /* PI( 1000) = 168, PRIME( 1000) = 7919 */ /* PI( 100) = 25, PRIME( 100) = 541 */ /* PI( 10) = 4, PRIME( 10) = 29 */ /* PI( 1) = 0, PRIME( 1) = 2 */ /* FACTOR( 2 ** (2 ** 5) +1 )= 641 *n, 5th Fermat number */ /* FACTOR( 2 ** (2 ** 6) +1 )= 274177 *n, 6th Fermat number */ /* FACTOR( 2**11 -1) = 23 *n, FACTOR( 2**23 -1) = 47 *n */ /* FACTOR( 2**37 -1) = 223 *n, FACTOR( 2**41 -1) = 13367 *n */ /* FACTOR( 2**43 -1) = 431 *n, FACTOR( 2**47 -1) = 2351 *n */ /* FACTOR( 2**53 -1) =6361 *n, FACTOR( 2**59 -1) =179951 *n */ /* Mersenne primes 2**k-1 for k= 2,3,5,7,13,17,19,31,61,89, */ /* 107,127,521,607,1279,2203,2281,3217,4253,4423,9689,9941. */ /* 2**44497 - 1, 13395 digits, unknown source */ /* 2**11213 - 1, (Mersenne), 3376 digits, by Gilles 1964, */ /* 1+ 180 * ((2**127-1) ** 2), 79 digits, by Miller 1951, */ /* (2**148 + 1) / 17, 44 digits, by Ferrier 1951, */ /* 2**127 - 1 (127th Mersenne), 39 digits, by Lucas 1876, */ /* 2**107 - 1 (107th Mersenne), 33 digits, prime test 1644. */ FACTOR: procedure expose PRIME. /* first divisor > 1 */ X = trunc( abs( arg( 1 ))) ; if X <= 1 then return 0 do N = 1 until X // Y = 0 ; Y = PRIME( N ) ; end N return Y MU: procedure expose PRIME. /* Moebius function: */ X = trunc( abs( arg( 1 ))) ; if X <= 1 then return 1 Y = 0 ; R = 1 do forever Z = Y ; Y = FACTOR( X ) ; if Y = Z then return 0 R = -R ; X = X % Y ; if X = 1 then return R end PI: procedure expose PRIME. /* count primes <= X */ do N = 1 while PRIME( N ) <= arg( 1 ) ; end N return N - 1 PHI: procedure expose PRIME. /* Euler's phi( X ): */ X = trunc( abs( arg( 1 ))) ; Y = 2 ; Z = X do N = 2 while Y <= X if X // Y = 0 then Z = Z * ( Y-1 ) / Y ; Y = PRIME( N ) end N return Z / 1 PRIME: procedure expose PRIME. /* return X-th prime */ if value( 'PRIME.0' ) <> 0 then do /* initialize sieve: */ PRIME.0 = 0 ; PRIME.! = -4 PRIME.1 = 2 ; N = -1 ; PRIME.N = 15 PRIME.2 = 3 ; N = -2 ; PRIME.N = 53 PRIME.3 = 5 ; N = -5 ; PRIME.N = 121 PRIME.4 = 7 ; PRIME.10 = 29 PRIME.5 = 11 ; PRIME.11 = 31 PRIME.6 = 13 ; PRIME.12 = 37 PRIME.7 = 17 ; PRIME.13 = 41 PRIME.8 = 19 ; PRIME.14 = 43 PRIME.9 = 23 ; PRIME.15 = 47 PRIME.? = ' 0 6 8 14 18 20 26 30' PRIME.? = PRIME.? ' 36 44 48 50 54 56 60 68' PRIME.? = PRIME.? ' 74 78 84 86 90 96 98 104' PRIME.? = PRIME.? '110 114 116 120 126 128 134 138' PRIME.? = PRIME.? '140 144 146 156 158 168 170 174' PRIME.? = PRIME.? '176 180 186 188 194 198 200 204' PRIME.. = '' do N = 1 to 15 PRIME.. = PRIME.. PRIME.N end N end /* PRIME.0 .. PRIME.N: 0..N-th prime, init. 0,2,..,47 */ /* PRIME.. string of found primes, only a */ /* hack to bypass REXXSAA problems */ /* with "too many" variables using */ /* PRIME.N = word( PRIME.. , N ) */ /* PRIME.-1: found primes counter initial 15 */ /* PRIME.-2: next prime candidate initial 53 */ /* PRIME.?: 48 candidates added to PRIME.-2 */ /* PRIME.!: 4 - sieved primes, init. 4-4= 0 */ /* PRIME.-3: not used => no 5 * n candidates */ /* PRIME.-4: not used => no 7 * n candidates */ /* PRIME.S-1: next test square, S-1=PRIME.!-1 */ /* PRIME.S .. PRIME.-5: sieve (S=PRIME.!), initial void */ /* The implemented "dumb" sieve tests 48 of 210 numbers, */ /* a "smarter" sieve could use the following progression */ /* 1= 1 of 2=2 1/ 2 */ /* 2*1= 2 of 6=2*3 1/ 3 */ /* 4*2*1= 8 of 30=2*3*5 4/ 15 */ /* 6*4*2*1= 48 of 210=2*3*5*7 8/ 35 */ /* 10*6*4*2*1= 480 of 2310=2*3*5*7*11 16/ 77 */ /* 12*10*6*4*2*1= 5760 of 30030=2*3*5*7*11*13 192/ 1001 */ /* etc. 92160 of 510510 3072/17017 */ X = trunc( abs( arg( 1 ))) ; numeric digits 12 S = -1 ; T = -2 /* T: Test candidate */ do while PRIME.S <= X /* S: Sieved primes */ N = PRIME.T ; PRIME.T = N + 210 /* 210 = 7*5*3*2 */ do M = 1 to 48 /* 48 = 210-162 */ call PRIME.! PRIME.?( N + word( PRIME.?, M )) end M /* 105 = (210 )/ 2, 35 = (210-105)/3 */ end /* 14 = (210-140)/ 5, 8 = (210-154)/7 */ /* 162 = 105+35+14+8 */ if X < - PRIME.! then return PRIME.X /* fast sieve access */ else return word( PRIME.. , X ) /* slow prime string */ PRIME.!: procedure expose PRIME. /* update the sieve, */ arg X /* called by PRIME() */ do N = PRIME.! to -5 if PRIME.N <= X then do S = -N ; PRIME.N = PRIME.N + 2 * PRIME.S end /* "lazy" sieve, odd multiples generated, */ end N /* but multiples of 3, 5, 7 never tested */ return PRIME.?: procedure expose PRIME. /* Eratosthenes test */ arg X /* called by PRIME() */ do N = -5 to PRIME.! by -1 if PRIME.N = X then return X /* matching no prime */ end N if X = PRIME.N then do /* test square prime */ S = PRIME.! -2 ; PRIME.! = S +1 ; N = -S PRIME.. = . subword( PRIME.. , 2 ) PRIME.N = word( PRIME.. , N ) /* fast sieve access */ PRIME.S = PRIME.N * PRIME.N ; return X end S = -1 ; N = PRIME.S + 1 ; PRIME.S = N PRIME.. = PRIME.. X ; return X /* found a new prime */