Find adjacent primes which differ by a square integer

From Rosetta Code
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


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

ALGOL 68

Library: ALGOL 68-primes

<lang algol68>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</lang>

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

C

<lang c>#include<stdio.h>

  1. 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;

}</lang>

Fermat

<lang 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;</lang>

FreeBASIC

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

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

GW-BASIC

<lang gwbasic>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</lang>

Julia

<lang 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)

</lang>

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).

PARI/GP

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

Perl

<lang 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 + 1] - $primeref->[$i - 1] = $diff\n";
 }</lang>
Output:
89759 - 89689 = 64
107449 - 107377 = 64
288649 - 288583 = 64
368029 - 367957 = 64
381169 - 381103 = 64
396871 - 396733 = 100
400837 - 400759 = 64
445433 - 445363 = 64
623209 - 623107 = 64
625777 - 625699 = 64
637073 - 637003 = 64
710779 - 710713 = 64
725293 - 725209 = 64
779489 - 779413 = 64
801949 - 801883 = 64
803819 - 803749 = 64
821747 - 821677 = 64
832591 - 832519 = 64
838351 - 838249 = 100
844847 - 844777 = 64
883877 - 883807 = 64
912173 - 912103 = 64
919519 - 919447 = 64
954829 - 954763 = 64
981889 - 981823 = 64
997879 - 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

<lang python> import math print("working...") limit = 1000000 Primes = [] oldPrime = 0 newPrime = 0 x = 0

def isPrime(n): 

   for i in range(2,int(n**0.5)+1):
       if n%i==0:
           return False
   return True

def issquare(x): for n in range(x): if (x == n*n): return 1 return 0

for n in range(limit): 

   if isPrime(n):
      Primes.append(n)

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))

print("done...") </lang>

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...

Raku

<lang perl6>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(', ');

}</lang>

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

<lang 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

</lang>

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...

Wren

Library: Wren-math
Library: Wren-fmt

<lang ecmascript>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)
       }
   }

}</lang>

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
Retrieved from "https://rosettacode.org/wiki/Find_adjacent_primes_which_differ_by_a_square_integer?oldid=317820"