Talk:ALGOL 68-primes

Source code

Note the is probably prime PROC is based on code from The Miller–Rabin primality test task and the prelude.

<lang algol68># 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 prime;
           prime[ 0 ] := prime[ 1 ] := FALSE;
           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 = 2      THEN TRUE
        ELIF NOT ODD n  THEN FALSE
        ELIF n < 3      THEN FALSE
        ELIF n < 10 000 THEN
            # smallish number = use trial division #
            BOOL is prime := TRUE;
            FOR i FROM 3 BY 2 TO ENTIER sqrt( SHORTEN SHORTEN n ) WHILE is prime := n MOD i /= 0 DO SKIP OD;
            is 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 prime := TRUE;
           TO k WHILE is 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 prime := FALSE FI FI
               FI
           OD;
           is 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 primes = ENTIER ( ( result / ln( result ) ) * 1.2 );
                  primes < n
            DO
                result *:= 4
            OD;
            result
        END # APPROXIMATESIEVESIZEFOR # ;

1. END primes.incl.a68 #</lang>