Talk:ALGOL 68-primes
Source code
This is the source of the "prime related" library used for some ALGOL 68 samples on Rosetta Code.
Note the is probably prime PROC is based on code from The Miller–Rabin primality test task and the prelude.
# primes.incl.a68: prime related operators, procedure etc. #
# returns a sieve of primes up to n #
OP PRIMESIEVE = ( INT n )[]BOOL:
BEGIN
[ 0 : n ]BOOL ps prime;
ps prime[ 0 ] := ps prime[ 1 ] := FALSE;
ps prime[ 2 ] := TRUE;
FOR i FROM 3 BY 2 TO UPB ps prime DO ps prime[ i ] := TRUE OD;
FOR i FROM 4 BY 2 TO UPB ps prime DO ps prime[ i ] := FALSE OD;
FOR i FROM 3 BY 2 TO ENTIER sqrt( UPB ps prime ) DO
IF ps prime[ i ] THEN
FOR s FROM i * i BY i + i TO UPB ps prime DO ps prime[ s ] := FALSE OD
FI
OD;
ps prime
END; # PRIMESIEVE #
# modes and operators for extracting primes from a sieve #
INT first n primes extract op = 0;
INT primes up to extract op = 1;
MODE PRIMEEXTRACT = STRUCT( INT extract op, INT n );
# returns a PRIMEEXTRACT with the first n primces extract op #
OP EXTRACTFIRST = ( INT n )PRIMEEXTRACT: PRIMEEXTRACT( first n primes extract op, n );
# returns a PRIMEEXTRACT with primes up to extract op #
OP EXTRACTPRIMESUPTO = ( INT n )PRIMEEXTRACT: PRIMEEXTRACT( primes up to extract op, n );
# returns an array of primes extracted from prime #
# as specified by by the extract op and number of primes in ex #
PRIO PRIMESFROMSIEVE = 9, FROMPRIMESIEVE = 9;
OP PRIMESFROMSIEVE = ( PRIMEEXTRACT ex, []BOOL fs prime )[]INT:
IF extract op OF ex = first n primes extract op
THEN
# extract the first n OF ex primes #
# count the primes up n OF ex #
INT pfs count := 0;
FOR i TO UPB fs prime WHILE pfs count < n OF ex DO
IF fs prime[ i ] THEN pfs count +:= 1 FI
OD;
# construct a list of the first n OF ex primes or p count, #
# whichever is smaller #
[ 1 : pfs count ]INT pfs low prime;
INT pfs low pos := 0;
FOR i WHILE pfs low pos < UPB pfs low prime DO
IF fs prime[ i ] THEN pfs low prime[ pfs low pos +:= 1 ] := i FI
OD;
pfs low prime
ELSE
# extract primes up to n OF ex primes #
INT max fs prime = n OF ex;
# count the primes up to max prime #
INT pfs count := 0;
FOR i TO max fs prime DO IF fs prime[ i ] THEN pfs count +:= 1 FI OD;
# construct a list of the primes up to max prime #
[ 1 : pfs count ]INT pfs low prime;
INT pfs low pos := 0;
FOR i WHILE pfs low pos < UPB pfs low prime DO
IF fs prime[ i ] THEN pfs low prime[ pfs low pos +:= 1 ] := i FI
OD;
pfs low prime
FI # PRIMESFROMSIEVE # ;
# alternative spelling of the operator name #
OP FROMPRIMESIEVE = ( PRIMEEXTRACT ex, []BOOL ex prime )[]INT: ex PRIMESFROMSIEVE ex prime;
PROC is probably prime = ( LONG LONG INT n )BOOL:
IF n < 3 THEN n = 2
ELIF NOT ODD n THEN FALSE
ELIF n MOD 3 = 0 THEN n = 3
ELIF n MOD 5 = 0 THEN n = 5
ELIF n MOD 7 = 0 THEN n = 7
ELIF n MOD 11 = 0 THEN n = 11
ELIF n < 10 000 000 THEN
INT short n = SHORTEN SHORTEN n;
INT f := 13;
INT f2 := 169;
INT to next := 56;
BOOL is mr prime := TRUE;
WHILE f2 <= n AND is mr prime DO
is mr prime := short n MOD f /= 0;
f +:= 2;
f2 +:= to next;
to next +:= 8
OD;
is mr prime
ELSE
# miller rabin primality test #
PROC mr pow mod = (LONG LONG INT b, in e, modulus)LONG LONG INT:
BEGIN
LONG LONG INT sq := b, e := in e;
LONG LONG INT result := IF ODD e THEN b ELSE 1 FI;
e OVERAB 2;
WHILE e /= 0 DO
sq := sq * sq MOD modulus;
IF ODD e THEN result := result * sq MOD modulus FI ;
e OVERAB 2
OD;
result
END # mr pow mod # ;
INT k = 10;
LONG LONG INT d := n - 1;
INT s := 0;
WHILE NOT ODD d DO
d OVERAB 2;
s +:= 1
OD;
BOOL is mr prime := TRUE;
TO k WHILE is mr prime DO
LONG LONG INT a := 2 + ENTIER ( random * ( n - 3 ) );
LONG LONG INT x := mr pow mod( a, d, n );
IF x /= 1 THEN
BOOL done := FALSE;
TO s WHILE NOT done DO
IF x = n - 1
THEN done := TRUE
ELSE x := x * x MOD n
FI
OD;
IF NOT done THEN IF x /= n-1 THEN is mr prime := FALSE FI FI
FI
OD;
is mr prime
FI # is probably prime # ;
# returns an approximate size of sieve required for n primes #
# uses Chebyshev's formula of n/ln(n), increased by 20% #
# as n/ln(n) is only pi(n) when n appraoches infinity #
OP APPROXIMATESIEVESIZEFOR = ( INT n )INT:
BEGIN
INT result := 10;
WHILE INT asp count = ENTIER ( ( result / ln( result ) ) * 1.2 );
asp count < n
DO
result *:= 4
OD;
result
END # APPROXIMATESIEVESIZEFOR # ;
IF FALSE THEN # avoid "not used" warnings from ALGOL 68G #
[]BOOL qzqxz = PRIMESIEVE 20;
[]INT qzqyz = EXTRACTPRIMESUPTO 10 FROMPRIMESIEVE qzqxz;
print( ( APPROXIMATESIEVESIZEFOR qzqyz[ 1 ] ) );
print( ( is probably prime( 3 ) ) );
print( ( EXTRACTFIRST 3 ) )
FI;
# END primes.incl.a68 #