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.

# 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 < 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                                                        #