/* edit REXX TRAP handler as needed (compare label TRAP near EOF) */ signal on novalue name TRAP ; signal on syntax name TRAP signal on failure name TRAP ; signal on halt name TRAP numeric digits 11 /* persistent results stored in file 'd:/tmp/primes_l.ist', */ /* length 300,000,000 % 30 = 10,000,000 bytes */ /* length 1,006,632,960 % 30 = 33,554,432 bytes (32 MB) */ call PRIME.! 'd:\progra~1\bin\projects\primes_l.ist', 300000000 /* -------------------------------------------------------------- */ /* the main program is essentially an interactive demo of PHI(), */ /* PI(), FACTOR(), and MU() as defined in the second part. This */ /* old demo could also be used with a sieve created on the fly, */ /* i.e. extended as needed. */ 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 '+': */ T = charin( PRIME.?, 1 ) A = charin( PRIME.?, 2 ) ; X = 0 do N = 3 to PRIME.LEN B = T ; T = A ; if X then say ; X = 0 A = charin( PRIME.?, N ) if bitand( T, '7E'x ) > '00'x then iterate N F = N * 30 - 45 ; LINE = right( F, 11 ) ; X = 1 F = ( F + 105 ) % 210 LINE = LINE '=' right( F, 9 ) '* 210 - 105: sextet' if F * 210 - 105 = N * 30 - 45 then call charout /**/, LINE else exit TRAP( 'oops' ) if bitand( B, '80'x ) > '00'x then iterate N if bitand( A, '01'x ) > '00'x then iterate N call charout /**/, ' + octet' ; if BREAK() then leave if bitand( B, '88'x ) > '00'x then iterate N if bitand( A, '11'x ) > '00'x then iterate N call charout /**/, '+ decade' ; if BREAK() then leave end N end /* ----------------- */ when N = '-' then do /* undocumented '-': */ 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 say LINE ; LINE = '' 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 ) otherwise say 'prime' 0-N 'is' PRIME( -N ) end /* ----------------- */ exit 0 HALT: say ' *...' ; exit 1 /* line end divisors */ BREAK: procedure /* await pressed key */ parse version X Y . /* (depending on OS) */ select when X = 'REXXSAA' | 6 <= Y 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 */ /* -------------------------------------------------------------- */ 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 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 DRAIM: procedure expose PRIME. /* Draim's factorization */ X = trunc( abs( arg( 1 ))) ; if X <= 1 then return 0 M = X ; Y = 2 if X // 2 then do Y = 3 by 2 until R = 0 Q = X % Y ; R = X // Y ; M = M - 2 * Q ; X = M + R end Y return Y PI: procedure expose PRIME. /* count primes <= X */ arg X ; L = 1 parse var PRIME.OLD K B /* file access accelerator */ if X < B * 30 then do ; K = 0 ; B = 1 ; end do B = B to PRIME.LEN while B * 30 < X C = charin( PRIME.TAB, B ) ; L = K ; K = K + PRIME.C end B if K > 1 then PRIME.OLD = L B /* file access accelerator */ do K = L while PRIME( K ) <= X ; end K ; return K - 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 ISPRIME: procedure expose PRIME. /* return 0 resp. 1 */ X = trunc( abs( arg( 1 ))) ; signal on notready name TRAP if 1 < X & X < 6 then return X <> 4 /* 2, 3, and 5 are prime */ C = PRIME.?( X ) ; X = X % 30 ; if C == '00'x then return 0 B = charin( PRIME.TAB, X + 1 ) ; return bitand( B, C ) = '00'x PRIME: procedure expose PRIME. /* return N-th prime */ call trace 'O' ; arg X select /* select special cases: */ when X = 0 then return 0 ; when X = 1 then return 2 when X = 2 then return 3 ; when X = 3 then return 5 when X < 0 | 0 % X then exit TRAP( 'invalid argument' ) otherwise X = X - 3 /* wanted set sieve bit */ end parse var PRIME.OLD K B /* file access accelerator */ if X < K then do ; K = 0 ; B = 1 ; end do B = B to PRIME.LEN /* evaluate sieve byte B */ C = charin( PRIME.TAB, B ) ; L = K K = K + PRIME.C ; if X > K then iterate B PRIME.OLD = L B /* file access accelerator */ do K = 1 to 29 /* evaluate bits in last B */ J = PRIME.?( K ) ; if J == '00'x then iterate K L = L + ( bitand( C, J ) == '00'x ) if X = L then return 30 * B + K - 30 end K /* - 1,7,11,13,17,19,23,29 */ exit TRAP( 'assertion failed in procedure PRIME' ) end B exit TRAP( 'prime' X 'not supported in' PRIME.TAB ) /* -------------------------------------------------------------- */ PRIME.!: procedure expose PRIME. /* create new primes file */ PRIME.LEN = arg( 2 ) % 30 + 1 /* hardwired sieve length */ PRIME.MAX = PRIME.LEN * 30 /* hardwired upper limit */ PRIME.TAB = arg( 1 ) /* hardwired sieve path */ PRIME.OLD = 0 1 /* file access accelerator */ PRIME.BIT = copies( '00'x, 30 ) /* 30 = 2 * 3 * 5 */ PRIME.BIT = overlay( d2c( 2 ** 0 ), PRIME.BIT, 1 + 1 ) PRIME.BIT = overlay( d2c( 2 ** 1 ), PRIME.BIT, 1 + 7 ) PRIME.BIT = overlay( d2c( 2 ** 2 ), PRIME.BIT, 1 + 11 ) PRIME.BIT = overlay( d2c( 2 ** 3 ), PRIME.BIT, 1 + 13 ) PRIME.BIT = overlay( d2c( 2 ** 4 ), PRIME.BIT, 1 + 17 ) PRIME.BIT = overlay( d2c( 2 ** 5 ), PRIME.BIT, 1 + 19 ) PRIME.BIT = overlay( d2c( 2 ** 6 ), PRIME.BIT, 1 + 23 ) PRIME.BIT = overlay( d2c( 2 ** 7 ), PRIME.BIT, 1 + 29 ) do B = 0 to 255 /* zero bit counter: */ X = B ; C = 8 /* (for sieve bytes) */ do while X > 0 ; C = C - X // 2 ; X = X % 2 ; end X = d2c( B ) ; PRIME.X = C /* CAVEAT: all PRIME.? for */ end B /* any character ? in use */ DONE = stream( PRIME.TAB, 'command', 'query size' ) if PRIME.LEN <= DONE then do /* use existing file */ signal on notready name TRAP /* trap write access */ return stream( PRIME.TAB, 'command', 'open read' ) end /* allow read access */ if DONE < 1 then do /* create new sieve */ DONE = 1 ; call charout PRIME.TAB, d2c( 1 ) end /* initialize new "black sieve" based on 2 * 3 * 5 = 30 using... */ /* 1= 30 -29 , 31= 60 -29 (bit 8), 61= 90 -29 (bit 16) */ /* 7= 30 -23 (bit 1), 37= 60 -23 (bit 9), 67= 90 -23 (bit 17) */ /* 11= 30 -19 (bit 2), 41= 60 -19 (bit 10), 71= 90 -19 (bit 18) */ /* 13= 30 -17 (bit 3), 43= 60 -17 (bit 11), 73= 90 -17 (bit 19) */ /* 17= 30 -13 (bit 4), 47= 60 -13 (bit 12), 77= 90 -13 */ /* 19= 30 -11 (bit 5), 49= 60 -11 , 79= 90 -11 (bit 21) */ /* 23= 30 - 7 (bit 6), 53= 60 - 7 (bit 14), 83= 90 - 7 (bit 22) */ /* 29= 30 - 1 (bit 7), 59= 60 - 1 (bit 15). 89= 90 - 1 (bit 23) */ C = max( 1, ( PRIME.LEN - DONE ) % 100 ) /* 100 NUL strings */ do B = DONE + 1 to PRIME.LEN by C /* of zero bits for */ C = min( C, PRIME.LEN - B ) /* prime candidates */ call charout PRIME.TAB, copies( '00'x, C ) end B do NEXT = 7 by 2 while NEXT * NEXT <= PRIME.MAX C = PRIME.?( NEXT ) ; if C == '00'x then iterate NEXT B = 1 + NEXT % 30 /* byte B and mask C */ if bitand( C, charin( PRIME.TAB, B )) == '00'x then do call charout 'stderr', d2c( 13 ) left( NEXT, 78 ) TEST = 2 * NEXT /* skip if N * 2 */ COMP = (( DONE * 30 ) % TEST ) * TEST + NEXT NICE = COMP + 50000 do COMP = max( COMP, NEXT * NEXT ) to PRIME.MAX by TEST C = PRIME.?( COMP ) ; if C = '00'x then iterate COMP B = 1 + COMP % 30 /* byte B and mask C */ X = charin( PRIME.TAB, B ) ; C = bitor( C, X ) if X == C then iterate COMP /* got old composite */ call charout PRIME.TAB, C, B /* set new composite */ if B = PRIME.LEN then call charout PRIME.TAB if COMP > NICE then do call charout 'stderr', d2c( 13 ) NEXT COMP NICE = COMP + 50000 end end COMP end end NEXT return charout( 'stderr', d2c( 13 ) left( '', 76 ) d2c( 13 )) PRIME.?: return substr( PRIME.BIT, 1 + arg( 1 ) // 30, 1 ) /* see , (c) F. Ellermann */ TRAP: /* select REXX exception handler */ call trace 'O' ; trace N /* don't trace interactive */ parse source TRAP /* source on separate line */ TRAP = x2c( 0D ) || right( '+++', 10 ) TRAP || x2c( 0D0A ) TRAP = TRAP || right( '+++', 10 ) /* = standard trace prefix */ TRAP = TRAP strip( condition( 'c' ) 'trap:' condition( 'd' )) select when wordpos( condition( 'c' ), 'ERROR FAILURE' ) > 0 then do if condition( 'd' ) > '' /* need an additional line */ then TRAP = TRAP || x2c( 0D0A ) || right( '+++', 10 ) TRAP = TRAP '(RC' rc || ')' /* any system error codes */ if condition( 'c' ) = 'FAILURE' then rc = -3 end when wordpos( condition( 'c' ), 'HALT SYNTAX' ) > 0 then do if condition( 'c' ) = 'HALT' then rc = 4 if condition( 'd' ) > '' & condition( 'd' ) <> rc then do if condition( 'd' ) <> errortext( rc ) then do TRAP = TRAP || x2c( 0D0A ) || right( '+++', 10 ) TRAP = TRAP errortext( rc ) end /* future condition( 'd' ) */ end /* may use errortext( rc ) */ else TRAP = TRAP errortext( rc ) rc = -rc /* rc < 0: REXX error code */ end when condition( 'c' ) = 'NOVALUE' then rc = -2 /* dubious */ when condition( 'c' ) = 'NOTREADY' then rc = -1 /* dubious */ otherwise /* force non-zero whole rc */ if datatype( value( 'RC' ), 'W' ) = 0 then rc = 1 if rc = 0 then rc = 1 if condition() = '' then TRAP = TRAP arg( 1 ) end /* direct: TRAP( message ) */ TRAP = TRAP || x2c( 0D0A ) || format( sigl, 6 ) signal on syntax name TRAP.SIGL /* throw syntax error 3... */ if 0 < sigl & sigl <= sourceline() /* if no handle for source */ then TRAP = TRAP '*-*' strip( sourceline( sigl )) else TRAP = TRAP '+++ (source line unavailable)' TRAP.SIGL: /* ...catch syntax error 3 */ if abbrev( right( TRAP, 2 + 6 ), x2c( 0D0A )) then do TRAP = TRAP '+++ (source line unreadable)' ; rc = -rc end select when 0 then do /* in pipes STDERR: output */ parse version TRAP.REXX /* REXX/Personal: \dev\con */ if abbrev( TRAP.REXX, 'REXXSAA ' ) | /**/ , 6 <= word( TRAP.REXX, 2 ) then TRAP.REXX = 'STDERR' else TRAP.REXX = '\dev\con' signal on syntax name TRAP.FAIL call lineout TRAP.REXX , TRAP /* fails if no more handle */ end when 0 then do /* OS/2 PM or ooREXX on NT */ signal on syntax name TRAP.FAIL call RxMessageBox translate( TRAP, ' ', x2c( 0D )), /**/ , 'Trap' time(),, 'ERROR' end otherwise say TRAP ; trace ?L /* interactive Label trace */ end if condition() = 'SIGNAL' then signal TRAP.EXIT TRAP.CALL: return rc /* continue after CALL ON */ TRAP.FAIL: say TRAP ; rc = 0 - rc /* force TRAP error output */ TRAP.EXIT: exit rc /* exit for any SIGNAL ON */