# Penta-power prime seeds

*is a*

**Penta-power prime seeds****draft**programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Generate the sequence of penta-power prime seeds: positive integers **n** such that:

nand^{0}+ n + 1, n^{1}+ n + 1, n^{2}+ n + 1, n^{3}+ n + 1nare all prime.^{4}+ n + 1

- Task

- Find and display the first thirty penta-power prime seeds. (Or as many as are reasonably supported by your languages math capability if it is less.)

- Stretch

- Find and display the position and value of first with a value greater than ten million.

- See also

*I can find no mention or record of this sequence anywhere. Perhaps I've invented it.*

## ALGOL 68

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 n+2 and 2n+1 numbers - Miller Rabin is used for n^2+n+1 etc. that are larger than the sieve.

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

NB: The source of the ALGOL 68-primes library is on a Rosetta Code code page linked from the above.

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

.

```
BEGIN # find some Penta power prime seeds, numbers n such that: #
# n^p + n + 1 is prime for p = 0. 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, note n^0 + n + 1 = n + 2 so n must be odd #
BOOL finished := FALSE;
INT count := 0;
INT next limit := 1 000 000;
INT limit increment = next limit;
[ 10 ]INT first, limit, index;
INT l count := 0;
print( ( "First 30 Penta power prime seeds:", newline ) );
FOR n BY 2 WHILE NOT finished DO
IF prime[ n + 2 ] THEN
IF INT n1 = n + 1;
prime[ n + n1 ]
THEN
# n^0 + n + 1 and n^1 + n + 1 are prime #
LONG INT np := LENG n * LENG 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 #
first[ l count +:= 1 ] := n;
limit[ l count ] := next limit;
index[ l count ] := count;
next limit +:= limit increment;
finished := l count >= UPB first
FI;
IF count <= 30 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
FI
OD;
print( ( newline ) );
FOR i TO UPB first DO
print( ( "First element over ", whole( limit[ i ], -10 )
, ": ", whole( first[ i ], -10 )
, ", index: ", whole( index[ i ], -6 )
, newline
)
)
OD
END
```

- Output:

First 30 Penta power prime seeds: 1 5 69 1665 2129 25739 29631 62321 77685 80535 82655 126489 207285 211091 234359 256719 366675 407945 414099 628859 644399 770531 781109 782781 923405 1121189 1158975 1483691 1490475 1512321 First element over 1000000: 1121189, index: 26 First element over 2000000: 2066079, index: 39 First element over 3000000: 3127011, index: 47 First element over 4000000: 4059525, index: 51 First element over 5000000: 5279175, index: 59 First element over 6000000: 6320601, index: 63 First element over 7000000: 7291361, index: 68 First element over 8000000: 8334915, index: 69 First element over 9000000: 9100671, index: 71 First element over 10000000: 10347035, index: 72

## Factor

```
USING: grouping io kernel lists lists.lazy math math.functions
math.primes prettyprint tools.memory.private ;
: seed? ( n -- ? )
5 [ dupd ^ 1 + + prime? ] with all-integers? ;
: pentas ( -- list )
1 lfrom [ seed? ] lfilter [ commas ] lmap-lazy ;
"First thirty penta-power prime seeds:" print
30 pentas ltake list>array 5 group simple-table.
```

- Output:

First thirty penta-power prime seeds: 1 5 69 1,665 2,129 25,739 29,631 62,321 77,685 80,535 82,655 126,489 207,285 211,091 234,359 256,719 366,675 407,945 414,099 628,859 644,399 770,531 781,109 782,781 923,405 1,121,189 1,158,975 1,483,691 1,490,475 1,512,321

## jq

The specified tasks are beyond the capabilties of the current implementations of jq:

- The C implementation (jq) does not support unbounded-precision integer arithmetic, and so is in effect only capable of generating the first few penta-power prime (ppp) seeds.

- The Go implementation (gojq) does support unbounded-precision integer arithmetic and the following program has been used to generate the first 13 ppp seeds, but gojq takes a long while to do so; other variants of the program have been considered but they run into gojq's memory-management limitations.

The program given below may nevertheless be of some interest as it illustrates how a JSON dictionary can be used to cache the primes up to a certain limit, and how this can be done using a sieve economically:

```
# Create a dictionary for the primes up to .
def primeDictionary:
# The array we use for the sieve only stores information for the odd integers greater than 1:
# index integer
# 0 3
# k 2*k + 3
# So if we wish to mark m = 2*k + 3, the relevant index is: m - 3 / 2
def ix:
if . % 2 == 0 then null
else ((. - 3) / 2)
end;
# erase(i) sets .[i*j] to false for odd integral j > i, and assumes i is odd
def erase(i):
((i - 3) / 2) as $k
# Consider relevant multiples:
| (((length * 2 + 3) / i)) as $upper
# ... only consider odd multiples from i onwards
| reduce range(i; $upper; 2) as $j (.;
(((i * $j) - 3) / 2) as $m
| if .[$m] then .[$m] = false else . end);
if . < 2 then {}
else (. + 1) as $n
| (($n|sqrt) / 2) as $s
| [range(3; $n; 2)|true]
| reduce (1 + (2 * range(1; $s)) ) as $i (.; erase($i))
| . as $sieve
| reduce (2, (range(3; $n; 2) | select($sieve[ix]))) as $i ({}; .[$i|tostring]=$i)
end ;
```

