Find squares n where n+1 is prime
- Task
Find squares n where n+1 is prime and n<1.000
ALGOL 68
<lang algol68>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</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
AutoHotkey
<lang AutoHotkey>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
}</lang>
- Output:
1, 4, 16, 36, 100, 196, 256, 400, 576, 676
AWK
<lang 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)
} </lang>
- Output:
1 4 16 36 100 196 256 400 576 676 Find squares 1-999: 10
BASIC
<lang 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</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
BCPL
<lang 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')
$)</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
C
<lang 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;
}</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
CLU
<lang 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</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
Cowgol
<lang 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;</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
F#
This task uses Extensible Prime Generator (F#) <lang fsharp> // 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 "" </lang>
- Output:
1 4 16 36 100 196 256 400 576 676
Fermat
<lang fermat>!!1; i:=2; i2:=4; while i2<1000 do
if Isprime(i2+1) then !!i2 fi; i:+2; i2:=i^2;
od;</lang>
- Output:
14 16 36 100 196 256 400 576 676
FreeBASIC
<lang 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</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
GW-BASIC
<lang gwbasic>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</lang>
- Output:
14 16 36 100 196 256 400 576 676
J
<lang j>((<.=])@%:#+)@(i.&.(p:^:_1)-1:) 1000</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
jq
Works with gojq, the Go implementation of jq
See Erdős-primes#jq for a suitable definition of `is_prime` as used here.
<lang jq>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</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
Julia
<lang julia>using Primes
isintegersquarebeforeprime(n) = isqrt(n)^2 == n && isprime(n + 1)
foreach(p -> print(lpad(last(p), 5)), filter(isintegersquarebeforeprime, 1:1000))
</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
MAD
<lang 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</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
Mathematica / Wolfram Language
<lang Mathematica>Cases[Table[n^2, {n, 101}], _?(PrimeQ[# + 1] &)]</lang>
- Output:
{1,4,16,36,100,196,256,400,576,676,1296,1600,2916,3136,4356,5476,7056,8100,8836}
Modula-2
<lang modula2>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.</lang>
- 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.
<lang parigp>for(n = 1, 1000, if(issquare(n)&&isprime(n+1),print(n)))</lang>
- Output:
14 16 36 100 196 256 400 576
676
Perl
Simple and Clear
<lang perl>#!/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";</lang>
- Output:
1 4 16 36 100 196 256 400 576 676
More Than One Way
TMTOWTDI, right? So do it. <lang perl>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;</lang>
- 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
<lang 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...") </lang>
- 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... <lang perl6>say ({$++²}…^*>Ⅿ).grep: (*+1).is-prime</lang>
- 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... <lang perl6>put (^Ⅿ).grep(*.is-prime).map(*-1).batch(20)».fmt("%3d").join: "\n"</lang>
- 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
<lang 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
</lang>
- Output:
working... 1 4 16 36 100 196 256 400 576 676 Found 10 numbers done...
Tiny BASIC
<lang tinybasic> 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</lang>
- Output:
14 16 36 100 196 256 400 576
676
Wren
<lang ecmascript>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)</lang>
- Output:
There are 10 square numbers 'n' where 'n+1' is prime, viz: [1, 4, 16, 36, 100, 196, 256, 400, 576, 676]
XPL0
<lang 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, ^ )];
]</lang>
- Output:
1 4 16 36 100 196 256 400 576 676