Find adjacent primes which differ by a square integer


Find adjacent primes under 1,000,000 whose difference (> 36) is a square integer.

Find adjacent primes which differ by a square integer 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

11l

F primes_upto(limit)
   V is_prime = [0B] * 2 [+] [1B] * (limit - 1)
   L(n) 0 .< Int(limit ^ 0.5 + 1.5)
      I is_prime[n]
         L(i) (n * n .. limit).step(n)
            is_prime[i] = 0B
   R enumerate(is_prime).filter((i, prime) -> prime).map((i, prime) -> i)

V primes = primes_upto(1'000'000)

F is_square(x)
   R Int(sqrt(x)) ^ 2 == x

L(n) 2 .< primes.len
   V pr1 = primes[n]
   V pr2 = primes[n - 1]
   V diff = pr1 - pr2
   I (is_square(diff) & diff > 36)
      print(pr1‘ ’pr2‘ diff = ’diff)
Output:
89753 89689 diff = 64
107441 107377 diff = 64
288647 288583 diff = 64
368021 367957 diff = 64
381167 381103 diff = 64
396833 396733 diff = 100
400823 400759 diff = 64
445427 445363 diff = 64
623171 623107 diff = 64
625763 625699 diff = 64
637067 637003 diff = 64
710777 710713 diff = 64
725273 725209 diff = 64
779477 779413 diff = 64
801947 801883 diff = 64
803813 803749 diff = 64
821741 821677 diff = 64
832583 832519 diff = 64
838349 838249 diff = 100
844841 844777 diff = 64
883871 883807 diff = 64
912167 912103 diff = 64
919511 919447 diff = 64
954827 954763 diff = 64
981887 981823 diff = 64
997877 997813 diff = 64

ALGOL 68

BEGIN # find a adjacent primes where the primes differ by a square > 36 #
    INT min diff  = 37;
    INT max prime = 1 000 000;
    PR read "primes.incl.a68" PR
    # form a list of primes to max prime #
    []INT prime = EXTRACTPRIMESUPTO max prime FROMPRIMESIEVE PRIMESIEVE max prime;
    # construct a table of squares, we will need at most the square root of max prime #
    # but in reality much less than that - assume 1000 will be enough #
    [ 1 : 1000 ]BOOL is square;
    FOR i TO UPB is square DO is square[ i ] := FALSE OD;
    FOR i WHILE INT i2 = i * i;
                i2 <= UPB is square
    DO
        is square[ i2 ] := TRUE
    OD;
    # find the primes #
    FOR p TO UPB prime - 1 DO
        INT q    = p + 1;
        INT diff = prime[ q ] - prime[ p ];
        IF diff > min diff AND is square[ diff ] THEN
            print( ( whole( prime[ q ], -6 ), " - ", whole( prime[ p ], -6 ), " = ", whole( diff, 0 ), newline ) )
        FI
    OD
END
Output:
 89753 -  89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64

Arturo

squares: map 7..15 'x -> x*x
primes: select 1..1000000 => prime?

loop.with:'i primes\[0..(size primes)-2] 'p [
    next: primes\[i+1]
    if contains? squares next-p ->
        print [pad to :string next 6 "-" pad to :string p 6 "=" next-p]
]
Output:
 89753 -  89689 = 64 
107441 - 107377 = 64 
288647 - 288583 = 64 
368021 - 367957 = 64 
381167 - 381103 = 64 
396833 - 396733 = 100 
400823 - 400759 = 64 
445427 - 445363 = 64 
623171 - 623107 = 64 
625763 - 625699 = 64 
637067 - 637003 = 64 
710777 - 710713 = 64 
725273 - 725209 = 64 
779477 - 779413 = 64 
801947 - 801883 = 64 
803813 - 803749 = 64 
821741 - 821677 = 64 
832583 - 832519 = 64 
838349 - 838249 = 100 
844841 - 844777 = 64 
883871 - 883807 = 64 
912167 - 912103 = 64 
919511 - 919447 = 64 
954827 - 954763 = 64 
981887 - 981823 = 64 
997877 - 997813 = 64

AWK

# syntax: GAWK -f FIND_ADJACENTS_PRIMES_WHICH_DIFFERENCE_IS_SQUARE_INTEGER.AWK
# converted from FreeBASIC
BEGIN {
    start = i = 3
    stop =  999999
    while (j <= stop) {
      j = next_prime(i)
      if (j-i > 36 && is_square(j-i)) {
        printf("%9d %9d %9d\n",i,j,j-i)
        count++
      }
      i = j
    }
    printf("Adjacent primes which difference is square integer (>36) %d-%d: %d\n",start,stop,count)
    exit(0)
}
function is_prime(n,  d) {
    d = 5
    if (n < 2) { return(0) }
    if (n % 2 == 0) { return(n == 2) }
    if (n % 3 == 0) { return(n == 3) }
    while (d*d <= n) {
      if (n % d == 0) { return(0) }
      d += 2
      if (n % d == 0) { return(0) }
      d += 4
    }
    return(1)
}
function is_square(n) {
    return (int(sqrt(n))^2 == n)
}
function next_prime(n,  q) { # finds next prime after n
    if (n == 0) { return(2) }
    if (n < 3) { return(++n) }
    q = n + 2
    while (!is_prime(q)) {
      q += 2
    }
    return(q)
}
Output:
    89689     89753        64
   107377    107441        64
   288583    288647        64
   367957    368021        64
   381103    381167        64
   396733    396833       100
   400759    400823        64
   445363    445427        64
   623107    623171        64
   625699    625763        64
   637003    637067        64
   710713    710777        64
   725209    725273        64
   779413    779477        64
   801883    801947        64
   803749    803813        64
   821677    821741        64
   832519    832583        64
   838249    838349       100
   844777    844841        64
   883807    883871        64
   912103    912167        64
   919447    919511        64
   954763    954827        64
   981823    981887        64
   997813    997877        64
Adjacent primes which difference is square integer (>36) 3-999999: 26

C

#include<stdio.h>
#include<stdlib.h>

