Cuban primes: Difference between revisions

Content added Content deleted
(→‎{{header|ALGOL 68}}: Use Miller Rabin for primality tesing)
Line 38: Line 38:
# p = n^3 - ( n - 1 )^3 #
# p = n^3 - ( n - 1 )^3 #
# for some n > 0) #
# for some n > 0) #
PR read "primes.incl.a68" PR # include prime utilities #

# returns a string representation of n with commas #
# returns a string representation of n with commas #
PROC commatise = ( LONG INT n )STRING:
PROC commatise = ( LONG INT n )STRING:
Line 54: Line 54:
END # commatise # ;
END # commatise # ;


# sieve the primes #
INT sieve max = 2 000 000;
INT sieve max = 2 000 000;
[ sieve max ]BOOL sieve; FOR i TO UPB sieve DO sieve[ i ] := TRUE OD;
[]BOOL sieve = PRIMESIEVE sieve max; # sieve the primes to max sieve #
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 - 1;
FOR s FROM LWB sieve TO UPB sieve DO
IF sieve[ s ] THEN primes[ prime pos +:= 1 ] := s FI
OD;


# find the cuban primes #
# find the cuban primes #
Line 79: Line 62:
INT max cuban = 100 000; # mximum number of cubans to find #
INT max cuban = 100 000; # mximum number of cubans to find #
INT print limit = 200; # show all cubans up to this one #
INT print limit = 200; # show all cubans up to this one #
print( ( "First ", commatise( print limit ), " cuban primes: ", newline ) );
print( ( "First ", commatise( print limit ), " cuban primes:", newline ) );
LONG INT prev cube := 1;
LONG INT prev cube := 1;
FOR n FROM 2 WHILE
FOR n FROM 2 WHILE
Line 87: Line 70:
IF ODD p THEN
IF ODD p THEN
# 2 is not a cuban prime so we only test odd numbers #
# 2 is not a cuban prime so we only test odd numbers #
BOOL is prime := TRUE;
IF IF p <= UPB sieve THEN sieve[ SHORTEN p ] ELSE is probably prime( p ) FI
IF p <= UPB sieve THEN
THEN
is prime := sieve[ SHORTEN p ]
ELSE
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
FI;
IF is prime THEN
# have a cuban prime #
# have a cuban prime #
IF ( cuban count +:= 1 ) <= print limit THEN
IF ( cuban count +:= 1 ) <= print limit THEN
Line 115: Line 90:
{{out}}
{{out}}
<pre>
<pre>
First 200 cuban primes:
First 200 cuban primes:
7 19 37 61 127 271 331 397 547 631
7 19 37 61 127 271 331 397 547 631
919 1,657 1,801 1,951 2,269 2,437 2,791 3,169 3,571 4,219
919 1,657 1,801 1,951 2,269 2,437 2,791 3,169 3,571 4,219