Talk:ALGOL 68-primes: Difference between revisions
→Source code: primes.incl.a68 - horrendous bug fix
(Created a page with the source code for ALGOL 68-primes) |
(→Source code: primes.incl.a68 - horrendous bug fix) |
||
| (9 intermediate revisions by the same user not shown) | |||
Line 1:
===Source code===
Note the is probably prime PROC is based on code from [https://www.rosettacode.org/wiki/Miller%E2%80%93Rabin_primality_test The Miller–Rabin primality test task] and the [https://www.rosettacode.org/wiki/ALGOL_68/prelude prelude].
<syntaxhighlight lang=algol68>
# primes.incl.a68: prime related operators, procedure etc. #
# returns a sieve of primes up to n #
OP PRIMESIEVE = ( INT n )[]BOOL:
BEGIN
[
prime[
prime[ 2 ] := TRUE;
FOR i FROM 3 BY 2 TO UPB prime DO prime[ i ] := TRUE OD;
FOR i FROM 4 BY 2 TO UPB prime DO prime[ i ] := FALSE OD;
FOR i FROM 3 BY 2 TO ENTIER sqrt( UPB prime ) DO
IF prime[ i ] THEN
FOR s FROM i * i BY i + i TO UPB prime DO prime[ s ] := FALSE OD
FI
OD;
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 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 p count := 0;
FOR i TO UPB prime WHILE p count < n OF ex DO
IF prime[ i ] THEN p count +:= 1 FI
OD;
# construct a list of the first n OF ex primes or p count, #
# whichever is smaller #
[ 1 : p count ]INT low prime;
INT low pos := 0;
FOR i WHILE low pos < UPB low prime DO
IF prime[ i ] THEN low prime[ low pos +:= 1 ] := i FI
OD;
low prime
ELSE
# extract primes up to n OF ex primes #
INT max prime = n OF ex;
# count the primes up to max prime #
INT p count := 0;
FOR i TO max prime DO IF prime[ i ] THEN p count +:= 1 FI OD;
# construct a list of the primes up to max prime #
[ 1 : p count ]INT low prime;
INT low pos := 0;
FOR i WHILE low pos < UPB low prime DO
IF prime[ i ] THEN low prime[ low pos +:= 1 ] := i FI
OD;
low prime
FI # PRIMESFROMSIEVE # ;
# alternative spelling of the operator name #
OP FROMPRIMESIEVE = ( PRIMEEXTRACT ex, []BOOL prime )[]INT: ex PRIMESFROMSIEVE 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 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 # 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 := 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 prime count = ENTIER ( ( result / ln( result ) ) * 1.2 );
prime count < n
DO
result *:= 4
OD;
result
END # APPROXIMATESIEVESIZEFOR # ;
# END primes.incl.a68 #
</syntaxhighlight>
| |||