int isprime( int p ) {
    int i;
    if(p==2) return 1;
    if(!(p%2)) return 0;
    for(i=3; i*i<=p; i+=2) {
       if(!(p%i)) return 0;
    }
    return 1;
}

int nextprime( int p ) {
    int i=0;
    if(p==0) return 2;
    if(p<3) return p+1;
    while(!isprime(++i + p));
    return i+p;
}

int issquare( int p ) {
    int i;
    for(i=0;i*i<p;i++);
    return i*i==p;
}

int main(void) {
    int i=3, j=2;
    for(i=3;j<=1000000;i=j) {
        j=nextprime(i);
        if(j-i>36&&issquare(j-i)) printf( "%d %d %d\n", i, j, j-i );
    }
    return 0;
}

C++

#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>

std::vector<uint32_t> list_prime_numbers(const uint32_t& limit) {
	std::vector<uint32_t> primes{};
	primes.emplace_back(2);
	const uint32_t half_limit = ( limit + 1 ) / 2;
	std::vector<bool> composite(half_limit, false);
	for ( uint32_t i = 1, p = 3; i < half_limit; p += 2, ++i ) {
		if ( ! composite[i] ) {
			primes.emplace_back(p);
			for ( uint32_t a = i + p; a < half_limit; a += p ) {
				composite[a] = true;
			}
		}
	}
	return primes;
}

bool is_square(const uint32_t& number) {
	const uint32_t root = std::floor(std::sqrt(number));
	return root * root == number;
}

int main() {
	const std::vector<uint32_t> primes = list_prime_numbers(1'000'000);

	for ( uint32_t i = 2; i < primes.size(); ++i ) {
		const uint32_t prime2 = primes[i - 1];
		const uint32_t prime1 = primes[i];
		const uint32_t difference = prime1 - prime2;
		if ( difference > 36 && is_square(difference) ) {
			std::cout << std::setw(7) << prime2 << " and " << std::setw(6) << prime1
					  << " : difference = " << difference << std::endl;
		}
	}
}
Output:
  89689 and  89753 : difference = 64
 107377 and 107441 : difference = 64
 288583 and 288647 : difference = 64
 367957 and 368021 : difference = 64
 381103 and 381167 : difference = 64
 396733 and 396833 : difference = 100
 400759 and 400823 : difference = 64
 445363 and 445427 : difference = 64
 623107 and 623171 : difference = 64
 625699 and 625763 : difference = 64
 637003 and 637067 : difference = 64
 710713 and 710777 : difference = 64
 725209 and 725273 : difference = 64
 779413 and 779477 : difference = 64
 801883 and 801947 : difference = 64
 803749 and 803813 : difference = 64
 821677 and 821741 : difference = 64
 832519 and 832583 : difference = 64
 838249 and 838349 : difference = 100
 844777 and 844841 : difference = 64
 883807 and 883871 : difference = 64
 912103 and 912167 : difference = 64
 919447 and 919511 : difference = 64
 954763 and 954827 : difference = 64
 981823 and 981887 : difference = 64
 997813 and 997877 : difference = 64

CLU

% Integer square root
isqrt = proc (s: int) returns (int)
    x0: int := s/2
    if x0=0 then return(s) end
    x1: int := (x0 + s/x0)/2
    while x1 < x0 do
        x0 := x1
        x1 := (x0 + s/x0)/2
    end
    return(x0)
end isqrt
    
% See if a number is square
% Note that all squares are 0, 1, 4, or 9 mod 16.
is_square = proc (n: int) returns (bool)
    d: int := n//16
    if d=0 cor d=1 cor d=4 cor d=9 then
        return(n = isqrt(n)**2)
    else    
        return(false)
    end
end is_square

% Find all primes up to a given number
sieve = proc (top: int) returns (array[int])
    prime: array[bool] := array[bool]$fill(2,top-1,true)
    for p: int in int$from_to(2,isqrt(top)) do
        if prime[p] then
            for c: int in int$from_to_by(p*p,top,p) do
                prime[c] := false
            end
        end
    end
    list: array[int] := array[int]$predict(1,isqrt(top))
    for p: int in int$from_to(2,top) do
        if prime[p] then array[int]$addh(list,p) end
    end
    return(list)
end sieve

start_up = proc ()
    MAX = 1000000
    DIFF = 36
    
    po: stream := stream$primary_output()
    primes: array[int] := sieve(MAX)
    for i: int in int$from_to(array[int]$low(primes)+1, 
                              array[int]$high(primes)) do
        d: int := primes[i] - primes[i-1]
        if d>DIFF cand is_square(d) then
            stream$putright(po, int$unparse(primes[i]), 6)
            stream$puts(po, " - ")
            stream$putright(po, int$unparse(primes[i-1]), 6)
            stream$puts(po, " = ")
            stream$putright(po, int$unparse(d), 4)
            stream$puts(po, " = ")
            stream$putright(po, int$unparse(isqrt(d)), 4)
            stream$putl(po, "^2")
        end
    end
end start_up
Output:
 89753 -  89689 =   64 =    8^2
107441 - 107377 =   64 =    8^2
288647 - 288583 =   64 =    8^2
368021 - 367957 =   64 =    8^2
381167 - 381103 =   64 =    8^2
396833 - 396733 =  100 =   10^2
400823 - 400759 =   64 =    8^2
445427 - 445363 =   64 =    8^2
623171 - 623107 =   64 =    8^2
625763 - 625699 =   64 =    8^2
637067 - 637003 =   64 =    8^2
710777 - 710713 =   64 =    8^2
725273 - 725209 =   64 =    8^2
779477 - 779413 =   64 =    8^2
801947 - 801883 =   64 =    8^2
803813 - 803749 =   64 =    8^2
821741 - 821677 =   64 =    8^2
832583 - 832519 =   64 =    8^2
838349 - 838249 =  100 =   10^2
844841 - 844777 =   64 =    8^2
883871 - 883807 =   64 =    8^2
912167 - 912103 =   64 =    8^2
919511 - 919447 =   64 =    8^2
954827 - 954763 =   64 =    8^2
981887 - 981823 =   64 =    8^2
997877 - 997813 =   64 =    8^2

