Cuban primes: Difference between revisions

→‎{{header|ALGOL 68}}: Use Miller Rabin for primality tesing
(→‎{{header|ALGOL 68}}: Use Miller Rabin for primality tesing)
Line 38:
# p = n^3 - ( n - 1 )^3 #
# for some n > 0) #
PR read "primes.incl.a68" PR # include prime utilities #
 
# returns a string representation of n with commas #
PROC commatise = ( LONG INT n )STRING:
Line 54:
END # commatise # ;
 
# sieve the primes #
INT sieve max = 2 000 000;
[]BOOL sieve max= ]BOOLPRIMESIEVE sieve max; FOR i TO UPB sieve DO # sieve[ ithe ]primes :=to TRUEmax OD;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 #
Line 79 ⟶ 62:
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
Line 87 ⟶ 70:
IF ODD p THEN
# 2 is not a cuban prime so we only test odd numbers #
BOOLIF IF p <= UPB sieve THEN sieve[ SHORTEN p ] ELSE is probably prime( :=p ) TRUE;FI
IF p <= UPB sieve 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 #
IF ( cuban count +:= 1 ) <= print limit THEN
Line 115 ⟶ 90:
{{out}}
<pre>
First 200 cuban primes:
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
3,043

edits