```
# Input should be an integer
def isPrime:
. as $n
| if ($n < 2) then false
elif ($n % 2 == 0) then $n == 2
elif ($n % 3 == 0) then $n == 3
else 5
| until( . <= 0;
if .*. > $n then -1
elif ($n % . == 0) then 0
else . + 2
| if ($n % . == 0) then 0
else . + 4
end
end)
| . == -1
end;
# $primedictionary should be a dictionary of primes up to $small
def ispentapowerprime($primedictionary; $small):
def isp: if . <= $small then $primedictionary[tostring] else isPrime end;
. as $n
| (. * .) as $n2
| (. * $n2) as $n3
| all($n + 2, $n + $n + 1, $n2 + $n + 1, $n3 + $n + 1, $n3 * $n + $n + 1; isp);
# Output: a stream of the first $count penta-power prime-seeds
# The size of the dictionary has been chosen with gojq in mind.
def ppprimes($count):
# The size of primeDictionary has been chosen with gojq's limitations in mind
($count | .*.*. | primeDictionary) as $pd
| limit($count; 1, 2, range(3; infinite; 2) | select(ispentapowerprime($pd; $small)) );
ppprimes(30)
```

- Output:

Using gojq, progress effectively grinds to a halt after 207285.

1 5 69 1665 2129 25739 29631 62321 77685 80535 82655 126489 207285 <terminated>

## Perl

```
use v5.36;
use bigint;
use ntheory 'is_prime';
use List::Util 'max';
sub comma { reverse ((reverse shift) =~ s/(.{3})/$1,/gr) =~ s/^,//r }
sub table ($c, @V) { my $t = $c * (my $w = 2 + max map {length} @V); ( sprintf( ('%'.$w.'s')x@V, @V) ) =~ s/.{1,$t}\K/\n/gr }
my($i,@ppps);
while (@ppps < 30) {
my $k = 1 + (my $n = 2 * $i++ + 1);
push @ppps, comma $n if
is_prime( 1 + $k) and
is_prime($n + $k) and
is_prime($n**2 + $k) and
is_prime($n**3 + $k) and
is_prime($n**4 + $k);
}
say 'First thirty penta-power prime seeds:';
say table(10,@ppps);
```

- Output:

First thirty penta-power prime seeds: 1 5 69 1,665 2,129 25,739 29,631 62,321 77,685 80,535 82,655 126,489 207,285 211,091 234,359 256,719 366,675 407,945 414,099 628,859 644,399 770,531 781,109 782,781 923,405 1,121,189 1,158,975 1,483,691 1,490,475 1,512,321

## Phix

with javascript_semantics include mpfr.e mpz {p,q} = mpz_inits(2) function isPentaPowerPrimeSeed(integer n) -- (we already know n+2 is prime) mpz_set_si(p,n) for i=1 to 4 do if i>1 then mpz_mul_si(p,p,n) end if mpz_add_ui(q,p,n+1) if not mpz_prime(q) then return false end if end for return true end function sequence ppps = {} integer pn = 1, m = 1, c = 0 string l = "" while m<=10 do integer n = get_prime(pn)-2 if isPentaPowerPrimeSeed(n) then c += 1 if c<31 then ppps &= n if c=30 then printf(1,"First thirty penta-power prime seeds:\n%s\n", {join_by(ppps,1,10," ",fmt:="%,9d")}) end if end if if n > m * 1e6 then l &= sprintf(" %5s million is %,d (the %s)\n", {ordinal(m,true), n, ordinal(c)}) m += 1 end if end if pn += 1 end while printf(1,"First penta-power prime seed greater than:\n%s",l)

- Output:

First thirty penta-power prime seeds: 1 5 69 1,665 2,129 25,739 29,631 62,321 77,685 80,535 82,655 126,489 207,285 211,091 234,359 256,719 366,675 407,945 414,099 628,859 644,399 770,531 781,109 782,781 923,405 1,121,189 1,158,975 1,483,691 1,490,475 1,512,321 First penta-power prime seed greater than: one million is 1,121,189 (the twenty-sixth) two million is 2,066,079 (the thirty-ninth) three million is 3,127,011 (the forty-seventh) four million is 4,059,525 (the fifty-first) five million is 5,279,175 (the fifty-ninth) six million is 6,320,601 (the sixty-third) seven million is 7,291,361 (the sixty-eighth) eight million is 8,334,915 (the sixty-ninth) nine million is 9,100,671 (the seventy-first) ten million is 10,347,035 (the seventy-second)

## Python

```
from sympy import isprime
def ispentapowerprime(n):
return all(isprime(i) for i in [n + 2, n + n + 1, n**2 + n + 1, n**3 + n + 1, n**4 + n + 1])
ppprimes = [i for i in range(10_400_000) if ispentapowerprime(i)]
for i in range(50):
print(f'{ppprimes[i]: 11,}', end='\n' if (i + 1) % 10 == 0 else '')
for n in range(1_000_000, 10_000_001, 1_000_000):
proot = next(filter(lambda x: x > n, ppprimes))
print(f'The first penta-power prime seed over {n:,} is {proot:,}')
```