Delphi

Works with: Delphi version 6.0


function IsPrime(N: integer): boolean;
{Optimised prime test - about 40% faster than the naive approach}
var I,Stop: integer;
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));
	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;

function GetNextPrime(Start: integer): integer;
{Get the next prime number after Start}
begin
repeat Inc(Start)
until IsPrime(Start);
Result:=Start;
end;



procedure ShowPrimeDiffs(Memo: TMemo);
var P1,P2,D: integer;
begin
P1:=GetNextPrime(2);
repeat
	begin
	P2:=GetNextPrime(P1);
	D:=P2 - P1;
	if (D>36) and (Frac(sqrt(D))=0) then
		begin
		Memo.Lines.Add(IntToStr(P2)+' - '+IntToStr(P1)+' = '+IntToStr(D));
		end;
	P1:=P2;
	end
until P2>=1000000;
end;
Output:
89753 - 89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64


DuckDB

Works with: DuckDB version V1.0
## Preliminaries

create or replace function primep(nnumber) as (
  select
    case
    when nnumber < 2 then false
    when nnumber = 2 then true
    else NOT exists
         ( select * from 
           ( select (nnumber % anumber) as modNumber
             from (select unnest(range(2, 1 + sqrt(nnumber)::BIGINT)) as anumber)
           ) 
           where modNumber = 0 
         )
    end
);

# primes up to and possibly including mx
create or replace function primes(mx) as table (
  select if(mx>0,2,null) as p
  union all
  ( select n
    from range(3,mx+1,2) _(n)
    where primep(n))
);

create or replace function issquare(n) as trunc(sqrt(n)) ^ 2 = n;

### The task:

with primes_t as (
  select prime, lead(prime) over () as nextprime
  from primes(1000000) _(prime)
)
select prime, nextprime
from primes_t
where (nextprime - prime) > 36 and issquare(nextprime - prime)
order by prime ;
Output:
┌────────┬───────────┐
│ prime  │ nextprime │
│ int64  │   int64   │
├────────┼───────────┤
│  89689 │     89753 │
│ 107377 │    107441 │
│ 288583 │    288647 │
│ 367957 │    368021 │
│ 381103 │    381167 │
│ 396733 │    396833 │
│ 400759 │    400823 │
│ 445363 │    445427 │
│ 623107 │    623171 │
│ 625699 │    625763 │
│ 637003 │    637067 │
│ 710713 │    710777 │
│ 725209 │    725273 │
│ 779413 │    779477 │
│ 801883 │    801947 │
│ 803749 │    803813 │
│ 821677 │    821741 │
│ 832519 │    832583 │
│ 838249 │    838349 │
│ 844777 │    844841 │
│ 883807 │    883871 │
│ 912103 │    912167 │
│ 919447 │    919511 │
│ 954763 │    954827 │
│ 981823 │    981887 │
│ 997813 │    997877 │
├────────┴───────────┤
│ 26 rows  2 columns │
└────────────────────┘

EasyLang

Translation of: AWK
fastfunc isprim num .
   i = 2
   while i <= sqrt num
      if num mod i = 0
         return 0
      .
      i += 1
   .
   return 1
.
prim = 2
proc nextprim . .
   repeat
      prim += 1
      until isprim prim = 1
   .
.
func is_square n .
   h = floor sqrt n
   return if h * h = n
.
while prim < 1000000
   prev = prim
   nextprim
   if prim - prev > 36 and is_square (prim - prev) = 1
      print prim & " - " & prev & " = " & prim - prev
   .
.
Output:
89753 - 89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64

F#

