Quad-power prime seeds: Difference between revisions

=={{header|ALGOL 68}}==

{{works with|ALGOL 68G|Any - tested with release 3.0.3 under Windopws}}

This uses ALGOL 68G's LONG LONG INT during the Miller Rabin testing, under ALGOL 68G version three, the default precision of LONG LONG INT is 72 digits and LONG INT is 128 bit.

A sieve is used to (hopefully) quickly eliminate non-prime 2n+1 numbers - Miller Rabin is used for n^2+n+1 etc. that are larger than the sieve.

This is about 10 times slower than the equivalent Penta-powwr prime seed program, possibly because even numbers are included and the n+2 test in the Penta-powers eliminates more numbers before the higher powers must be calculated.

<br>

Note that to run this with ALGOL 68G under Windows (and probably Linux) a large heap sie must be specified on the command line, e.g. <code>-heap 1024M</code>.

<lang algol68>BEGIN # find some Quad power prime seeds, numbers n such that: #

# n^p + n + 1 is prime for p = 1, 2, 3, 4 #

PR read "primes.incl.a68" PR # include prime utilities #

INT max prime = 22 000 000;

# sieve the primes to max prime - 22 000 000 should be enough to allow #

# checking primality of 2n+1 up to a little under 11 000 000 which should #

# be enough for the task #

[ 0 : max prime ]BOOL prime;

FOR i FROM LWB prime TO UPB prime DO prime[ i ] := FALSE OD;

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;

# returns TRUE if p is (probably) prime, FALSE otherwise #

# uses the sieve if possible, Miller Rabin otherwise #

PROC is prime = ( LONG INT p )BOOL:

IF p <= max prime

THEN prime[ SHORTEN p ]

ELSE is probably prime( p )

FI;

# attempt to find the numbers #

BOOL finished := FALSE;

INT count := 0;

INT next limit := 1 000 000;

INT last limit = 10 000 000;

INT limit increment = next limit;

print( ( "First 50 Quad power prime seeds:", newline ) );

FOR n WHILE NOT finished DO

IF INT n1 = n + 1;

prime[ n + n1 ]

THEN

# n^0 + n + 1 and n^1 + n + 1 are prime #

LONG INT np := n * n;

IF is prime( np + n1 ) THEN

# n^2 + n + 1 is prime #

IF is prime( ( np *:= n ) + n1 ) THEN

# n^3 + n + 1 is prime #

IF is prime( ( np * n ) + n1 ) THEN

# n^4 + n + 1 is prime - have a suitable number #

count +:= 1;

IF n > next limit THEN

# found the first element over the next limit #

print( ( newline

, "First element over ", whole( next limit, -10 )

, ": ", whole( n, -10 )

, ", index: ", whole( count, -6 )

)

);

next limit +:= limit increment;

finished := next limit > last limit

FI;

IF count <= 50 THEN

# found one of the first 30 numbers #

print( ( " ", whole( n, -8 ) ) );

IF count MOD 10 = 0 THEN print( ( newline ) ) FI

FI

FI

FI

FI

FI

OD

END</lang>

{{out}}

<pre>

First 50 Quad power prime seeds:

1 2 5 6 69 131 426 1665 2129 2696

5250 7929 9689 13545 14154 14286 16434 19760 25739 27809

29631 36821 41819 46619 48321 59030 60500 61955 62321 73610

77685 79646 80535 82655 85251 86996 91014 96566 97739 105939

108240 108681 119754 122436 123164 126489 140636 150480 153179 163070

First element over 1000000: 1009286, index: 141

First element over 2000000: 2015496, index: 234

First element over 3000000: 3005316, index: 319

First element over 4000000: 4004726, index: 383

First element over 5000000: 5023880, index: 452

First element over 6000000: 6000554, index: 514

First element over 7000000: 7047129, index: 567

First element over 8000000: 8005710, index: 601

First element over 9000000: 9055151, index: 645

First element over 10000000: 10023600, index: 701

</pre>