Find squares n where n+1 is prime

Find squares n where n+1 is prime 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 squares n where n+1 is prime and n<1.000

ALGOL 68

Library: ALGOL 68-primes

BEGIN # find squares n where n + 1 is prime #
    PR read "primes.incl.a68" PR
    []BOOL prime = PRIMESIEVE 1 000; # construct a sieve of primes up to 1000 #
    # find the squares 1 less than a prime (ignoring squares of non-integers) #
    # other than 1, the numbers must be even                                  #
    IF prime[ 2 # i.e.: ( 1 * 1 ) + 1 # ] THEN print( ( " 1" ) ) FI;
    FOR i FROM 2 BY 2 TO UPB prime WHILE INT i2 = i * i;
                                         i2 < UPB prime
    DO
        IF prime[ i2 + 1 ] THEN
            print( ( " ", whole( i2, 0 ) ) )
        FI
    OD
END

Output:

 1 4 16 36 100 196 256 400 576 676

AutoHotkey

Translation of: FreeBASIC

n := 0
while ((n2 := (n+=2)**2) < 1000)
    if isPrime(n2+1)
        result .= (result ? ", ":"" ) n2
MsgBox % result := 1 ", " result
return

isPrime(n, i:=2){
    while (i < Sqrt(n)+1)
        if !Mod(n, i++)
            return False
    return True
}

Output:

1, 4, 16, 36, 100, 196, 256, 400, 576, 676

AWK

# syntax: GAWK -f FIND_SQUARES_N_WHERE_N+1_IS_PRIME.AWK
BEGIN {
    start = 1
    stop = 999
    n = 2
    n2 = n^2
    printf("1")
    count++
    while (n2 < stop) {
      if (is_prime(n2+1)) {
        printf(" %d",n2)
        count++
      }
      n += 2
      n2 = n^2
    }
    printf("\nFind squares %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)
}

Output:

1 4 16 36 100 196 256 400 576 676
Find squares 1-999: 10

BASIC

10 DEFINT A-Z: N=1000
20 DIM C(N)
30 FOR P=2 TO SQR(N)
40 IF NOT C(P) THEN FOR C=P*P TO N STEP P: C(C)=1=1: NEXT
50 NEXT
60 FOR I=2 TO N
70 IF C(I) THEN 100
80 X=I-1: R=SQR(X)
90 IF R*R=X THEN PRINT X;
100 NEXT

Output:

 1  4  16  36  100  196  256  400  576  676

BCPL

get "libhdr"
manifest $( MAX = 1000 $)

let isqrt(s) = valof
$(  let x0 = s>>1 and x1 = ?
    if x0 = 0 resultis s
    x1 := (x0 + s/x0)>>1
    while x1<x0
    $(  x0 := x1
        x1 := (x0 + s/x0)>>1
    $)
    resultis x0
$)

let sieve(prime, n) be
$(  0!prime := false
    1!prime := false
    for i = 2 to n do i!prime := true
    for p = 2 to isqrt(n) if p!prime
    $(  let c = p*p
        while c<n
        $(  c!prime := false
            c := c + p
        $)
    $)
$)

let square(n) = valof
$(  let sq = isqrt(n)
    resultis sq*sq = n
$)

let start() be
$(  let prime = vec MAX
    sieve(prime, MAX)
    
    for i=2 to MAX if i!prime
    $(  let sq = i-1
        if square(sq) then writef("%N ",sq)
    $)
    wrch('*N')
$)

Output:

1 4 16 36 100 196 256 400 576 676

C

#include <stdio.h>
#include <stdbool.h>
#include <math.h>

#define MAX 1000

void sieve(int n, bool *prime) {
    prime[0] = prime[1] = false;
    for (int i=2; i<=n; i++) prime[i] = true;
    for (int p=2; p*p<=n; p++) 
        if (prime[p])
            for (int c=p*p; c<=n; c+=p) prime[c] = false;
}

bool square(int n) {
    int sq = sqrt(n);
    return (sq * sq == n);
}

int main() {
    bool prime[MAX + 1];
    sieve(MAX, prime);
    for (int i=2; i<=MAX; i++) if (prime[i]) {
        int sq = i-1;
        if (square(sq)) printf("%d ", sq);
    }
    printf("\n");
    return 0;
}

Output:

1 4 16 36 100 196 256 400 576 676

CLU

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

sieve = proc (n: int) returns (array[int])
    prime: array[bool] := array[bool]$fill(2,n-1,true)
    primes: array[int] := array[int]$predict(1,isqrt(n))
    for p: int in int$from_to(2,isqrt(n)) do
        if prime[p] then
            for c: int in int$from_to_by(p*p,n,p) do
                prime[c] := false
            end
        end
    end
    for p: int in array[bool]$indexes(prime) do
        if prime[p] then array[int]$addh(primes,p) end
    end
    return(primes)
end sieve

is_square = proc (n: int) returns (bool)
    return(isqrt(n) ** 2 = n)
end is_square

start_up = proc ()
    po: stream := stream$primary_output()
    primes: array[int] := sieve(1000)
    
    for prime: int in array[int]$elements(primes) do
        n: int := prime-1
        if is_square(n) then stream$puts(po, int$unparse(n) || " ") end
    end
end start_up

Output:

1 4 16 36 100 196 256 400 576 676

COBOL

       IDENTIFICATION DIVISION.
       PROGRAM-ID. SQUARE-PLUS-1-PRIME.
       
       DATA DIVISION.
       WORKING-STORAGE SECTION.
       01 N                   PIC 999.
       01 P                   PIC 9999 VALUE ZERO.
       01 PRIMETEST.
          03 DSOR             PIC 9999.
          03 PRIME-FLAG       PIC X.
             88 PRIME         VALUE '*'.
          03 DIVTEST          PIC 9999V999.
          03 FILLER           REDEFINES DIVTEST.
             05 FILLER        PIC 9999.
             05 FILLER        PIC 999.
                88 DIVISIBLE  VALUE ZERO.
           
       PROCEDURE DIVISION.
       BEGIN.
           PERFORM CHECK-N VARYING N FROM 1 BY 1 
                UNTIL P IS GREATER THAN 1000.
           STOP RUN.
       
       CHECK-N.
           MULTIPLY N BY N GIVING P.
           ADD 1 TO P.
           PERFORM CHECK-PRIME.
           SUBTRACT 1 FROM P.
           IF PRIME, DISPLAY P.
           
       CHECK-PRIME.
           IF P IS LESS THAN 2, MOVE SPACE TO PRIME-FLAG,
           ELSE, MOVE '*' TO PRIME-FLAG.
           PERFORM CHECK-DSOR VARYING DSOR FROM 2 BY 1
               UNTIL NOT PRIME OR DSOR IS GREATER THAN N.
       
       CHECK-DSOR.
           DIVIDE P BY DSOR GIVING DIVTEST.
           IF DIVISIBLE, MOVE SPACE TO PRIME-FLAG.

Output:

0001
0004
0016
0036
0100
0196
0256
0400
0576
0676

Cowgol

include "cowgol.coh";

const MAX := 1000;

sub isqrt(s: uint16): (x0: uint16) is
    x0 := s>>1;
    if x0 == 0 then
        x0 := s;
    else
        loop
            var x1: uint16 := (x0 + s/x0) >> 1;
            if x1 >= x0 then return; end if;
            x0 := x1;
        end loop;
    end if;
end sub;

var prime: uint8[MAX + 1];
MemSet(&prime[0], 1, @bytesof prime);

var p: uint16 := 2;
while p*p <= MAX loop
    if prime[p] != 0 then
        var c := p*p;
        while c <= MAX loop
            prime[c] := 0;
            c := c + p;
        end loop;
    end if;
    p := p + 1;
end loop;

var i: uint16 := 2;
while i <= MAX loop
    if prime[i] != 0 then
        var sq := i - 1;
        var sqr := isqrt(sq);
        if sqr*sqr == sq then
            print_i16(sq);
            print_nl();
        end if;
    end if;
    i := i + 1;
end loop;

Output:

1
4
16
36
100
196
256
400
576
676

F#

This task uses Extensible Prime Generator (F#)

// Find squares n where n+1 is prime. Nigel Galloway: December 17th., 2021
seq{yield 1; for g in 2..2..30 do let n=g*g in if isPrime(n+1) then yield n}|>Seq.iter(printf "%d "); printfn ""

Output:

1 4 16 36 100 196 256 400 576 676

Fermat

!!1;
i:=2;
i2:=4;
while i2<1000 do
    if Isprime(i2+1) then !!i2 fi;
    i:+2;
    i2:=i^2;
od;

Output:

1
4
16
36
100
196
256
400
576
676

FreeBASIC

function isprime(n as integer) as boolean
    if n<0 then return isprime(-n)
    if n<2 then return false
    if n<4 then return true
    dim as uinteger i=3
    while i*i<n
        if n mod i = 0 then return false
        i+=2
    wend
    return true
end function


print 1;"  ";
dim as integer n=2, n2=4
while n2<1000
    if isprime(1+n2) then print n2;"  ";
    n+=2
    n2=n^2
wend

Output:

 1   4   16   36   100   196   256   400   576   676

Go

Translation of: Wren

Library: Go-rcu

package main

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

func main() {
    var squares []int
    limit := int(math.Sqrt(1000))
    i := 1
    for i <= limit {
        n := i * i
        if rcu.IsPrime(n + 1) {
            squares = append(squares, n)
        }
        if i == 1 {
            i = 2
        } else {
            i += 2
        }
    }
    fmt.Println("There are", len(squares), "square numbers 'n' where 'n+1' is prime, viz:")
    fmt.Println(squares)
}

Output:

There are 10 square numbers 'n' where 'n+1' is prime, viz:
[1 4 16 36 100 196 256 400 576 676]

GW-BASIC

10 PRINT 1
20 N = 2 : N2 = 4
30 WHILE N2 < 1000
40 J = N2+1
50 GOSUB 110
60 IF PRIME = 1 THEN PRINT N2
70 N = N + 2
80 N2 = N*N
90 WEND
100 END
110 PRIME = 0
120 IF J < 2 THEN RETURN
130 PRIME = 1
140 IF J<4 THEN RETURN
150 I=5
160 WHILE I*I<J
170 IF J MOD I = 0 THEN PRIME = 0 : RETURN
180 I=I +2
190 WEND
200 RETURN

Output:

1
4
16
36
100
196
256
400
576
676

J

((<.=])@%:#+)@(i.&.(p:^:_1)-1:) 1000

Output:

1 4 16 36 100 196 256 400 576 676

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.

def squares_for_which_successor_is_prime:
  (. // infinite) as $limit
  | {i:1, sq: 1}
  | while( .sq < $limit; .i += 1 | .sq = .i*.i)
  | .sq
  | select((. + 1)|is_prime) ;

1000 | squares_for_which_successor_is_prime

Output:

1
4
16
36
100
196
256
400
576
676

Julia

using Primes

isintegersquarebeforeprime(n) = isqrt(n)^2 == n && isprime(n + 1)

foreach(p -> print(lpad(last(p), 5)), filter(isintegersquarebeforeprime, 1:1000))

Output:

    1    4   16   36  100  196  256  400  576  676 

MAD

            NORMAL MODE IS INTEGER
            BOOLEAN PRIME
            DIMENSION PRIME(1000)
            
            INTERNAL FUNCTION(S)
            ENTRY TO ISQRT.
            X0 = S/2
            WHENEVER X0.E.0, FUNCTION RETURN S
FNDRT       X1 = (X0 + S/X0)/2
            WHENEVER X1.GE.X0, FUNCTION RETURN X0
            X0 = X1
            TRANSFER TO FNDRT
            END OF FUNCTION
            
            THROUGH INIT, FOR P=2, 1, P.G.1000
INIT        PRIME(P) = 1B
            
            THROUGH SIEVE, FOR P=2, 1, P*P.G.1000
            THROUGH SIEVE, FOR C=P*P, P, C.G.1000
SIEVE       PRIME(C) = 0B

            THROUGH TEST, FOR P=2, 1, P.G.1000
            WHENEVER PRIME(P)
                SQ = P-1
                SQR = ISQRT.(SQ)
                WHENEVER SQR*SQR.E.SQ
                    PRINT FORMAT FMT, SQ
                END OF CONDITIONAL
            END OF CONDITIONAL
TEST        CONTINUE

            VECTOR VALUES FMT = $I4*$
            END OF PROGRAM

Output:

   1
   4
  16
  36
 100
 196
 256
 400
 576
 676

Mathematica / Wolfram Language

Cases[Table[n^2, {n, 101}], _?(PrimeQ[# + 1] &)]

Output:

{1,4,16,36,100,196,256,400,576,676,1296,1600,2916,3136,4356,5476,7056,8100,8836}

Modula-2

MODULE SquareAlmostPrime;
FROM InOut IMPORT WriteCard, WriteLn;
FROM MathLib IMPORT sqrt;

CONST Max = 1000;

VAR prime: ARRAY [0..Max] OF BOOLEAN;
    i, sq: CARDINAL;
    
PROCEDURE Sieve;
    VAR i, j, sqmax: CARDINAL;
BEGIN
    sqmax := TRUNC(sqrt(FLOAT(Max)));
    FOR i := 2 TO Max DO prime[i] := TRUE; END;
    FOR i := 2 TO sqmax DO
        IF prime[i] THEN
            j := i * i;
            WHILE j <= Max DO
                prime[j] := FALSE;
                j := j + i;
            END;
        END;
    END;
END Sieve;

PROCEDURE isSquare(n: CARDINAL): BOOLEAN;
    VAR sq: CARDINAL;
BEGIN
    sq := TRUNC(sqrt(FLOAT(n)));
    RETURN sq * sq = n;
END isSquare;

BEGIN
    Sieve;
    FOR i := 2 TO Max DO
        IF prime[i] THEN
            sq := i-1;
            IF isSquare(sq) THEN
                WriteCard(sq, 4);
                WriteLn;
            END;
        END;
    END;
END SquareAlmostPrime.

Output:

   1
   4
  16
  36
 100
 196
 256
 400
 576
 676

PARI/GP

This is not terribly efficient, but it does show off the issquare and isprime functions.

for(n = 1, 1000, if(issquare(n)&&isprime(n+1),print(n)))

Output:

1
4  
16  
36  
100  
196  
256  
400  
576  

676

Perl

Simple and Clear

#!/usr/bin/perl
use strict; # https://rosettacode.org/wiki/Find_squares_n_where_n%2B1_is_prime
use warnings;
use ntheory qw( primes is_square );

my @answer = grep is_square($_), map $_ - 1, @{ primes(1000) };
print "@answer\n";

Output:

1 4 16 36 100 196 256 400 576 676

More Than One Way

TMTOWTDI, right? So do it.

use strict;
use warnings;
use feature 'say';
use ntheory 'is_prime';

my $a; is_prime $_ and $a = sqrt $_-1 and $a == int $a and say $_-1 for 1..1000; # backwards approach
my $b; do { say $b**2 if is_prime 1 + ++$b**2 } until $b > int sqrt 1000;        # do/until
my $c; while (++$c < int sqrt 1000) { say $c**2 if is_prime 1 + $c**2 }          # while/if
say for map $_**2, grep is_prime 1 + $_**2, 1 .. int sqrt 1000;                  # for/map/grep
for (1 .. int sqrt 1000) { say $_**2 if is_prime 1 + $_**2 }                     # for/if
say $_**2 for grep is_prime 1 + $_**2, 1 .. int sqrt 1000;                       # for/grep
is_prime 1 + $_**2 and say $_**2 for 1 .. int sqrt 1000;                         # and/for
is_prime 1+$_**2&&say$_**2for 1..31;                                             # and/for golf, FTW

# or dispense with the module and find primes the slowest way possible
(1 x (1+$_**2)) !~ /^(11+)\1+$/ and say $_**2 for 1 .. int sqrt 1000;

Output:

In all cases:

1
4
16
36
100
196
256
400
576
676

Phix

with javascript_semantics
sequence res = {1}
integer sq = 4, d = 2
while sq<1000 do
    if is_prime(sq+1) then
        res &= sq
    end if
    d += 4
    sq += 2*d
end while
printf(1,"%V\n",{res})

Output:

{1,4,16,36,100,196,256,400,576,676}

Alternative, same output, but 168 iterations/tests compared to just 16 by the above:

with javascript_semantics
function sq(integer n) return integer(sqrt(n)) end function
pp(filter(sq_sub(get_primes_le(1000),1),sq))

Drop the filter to get the 168 (cheekily humorous) squares-of-integers-and-non-integers result of Raku (and format/arrange them identically):

puts(1,join_by(apply(true,sprintf,{{"%3d"},sq_sub(get_primes_le(1000),1)}),1,20," "))

Python

limit = 1000
print("working...")

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(1,x+1):
		if (x == n*n):
			return 1
	return 0

for n in range(limit-1):
	if issquare(n) and isprime(n+1):
		print(n,end=" ")

print()
print("done...")

Output:

working...
1 4 16 36 100 196 256 400 576 676 
done...

Raku

Use up to to one thousand (1,000) rather than up to one (1.000) as otherwise it would be a pretty short list...

say ({$++²}…^*>Ⅿ).grep: (*+1).is-prime

Output:

(1 4 16 36 100 196 256 400 576 676)

Although, technically, there is absolutely nothing in the task directions specifying that n needs to be the square of an integer. So, more accurately...

put (^Ⅿ).grep(*.is-prime).map(*-1).batch(20)».fmt("%3d").join: "\n"

Output:

  1   2   4   6  10  12  16  18  22  28  30  36  40  42  46  52  58  60  66  70
 72  78  82  88  96 100 102 106 108 112 126 130 136 138 148 150 156 162 166 172
178 180 190 192 196 198 210 222 226 228 232 238 240 250 256 262 268 270 276 280
282 292 306 310 312 316 330 336 346 348 352 358 366 372 378 382 388 396 400 408
418 420 430 432 438 442 448 456 460 462 466 478 486 490 498 502 508 520 522 540
546 556 562 568 570 576 586 592 598 600 606 612 616 618 630 640 642 646 652 658
660 672 676 682 690 700 708 718 726 732 738 742 750 756 760 768 772 786 796 808
810 820 822 826 828 838 852 856 858 862 876 880 882 886 906 910 918 928 936 940
946 952 966 970 976 982 990 996

Ring

load "stdlib.ring"
row = 0
limit = 1000
see "working..." + nl

for n = 1 to limit-1
    if issquare(n) and isprime(n+1)
       row++
       see "" + n +nl
    ok
next

see "Found " + row + " numbers" + nl
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...
1
4
16
36
100
196
256
400
576
676
Found 10 numbers
done...

Sidef

1..1000.isqrt -> map { _**2 }.grep { is_prime(_+1) }.say

Output:

[1, 4, 16, 36, 100, 196, 256, 400, 576, 676]

Tiny BASIC

      PRINT 1
      LET N = 2
      LET M = 4
   10 LET J = M + 1
      GOSUB 20
      IF P = 1 THEN PRINT M
      LET N = N + 2
      LET M = N*N
      IF M < 1000 THEN GOTO 10
      END
   20 LET P = 0
      LET I = 3
   30 IF (J/I)*I = J THEN RETURN
      LET I = I + 2
      IF I*I < J THEN GOTO 30
      LET P = 1
      RETURN

Output:

1
4
16
36
100
196
256
400
576

676

Wren

Library: Wren-math

import "./math" for Int

var squares = []
var limit = 1000.sqrt.floor
var i = 1
while (i <= limit) {
    var n = i * i
    if (Int.isPrime(n+1)) squares.add(n)
    i = (i == 1) ? 2 : i + 2
}
System.print("There are %(squares.count) square numbers 'n' where 'n+1' is prime, viz:")
System.print(squares)

Output:

There are 10 square numbers 'n' where 'n+1' is prime, viz:
[1, 4, 16, 36, 100, 196, 256, 400, 576, 676]

XPL0

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

int  N;
[for N:= 1 to sqrt(1000-1) do
    if IsPrime(N*N+1) then
        [IntOut(0, N*N);  ChOut(0, ^ )];
]

Output:

1 4 16 36 100 196 256 400 576 676