Quadrat special primes
- Task
Find the sequence of increasing primes, q, from 2 up to but excluding 16,000,
where the successor of q is the least prime, p, such that p - q is a perfect square.
11l
F is_prime(n)
I n == 2
R 1B
I n < 2 | n % 2 == 0
R 0B
L(i) (3 .. Int(sqrt(n))).step(2)
I n % i == 0
R 0B
R 1B
V Max = 16'000
V quadraPrimes = [2, 3]
V n = 3
L
L(i) (2 .. Int(sqrt(Max))).step(2)
V next = n + i * i
I next >= Max
^L.break
I is_prime(next)
n = next
quadraPrimes [+]= n
L.break
print(‘Quadrat special primes < 16000:’)
L(qp) quadraPrimes
print(‘#5’.format(qp), end' I (L.index + 1) % 7 == 0 {"\n"} E ‘ ’)
- Output:
Quadrat special primes < 16000: 2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Action!
INCLUDE "H6:SIEVE.ACT"
DEFINE MAX="15999"
DEFINE MAXSQUARES="126"
BYTE ARRAY primes(MAX+1)
INT ARRAY squares(MAXSQUARES)
PROC CalcSquares()
INT i
FOR i=1 TO MAXSQUARES
DO
squares(i-1)=i*i
OD
RETURN
INT FUNC FindNextQuadraticPrime(INT x)
INT i,next
FOR i=0 TO MAXSQUARES-1
DO
next=x+squares(i)
IF next>MAX THEN
RETURN (-1)
FI
IF primes(next) THEN
RETURN (next)
FI
OD
RETURN (-1)
PROC Main()
INT p=[2]
Put(125) PutE() ;clear the screen
Sieve(primes,MAX+1)
CalcSquares()
WHILE p>0
DO
PrintI(p) Put(32)
p=FindNextQuadraticPrime(p)
OD
RETURN
- Output:
Screenshot from Atari 8-bit computer
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
ALGOL 68
BEGIN # find some primes where the gap between the current prime and the next is a square #
# an array of squares #
PROC get squares = ( INT n )[]INT:
BEGIN
[ 1 : n ]INT s;
FOR i TO n DO s[ i ] := i * i OD;
s
END # get squares # ;
INT max number = 16000;
PR read "primes.incl.a68" PR
[]BOOL prime = PRIMESIEVE max number;
[]INT square = get squares( max number );
INT p count, this prime, next prime;
# the first square gap is 1 (between 2 and 3) the gap between all other primes is even #
# so we treat 2-3 as a special case #
p count := 1; this prime := 2; next prime := 3;
print( ( " ", whole( this prime, -5 ) ) );
WHILE next prime < max number DO
this prime := next prime;
p count +:= 1;
print( ( " ", whole( this prime, -5 ) ) );
IF p count MOD 12 = 0 THEN print( ( newline ) ) FI;
INT sq pos := 2;
WHILE next prime := this prime + square[ sq pos ];
IF next prime < max number THEN NOT prime[ next prime ] ELSE FALSE FI
DO sq pos +:= 2 OD
OD
END
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
ALGOL W
begin % find some primes where the gap between the current prime and the next is a square %
% sets p( 1 :: n ) to a sieve of primes up to n %
procedure Eratosthenes ( logical array p( * ) ; integer value n ) ;
begin
p( 1 ) := false; p( 2 ) := true;
for i := 3 step 2 until n do p( i ) := true;
for i := 4 step 2 until n do p( i ) := false;
for i := 3 step 2 until truncate( sqrt( n ) ) do begin
integer ii; ii := i + i;
if p( i ) then for pr := i * i step ii until n do p( pr ) := false
end for_i ;
end Eratosthenes ;
% sets s( 1 :: n ) to the squares %
procedure getSquares ( integer array s ( * ) ; integer value n ) ;
for i := 1 until n do s( i ) := i * i;
integer MAX_NUMBER;
MAX_NUMBER := 16000;
begin
logical array prime( 1 :: MAX_NUMBER );
integer array square( 1 :: MAX_NUMBER );
integer pCount, thisPrime, nextPrime;
% sieve the primes to MAX_NUMBER %
Eratosthenes( prime, MAX_NUMBER );
% calculate the squares to MAX_NUMBER %
getSquares( square, MAX_NUMBER );
% the first gap is 1 (between 2 and 3) the gap between all other primes is even %
% so we treat 2-3 as a special case %
pCount := 1; thisPrime := 2; nextPrime := 3;
write( i_w := 6, s_w := 0, " ", thisPrime );
while nextPrime < MAX_NUMBER do begin
integer sqPos;
thisPrime := nextPrime;
pCount := pCount + 1;
writeon( i_w := 6, s_w := 0, " ", thisPrime );
if pCount rem 12 = 0 then write();
sqPos := 2;
while begin
nextPrime := thisPrime + square( sqPos );
nextPrime < MAX_NUMBER and not prime( nextPrime )
end do sqPos := sqPos + 2;
end while_thisPrime_lt_MAX_NUMBER
end
end.
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
AWK
# syntax: GAWK -f QUADRAT_SPECIAL_PRIMES.AWK
# converted from FreeBASIC
BEGIN {
stop = 15999
p = 2
j = 1
printf("%5d ",p)
count++
while (1) {
while (1) {
if (is_prime(p+j*j)) { break }
j++
}
p += j*j
if (p > stop) { break }
printf("%5d%1s",p,++count%10?"":"\n")
j = 1
}
printf("\nQuadrat special primes 1-%d: %d\n",stop,count)
exit(0)
}
function is_prime(x, i) {
if (x <= 1) {
return(0)
}
for (i=2; i<=int(sqrt(x)); i++) {
if (x % i == 0) {
return(0)
}
}
return(1)
}
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671 Quadrat special primes 1-15999: 49
BASIC
Applesoft BASIC
100 FOR P = 2 TO 16000
110 PRINT S$P;
120 LET S$ = " "
130 FOR J = 1 TO 1E9
140 LET V = P + J * J
150 GOSUB 200"IS PRIME"
160 IF NOT ISPRIME THEN NEXT J
170 LET P = V - 1
180 NEXT P
190 END
200 DEF FN MOD(DIVISR) = V - INT (V / DIVISR) * DIVISR
210 ISPRIME = FALSE
220 IF V < 2 THEN RETURN
230 ISPRIME = V = 2
240 IF FN MOD(2) = 0 THEN RETURN
250 ISPRIME = V = 3
260 IF FN MOD(3) = 0 THEN RETURN
270 FOR D = 5 TO SQR (V) STEP 2
280 LET ISPRIME = FN MOD(D)
290 IF NOT ISPRIME THEN RETURN
300 NEXT D
310 RETURN
BASIC256
function isPrime(v)
if v < 2 then return False
if v mod 2 = 0 then return v = 2
if v mod 3 = 0 then return v = 3
d = 5
while d * d <= v
if v mod d = 0 then return False else d += 2
end while
return True
end function
p = 2
j = 1
print 2; " ";
while True
while True
if isPrime(p + j*j) then exit while
j += 1
end while
p += j*j
if p > 16000 then exit while
print p; " ";
j = 1
end while
end
PureBasic
Procedure isPrime(v.i)
If v <= 1 : ProcedureReturn #False
ElseIf v < 4 : ProcedureReturn #True
ElseIf v % 2 = 0 : ProcedureReturn #False
ElseIf v < 9 : ProcedureReturn #True
ElseIf v % 3 = 0 : ProcedureReturn #False
Else
Protected r = Round(Sqr(v), #PB_Round_Down)
Protected f = 5
While f <= r
If v % f = 0 Or v % (f + 2) = 0
ProcedureReturn #False
EndIf
f + 6
Wend
EndIf
ProcedureReturn #True
EndProcedure
OpenConsole()
p.i = 2
j.i = 1
Print("2" + #TAB$)
Repeat
Repeat
If isPrime(p + j*j)
Break
EndIf
j + 1
ForEver
p + j*j
If p > 16000
Break
EndIf
Print(Str(p) + #TAB$)
j = 1
ForEver
Input()
CloseConsole()
Yabasic
sub isPrime(v)
if v < 2 then return False : fi
if mod(v, 2) = 0 then return v = 2 : fi
if mod(v, 3) = 0 then return v = 3 : fi
d = 5
while d * d <= v
if mod(v, d) = 0 then return False else d = d + 2 : fi
wend
return True
end sub
p = 2
j = 1
print 2, " ";
do
do
if isPrime(p + j*j) then break : fi
j = j + 1
loop
p = p + j*j
if p > 16000 then break : fi
print p, " ";
j = 1
loop
end
C++
#include <cstdint>
#include <iostream>
bool is_prime(const uint32_t& number) {
if ( number % 2 == 0 ) {
return number == 2;
}
int k = 3;
while ( k * k <= number ) {
if ( number % k == 0 ) {
return false;
}
k += 2;
}
return true;
}
int main() {
uint32_t i = 1;
uint32_t prime = 2;
std::cout << prime << " ";
while ( true ) {
while ( ! is_prime(prime + i * i) ) {
i += 1;
}
prime += i * i;
if ( prime > 16'000 ) {
break;
}
std::cout << prime << " ";
i = 1;
}
std::cout << std::endl;
}
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Delphi
function IsPrime(N: int64): boolean;
{Fast, optimised prime test}
var I,Stop: int64;
begin
if (N = 2) or (N=3) then Result:=true
else if (n <= 1) or ((n mod 2) = 0) or ((n mod 3) = 0) then Result:= false
else
begin
I:=5;
Stop:=Trunc(sqrt(N+0.0));
Result:=False;
while I<=Stop do
begin
if ((N mod I) = 0) or ((N mod (I + 2)) = 0) then exit;
Inc(I,6);
end;
Result:=True;
end;
end;
procedure QuadratSpecialPrimes(Memo: TMemo);
var Q,P,Cnt: integer;
var IA: TIntegerDynArray;
begin
Memo.Lines.Add('Count Prime1 Prime2 Gap Sqrt');
Memo.Lines.Add('---------------------------------');
Cnt:=0;
Q:=2;
for P:=3 to 16000-1 do
if IsPrime(P) then
begin
if frac(sqrt(P - Q))=0 then
begin
Inc(Cnt);
Memo.Lines.Add(Format('%5D%7D%8D%7D%6.1f',[Cnt,Q,P,P-Q,Sqrt(P-Q)]));
Q:=P;
end;
end;
Memo.Lines.Add('Count = '+IntToStr(Cnt));
end;
- Output:
Count Prime1 Prime2 Gap Sqrt --------------------------------- 1 2 3 1 1.0 2 3 7 4 2.0 3 7 11 4 2.0 4 11 47 36 6.0 5 47 83 36 6.0 6 83 227 144 12.0 7 227 263 36 6.0 8 263 587 324 18.0 9 587 911 324 18.0 10 911 947 36 6.0 11 947 983 36 6.0 12 983 1019 36 6.0 13 1019 1163 144 12.0 14 1163 1307 144 12.0 15 1307 1451 144 12.0 16 1451 1487 36 6.0 17 1487 1523 36 6.0 18 1523 1559 36 6.0 19 1559 2459 900 30.0 20 2459 3359 900 30.0 21 3359 4259 900 30.0 22 4259 4583 324 18.0 23 4583 5483 900 30.0 24 5483 5519 36 6.0 25 5519 5843 324 18.0 26 5843 5879 36 6.0 27 5879 6203 324 18.0 28 6203 6779 576 24.0 29 6779 7103 324 18.0 30 7103 7247 144 12.0 31 7247 7283 36 6.0 32 7283 7607 324 18.0 33 7607 7643 36 6.0 34 7643 8219 576 24.0 35 8219 8363 144 12.0 36 8363 10667 2304 48.0 37 10667 11243 576 24.0 38 11243 11279 36 6.0 39 11279 11423 144 12.0 40 11423 12323 900 30.0 41 12323 12647 324 18.0 42 12647 12791 144 12.0 43 12791 13367 576 24.0 44 13367 13691 324 18.0 45 13691 14591 900 30.0 46 14591 14627 36 6.0 47 14627 14771 144 12.0 48 14771 15671 900 30.0 Count = 48 Elapsed Time: 111.233 ms.
EasyLang
fastfunc isprim num .
i = 2
while i <= sqrt num
if num mod i = 0
return 0
.
i += 1
.
return 1
.
write 2 & " "
p = 2
j = 1
repeat
while isprim (p + j * j) = 0
j += 1
.
p += j * j
until p > 16000
write p & " "
j = 1
.
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
F#
This task uses Extensible Prime Generator (F#)
//Quadrat special primes. Nigel Galloway: January 16th., 2023
let isPs(n:int)=MathNet.Numerics.Euclid.IsPerfectSquare n
let rec fG n g=seq{match Seq.tryHead g with Some h when isPs(h-n)->yield h; yield! fG h g |Some _->yield! fG n g |_->()}
fG 2 (primes32()|>Seq.takeWhile((>)16000))|>Seq.iter(printf "%d "); printfn ""
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Factor
USING: fry io kernel lists lists.lazy math math.primes prettyprint ;
2 [ 1 lfrom swap '[ sq _ + ] lmap-lazy [ prime? ] lfilter car ]
lfrom-by [ 16000 < ] lwhile [ pprint bl ] leach nl
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
FreeBASIC
#include "isprime.bas"
dim as integer p = 2, j = 1
print 2;" ";
do
do
if isprime(p + j*j) then exit do
j += 1
loop
p += j*j
if p > 16000 then exit do
print p;" ";
j = 1
loop
print
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
FutureBasic
local fn IsPrime( n as NSUInteger ) as BOOL
BOOL isPrime = YES
NSUInteger i
if n < 2 then exit fn = NO
if n = 2 then exit fn = YES
if n mod 2 == 0 then exit fn = NO
for i = 3 to int(n^.5) step 2
if n mod i == 0 then exit fn = NO
next
end fn = isPrime
local fn QuadratSpecialPrimes
NSUInteger p = 2, j = 1, count = 1
printf @"%6lu \b", 2
while (1)
count++
while (1)
if fn IsPrime( p + j*j ) then exit while
j += 1
wend
p += j*j
if p > 16000 then exit while
printf @"%6lu \b", p
if count == 7 then count = 0 : print
j = 1
wend
print
end fn
fn QuadratSpecialPrimes
HandleEvents
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Go
package main
import (
"fmt"
"math"
)
func sieve(limit int) []bool {
limit++
// True denotes composite, false denotes prime.
c := make([]bool, limit) // all false by default
c[0] = true
c[1] = true
// no need to bother with even numbers over 2 for this task
p := 3 // Start from 3.
for {
p2 := p * p
if p2 >= limit {
break
}
for i := p2; i < limit; i += 2 * p {
c[i] = true
}
for {
p += 2
if !c[p] {
break
}
}
}
return c
}
func isSquare(n int) bool {
s := int(math.Sqrt(float64(n)))
return s*s == n
}
func commas(n int) string {
s := fmt.Sprintf("%d", n)
if n < 0 {
s = s[1:]
}
le := len(s)
for i := le - 3; i >= 1; i -= 3 {
s = s[0:i] + "," + s[i:]
}
if n >= 0 {
return s
}
return "-" + s
}
func main() {
c := sieve(15999)
fmt.Println("Quadrat special primes under 16,000:")
fmt.Println(" Prime1 Prime2 Gap Sqrt")
lastQuadSpecial := 3
gap := 1
count := 1
fmt.Printf("%7d %7d %6d %4d\n", 2, 3, 1, 1)
for i := 5; i < 16000; i += 2 {
if c[i] {
continue
}
gap = i - lastQuadSpecial
if isSquare(gap) {
sqrt := int(math.Sqrt(float64(gap)))
fmt.Printf("%7s %7s %6s %4d\n", commas(lastQuadSpecial), commas(i), commas(gap), sqrt)
lastQuadSpecial = i
count++
}
}
fmt.Println("\n", count+1, "such primes found.")
}
- Output:
Same as Wren example.
J
{{
j=. 0
r=. 2
while. (j=.j+1)<#y do.
if. (=<.)%:(j{y)-{:r do. r=. r, j{y end.
end.
}} p:i.p:inv 16e3
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Java
public final class QuadratSpecialPrimes {
public static void main(String[] args) {
int i = 1;
int prime = 2;
System.out.print(prime + " ");
while ( true ) {
while ( ! isPrime(prime + i * i) ) {
i += 1;
}
prime += i * i;
if ( prime > 16_000 ) {
break;
}
System.out.print(prime + " ");
i = 1;
}
System.out.println();
}
private static boolean isPrime(int number) {
if ( number % 2 == 0 ) {
return number == 2;
}
int k = 3;
while ( k * k <= number ) {
if ( number % k == 0 ) {
return false;
}
k += 2;
}
return true;
}
}
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
jq
Adaptation of Julia
Works with gojq, the Go implementation of jq
For the definition of `is_prime` used here, see https://rosettacode.org/wiki/Additive_primes
# Input: a number > 2
# Output: an array of the quadrat primes less than `.`
def quadrat:
. as $N
| ($N|sqrt) as $lastn
| { qprimes: [2], q: 2 }
| until ( .qprimes[-1] >= $N or .q >= $N;
label $out
| foreach range(1; $lastn + 1) as $i (.;
.q = .qprimes[-1] + $i * $i
| if .q >= $N then .emit = true
elif .q|is_prime then .qprimes = .qprimes + [.q]
| .emit = true
else .
end;
select(.emit)) | {qprimes, q}, break $out )
| .qprimes ;
"Quadrat special primes < 16000:",
(16000 | quadrat[])
- Output:
Quadrat special primes < 16000: 2 3 7 ... 14627 14771 15671
Julia
using Primes
function quadrat(N = 16000)
pmask = primesmask(1, N)
qprimes, lastn = [2], isqrt(N)
while (n = qprimes[end]) < N
for i in 1:lastn
q = n + i * i
if q > N
return qprimes
elseif pmask[q] # got next qprime
push!(qprimes, q)
break
end
end
end
end
println("Quadrat special primes < 16000:")
foreach(p -> print(rpad(p[2], 6), p[1] % 10 == 0 ? "\n" : ""), enumerate(quadrat()))
- Output:
Quadrat special primes < 16000: 2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Ksh
#!/bin/ksh
# Quadrat Special Primes
# # Variables:
#
alias SHORTINT='typeset -si'
SHORTINT MAXN=16000
# # Functions:
#
# # Function _isquadrat(n, m) return 1 when (m-n) is a perfect square
#
function _isquadrat {
typeset _n ; SHORTINT _n=$1
typeset _m ; SHORTINT _m=$2
[[ $(( sqrt(_m - _n) )) == +(\d).+(\d) ]] && return 0
return 1
}
# # Function _isprime(n) return 1 for prime, 0 for not prime
#
function _isprime {
typeset _n ; integer _n=$1
typeset _i ; integer _i
(( _n < 2 )) && return 0
for (( _i=2 ; _i*_i<=_n ; _i++ )); do
(( ! ( _n % _i ) )) && return 0
done
return 1
}
######
# main #
######
SHORTINT i prev_pr=2
for ((i=2; i<MAXN; i++)); do
_isprime ${i}
if (( $? )); then
_isquadrat ${prev_pr} ${i}
if (( $? )); then
printf "%d " ${i}
prev_pr=${i}
fi
fi
done
printf "\n"
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Mathematica /Wolfram Language
ps = {2};
Do[
Do[
q = Last[ps] + i^2;
If[PrimeQ[q],
AppendTo[ps, q];
Break[]
]
,
{i, 1, \[Infinity]}
];
If[Last[ps] >= 16000,
Break[]]
,
{\[Infinity]}
];
ps //= Most;
Multicolumn[ps, {Automatic, 7}, Appearance -> "Horizontal"]
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Maxima
quadrat(n):=block(aux:next_prime(n),while not integerp(sqrt(aux-n)) do aux:next_prime(aux),aux)$
block(a:2,append([a],makelist(a:quadrat(a),48)));
- Output:
[2,3,7,11,47,83,227,263,587,911,947,983,1019,1163,1307,1451,1487,1523,1559,2459,3359,4259,4583,5483,5519,5843,5879,6203,6779,7103,7247,7283,7607,7643,8219,8363,10667,11243,11279,11423,12323,12647,12791,13367,13691,14591,14627,14771,15671]
Nim
import math, strutils, sugar
func isPrime(n: Natural): bool =
if n < 2: return false
if n mod 2 == 0: return n == 2
if n mod 3 == 0: return n == 3
var d = 5
while d * d <= n:
if n mod d == 0: return false
inc d, 2
if n mod d == 0: return false
inc d, 4
result = true
const
Max = 16_000
Squares = collect(newSeq):
for n in countup(2, sqrt(Max.float).int, 2): n * n
iterator quadraPrimes(lim: Positive): int =
assert lim >= 3
yield 2
yield 3
var n = 3
block mainloop:
while true:
for square in Squares:
let next = n + square
if next > lim: break mainloop
if next.isPrime:
n = next
yield n
break
echo "Quadrat special primes < 16000:"
var count = 0
for qp in quadraPrimes(Max):
inc count
stdout.write ($qp).align(5), if count mod 7 == 0: '\n' else: ' '
- Output:
Quadrat special primes < 16000: 2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Perl
use strict;
use warnings;
use feature 'say';
use ntheory 'is_prime';
my @sp = my $previous = 2;
do {
my($next,$n);
while () { last if is_prime( $next = $previous + ++$n**2 ) }
push @sp, $next;
$previous = $next;
} until $sp[-1] >= 16000;
pop @sp and say ((sprintf '%-7d'x@sp, @sp) =~ s/.{1,$#sp}\K\s/\n/gr);
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Phix
constant desc = split("linear quadratic cubic quartic quintic sextic septic octic nonic decic"), limits = {1, 16000, 15000, 14e9, 8025e5, 25e12, 195e12,75e11, 3e9, 11e8} for p=2 to length(desc) do atom N = limits[p], lastn = ceil(power(N,1/p)) sequence res = {2} bool done = false while not done do for n=1 to lastn do atom m = res[$] + power(n,p) if m>N then done = true exit elsif is_prime(m) then res &= m exit end if end for end while string r = join_by(apply(true,sprintf,{{"%,6d"},res}),1,p+5) printf(1,"Found %d %s special primes < %g:\n%s\n",{length(res),desc[p],N,r}) end for
- Output:
Found 49 quadratic special primes < 16000: 2 3 7 11 47 83 227 263 587 911 947 983 1,019 1,163 1,307 1,451 1,487 1,523 1,559 2,459 3,359 4,259 4,583 5,483 5,519 5,843 5,879 6,203 6,779 7,103 7,247 7,283 7,607 7,643 8,219 8,363 10,667 11,243 11,279 11,423 12,323 12,647 12,791 13,367 13,691 14,591 14,627 14,771 15,671 Found 23 cubic special primes < 15000: 2 3 11 19 83 1,811 2,027 2,243 2,251 2,467 2,531 2,539 3,539 3,547 4,547 5,059 10,891 12,619 13,619 13,627 13,691 13,907 14,419 Found 9 quartic special primes < 1.4e+10: 2 3 19 160,019 1,049,920,019 1,050,730,019 1,051,540,019 12,910,750,019 13,960,510,019 Found 9 quintic special primes < 8.025e+8: 2 3 32,771 32,803 33,827 41,603 579,427 778,179,427 802,479,427 Found 6 sextic special primes < 2.5e+13: 2 3 67 131 2,176,782,467 22,485,250,805,891 Found 7 septic special primes < 1.95e+14: 2 3 131 194,871,710,000,131 194,893,580,000,131 194,893,580,280,067 194,971,944,444,163 Found 4 octic special primes < 7.5e+12: 2 3 65,539 6,553,600,065,539 Found 6 nonic special primes < 3e+9: 2 3 262,147 10,339,843 20,417,539 1,020,417,539 Found 4 decic special primes < 1.1e+9: 2 3 1,073,741,827 1,073,742,851
Python
#!/usr/bin/python
def isPrime(n):
for i in range(2, int(n**0.5) + 1):
if n % i == 0:
return False
return True
if __name__ == '__main__':
p = 2
j = 1
print(2, end = " ");
while True:
while True:
if isPrime(p + j*j):
break
j += 1
p += j*j
if p > 16000:
break
print(p, end = " ");
j = 1
Raku
my @sqp = 2, -> $previous {
my $next;
for (1..∞).map: *² {
$next = $previous + $_;
last if $next.is-prime;
}
$next
} … *;
say "{+$_} matching numbers:\n", $_».fmt('%5d').batch(7).join: "\n" given
@sqp[^(@sqp.first: * > 16000, :k)];
- Output:
49 matching numbers: 2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
REXX
Version 1 Inline code
/*REXX program finds the smallest primes such that the difference of successive terms */
/*─────────────────────────────────────────────────── are the smallest quadrat numbers. */
parse arg hi cols . /*obtain optional argument from the CL.*/
if hi=='' | hi=="," then hi= 16000 /* " " " " " " */
if cols=='' | cols=="," then cols= 10 /* " " " " " " */
call genP /*build array of semaphores for primes.*/
w= 10 /*width of a number in any column. */
title= 'the smallest primes < ' commas(hi) " such that the" ,
'difference of successive terms are the smallest quadrat numbers'
if cols>0 then say ' index │'center(title, 1 + cols*(w+1) )
if cols>0 then say '───────┼'center("" , 1 + cols*(w+1), '─')
sqp= 0; idx= 1 /*initialize number of sqp and index.*/
op= 1
$= /*list of quad─special primes (so far).*/
do j=0 by 0
do k=1 until !.np; np= op + k*k /*find the next square + oldPrime.*/
if np>= hi then leave j /*Is newPrime too big? Then quit.*/
end /*k*/
sqp= sqp + 1 /*bump the number of sqp's. */
op= np /*assign the newPrime to the oldPrime*/
if cols<0 then iterate /*Build the list (to be shown later)? */
c= commas(np) /*maybe add commas to the number. */
$= $ right(c, max(w, length(c) ) ) /*add quadratic─special prime ──► list.*/
if sqp//cols\==0 then iterate /*have we populated a line of output? */
say center(idx, 7)'│' substr($, 2); $= /*display what we have so far (cols). */
idx= idx + cols /*bump the index count for the output*/
end /*j*/
if $\=='' then say center(idx, 7)"│" substr($, 2) /*possible display residual output.*/
if cols>0 then say '───────┴'center("" , 1 + cols*(w+1), '─')
say
say 'Found ' commas(sqp) " of " title
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ?
/*──────────────────────────────────────────────────────────────────────────────────────*/
genP: !.= 0 /*placeholders for primes (semaphores).*/
@.1=2; @.2=3; @.3=5; @.4=7; @.5=11 /*define some low primes. */
!.2=1; !.3=1; !.5=1; !.7=1; !.11=1 /* " " " " flags. */
#=5; s.#= @.# **2 /*number of primes so far; prime². */
/* [↓] generate more primes ≤ high.*/
do j=@.#+2 by 2 to hi /*find odd primes from here on. */
parse var j '' -1 _; if _==5 then iterate /*J divisible by 5? (right dig)*/
if j// 3==0 then iterate /*" " " 3? */
if j// 7==0 then iterate /*" " " 7? */
do k=5 while s.k<=j /* [↓] divide by the known odd primes.*/
if j // @.k == 0 then iterate j /*Is J ÷ X? Then not prime. ___ */
end /*k*/ /* [↑] only process numbers ≤ √ J */
#= #+1; @.#= j; s.#= j*j; !.j= 1 /*bump # of Ps; assign next P; P²; P# */
end /*j*/; return
- output when using the default inputs:
index │ the smallest primes < 16,000 such that the difference of successive terma are the smallest quadrat numbers ───────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ 2 3 7 11 47 83 227 263 587 911 11 │ 947 983 1,019 1,163 1,307 1,451 1,487 1,523 1,559 2,459 21 │ 3,359 4,259 4,583 5,483 5,519 5,843 5,879 6,203 6,779 7,103 31 │ 7,247 7,283 7,607 7,643 8,219 8,363 10,667 11,243 11,279 11,423 41 │ 12,323 12,647 12,791 13,367 13,691 14,591 14,627 14,771 15,671 ───────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────── Found 49 of the smallest primes < 16,000 such that the difference of successive terma are the smallest quadrat numbers
Version 2 Using REXX libraries and simplified
Libraries: How to use
Library: Sequences
Library: Functions
See Cubic special primes. The Versions 2 are almost the same, only line 'b = a + k*k*k' is changed in 'b = a + k*k'.
call Time('r')
say 'Quadrat special primes - Using REXX libraries'
parse version version; say version; say
glob. = ''
t = 16000; d = Xpon(t)+2; numeric digits d
call Primes t
call Quadrats t,d
say
say Format(Time('e'),,3) 'seconds'
exit
Quadrats:
procedure expose flag.
arg x,y
call Charout ,Right(2,y)
a = 2; b = 3; k = 1; n = 1
do while b <= x
b = a + k*k
if flag.prime.b = 1 then do
n = n+1
call Charout ,Right(b,y)
if n//10 = 0 then
say
a = b; k = 1
end
k = k+1
end
say; say
say n 'quadratic special primes found <' x
return
include Sequences
include Functions
- Output:
Quadrat special primes - Using REXX libraries REXX-Regina_3.9.6(MT) 5.00 29 Apr 2024 2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671 49 quadratic special primes found below 16000 0.005 seconds
Ring
load "stdlib.ring"
/* Searching for the smallest prime gaps under a limit,
such that the difference of successive terms (gaps)
is of the smallest degree. */
? "working..."
desc = split("na quadratic cubic quartic quintic sextic septic octic nonic decic"," ")
limits = [1, 16000, 15000, 30000, 50000, 50000, 50000, 75000, 300000, 500000]
for deg = 2 to len(desc)
Primes = []
limit = limits[deg]
oldPrime = 2
add(Primes, 2)
for n = 1 to sqrt(limit)
nextPrime = oldPrime + pow(n, deg)
if isprime(nextPrime)
n = 1
if nextPrime < limit add(Primes, nextPrime) ok
oldPrime = nextPrime
else
nextPrime = nextPrime - oldPrime
ok
if nextPrime > limit exit ok
next
? nl + desc[deg] + ":" + nl + " prime1 prime2 Gap Rt"
for n = 1 to Len(Primes) - 1
diff = Primes[n + 1] - Primes[n]
? sf(Primes[n], 7) + " " + sf(Primes[n+1], 7) + " " + sf(diff, 6) + " " + sf(floor(0.49 + pow(diff, 1 / deg)), 4)
next
? "Found " + Len(Primes) + " primes under " + limit + " for " + desc[deg] + " gaps."
next
? nl + "done..."
# a very plain string formatter, intended to even up columnar outputs
def sf x, y
s = string(x) l = len(s)
if l > y y = l ok
return substr(" ", 11 - y + l) + s
- Output:
working... quadratic: prime1 prime2 Gap Rt 2 3 1 1 3 7 4 2 7 11 4 2 11 47 36 6 47 83 36 6 83 227 144 12 227 263 36 6 263 587 324 18 587 911 324 18 911 947 36 6 947 983 36 6 983 1019 36 6 1019 1163 144 12 1163 1307 144 12 1307 1451 144 12 1451 1487 36 6 1487 1523 36 6 1523 1559 36 6 1559 2459 900 30 2459 3359 900 30 3359 4259 900 30 4259 4583 324 18 4583 5483 900 30 5483 5519 36 6 5519 5843 324 18 5843 5879 36 6 5879 6203 324 18 6203 6779 576 24 6779 7103 324 18 7103 7247 144 12 7247 7283 36 6 7283 7607 324 18 7607 7643 36 6 7643 8219 576 24 8219 8363 144 12 8363 10667 2304 48 10667 11243 576 24 11243 11279 36 6 11279 11423 144 12 11423 12323 900 30 12323 12647 324 18 12647 12791 144 12 12791 13367 576 24 13367 13691 324 18 13691 14591 900 30 14591 14627 36 6 14627 14771 144 12 14771 15671 900 30 Found 49 primes under 16000 for quadratic gaps. cubic: prime1 prime2 Gap Rt 2 3 1 1 3 11 8 2 11 19 8 2 19 83 64 4 83 1811 1728 12 1811 2027 216 6 2027 2243 216 6 2243 2251 8 2 2251 2467 216 6 2467 2531 64 4 2531 2539 8 2 2539 3539 1000 10 3539 3547 8 2 3547 4547 1000 10 4547 5059 512 8 5059 10891 5832 18 10891 12619 1728 12 12619 13619 1000 10 13619 13627 8 2 13627 13691 64 4 13691 13907 216 6 13907 14419 512 8 Found 23 primes under 15000 for cubic gaps. quartic: prime1 prime2 Gap Rt 2 3 1 1 3 19 16 2 Found 3 primes under 30000 for quartic gaps. quintic: prime1 prime2 Gap Rt 2 3 1 1 3 32771 32768 8 32771 32803 32 2 32803 33827 1024 4 33827 41603 7776 6 Found 6 primes under 50000 for quintic gaps. sextic: prime1 prime2 Gap Rt 2 3 1 1 3 67 64 2 67 131 64 2 Found 4 primes under 50000 for sextic gaps. septic: prime1 prime2 Gap Rt 2 3 1 1 3 131 128 2 Found 3 primes under 50000 for septic gaps. octic: prime1 prime2 Gap Rt 2 3 1 1 3 65539 65536 4 Found 3 primes under 75000 for octic gaps. nonic: prime1 prime2 Gap Rt 2 3 1 1 3 262147 262144 4 Found 3 primes under 300000 for nonic gaps. decic: prime1 prime2 Gap Rt 2 3 1 1 Found 2 primes under 500000 for decic gaps. done...
RPL
≪ { 2 } 2 DUP DO DUP NEXTPRIME IF DUP2 SWAP - √ FP NOT THEN NIP SWAP OVER + SWAP DUP END UNTIL DUP 16000 ≥ END DROP2 ≫ 'TASK' STO
- Output:
1: {2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671}
Runs in 6 minutes 25 seconds on a HP-50g.
Ruby
require 'prime'
res = [2]
until res.last > 16000 do
res << (1..).detect{|n| (res.last + n**2).prime? } ** 2 + res.last
end
puts res[..-2].join(" ")
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671
Sidef
func quadrat_primes(callback) {
var prev = 2
callback(prev)
loop {
var curr = (1..Inf -> lazy.map { prev + _**2 }.first { .is_prime })
callback(curr)
prev = curr
}
}
say gather {
quadrat_primes({|k|
break if (k >= 16000)
take(k)
})
}
- Output:
[2, 3, 7, 11, 47, 83, 227, 263, 587, 911, 947, 983, 1019, 1163, 1307, 1451, 1487, 1523, 1559, 2459, 3359, 4259, 4583, 5483, 5519, 5843, 5879, 6203, 6779, 7103, 7247, 7283, 7607, 7643, 8219, 8363, 10667, 11243, 11279, 11423, 12323, 12647, 12791, 13367, 13691, 14591, 14627, 14771, 15671]
Wren
import "./math" for Int
import "./fmt" for Fmt
var isSquare = Fn.new { |n|
var s = n.sqrt.floor
return s*s == n
}
var primes = Int.primeSieve(15999)
System.print("Quadrat special primes under 16,000:")
System.print(" Prime1 Prime2 Gap Sqrt")
var lastQuadSpecial = 3
var gap = 1
var count = 1
Fmt.print("$,7d $,7d $,6d $4d", 2, 3, 1, 1)
for (p in primes.skip(2)) {
gap = p - lastQuadSpecial
if (isSquare.call(gap)) {
Fmt.print("$,7d $,7d $,6d $4d", lastQuadSpecial, p, gap, gap.sqrt)
lastQuadSpecial = p
count = count + 1
}
}
System.print("\n%(count+1) such primes found.")
- Output:
Quadrat special primes under 16,000: Prime1 Prime2 Gap Sqrt 2 3 1 1 3 7 4 2 7 11 4 2 11 47 36 6 47 83 36 6 83 227 144 12 227 263 36 6 263 587 324 18 587 911 324 18 911 947 36 6 947 983 36 6 983 1,019 36 6 1,019 1,163 144 12 1,163 1,307 144 12 1,307 1,451 144 12 1,451 1,487 36 6 1,487 1,523 36 6 1,523 1,559 36 6 1,559 2,459 900 30 2,459 3,359 900 30 3,359 4,259 900 30 4,259 4,583 324 18 4,583 5,483 900 30 5,483 5,519 36 6 5,519 5,843 324 18 5,843 5,879 36 6 5,879 6,203 324 18 6,203 6,779 576 24 6,779 7,103 324 18 7,103 7,247 144 12 7,247 7,283 36 6 7,283 7,607 324 18 7,607 7,643 36 6 7,643 8,219 576 24 8,219 8,363 144 12 8,363 10,667 2,304 48 10,667 11,243 576 24 11,243 11,279 36 6 11,279 11,423 144 12 11,423 12,323 900 30 12,323 12,647 324 18 12,647 12,791 144 12 12,791 13,367 576 24 13,367 13,691 324 18 13,691 14,591 900 30 14,591 14,627 36 6 14,627 14,771 144 12 14,771 15,671 900 30 49 such primes found.
XPL0
Find primes where the difference between the current one and a following one is a perfect square.
func IsPrime(N); \Return 'true' if N is a prime number
int N, I;
[if N <= 1 then return false;
for I:= 2 to sqrt(N) do
if rem(N/I) = 0 then return false;
return true;
];
int Count, P, Q;
[Count:= 0;
P:= 2; Q:= 3;
repeat if IsPrime(Q) then
[if sq(sqrt(Q-P)) = Q-P then
[IntOut(0, P);
P:= Q;
Count:= Count+1;
if rem(Count/10) then ChOut(0, 9\tab\) else CrLf(0);
];
];
Q:= Q+2;
until P >= 16000;
CrLf(0);
IntOut(0, Count);
Text(0, " such primes found below 16000.
");
]
- Output:
2 3 7 11 47 83 227 263 587 911 947 983 1019 1163 1307 1451 1487 1523 1559 2459 3359 4259 4583 5483 5519 5843 5879 6203 6779 7103 7247 7283 7607 7643 8219 8363 10667 11243 11279 11423 12323 12647 12791 13367 13691 14591 14627 14771 15671 49 such primes found below 16000.
- Draft Programming Tasks
- Prime Numbers
- 11l
- Action!
- Action! Sieve of Eratosthenes
- ALGOL 68
- ALGOL 68-primes
- ALGOL W
- AWK
- BASIC
- Applesoft BASIC
- BASIC256
- PureBasic
- Yabasic
- C++
- Delphi
- SysUtils,StdCtrls
- EasyLang
- F Sharp
- Factor
- FreeBASIC
- FutureBasic
- Go
- J
- Java
- Jq
- Julia
- Ksh
- Mathematica
- Wolfram Language
- Maxima
- Nim
- Perl
- Ntheory
- Phix
- Python
- Raku
- REXX
- Ring
- RPL
- Ruby
- Sidef
- Wren
- Wren-math
- Wren-fmt
- XPL0