Find squares n where n+1 is prime
Find squares n where n+1 is prime and n<1.000
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
Ada
-- Rosetta Code Task written in Ada
-- Find squares n where n+1 is prime
-- https://rosettacode.org/wiki/Find_squares_n_where_n%2B1_is_prime
-- Loosely translated from the AWK solution
-- September 2024, R. B. E.
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Integer_Text_IO; use Ada.Integer_Text_IO;
procedure N_Squared_Where_N_Plus_One_Is_Prime is
function Is_Prime (N : in Natural) return Boolean is
Temp : Natural := 5;
begin
if N < 2 then
return False;
end if;
if N mod 2 = 0 then
return N = 2;
end if;
if N mod 3 = 00 then
return N = 3;
end if;
while Temp * Temp <= N loop
if N mod Temp = 0 then
return False;
end if;
Temp := Temp + 2;
if N mod Temp = 0 then
return False;
end if;
Temp := Temp + 4;
end loop;
return True;
end Is_Prime;
procedure Display (Start, Stop, Count : Positive) is
begin
New_Line;
Put ("Find squares: ");
Put (Start, 0);
Put ("-");
Put (Stop, 0);
Put (": ");
Put (Count, 3);
New_Line;
end Display;
Count : Natural := 0;
Start : Positive := 1;
Stop : Positive := 999;
N : Positive := 2;
N2 : Positive := N ** 2;
begin
Put (1, 4);
Count := Count + 1;
while (N2 < Stop) loop
if (Is_Prime (N2 + 1)) then
Put (N2, 4);
Count := Count + 1;
end if;
N := N + 2;
N2 := N ** 2;
end loop;
Display (Start, Stop, Count);
end N_Squared_Where_N_Plus_One_Is_Prime;
- Output:
1 4 16 36 100 196 256 400 576 676 Find squares: 1-999: 10
ALGOL 68
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
ALGOL W
Using the difference of two squares to optimise the primality tests and the square loop (similar to the Phix and Applescript samples), though with the small number of values to test, that probably doesn't affect runtime much...
begin % find squares n where n + 1 is prime %
% returns true if n is prime, false otherwise, uses trial division %
logical procedure isPrime ( integer value n ) ;
if n < 3 then n = 2
else if n rem 3 = 0 then n = 3
else if not odd( n ) then false
else begin
logical prime;
integer f, f2, toNext;
prime := true;
f := 5;
f2 := 25;
toNext := 24; % note: ( 2n + 1 )^2 - ( 2n - 1 )^2 = 8n %
while f2 <= n and prime do begin
prime := n rem f not = 0;
f := f + 2;
f2 := toNext;
toNext := toNext + 8
end while_f2_le_n_and_prime ;
prime
end isPrime ;
% other than 1, the numbers must be even %
if isPrime( 2 % i.e.: ( 1 * 1 ) + 1 % ) then write( ( " 1" ) );
begin
integer i2, toNext;
toNext := i2 := 4; % note: ( 2n + 2 )^2 - 2n^2 = 8n + 4 %
while i2 < 1000 do begin
if isPrime( i2 + 1 ) then writeon( i_w := 1, s_w := 0, " ", i2 );
toNext := toNext + 8;
i2 := i2 + toNext
end while_i2_lt_1000
end
end.- Output:
1 4 16 36 100 196 256 400 576 676
AppleScript
on isPrime(n)
if (n < 4) then return (n > 1)
if ((n mod 2 is 0) or (n mod 3 is 0)) then return false
repeat with i from 5 to (n ^ 0.5) div 1 by 6
if ((n mod i is 0) or (n mod (i + 2) is 0)) then return false
end repeat
return true
end isPrime
on task()
set output to {}
if (isPrime(1 * 1 + 1)) then set end of output to 1 * 1
repeat with sqrt from 2 to (1000 ^ 0.5) by 2
set n to sqrt * sqrt
if (isPrime(n + 1)) then set end of output to n
end repeat
return output
end task
task()
- Output:
{1, 4, 16, 36, 100, 196, 256, 400, 576, 676}
The first Phix solution's method of incrementing the square is fun, but not more efficient in AppleScript than the more straightforward method above. It can be optimised slightly by incrementing the d variable by 8 instead of incrementing by 4 and multiplying the result by 2. Also by incrementing the square + 1 each time instead of the square itself, so that 1 only has to be subtracted from the hits instead of being added to every square.
on task()
set output to {1}
set nPlus1 to 5
repeat with d from 12 to (1000 ^ 0.5 div 0.25) by 8
if (isPrime(nPlus1)) then set end of output to nPlus1 - 1
set nPlus1 to nPlus1 + d
end repeat
return output
end task
Arturo
1..31 | select 'x -> prime? 1 + x^2
| map 'x -> x^2
| print
- Output:
1 4 16 36 100 196 256 400 576 676
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
}
- 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
Delphi
function IsPrime(N: integer): boolean;
{Fast, optimised prime test}
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;
procedure ShowPrimeSquares(Memo: TMemo);
var N,S2: integer;
begin
for N:= 1 to Trunc(sqrt(1000-1)) do
begin
S2:=N*N;
if IsPrime(S2+1) then Memo.Text:=Memo.Text+' '+IntToStr(S2);
end;
end;
- Output:
1 4 16 36 100 196 256 400 576 676
EasyLang
fastfunc isprim num .
i = 2
while i <= sqrt num
if num mod i = 0
return 0
.
i += 1
.
return 1
.
n0 = 1
repeat
n = n0 * n0
until n >= 1000
if isprim (n + 1) = 1
write n & " "
.
n0 += 1
.
- 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:
14 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
FutureBasic
local fn IsPrime( n as NSUInteger ) as BOOL
NSUInteger i
if n = 2 then return YES
if n < 2 || n % 2 == 0 then return NO
for i = 3 to int(n^.5) step 2
if n % i == 0 then return NO
next
end fn = YES
printf @"1 \b"
int n = 2, n2 = 4
while ( n2 < 1000 )
if fn IsPrime( 1 + n2 ) then printf @"%d \b", n2
n += 2 : n2 = n^2
wend
HandleEvents- Output:
1 4 16 36 100 196 256 400 576 676
Go
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:
14 16 36 100 196 256 400 576 676
Haskell
module Squares
where
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
isSquare :: Int -> Bool
isSquare n = theFloor * theFloor == n
where
theFloor :: Int
theFloor = floor $ sqrt $ fromIntegral n
solution :: [Int]
solution = [d | d <- [1..999] , isSquare d && isPrime ( d + 1 )]
- 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 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
Nim
import std/strutils
func isPrime(n: Positive): bool =
if n < 2: return false
if (n and 1) == 0: return n == 2
var d = 3
while d * d <= n:
if n mod d == 0:
return false
inc d, 2
result = true
var list = @[1]
var n = 2
var n2 = 4
while n2 < 1000:
if isPrime(n2 + 1):
list.add n2
inc n, 2
n2 = n * n
echo list.join(" ")
- Output:
1 4 16 36 100 196 256 400 576 676
OCaml
let is_prime n =
let rec test x =
x * x > n || n mod x <> 0 && n mod (x + 2) <> 0 && test (x + 6)
in if n < 5 then n lor 1 = 3 else n land 1 <> 0 && n mod 3 <> 0 && test 5
let seq_squares =
let rec next n a () = Seq.Cons (n, next (n + a) (a + 2)) in
next 0 1
let () =
let cond n = is_prime (succ n) in
seq_squares |> Seq.take_while ((>) 1000) |> Seq.filter cond
|> Seq.iter (Printf.printf " %u") |> print_newline
- 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:
14 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," "))
PL/M
... under CP/M (or an emulator)
100H: /* FIND SQUARES N WHERE N + ! IS PRIME */
/* CP/M BDOS SYSTEM CALL AND I/O ROUTINES */
BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
PR$CHAR: PROCEDURE( C ); DECLARE C BYTE; CALL BDOS( 2, C ); END;
PR$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
PR$NL: PROCEDURE; CALL PR$CHAR( 0DH ); CALL PR$CHAR( 0AH ); END;
PR$NUMBER: PROCEDURE( N ); /* PRINTS A NUMBER IN THE MINIMUN FIELD WIDTH */
DECLARE N ADDRESS;
DECLARE V ADDRESS, N$STR ( 6 )BYTE, W BYTE;
V = N;
W = LAST( N$STR );
N$STR( W ) = '$';
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
DO WHILE( ( V := V / 10 ) > 0 );
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
END;
CALL PR$STRING( .N$STR( W ) );
END PR$NUMBER;
/* RETURNS TRUE IF N IS PRIME, FALSE OTHERWISE, USES TRIAL DIVISION */
IS$PRIME: PROCEDURE( N )BYTE;
DECLARE N ADDRESS;
DECLARE PRIME BYTE;
IF N < 3 THEN PRIME = N = 2;
ELSE IF N MOD 3 = 0 THEN PRIME = N = 3;
ELSE IF N MOD 2 = 0 THEN PRIME = 0;
ELSE DO;
DECLARE ( F, F2, TO$NEXT ) ADDRESS;
PRIME = 1;
F = 5;
F2 = 25;
TO$NEXT = 24; /* NOTE: ( 2N + 1 )^2 - ( 2N - 1 )^2 = 8N */
DO WHILE F2 <= N AND PRIME;
PRIME = N MOD F <> 0;
F = F + 2;
F2 = F2 + TO$NEXT;
TO$NEXT = TO$NEXT + 8;
END;
END;
RETURN PRIME;
END IS$PRIME;
/* TASK */
/* OTHER THAN 1, THE NUMBERS MUST BE EVEN */
IF IS$PRIME( 2 /* I.E.: ( 1 * 1 ) + 1 */ ) THEN DO;
CALL PR$CHAR( ' ' );
CALL PR$CHAR( '1' );
END;
DECLARE ( I2, TO$NEXT ) ADDRESS;
TO$NEXT, I2 = 4; /* NOTE: ( 2N + 2 )^2 - 2N^2 = 8N + 4 */
DO WHILE I2 < 1000;
IF IS$PRIME( I2 + 1 ) THEN DO;
CALL PR$CHAR( ' ' );
CALL PR$NUMBER( I2 );
END;
TO$NEXT = TO$NEXT + 8;
I2 = I2 + TO$NEXT;
END;
EOF- Output:
1 4 16 36 100 196 256 400 576 676
PROMAL
;;; Find squares n where n + 1 is prime
PROGRAM primesq
INCLUDE library
;;; returns TRUE(1) if p is prime, FALSE(0) otherwise
FUNC BYTE isPrime
ARG WORD n
WORD i
WORD f
WORD f2
WORD toNext
BYTE prime
BEGIN
IF n < 3
prime = n = 2
ELSE IF n % 3 = 0
prime = n = 3
ELSE IF n % 2 = 0
prime = 0
ELSE
prime = 1
f = 5
f2 = 25
toNext = 24 ; note: ( 2n + 1 )^2 - ( 2n - 1 )^2 = 8n
WHILE f2 <= n AND prime
prime = n % f <> 0
f = f + 2
f2 = toNext
toNext = toNext + 8
RETURN prime
END
WORD i2
WORD toNext
BEGIN
IF isPrime( ( 1 * 1 ) + 1 ) ; 1 is the only possible odd number
OUTPUT " 1"
i2 = 4
toNext = 4 ; note: ( 2n + 2 )^2 - 2n^2 = 8n + 4
WHILE i2 < 1000
IF isPrime( i2 + 1 )
OUTPUT " #W", i2
toNext = toNext + 8
i2 = i2 + toNext
END- Output:
1 4 16 36 100 196 256 400 576 676
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...
Quackery
isprime is defined at Primality by trial division#Quackery.
[] [] 0
[ 1+ dup 2 **
dup 1000 < while
1+ isprime if
[ dup dip join ]
again ]
2drop
witheach [ 2 ** join ]
echo- Output:
[ 1 4 16 36 100 196 256 400 576 676 ]
Racket
#lang racket
(define (find-subprime-squares up-to)
(let rc ([curr-num 1]
[found '()])
(let ([n-sq (* curr-num curr-num)])
(cond [(>= n-sq up-to) (reverse found)]
[(prime? (add1 n-sq)) (rc (add1 curr-num) (cons n-sq found))]
[else (rc (add1 curr-num) found)]))))
(define (prime? n)
(let iter ([counter 2])
(cond [(eq? n 1) #f]
[(<= (expt counter 2) n)
(if (zero? (remainder n counter))
#f
(iter (add1 counter)))]
[else #t])))
(find-subprime-squares 1000)
- Output:
'(1 4 16 36 100 196 256 400 576 676)
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...
RPL
≪ { }
1 1000 √ FOR j
IF j SQ 1 + ISPRIME? THEN j SQ + END
NEXT
≫ 'TASK' STO
- Output:
1: { 1 4 16 36 100 196 256 400 576 676 }
Ruby
require 'prime'
p (1..Integer.sqrt(1000)).filter_map{|n| sqr = n*n; sqr if (sqr+1).prime? }
- Output:
[1, 4, 16, 36, 100, 196, 256, 400, 576, 676]
Rust
use primes::is_prime ;
fn is_square( number : u64 ) -> bool {
let floor : u64 = (number as f64).sqrt( ).floor( ) as u64 ;
floor * floor == number
}
fn main() {
let solution : Vec<u64> = (1..1000).into_iter( ).
filter( | d | is_square( *d ) && is_prime( *d + 1 )).collect( ) ;
println!("{:?}" , solution);
}
- Output:
[1, 4, 16, 36, 100, 196, 256, 400, 576, 676]
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:
14 16 36 100 196 256 400 576
676
VTL-2
1000 ?=1
1010 N=2
1020 M=4
1030 J=M+1
1040 #=2000
1050 #=P=1=0*1080
1060 $=32
1070 ?=M
1080 N=N+2
1090 M=N*N
1100 #=M<1000*1030
1110 #=9999
2000 R=!
2010 P=0
2020 I=3
2030 #=J/I*0+%=0*R
2040 I=I+2
2050 #=I*I<J*2030
2060 P=1
2070 #=R- Output:
1 4 16 36 100 196 256 400 576 676
Wren
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