Quadrat special primes

Quadrat special primes 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.

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

Translation of: Nim

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!

Library: Action! Sieve of Eratosthenes

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

Translation of: ALGOL W

Library: ALGOL 68-primes

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

Translation of: BASIC256

 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

Delphi

Works with: Delphi version 6.0

Library: SysUtils,StdCtrls

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.

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

Works with: Factor version 0.98

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

Translation of: Wren

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

jq

Adaptation of Julia

Works with: jq

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

Library: ntheory

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

/*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

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

Works with: HP version 49

≪ { 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

Library: Wren-math

Library: Wren-fmt

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.