Quad-power prime seeds
Generate the sequence of quad-power prime seeds: positive integers n such that:
Quad-power prime seeds is a 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.
n1 + n + 1, n2 + n + 1, n3 + n + 1 and n4 + n + 1 are all prime.
- Task
- Find and display the first fifty quad-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 one million, two million, three million.
- See also
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 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.
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 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^1 + n + 1 is 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 #
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- Output:
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
Arturo
quadPowerPrime?: function [n]->
every? [[n+n+1] [1+n+n^2] [1+n+n^3] [1+n+n^4]] 'x ->
prime? do x
first50qpps: select.first:50 1..∞ => quadPowerPrime?
loop split.every: 10 first50qpps 'x ->
print map x 's -> pad to :string s 7
- Output:
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
C
#include <stdio.h>
#include <stdbool.h>
#include <locale.h>
#include <gmp.h>
mpz_t p, p2;
bool isQuadPowerPrimeSeed(unsigned int n) {
int i;
mpz_set_ui(p, n);
unsigned int k = n + 1;
mpz_add_ui(p2, p, k);
if (!mpz_probab_prime_p(p2, 15)) return false;
for (i = 0; i < 3; ++i) {
mpz_mul_ui(p, p, n);
mpz_set(p2, p);
mpz_add_ui(p2, p2, k);
if (!mpz_probab_prime_p(p2, 15)) return false;
}
return true;
}
const char *ord(int c) {
int m = c % 100;
if (m >= 4 && m <= 20) return "th";
m %= 10;
return (m == 1) ? "st" :
(m == 2) ? "nd" :
(m == 3) ? "rd" : "th";
}
int main() {
unsigned int n;
int c = 0, m = 1;
mpz_init(p);
mpz_init(p2);
setlocale(LC_NUMERIC, "");
printf("First fifty quad-power prime seeds:\n");
for (n = 1; c < 50; ++n) {
if (isQuadPowerPrimeSeed(n)) {
printf("%'7u ", n);
if (!((++c) % 10)) printf("\n");
}
}
printf("\nFirst quad-power prime seed greater than:\n");
while (1) {
if (isQuadPowerPrimeSeed(n)) {
++c;
if (n > 1000000 * m) {
printf(" %2d million is the %d%s: %'10u\n", m, c, ord(c), n);
if (++m == 11) break;
}
}
++n;
}
return 0;
}
- Output:
First fifty quad-power prime seeds:
1 2 5 6 69 131 426 1,665 2,129 2,696
5,250 7,929 9,689 13,545 14,154 14,286 16,434 19,760 25,739 27,809
29,631 36,821 41,819 46,619 48,321 59,030 60,500 61,955 62,321 73,610
77,685 79,646 80,535 82,655 85,251 86,996 91,014 96,566 97,739 105,939
108,240 108,681 119,754 122,436 123,164 126,489 140,636 150,480 153,179 163,070
First quad-power prime seed greater than:
1 million is the 141st: 1,009,286
2 million is the 234th: 2,015,496
3 million is the 319th: 3,005,316
4 million is the 383rd: 4,004,726
5 million is the 452nd: 5,023,880
6 million is the 514th: 6,000,554
7 million is the 567th: 7,047,129
8 million is the 601st: 8,005,710
9 million is the 645th: 9,055,151
10 million is the 701st: 10,023,600
F#
// Quad-power prime seeds. Nigel Galloway: August 22nd., 2022
let fG n g=let n=bigint(n:int) in let n=n**g+n+1I in Open.Numeric.Primes.MillerRabin.IsProbablePrime &n
let fN(n,g)=Seq.initInfinite((+)n)|>Seq.filter(fun n->let g=fG n in g 1&&g 2&&g 3&&g 4)|>Seq.mapi(fun n g->(n,g))|>Seq.find(snd>>(<)g)
Seq.initInfinite((+)1)|>Seq.filter(fun n->let g=fG n in g 1&&g 2&&g 3&&g 4)|>Seq.take 50|>Seq.iter(printf "%d "); printfn "\n"
[1000000..1000000..10000000]|>Seq.scan(fun(n,g,x) l->let i,e=fN(g,l) in (n+i,e,l))(0,0,0)|>Seq.skip 1|>Seq.iter(fun(n,g,l)->printfn $"First element over %8d{l} is %9d{g} at index %3d{n}")
- Output:
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 is 1009286 at index 140 First element over 2000000 is 2015496 at index 233 First element over 3000000 is 3005316 at index 318 First element over 4000000 is 4004726 at index 382 First element over 5000000 is 5023880 at index 451 First element over 6000000 is 6000554 at index 513 First element over 7000000 is 7047129 at index 566 First element over 8000000 is 8005710 at index 600 First element over 9000000 is 9055151 at index 644 First element over 10000000 is 10023600 at index 700
Factor
USING: grouping io kernel lists lists.lazy math math.functions
math.primes prettyprint sequences tools.memory.private ;
: seed? ( n -- ? )
{ 1 2 3 4 } [ dupd ^ 1 + + prime? ] with all? ;
: quads ( -- list )
1 lfrom [ seed? ] lfilter [ commas ] lmap-lazy ;
"First fifty quad-power prime seeds:" print
50 quads ltake list>array 10 group simple-table.
- Output:
First fifty quad-power prime seeds: 1 2 5 6 69 131 426 1,665 2,129 2,696 5,250 7,929 9,689 13,545 14,154 14,286 16,434 19,760 25,739 27,809 29,631 36,821 41,819 46,619 48,321 59,030 60,500 61,955 62,321 73,610 77,685 79,646 80,535 82,655 85,251 86,996 91,014 96,566 97,739 105,939 108,240 108,681 119,754 122,436 123,164 126,489 140,636 150,480 153,179 163,070
FreeBASIC
#include "isprime.bas"
Function isQuadPowerPrime(n As Uinteger) As Boolean
Dim As Uinteger p(n)
Dim As Uinteger k = n + 1
Return isPrime(n+n+1) And isPrime((n*n)+(n+1)) And _
isPrime((n*n*n)+(n+1)) And isPrime((n^4)+(n+1))
End Function
Dim As Uinteger n = 1, cont = 1
Dim As Double t0 = Timer
Print "First fifty quad-power prime seeds:"
While cont < 50
If isQuadPowerPrime(n) Then
Print Using "###,###"; n;
If cont Mod 10 = 0 Then Print
cont += 1
End If
n += 1
Wend
Print: Print Timer - t0
Print !"\nFirst quad-power prime seed greater than:"
Dim As Uinteger m = 1, c = 50
Do
If isQuadPowerPrime(n) Then
c += 1
If n > m * 1e6 Then
Print Using " ## million is the ###: ##########"; m; c; n
m += 1
If m = 11 Then Exit Do
End If
End If
n += 1
Loop
SleepPerl
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($n,@qpps);
while (@qpps < 50) {
my $k = 1 + ++$n;
push @qpps, comma $n if
is_prime($n + $k) and
is_prime($n**2 + $k) and
is_prime($n**3 + $k) and
is_prime($n**4 + $k);
}
say 'First fifty quad-power prime seeds:';
say table(10,@qpps);
- Output:
First fifty quad-power prime seeds:
1 2 5 6 69 131 426 1,665 2,129 2,696
5,250 7,929 9,689 13,545 14,154 14,286 16,434 19,760 25,739 27,809
29,631 36,821 41,819 46,619 48,321 59,030 60,500 61,955 62,321 73,610
77,685 79,646 80,535 82,655 85,251 86,996 91,014 96,566 97,739 105,939
108,240 108,681 119,754 122,436 123,164 126,489 140,636 150,480 153,179 163,070
Phix
with javascript_semantics include mpfr.e mpz {p,q} = mpz_inits(2) function isQuadPowerPrimeSeed(integer n) 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 qpps = {} integer n = 1 while length(qpps)<50 do if isQuadPowerPrimeSeed(n) then qpps &= n end if n += 1 end while printf(1,"First fifty quad-power prime seeds:\n%s\n", {join_by(qpps,1,10," ",fmt:="%,7d")}) printf(1,"First quad-power prime seed greater than:\n") integer m = 1, c = 50 while m<=10 do if isQuadPowerPrimeSeed(n) then c += 1 if n > m * 1e6 then printf(1," %5s million is %,d (the %s)\n", {ordinal(m,true), n, ordinal(c)}) m += 1 end if end if n += 1 end while
- Output:
First fifty quad-power prime seeds:
1 2 5 6 69 131 426 1,665 2,129 2,696
5,250 7,929 9,689 13,545 14,154 14,286 16,434 19,760 25,739 27,809
29,631 36,821 41,819 46,619 48,321 59,030 60,500 61,955 62,321 73,610
77,685 79,646 80,535 82,655 85,251 86,996 91,014 96,566 97,739 105,939
108,240 108,681 119,754 122,436 123,164 126,489 140,636 150,480 153,179 163,070
First quad-power prime seed greater than:
one million is 1,009,286 (the one hundred and forty-first)
two million is 2,015,496 (the two hundred and thirty-fourth)
three million is 3,005,316 (the three hundred and nineteenth)
four million is 4,004,726 (the three hundred and eighty-third)
five million is 5,023,880 (the four hundred and fifty-second)
six million is 6,000,554 (the five hundred and fourteenth)
seven million is 7,047,129 (the five hundred and sixty-seventh)
eight million is 8,005,710 (the six hundred and first)
nine million is 9,055,151 (the six hundred and forty-fifth)
ten million is 10,023,600 (the seven hundred and first)
Python
""" quad-power prime root numbers """
from sympy import isprime
def isquadpowerprime(cand):
""" return if is a quad power prime root number """
return all(isprime(i) for i in
[cand + cand + 1, cand**2 + cand + 1, cand**3 + cand + 1, cand**4 + cand + 1])
qpprimes = [k for k in range(10_100_000) if isquadpowerprime(k)]
for i in range(50):
print(f'{qpprimes[i]: 9,}', end='\n' if (i + 1) % 10 == 0 else '')
for j in range(1_000_000, 10_000_001, 1_000_000):
for p in qpprimes:
if p > j:
print(f'The first quad-power prime seed over {j:,} is {p:,}')
break
- Output:
1 2 5 6 69 131 426 1,665 2,129 2,696
5,250 7,929 9,689 13,545 14,154 14,286 16,434 19,760 25,739 27,809
29,631 36,821 41,819 46,619 48,321 59,030 60,500 61,955 62,321 73,610
77,685 79,646 80,535 82,655 85,251 86,996 91,014 96,566 97,739 105,939
108,240 108,681 119,754 122,436 123,164 126,489 140,636 150,480 153,179 163,070
The first quad-power prime seed over 1,000,000 is 1,009,286
The first quad-power prime seed over 2,000,000 is 2,015,496
The first quad-power prime seed over 3,000,000 is 3,005,316
The first quad-power prime seed over 4,000,000 is 4,004,726
The first quad-power prime seed over 5,000,000 is 5,023,880
The first quad-power prime seed over 6,000,000 is 6,000,554
The first quad-power prime seed over 7,000,000 is 7,047,129
The first quad-power prime seed over 8,000,000 is 8,005,710
The first quad-power prime seed over 9,000,000 is 9,055,151
The first quad-power prime seed over 10,000,000 is 10,023,600
Raku
use Lingua::EN::Numbers;
my @qpps = lazy (1..*).hyper(:5000batch).grep: -> \n { my \k = n + 1; (n+k).is-prime && (n²+k).is-prime && (n³+k).is-prime && (n⁴+k).is-prime }
say "First fifty quad-power prime seeds:\n" ~ @qpps[^50].batch(10)».&comma».fmt("%7s").join: "\n";
say "\nFirst quad-power prime seed greater than:";
for 1..10 {
my $threshold = Int(1e6 × $_);
my $key = @qpps.first: * > $threshold, :k;
say "{$threshold.&cardinal.fmt: '%13s'} is the {ordinal-digit $key + 1}: {@qpps[$key].&comma}";
}
- Output:
First fifty quad-power prime seeds:
1 2 5 6 69 131 426 1,665 2,129 2,696
5,250 7,929 9,689 13,545 14,154 14,286 16,434 19,760 25,739 27,809
29,631 36,821 41,819 46,619 48,321 59,030 60,500 61,955 62,321 73,610
77,685 79,646 80,535 82,655 85,251 86,996 91,014 96,566 97,739 105,939
108,240 108,681 119,754 122,436 123,164 126,489 140,636 150,480 153,179 163,070
First quad-power prime seed greater than:
one million is the 141st: 1,009,286
two million is the 234th: 2,015,496
three million is the 319th: 3,005,316
four million is the 383rd: 4,004,726
five million is the 452nd: 5,023,880
six million is the 514th: 6,000,554
seven million is the 567th: 7,047,129
eight million is the 601st: 8,005,710
nine million is the 645th: 9,055,151
ten million is the 701st: 10,023,600
Ruby
require 'openssl'
quad_pow_primes = (1..).lazy.select{|n| (1..4).all?{|exp| OpenSSL::BN.new(n**exp + n + 1).prime?} }
n = 50
puts "The first #{n} quad-power prime seeds:"
quad_pow_primes.take(n).each_slice(10){|s| puts "%8s"*s.size % s}
- Output:
The 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
Wren
GMP allows us to stretch a little more.
import "./gmp" for Mpz
import "./fmt" for Fmt
var p = Mpz.new()
var isQuadPowerPrimeSeed = Fn.new { |n|
p.setUi(n)
var k = n + 1
return (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 qpps = []
var n = 1
while (qpps.count < 50) {
if (isQuadPowerPrimeSeed.call(n)) qpps.add(n)
n = n + 1
}
System.print("First fifty quad-power prime seeds:")
Fmt.tprint("$,7d", qpps, 10)
System.print("\nFirst quad-power prime seed greater than:")
var m = 1
var c = 50
while (true) {
if (isQuadPowerPrimeSeed.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 + 1
}- Output:
First fifty quad-power prime seeds:
1 2 5 6 69 131 426 1,665 2,129 2,696
5,250 7,929 9,689 13,545 14,154 14,286 16,434 19,760 25,739 27,809
29,631 36,821 41,819 46,619 48,321 59,030 60,500 61,955 62,321 73,610
77,685 79,646 80,535 82,655 85,251 86,996 91,014 96,566 97,739 105,939
108,240 108,681 119,754 122,436 123,164 126,489 140,636 150,480 153,179 163,070
First quad-power prime seed greater than:
1 million is the 141st: 1,009,286
2 million is the 234th: 2,015,496
3 million is the 319th: 3,005,316
4 million is the 383rd: 4,004,726
5 million is the 452nd: 5,023,880
6 million is the 514th: 6,000,554
7 million is the 567th: 7,047,129
8 million is the 601st: 8,005,710
9 million is the 645th: 9,055,151
10 million is the 701st: 10,023,600