This task uses Extensible Prime Generator (F#)

// Find adjacents primes which difference is square integer . Nigel Galloway: November 23rd., 2021
primes32()|>Seq.takeWhile((>)1000000)|>Seq.pairwise|>Seq.filter(fun(n,g)->let n=g-n in let g=(float>>sqrt>>int)n in g>6 && n=g*g)|>Seq.iter(printfn "%A")
Output:
(89689, 89753)
(107377, 107441)
(288583, 288647)
(367957, 368021)
(381103, 381167)
(396733, 396833)
(400759, 400823)
(445363, 445427)
(623107, 623171)
(625699, 625763)
(637003, 637067)
(710713, 710777)
(725209, 725273)
(779413, 779477)
(801883, 801947)
(803749, 803813)
(821677, 821741)
(832519, 832583)
(838249, 838349)
(844777, 844841)
(883807, 883871)
(912103, 912167)
(919447, 919511)
(954763, 954827)
(981823, 981887)
(997813, 997877)

Factor

Works with: Factor version 0.99 2021-06-02
USING: formatting io kernel lists lists.lazy math math.functions
math.primes.lists sequences ;

: adj-primes ( -- list ) lprimes dup cdr lzip ;

: diff ( pair -- n ) first2 swap - ;

: adj-primes-diff ( -- list )
    adj-primes [ dup diff suffix ] lmap-lazy ;

: big-adj-primes-diff ( -- list )
    adj-primes-diff [ last 36 > ] lfilter ;

: square? ( n -- ? ) sqrt dup >integer number= ;

: big-sq-adj-primes-diff ( -- list )
    big-adj-primes-diff [ last square? ] lfilter ;

"Adjacent primes under a million whose difference is a square > 36:" print nl
"p1       p2       difference" print
"============================" print
big-sq-adj-primes-diff [ second 1,000,000 < ] lwhile
[ "%-6d   %-6d   %d\n" vprintf ] leach
Output:
Adjacent primes under a million whose difference is a square > 36:

p1       p2       difference
============================
89689    89753    64
107377   107441   64
288583   288647   64
367957   368021   64
381103   381167   64
396733   396833   100
400759   400823   64
445363   445427   64
623107   623171   64
625699   625763   64
637003   637067   64
710713   710777   64
725209   725273   64
779413   779477   64
801883   801947   64
803749   803813   64
821677   821741   64
832519   832583   64
838249   838349   100
844777   844841   64
883807   883871   64
912103   912167   64
919447   919511   64
954763   954827   64
981823   981887   64
997813   997877   64

Fermat

Func Issqr( n ) = if (Sqrt(n))^2=n then 1 else 0 fi.;
i:=3;
j:=3;
while j<1000000 do
    j:=i+2;
    while j < 1000000 do
        if Isprime(j) then
            if j-i>36 and Issqr(j-i) then !!(i,j,j-i) fi;
            i:=j;
        fi;
        j:=j+2;
    od;
od;

FreeBASIC

#include "isprime.bas"

function nextprime( n as uinteger ) as uinteger
    'finds the next prime after n
    if n = 0 then return 2
    if n < 3 then return n + 1
    dim as integer q = n + 2
    while not isprime(q)
        q+=2
    wend
    return q
end function

function issquare( n as uinteger ) as boolean
    if int(sqr(n))^2 = n then return true else return false
end function

dim as uinteger i=3, j=0
while j<1000000
    j = nextprime(i)
    if j-i > 36 and issquare(j-i) then print i, j, j-i
    i = j
wend
Output:

89689 89753 64 107377 107441 64 288583 288647 64 367957 368021 64 381103 381167 64 396733 396833 100 400759 400823 64 445363 445427 64 623107 623171 64 625699 625763 64 637003 637067 64 710713 710777 64 725209 725273 64 779413 779477 64 801883 801947 64 803749 803813 64 821677 821741 64 832519 832583 64 838249 838349 100 844777 844841 64 883807 883871 64 912103 912167 64 919447 919511 64 954763 954827 64 981823 981887 64 997813 997877 64

FutureBasic

Translation of: FreeBasic
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 NextPrime( n as UInt32 ) as UInt32
  if n = 0 then return 2
  if n < 3 then return n + 1
  NSInteger q = n + 2
  while ( fn IsPrime(q) == NO )
    q += 2
  wend
end fn = q

local fn IsSquare( n as UInt32 ) as BOOL
  if int(sqr(n))^2 == n then return YES else return NO
end fn = NO

UInt32 i, j
i = 3 : j = 0

while ( j < 1000000 )
  j = fn NextPrime(i)
  if j - i > 36 && fn IsSquare( j-i ) == YES then printf @"%6lu - %6lu = %2lu", j, i, j-i
  i = j
wend

HandleEvents
Output:
 89753 -  89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64


Go

Translation of: Wren
Library: Go-rcu
package main

import (
    "fmt"
    "math"
    "rcu"
)

func main() {
    limit := 999999
    primes := rcu.Primes(limit)
    fmt.Println("Adjacent primes under 1,000,000 whose difference is a square > 36:")
    for i := 1; i < len(primes); i++ {
        diff := primes[i] - primes[i-1]
        if diff > 36 {
            s := int(math.Sqrt(float64(diff)))
            if diff == s*s {
                cp1 := rcu.Commatize(primes[i])
                cp2 := rcu.Commatize(primes[i-1])
                fmt.Printf("%7s - %7s = %3d = %2d x %2d\n", cp1, cp2, diff, s, s)
            }
        }
    }
}
Output:
Same as Wren example.

GW-BASIC

10 P=3 : P2=0
20 GOSUB 180
30 IF P2>1000000! THEN END
40 R = P2-P
50 IF R > 36 AND INT(SQR(R))^2=R THEN PRINT P,P2,R
60 P=P2
70 GOTO 20
80 REM tests if a number is prime
90 Q=0
100 IF P = 2 THEN Q = 1:RETURN
110 IF P=3 THEN Q=1:RETURN
120 I=1
130 I=I+1
140 IF INT(P/I)*I = P THEN RETURN
150 IF I*I<=P THEN GOTO 130
160 Q = 1
170 RETURN
180 REM finds the next prime after P, result in P2
190 IF P = 0 THEN P2 = 2: RETURN
200 IF P<3 THEN P2 = P + 1: RETURN
210 T = P
220 P = P + 1
230 GOSUB 80
240 IF Q = 1 THEN P2 = P: P = T: RETURN
250 GOTO 220

Haskell

import Data.List.Split ( divvy )

isSquare :: Int -> Bool
isSquare n = (snd $ properFraction $ sqrt $ fromIntegral n) == 0.0

isPrime :: Int -> Bool
isPrime n 
   |n == 2 = True
   |n == 1 = False
   |otherwise = null $ filter (\i -> mod n i == 0 ) [2 .. root]
   where
      root :: Int
      root = floor $ sqrt $ fromIntegral n

solution :: [[Int]]
solution = filter (\li -> isSquare (last li - head li ) &&
 ( last li - head li ) > 36 ) $ divvy 2 1 $ filter isPrime [2..1000000]

printResultLine :: [Int] -> String
printResultLine list = show ( last list ) ++ " - " ++ ( show $ head list )
 ++ " = " ++ ( show ( last list - head list ))

main :: IO ( )
main = do
   let resultPairs = solution
   mapM_ (\li -> putStrLn $ printResultLine li ) resultPairs
Output:
89753 - 89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64


J

   #(,.-~/"1) p:0 1+/~I.(= <.)6.5>.%:2-~/\p:i.p:inv 1e6  NB. count them
26
   (,.-~/"1) p:0 1+/~I.(= <.)6.5>.%:2-~/\p:i.p:inv 1e6   NB. show them
 89689  89753  64
107377 107441  64
288583 288647  64
367957 368021  64
381103 381167  64
396733 396833 100
400759 400823  64
445363 445427  64
623107 623171  64
625699 625763  64
637003 637067  64
710713 710777  64
725209 725273  64
779413 779477  64
801883 801947  64
803749 803813  64
821677 821741  64
832519 832583  64
838249 838349 100
844777 844841  64
883807 883871  64
912103 912167  64
919447 919511  64
954763 954827  64
981823 981887  64
997813 997877  64

In other words: enumerate primes less than 1e6, find the pairwise differences, find where the prime pairs where maximum of their square root and 6.5 is an integer, and list those pairs with their differences.

Java

import java.util.ArrayList;
import java.util.List;

public final class FindAdjacentPrimesWhichDifferByASquareInteger {

	public static void main(String[] args) {
		List<Integer> primes = listPrimeNumbers(1_000_000);
		
		for ( int i = 2; i < primes.size(); i++ ) {
			final int prime2 = primes.get(i - 1);
			final int prime1 = primes.get(i); 
		    final int difference = prime1 - prime2;
		    if ( difference > 36 && isSquare(difference) ) {
		    	System.out.println(String.format("%12s%9s%s",
		    		prime2 + " and ", prime1 + " : ", "difference = " + difference));
		    }
		}
	}
	
	private static boolean isSquare(int number) {
		return Math.pow((int) Math.sqrt(number), 2) == number;
	}
	
	private static List<Integer> listPrimeNumbers(int limit) {
		List<Integer> primes = new ArrayList<Integer>();
		primes.add(2);
		final int halfLimit = ( limit + 1 ) / 2;
		boolean[] composite = new boolean[halfLimit];
		for ( int i = 1, p = 3; i < halfLimit; p += 2, i++ ) {
			if ( ! composite[i] ) {
				primes.add(p);
				for ( int a = i + p; a < halfLimit; a += p ) {
					composite[a] = true;
				}
			}
		}
		return primes;
	}	

}
Output:
  89689 and  89753 : difference = 64
 107377 and 107441 : difference = 64
 288583 and 288647 : difference = 64
 367957 and 368021 : difference = 64
 381103 and 381167 : difference = 64
 396733 and 396833 : difference = 100
 400759 and 400823 : difference = 64
 445363 and 445427 : difference = 64
 623107 and 623171 : difference = 64
 625699 and 625763 : difference = 64
 637003 and 637067 : difference = 64
 710713 and 710777 : difference = 64
 725209 and 725273 : difference = 64
 779413 and 779477 : difference = 64
 801883 and 801947 : difference = 64
 803749 and 803813 : difference = 64
 821677 and 821741 : difference = 64
 832519 and 832583 : difference = 64
 838249 and 838349 : difference = 100
 844777 and 844841 : difference = 64
 883807 and 883871 : difference = 64
 912103 and 912167 : difference = 64
 919447 and 919511 : difference = 64
 954763 and 954827 : difference = 64
 981823 and 981887 : difference = 64
 997813 and 997877 : difference = 64

jq

Works with: jq

Works with gojq, the Go implementation of jq

See Erdős-primes#jq for a suitable definition of `is_prime` as used here.

Preliminaries

def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;

# Primes less than . // infinite
def primes:
  (. // infinite) as $n
  | if $n < 3 then empty
    else 2, (range(3; $n) | select(is_prime))
    end;

The task

# Input is given to primes/0  - to determine the maximum prime to consider
# Output: stream of [$prime, $nextPrime]
def adjacentPrimesWhichDifferBySquare:
  def isSquare: sqrt | . == floor;
  
  foreach primes as $p ( {previous: null};
    .emit = null
    | if .previous != null
         and (($p - .previous) | isSquare)
      then .emit = [.previous, $p]
      else .
      end
      | .previous = $p;
      select(.emit).emit);

# Input is given to primes/0 to determine the maximum prime to consider.
# Gap must be greater than $gap
def task($gap):
  def l: lpad(6);
  "Adjacent primes under \(.) whose difference is a square > \($gap):",
   (adjacentPrimesWhichDifferBySquare
    | (.[1] - .[0]) as $diff
    | select($diff > $gap)
    | "\(.[1]|l) - \(.[0]|l) = \($diff|lpad(4))" ) ;

1E6 | task(36)
Output:

As for #ALGOL_68.

Julia

using Primes

function squareprimegaps(limit)
    pri = primes(limit)
    squares = Set([1; [x * x for x in 2:2:100]])
    diffs = [pri[i] - pri[i - 1] for i in 2:length(pri)]
    squarediffs = sort(unique(filter(n -> n in squares, diffs)))
    println("\n\nSquare prime gaps to $limit:")
    for sq in squarediffs
        i = findfirst(x -> x == sq, diffs)
        n = count(x -> x == sq, diffs)
        if limit == 1000000 && sq > 36
            println("Showing all $n with square difference $sq:")
            pairs = [(pri[i], pri[i + 1]) for i in findall(x -> x == sq, diffs)]
            foreach(p -> print(last(p), first(p) % 4 == 0 ? "\n" : " "), enumerate(pairs))
        else
            println("Square difference $sq: $n found. Example: ($(pri[i]), $(pri[i + 1])).")
        end
    end
end

squareprimegaps(1_000_000)
squareprimegaps(10_000_000_000)
Output:

Square prime gaps to 1000000:
Square difference 1: 1 found. Example: (2, 3).
Square difference 4: 8143 found. Example: (7, 11).
Square difference 16: 2881 found. Example: (1831, 1847).
Square difference 36: 767 found. Example: (9551, 9587).
Showing all 24 with square difference 64:
(89689, 89753) (107377, 107441) (288583, 288647) (367957, 368021)
(381103, 381167) (400759, 400823) (445363, 445427) (623107, 623171)
(625699, 625763) (637003, 637067) (710713, 710777) (725209, 725273)
(779413, 779477) (801883, 801947) (803749, 803813) (821677, 821741)
(832519, 832583) (844777, 844841) (883807, 883871) (912103, 912167)
(919447, 919511) (954763, 954827) (981823, 981887) (997813, 997877)
Showing all 2 with square difference 100:
(396733, 396833) (838249, 838349)

Square prime gaps to 10000000000:
Square difference 1: 1 found. Example: (2, 3).
Square difference 4: 27409998 found. Example: (7, 11).
Square difference 16: 15888305 found. Example: (1831, 1847).
Square difference 36: 11593345 found. Example: (9551, 9587).
Square difference 64: 1434957 found. Example: (89689, 89753).
Square difference 100: 268933 found. Example: (396733, 396833).
Square difference 144: 35563 found. Example: (11981443, 11981587).
Square difference 196: 1254 found. Example: (70396393, 70396589).
Square difference 256: 41 found. Example: (1872851947, 1872852203).

Mathematica /Wolfram Language

ps = Prime[Range[PrimePi[10^6]]];
ps = Partition[ps, 2, 1];
ps = {#1, #2, #2 - #1} & @@@ ps;
ps //= Select[Extract[{3}]/*GreaterThan[36]];
ps //= Select[Extract[{3}]/*Sqrt/*IntegerQ];
ps // Grid
Output:
89689	89753	64
107377	107441	64
288583	288647	64
367957	368021	64
381103	381167	64
396733	396833	100
400759	400823	64
445363	445427	64
623107	623171	64
625699	625763	64
637003	637067	64
710713	710777	64
725209	725273	64
779413	779477	64
801883	801947	64
803749	803813	64
821677	821741	64
832519	832583	64
838249	838349	100
844777	844841	64
883807	883871	64
912103	912167	64
919447	919511	64
954763	954827	64
981823	981887	64
997813	997877	64

PARI/GP

for(i=3,1000000,j=nextprime(i+1);if(isprime(i)&&j-i>36&&issquare(j-i),print(i," ",j," ",j-i)))

Perl

#!/usr/bin/perl

use strict;    # https://rosettacode.org/wiki/Find_adjacents_primes_which_difference_is_square_integer
use warnings;
use ntheory qw( primes is_square );

my $primeref = primes(1e6);
for my $i (1 .. $#$primeref) {
    (my $diff = $primeref->[$i] - $primeref->[$i - 1]) > 36 or next;
    is_square($diff) and print "$primeref->[$i] - $primeref->[$i - 1] = $diff\n";
}
Output:
89753 - 89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64

Phix

with javascript_semantics
constant limit = 1_000_000
sequence primes = get_primes_le(limit),
         square = repeat(false,floor(sqrt(limit)))
integer sq = 7
while sq*sq<=length(square) do
    square[sq*sq] = true
    sq += 1
end while
for i=2 to length(primes) do
    integer p = primes[i],
            q = primes[i-1],
            d = p-q
    if square[d] then
        printf(1,"%6d - %6d = %d\n",{p,q,d})
    end if
end for
Output:
 89753 -  89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64

Python

import math

limit = 1000000
Primes = []
oldPrime = 0
newPrime = 0
x = 0

def sieve(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False

    for i in range(2, int(n ** 0.5) + 1):
        if is_prime[i]:
            for j in range(i*i, n + 1, i):
                is_prime[j] = False

    return [i for i in range(n+1) if is_prime[i]]

def issquare(x):
	n = math.isqrt(x)
	return n * n == x

Primes = sieve(limit)

for n in range(2, len(Primes)):
    pr1 = Primes[n]
    pr2 = Primes[n-1]
    diff = pr1 - pr2
    flag = issquare(diff)
    if (flag == 1 and diff > 36):
       print(str(pr1) + " " + str(pr2) + " diff = " + str(diff))
Output:
89753 89689 diff = 64
107441 107377 diff = 64
288647 288583 diff = 64
368021 367957 diff = 64
381167 381103 diff = 64
396833 396733 diff = 100
400823 400759 diff = 64
445427 445363 diff = 64
623171 623107 diff = 64
625763 625699 diff = 64
637067 637003 diff = 64
710777 710713 diff = 64
725273 725209 diff = 64
779477 779413 diff = 64
801947 801883 diff = 64
803813 803749 diff = 64
821741 821677 diff = 64
832583 832519 diff = 64
838349 838249 diff = 100
844841 844777 diff = 64
883871 883807 diff = 64
912167 912103 diff = 64
919511 919447 diff = 64
954827 954763 diff = 64
981887 981823 diff = 64
997877 997813 diff = 64

Quackery

eratosthenes, isprime, and sqrt are defined at Sieve of Eratosthenes#Quackery.

  1000000 eratosthenes

  0 0
  1000000 times
    [ i^ isprime if
        [ nip i^ 2dup swap -
          dup 36 > iff
            [ dup sqrt dup * = if
                [ 2dup swap
                  2dup - unrot
                  echo say " - "
                  echo say " = "
                  echo cr ] ]
          else drop ] ]
  2drop
Output:
89689 - 89753 = 64
107377 - 107441 = 64
288583 - 288647 = 64
367957 - 368021 = 64
381103 - 381167 = 64
396733 - 396833 = 100
400759 - 400823 = 64
445363 - 445427 = 64
623107 - 623171 = 64
625699 - 625763 = 64
637003 - 637067 = 64
710713 - 710777 = 64
725209 - 725273 = 64
779413 - 779477 = 64
801883 - 801947 = 64
803749 - 803813 = 64
821677 - 821741 = 64
832519 - 832583 = 64
838249 - 838349 = 100
844777 - 844841 = 64
883807 - 883871 = 64
912103 - 912167 = 64
919447 - 919511 = 64
954763 - 954827 = 64
981823 - 981887 = 64
997813 - 997877 = 64

Raku

use Lingua::EN::Numbers;
use Math::Primesieve;

my $iterator = Math::Primesieve::iterator.new;
my $limit    = 1e10;
my @squares  = (1..30).map: *²;
my $last     = 2;
my @gaps;
my @counts;

loop {
    my $this = (my $p = $iterator.next) - $last;
    quietly @gaps[$this].push($last) if +@gaps[$this] < 10;
    @counts[$this]++;
    last if $p > $limit;
    $last = $p;
}

print "Adjacent primes up to {comma $limit.Int} with a gap value that is a perfect square:";
for @gaps.pairs.grep: { (.key@squares) && .value.defined} -> $p {
    my $ten = (@counts[$p.key] > 10) ?? ', (first ten)' !! '';
    say "\nGap {$p.key}: {comma @counts[$p.key]} found$ten:";
    put join "\n", $p.value.batch(5)».map({"($_, {$_+ $p.key})"})».join(', ');
}
Output:
Adjacent primes up to 10,000,000,000 with a gap value that is a perfect square:
Gap 1: 1 found:
(2, 3)

Gap 4: 27,409,998 found, (first ten):
(7, 11), (13, 17), (19, 23), (37, 41), (43, 47)
(67, 71), (79, 83), (97, 101), (103, 107), (109, 113)

Gap 16: 15,888,305 found, (first ten):
(1831, 1847), (1933, 1949), (2113, 2129), (2221, 2237), (2251, 2267)
(2593, 2609), (2803, 2819), (3121, 3137), (3373, 3389), (3391, 3407)

Gap 36: 11,593,345 found, (first ten):
(9551, 9587), (12853, 12889), (14107, 14143), (15823, 15859), (18803, 18839)
(22193, 22229), (22307, 22343), (22817, 22853), (24281, 24317), (27143, 27179)

Gap 64: 1,434,957 found, (first ten):
(89689, 89753), (107377, 107441), (288583, 288647), (367957, 368021), (381103, 381167)
(400759, 400823), (445363, 445427), (623107, 623171), (625699, 625763), (637003, 637067)

Gap 100: 268,933 found, (first ten):
(396733, 396833), (838249, 838349), (1313467, 1313567), (1648081, 1648181), (1655707, 1655807)
(2345989, 2346089), (2784373, 2784473), (3254959, 3255059), (3595489, 3595589), (4047157, 4047257)

Gap 144: 35,563 found, (first ten):
(11981443, 11981587), (18687587, 18687731), (20024339, 20024483), (20388583, 20388727), (21782503, 21782647)
(25507423, 25507567), (27010003, 27010147), (28716287, 28716431), (31515413, 31515557), (32817493, 32817637)

Gap 196: 1,254 found, (first ten):
(70396393, 70396589), (191186251, 191186447), (208744777, 208744973), (233987851, 233988047), (288568771, 288568967)
(319183093, 319183289), (336075937, 336076133), (339408151, 339408347), (345247753, 345247949), (362956201, 362956397)

Gap 256: 41 found, (first ten):
(1872851947, 1872852203), (2362150363, 2362150619), (2394261637, 2394261893), (2880755131, 2880755387), (2891509333, 2891509589)
(3353981623, 3353981879), (3512569873, 3512570129), (3727051753, 3727052009), (3847458487, 3847458743), (4008610423, 4008610679)

Ring

load "stdlib.ring"
see "working..." + nl
limit = 1000000
Primes = []
oldPrime = 0
newPrime = 0
x = 0

for n = 1 to limit
    if isprime(n)
       add(Primes,n)
    ok
next

for n = 2 to len(Primes)
    pr1 = Primes[n]
    pr2 = Primes[n-1]
    diff = pr1 - pr2
    flag = issquare(diff)
    if flag = 1 and diff > 36
       see "" + pr1 + " " + pr2 + " diff = " + diff + nl
    ok
next

see "done..." + nl

func issquare(x)
     for n = 1 to sqrt(x)
         if x = pow(n,2)
            return 1
         ok
     next
     return 0
Output:
working...
89753 89689 diff = 64
107441 107377 diff = 64
288647 288583 diff = 64
368021 367957 diff = 64
381167 381103 diff = 64
396833 396733 diff = 100
400823 400759 diff = 64
445427 445363 diff = 64
623171 623107 diff = 64
625763 625699 diff = 64
637067 637003 diff = 64
710777 710713 diff = 64
725273 725209 diff = 64
779477 779413 diff = 64
801947 801883 diff = 64
803813 803749 diff = 64
821741 821677 diff = 64
832583 832519 diff = 64
838349 838249 diff = 100
844841 844777 diff = 64
883871 883807 diff = 64
912167 912103 diff = 64
919511 919447 diff = 64
954827 954763 diff = 64
981887 981823 diff = 64
997877 997813 diff = 64
done...

Ruby

require "prime"

Prime.each(1_000_000).each_cons(2) do |a, b|
  diff = b - a
  next unless diff > 36
  isqrt = Integer.sqrt(diff)
  puts "#{b} - #{a} = #{diff}" if isqrt*isqrt == diff
end
Output:
89753 - 89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64

Rust

use prime_tools ;

fn is_square_number( num : u32 ) -> bool {
   let comp_num : f32 = num as f32 ;
   let root = comp_num.sqrt( ) ;
   return root == root.floor( ) ;
}

fn main() {
   let primes: Vec<u32> = prime_tools::get_primes_less_than_x(1000000_u32) ;
   let len = primes.len( ) ;
   let mut i : usize = 0 ;
   while  i < len - 1  {
      let diff : u32 = primes[ i + 1 ] - primes[ i ] ;
      if diff > 36 && is_square_number( diff ) {
         println!("{} - {} = {}" , primes[ i + 1 ] , primes[ i ] , diff) ;
      }
      i += 1 ;
   }
}
Output:
89753 - 89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64

Sidef

var p = 2
var upto = 1e6

each_prime(p.next_prime, upto, {|q|
    if (q-p > 36 && is_square(q-p)) {
        say "#{'%6s' % q} - #{'%6s' % p} = #{'%2s' % isqrt(q-p)}^2"
    }
    p = q
})
Output:
 89753 -  89689 =  8^2
107441 - 107377 =  8^2
288647 - 288583 =  8^2
368021 - 367957 =  8^2
381167 - 381103 =  8^2
396833 - 396733 = 10^2
400823 - 400759 =  8^2
445427 - 445363 =  8^2
623171 - 623107 =  8^2
625763 - 625699 =  8^2
637067 - 637003 =  8^2
710777 - 710713 =  8^2
725273 - 725209 =  8^2
779477 - 779413 =  8^2
801947 - 801883 =  8^2
803813 - 803749 =  8^2
821741 - 821677 =  8^2
832583 - 832519 =  8^2
838349 - 838249 = 10^2
844841 - 844777 =  8^2
883871 - 883807 =  8^2
912167 - 912103 =  8^2
919511 - 919447 =  8^2
954827 - 954763 =  8^2
981887 - 981823 =  8^2
997877 - 997813 =  8^2

Wren

Library: Wren-math
Library: Wren-fmt
import "./math" for Int
import "./fmt" for Fmt

var limit = 1e6 - 1
var primes = Int.primeSieve(limit)
System.print("Adjacent primes under 1,000,000 whose difference is a square > 36:")
for (i in 1...primes.count) {
    var diff = primes[i] - primes[i-1]
    if (diff > 36) {
        var s = diff.sqrt.floor
        if (diff == s * s) {
            Fmt.print ("$,7d - $,7d = $3d = $2d x $2d", primes[i], primes[i-1], diff, s, s)
        }
    }
}
Output:
Adjacent primes under 1,000,000 whose difference is a square > 36:
 89,753 -  89,689 =  64 =  8 x  8
107,441 - 107,377 =  64 =  8 x  8
288,647 - 288,583 =  64 =  8 x  8
368,021 - 367,957 =  64 =  8 x  8
381,167 - 381,103 =  64 =  8 x  8
396,833 - 396,733 = 100 = 10 x 10
400,823 - 400,759 =  64 =  8 x  8
445,427 - 445,363 =  64 =  8 x  8
623,171 - 623,107 =  64 =  8 x  8
625,763 - 625,699 =  64 =  8 x  8
637,067 - 637,003 =  64 =  8 x  8
710,777 - 710,713 =  64 =  8 x  8
725,273 - 725,209 =  64 =  8 x  8
779,477 - 779,413 =  64 =  8 x  8
801,947 - 801,883 =  64 =  8 x  8
803,813 - 803,749 =  64 =  8 x  8
821,741 - 821,677 =  64 =  8 x  8
832,583 - 832,519 =  64 =  8 x  8
838,349 - 838,249 = 100 = 10 x 10
844,841 - 844,777 =  64 =  8 x  8
883,871 - 883,807 =  64 =  8 x  8
912,167 - 912,103 =  64 =  8 x  8
919,511 - 919,447 =  64 =  8 x  8
954,827 - 954,763 =  64 =  8 x  8
981,887 - 981,823 =  64 =  8 x  8
997,877 - 997,813 =  64 =  8 x  8

XPL0

func IsPrime(N);        \Return 'true' if odd N > 2 is prime
int  N, I;
[for I:= 3 to sqrt(N) do
    [if rem(N/I) = 0 then return false;
    I:= I+1;
    ];
return true;
];

int N, P0, P1, D, RD;
[P0:= 2;
for N:= 3 to 1_000_000-1 do
    [if IsPrime(N) then
        [P1:= N;
        D:= P1 - P0;            \D is even because odd - odd = even
        if D >= 64 then         \the next even square > 36 is 64
            [RD:= sqrt(D);
            if RD*RD = D then
                [IntOut(0, P1);  Text(0, " - ");
                 IntOut(0, P0);  Text(0, " = ");
                 IntOut(0, D);  CrLf(0);
                ];
            ];
        P0:= P1;
        ];
    N:= N+1;                    \step by 1+1 = 2 (for odd numbers)
    ];
]
Output:
89753 - 89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64

Zig

Works with: Zig version 0.14dev
const std = @import("std");

pub fn main() !void {
    var bw = std.io.bufferedWriter(std.io.getStdOut().writer());
    const writer = bw.writer();

    const limit = 1_000_000;
    try writer.print("Adjacent primes under {} whose difference is a square > 36:\n", .{limit});

    var a: u32 = undefined;
    var b: u32 = 3;
    while (b < limit) : (b = a) {
        a = nextPrime(b);
        const diff = a - b;
        if (diff > 36 and isSquare(diff))
            try writer.print("{} - {} = {}\n", .{ a, b, diff });
    }
    try bw.flush();
}

fn nextPrime(n_: anytype) @TypeOf(n_) {
    const T = @TypeOf(n_);
    if (@typeInfo(T) != .int or @typeInfo(T).int.signedness != .unsigned)
        @compileError("nextPrime() expected unsigned integer argument, found " ++ @typeName(T));

    if (n_ < 2) return 2;
    if (n_ == 2) return 3;
    if (n_ % 2 == 0) return n_ + 1;

    var n = n_ + 2;
    while (!isPrime(n))
        n += 2;
    return n;
}

fn isPrime(n: anytype) bool {
    const T = @TypeOf(n);
    if (@typeInfo(T) != .int or @typeInfo(T).int.signedness != .unsigned)
        @compileError("isPrime() expected unsigned integer argument, found " ++ @typeName(T));

    if (n < 2) return false;

    inline for ([3]u3{ 2, 3, 5 }) |p| if (n % p == 0) return n == p;

    const wheel = comptime [_]u3{ 4, 2, 4, 2, 4, 6, 2, 6 };

    var p: T = 7;
    while (true)
        for (wheel) |w| {
            if (p * p > n) return true;
            if (n % p == 0) return false;
            p += w;
        };
}

fn isSquare(n: anytype) bool {
    const T = @TypeOf(n);
    if (@typeInfo(T) != .int or @typeInfo(T).int.signedness != .unsigned)
        @compileError("isSquare() expected unsigned integer argument, found " ++ @typeName(T));

    const root: T = std.math.sqrt(n);
    return root * root == n;
}
Output:
Adjacent primes under 1000000 whose difference is a square > 36:
89753 - 89689 = 64
107441 - 107377 = 64
288647 - 288583 = 64
368021 - 367957 = 64
381167 - 381103 = 64
396833 - 396733 = 100
400823 - 400759 = 64
445427 - 445363 = 64
623171 - 623107 = 64
625763 - 625699 = 64
637067 - 637003 = 64
710777 - 710713 = 64
725273 - 725209 = 64
779477 - 779413 = 64
801947 - 801883 = 64
803813 - 803749 = 64
821741 - 821677 = 64
832583 - 832519 = 64
838349 - 838249 = 100
844841 - 844777 = 64
883871 - 883807 = 64
912167 - 912103 = 64
919511 - 919447 = 64
954827 - 954763 = 64
981887 - 981823 = 64
997877 - 997813 = 64