The name   cuban   has nothing to do with   Cuba  (the country),   but has to do with the
fact that cubes   (3rd powers)   play a role in its definition.



Some definitions of cuban primes

   primes which are the difference of two consecutive cubes.
   primes of the form:   (n+1)3 - n3.
   primes of the form:   n3 - (n-1)3.
   primes   p   such that   n2(p+n)   is a cube for some   n>0.
   primes   p   such that   4p = 1 + 3n2.


Cuban primes were named in 1923 by Allan Joseph Champneys Cunningham.



Task requirements

   show the first   200   cuban primes   (in a multi─line horizontal format).
   show the   100,000th   cuban prime.
   show all cuban primes with commas   (if appropriate).
   show all output here.


Note that   cuban prime   isn't capitalized   (as it doesn't refer to the nation of Cuba).



Also see

   Wikipedia entry:     cuban prime.
   MathWorld entry:   cuban prime.
   The OEIS entry:     A002407.     The   100,000th   cuban prime can be verified in the   2nd   example   on this OEIS web page.


BEGIN
    # find some cuban primes (using the definition: a prime p is a cuban prime if   #
    #                                                   p = n^3 - ( n - 1 )^3       #
    #                                                   for some n > 0)             #
 
    # returns a string representation of n with commas                              #
    PROC commatise = ( LONG INT n )STRING:
         BEGIN
            STRING result      := "";
            STRING unformatted  = whole( n, 0 );
            INT    ch count    := 0;
            FOR c FROM UPB unformatted BY -1 TO LWB unformatted DO
                IF   ch count <= 2 THEN ch count +:= 1
                ELSE                    ch count  := 1; "," +=: result
                FI;
                unformatted[ c ] +=: result
            OD;
            result
         END # commatise # ;
 
    # sieve the primes                                                              #
    INT sieve max = 2 000 000;
    [ sieve max ]BOOL sieve; FOR i TO UPB sieve DO sieve[ i ] := TRUE OD;
    sieve[ 1 ] := FALSE;
    FOR s FROM 2 TO ENTIER sqrt( sieve max ) DO
        IF sieve[ s ] THEN
            FOR p FROM s * s BY s TO sieve max DO sieve[ p ] := FALSE OD
        FI
    OD;
    # count the primes, we can ignore 2, as we know it isn't a cuban prime          #
    sieve[ 2 ] := FALSE;
    INT prime count := 0;
    FOR s TO UPB sieve DO IF sieve[ s ] THEN prime count +:= 1 FI OD;
    # construct a list of the primes                                                #
    [ 1 : prime count ]INT primes;
    INT prime pos := LWB primes;
    FOR s FROM LWB sieve TO UPB sieve DO
        IF sieve[ s ] THEN primes[ prime pos ] := s; prime pos +:= 1 FI
    OD;
 
    # find the cuban primes                                                         #
    INT       cuban count   := 0;
    LONG INT  final cuban   := 0;
    INT       max cuban      = 100 000; # mximum number of cubans to find           #
    INT       print limit    =     200; # show all cubans up to this one            #
    print( ( "First ", commatise( print limit ), " cuban primes: ", newline ) );
    LONG INT prev cube      := 1;
    FOR n FROM 2 WHILE
        LONG INT this cube   = ( LENG n * n ) * n;
        LONG INT p           = this cube - prev cube;
        prev cube           := this cube;
        IF ODD p THEN
            # 2 is not a cuban prime so we only test odd numbers                    #
            BOOL is prime   := TRUE;
            INT max factor   = SHORTEN ENTIER long sqrt( p );
            FOR f FROM LWB primes WHILE is prime AND primes[ f ] <= max factor DO
                is prime    := p MOD primes[ f ] /= 0
            OD;
            IF is prime THEN
                # have a cuban prime                                                #
                cuban count +:= 1;
                IF cuban count <= print limit THEN
                    # must show this cuban                                          #
                    STRING p formatted = commatise( p );
                    print( ( "          "[ UPB p formatted : ], p formatted ) );
                    IF cuban count MOD 10 = 0 THEN print( ( newline ) ) FI
                FI;
                final cuban := p
            FI
        FI;
        cuban count < max cuban
    DO SKIP OD;
    IF cuban count MOD 10 /= 0 THEN print( ( newline ) ) FI;
    print( ( "The ", commatise( max cuban ), " cuban prime is: ", commatise( final cuban ), newline ) )
END