- Output:

1 5 69 1,665 2,129 25,739 29,631 62,321 77,685 80,535 82,655 126,489 207,285 211,091 234,359 256,719 366,675 407,945 414,099 628,859 644,399 770,531 781,109 782,781 923,405 1,121,189 1,158,975 1,483,691 1,490,475 1,512,321 1,711,991 1,716,989 1,780,485 1,791,041 1,835,589 1,860,011 1,861,259 1,980,441 2,066,079 2,211,705 2,215,529 2,271,009 2,413,265 2,514,161 2,915,109 2,940,405 3,127,011 3,587,319 3,890,769 3,992,379 The first penta-power prime seed over 1,000,000 is 1,121,189 The first penta-power prime seed over 2,000,000 is 2,066,079 The first penta-power prime seed over 3,000,000 is 3,127,011 The first penta-power prime seed over 4,000,000 is 4,059,525 The first penta-power prime seed over 5,000,000 is 5,279,175 The first penta-power prime seed over 6,000,000 is 6,320,601 The first penta-power prime seed over 7,000,000 is 7,291,361 The first penta-power prime seed over 8,000,000 is 8,334,915 The first penta-power prime seed over 9,000,000 is 9,100,671 The first penta-power prime seed over 10,000,000 is 10,347,035

## Raku

```
use Lingua::EN::Numbers;
my @ppps = lazy (^∞).hyper(:5000batch).map(* × 2 + 1).grep: -> \n { my \k = n + 1; (1+k).is-prime && (n+k).is-prime && (n²+k).is-prime && (n³+k).is-prime && (n⁴+k).is-prime }
say "First thirty penta-power prime seeds:\n" ~ @ppps[^30].batch(10)».&comma».fmt("%9s").join: "\n";
say "\nFirst penta-power prime seed greater than:";
for 1..10 {
my $threshold = Int(1e6 × $_);
my $key = @ppps.first: * > $threshold, :k;
say "{$threshold.&cardinal.fmt: '%13s'} is the {ordinal-digit $key + 1}: {@ppps[$key].&comma}";
}
```

- Output:

First thirty penta-power prime seeds: 1 5 69 1,665 2,129 25,739 29,631 62,321 77,685 80,535 82,655 126,489 207,285 211,091 234,359 256,719 366,675 407,945 414,099 628,859 644,399 770,531 781,109 782,781 923,405 1,121,189 1,158,975 1,483,691 1,490,475 1,512,321 First penta-power prime seed greater than: one million is the 26th: 1,121,189 two million is the 39th: 2,066,079 three million is the 47th: 3,127,011 four million is the 51st: 4,059,525 five million is the 59th: 5,279,175 six million is the 63rd: 6,320,601 seven million is the 68th: 7,291,361 eight million is the 69th: 8,334,915 nine million is the 71st: 9,100,671 ten million is the 72nd: 10,347,035

## Wren

```
import "./gmp" for Mpz
import "./fmt" for Fmt
var p = Mpz.new()
var q = Mpz.one
var isPentaPowerPrimeSeed = Fn.new { |n|
p.setUi(n)
var k = n + 1
return (q + k).probPrime(15) > 0 &&
(p + k).probPrime(15) > 0 &&
(p.mul(n) + k).probPrime(15) > 0 &&
(p.mul(n) + k).probPrime(15) > 0 &&
(p.mul(n) + k).probPrime(15) > 0
}
var ppps = []
var n = 1
while (ppps.count < 30) {
if (isPentaPowerPrimeSeed.call(n)) ppps.add(n)
n = n + 2 // n must be odd
}
System.print("First thirty penta-power prime seeds:")
Fmt.tprint("$,9d", ppps, 10)
System.print("\nFirst penta-power prime seed greater than:")
n = 1
var m = 1
var c = 0
while (true) {
if (isPentaPowerPrimeSeed.call(n)) {
c = c + 1
if (n > m * 1e6) {
Fmt.print(" $2d million is the $r: $,10d", m, c, n)
m = m + 1
if (m == 11) return
}
}
n = n + 2
}
```

- Output:

First thirty penta-power prime seeds: 1 5 69 1,665 2,129 25,739 29,631 62,321 77,685 80,535 82,655 126,489 207,285 211,091 234,359 256,719 366,675 407,945 414,099 628,859 644,399 770,531 781,109 782,781 923,405 1,121,189 1,158,975 1,483,691 1,490,475 1,512,321 First penta-power prime seed greater than: 1 million is the 26th: 1,121,189 2 million is the 39th: 2,066,079 3 million is the 47th: 3,127,011 4 million is the 51st: 4,059,525 5 million is the 59th: 5,279,175 6 million is the 63rd: 6,320,601 7 million is the 68th: 7,291,361 8 million is the 69th: 8,334,915 9 million is the 71st: 9,100,671 10 million is the 72nd: 10,347,035