Safe and Sophie Germain primes
Safe and Sophie Germain primes 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.
A prime number p is a Sophie Germain prime if 2p + 1 is also prime.
The number 2p + 1 associated with a Sophie Germain prime is called a safe prime.
- Task
Generate the first 50 Sophie Germain prime numbers.
- See also
11l
F is_prime(a)
I a == 2
R 1B
I a < 2 | a % 2 == 0
R 0B
L(i) (3 .. Int(sqrt(a))).step(2)
I a % i == 0
R 0B
R 1B
V cnt = 0
L(n) 1..
I is_prime(n) & is_prime(2 * n + 1)
print(n, end' ‘ ’)
I ++cnt == 50
L.break
print()- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
8080 Assembly
NPRIMES equ 3000
org 100h
lxi b,NPRIMES
lxi h,primes
mvi e,0
zerolp: mov m,e
inx h
dcx b
mov a,b
ora c
jnz zerolp
lxi b,-NPRIMES
lxi d,2
sieve: mov h,d
mov l,e
sieve2: dad d
push h
lxi b,primes
dad b
inr m
pop h
push h
lxi b,-NPRIMES
dad b
pop h
jnc sieve2
inx d
mov h,d
mov l,e
dad b
jnc sieve
lxi b,50*256+5
lxi d,1
search: inx d
lxi h,primes
dad d
mov a,m
ora a
jnz search
mov h,d
mov l,e
dad d
inx h
push b
lxi b,primes
dad b
pop b
mov a,m
ora a
jnz search
push d
push b
lxi h,buf
push h
lxi b,-10
itoa: xchg
lxi d,-1
digit: inx d
dad b
jc digit
mvi a,'0'+10
add l
pop h
dcx h
mov m,a
push h
mov a,d
ora e
jnz itoa
pop d
mvi c,9
call 5
pop b
dcr c
push b
jz prnl
mvi e,9
mvi c,2
call 5
pop b
jmp next
prnl: lxi d,nl
mvi c,9
call 5
pop b
mvi c,5
next: pop d
dcr b
jnz search
ret
db '.....'
buf: db '$'
nl: db 13,10,'$'
primes: equ $
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
ABC
HOW TO REPORT prime n:
IF n<2: FAIL
REPORT NO d IN {2..floor (root n)} HAS n mod d = 0
HOW TO REPORT sophie.germain.prime n:
REPORT prime n AND prime (2*n + 1)
PUT 50 IN amount
PUT 1 IN n
WHILE amount > 0:
WHILE NOT sophie.germain.prime n: PUT n+1 IN n
WRITE n >> 6
IF amount mod 10 = 1: WRITE /
PUT amount-1 IN amount
PUT n+1 IN n- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
ALGOL 68
BEGIN # find some Sophie Germain primes: primes p such that 2p + 1 is prime #
PR read "primes.incl.a68" PR
[]BOOL prime = PRIMESIEVE 10 000; # hopefully, enough primes #
INT sg count := 0;
FOR p WHILE sg count < 50 DO # find the first 50 Sophie Germain primes #
IF prime[ p ] THEN
IF prime[ p + p + 1 ] THEN
print( ( " ", whole( p, -6 ) ) );
IF ( sg count +:= 1 ) MOD 12 = 0 THEN print( ( newline ) ) FI
FI
FI
OD
END- Output:
2 3 5 11 23 29 41 53 83 89 113 131
173 179 191 233 239 251 281 293 359 419 431 443
491 509 593 641 653 659 683 719 743 761 809 911
953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451
1481 1499
APL
5 10⍴(⊢(/⍨)⊢∊⍨1+2×⊢) ((⊢(/⍨)⊢(~∊)∘.×⍨) 1↓⍳) 3000
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Arturo
sophieG?: function [p][
and? [prime? p][prime? 1 + 2*p]
]
sophieGermaines: new [2]
i: 3
while [50 > size sophieGermaines][
if sophieG? i ->
'sophieGermaines ++ i
i: i + 2
]
loop split.every:10 sophieGermaines 'a ->
print map a => [pad to :string & 4]
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
AWK
# syntax: GAWK -f SAFE_AND_SOPHIE_GERMAIN_PRIMES.AWK
BEGIN {
limit = 50
printf("The first %d Sophie Germain primes:\n",limit)
while (count < limit) {
if (is_prime(++i)) {
if (is_prime(i+i+1)) {
printf("%5d%1s",i,++count%10?"":"\n")
}
}
}
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:
The first 50 Sophie Germain primes:
2 3 5 11 23 29 41 53 83 89
113 131 173 179 191 233 239 251 281 293
359 419 431 443 491 509 593 641 653 659
683 719 743 761 809 911 953 1013 1019 1031
1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
BQN
Prime ← (2≤⊢)∧(¬0⊑∘∊(1↓1+↕∘⌊∘√)⊸|)
Sophie ← Prime ∧ (Prime 1+2×⊢)
•Show 5‿10⥊Sophie¨⊸/↕3000
- Output:
┌─
╵ 2 3 5 11 23 29 41 53 83 89
113 131 173 179 191 233 239 251 281 293
359 419 431 443 491 509 593 641 653 659
683 719 743 761 809 911 953 1013 1019 1031
1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
┘
C
#include <stdbool.h>
#include <stdio.h>
bool prime(unsigned n) {
if (n < 2) return false;
if (n % 2 == 0) return n == 2;
if (n % 3 == 0) return n == 3;
for (unsigned d = 5; d*d <= n; d += 4) {
if (n % d == 0) return false;
d += 2;
if (n % d == 0) return false;
}
return true;
}
bool sophie_germain(unsigned n) {
return prime(n) && prime(2*n + 1);
}
int main(void) {
unsigned n = 0;
for (int amount = 50; amount > 0; amount--) {
do { n++; } while (!sophie_germain(n));
printf("%-6u", n);
if (amount % 10 == 1) printf("\n");
}
return 0;
}
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Cowgol
include "cowgol.coh";
sub prime(n: uint16): (r: uint8) is
if n<2 then r := 0; return; end if;
if n==2 then r := 1; return; end if;
if n%2==0 then r := 0; return; end if;
var d: uint16 := 3;
while d*d <= n loop
if n%d == 0 then r := 0; return; end if
d := d + 2;
end loop;
r := 1;
end sub;
sub sophie_germain(n: uint16): (r: uint8) is
r := 0;
if prime(n) != 0 and prime(2*n+1) != 0 then
r := 1;
end if;
end sub;
var n: uint16 := 1;
var amount: uint16 := 50;
while amount > 0 loop
if sophie_germain(n) != 0 then
print_i16(n);
if amount % 10 == 1
then print_nl();
else print_char('\t');
end if;
amount := amount - 1;
end if;
n := n + 1;
end loop;
print_nl();- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Delphi
function IsPrime(N: int64): boolean;
{Fast, optimised prime test}
var I,Stop: int64;
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+0.0));
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 SophieGermainPrimes(Memo: TMemo);
var I,Cnt: integer;
var S: string;
begin
Cnt:=0;
S:='';
for I:=0 to high(integer) do
if IsPrime(I) then
if IsPrime(2 * I + 1) then
begin
Inc(Cnt);
S:=S+Format('%5D',[I]);
if Cnt>=50 then break;
If (Cnt mod 5)=0 then S:=S+CRLF;
end;
Memo.Lines.Add(S);
Memo.Lines.Add('Count = '+IntToStr(Cnt));
end;
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499 Count = 50 Elapsed Time: 2.520 ms.
Draco
proc prime(word n) bool:
word d;
bool is_prime;
if n<2 then false
elif n%2 = 0 then n=2
elif n%3 = 0 then n=3
else
d := 5;
is_prime := true;
while
is_prime and d*d <= n
do
is_prime := is_prime and n % d /= 0;
d := d + 2;
is_prime := is_prime and n % d /= 0;
d := d + 4
od;
is_prime
fi
corp
proc sophie(word n) bool:
prime(n) and prime(2*n + 1)
corp
proc main() void:
word n, i;
n := 0;
for i from 1 upto 50 do
while n := n + 1; not sophie(n) do od;
write(n:6);
if i % 10 = 0 then writeln() fi
od
corp- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
EasyLang
fastfunc isprim num .
i = 2
while i <= sqrt num
if num mod i = 0
return 0
.
i += 1
.
return 1
.
n = 2
repeat
if isprim n = 1 and isprim (2 * n + 1) = 1
write n & " "
cnt += 1
.
until cnt = 50
n += 1
.
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Factor
USING: lists lists.lazy math math.primes math.primes.lists prettyprint ;
50 lprimes [ 2 * 1 + prime? ] lfilter ltake [ . ] leach
- Output:
2 3 5 ... 1451 1481 1499
Fermat
c:=1;
n:=3;
!!2;
while c<50 do
if Isprime(n) and Isprime(2*n+1) then
c:+;
!!n;
fi;
n:+2;
od;BASIC
FreeBASIC
function isprime(n as integer) as boolean
if n < 2 then return false
if n < 4 then return true
if n mod 2 = 0 then return false
dim as uinteger i = 1
while i*i<=n
i+=2
if n mod i = 0 then return false
wend
return true
end function
function is_sg( n as integer ) as boolean
if not isprime(n) then return false
return isprime(2*n+1)
end function
dim as uinteger c = 1, i = 1
print "2 ";
while c<50
i+=2
if is_sg(i) then
print i;" ";
c+=1
if c mod 10 = 0 then print
end if
wend
- Output:
2 3 5 11 23 29 41 53 83 89113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031
1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
GW-BASIC
10 PRINT "2 ";
20 C = 1
30 N = 3
40 WHILE C < 51
50 P = N
60 GOSUB 170
70 IF Z = 0 THEN GOTO 140
80 P = 2 * N + 1
90 GOSUB 170
100 IF Z = 0 THEN GOTO 140
110 C = C + 1
120 PRINT N;" ";
130 IF C MOD 10 = 0 THEN PRINT
140 N = N + 2
150 WEND
160 END
170 Z = 0
180 IF P < 2 THEN RETURN
190 Z = 1
200 IF P < 4 THEN RETURN
210 Z = 0
220 IF P MOD 2 = 0 THEN RETURN
230 I = 3
240 WHILE I*I<P
250 IF P MOD I = 0 THEN RETURN
260 I = I + 1
270 WEND
280 Z = 1
290 RETURN
BASIC256
function isPrime(v)
if v < 2 then return False
if v mod 2 = 0 then return v = 2
if v mod 3 = 0 then return v = 3
d = 5
while d * d <= v
if v mod d = 0 then return False else d += 2
end while
return True
end function
function isSG(n)
if not isPrime(n) then return False
return isPrime(2*n+1)
end function
c = 1
i = 1
print "2 ";
while c < 50
i += 2
if isSG(i) then
print i; chr(9);
c += 1
if c mod 10 = 0 then print
end if
end while
end
PureBasic
Procedure isPrime(v.i)
If v <= 1 : ProcedureReturn #False
ElseIf v < 4 : ProcedureReturn #True
ElseIf v % 2 = 0 : ProcedureReturn #False
ElseIf v < 9 : ProcedureReturn #True
ElseIf v % 3 = 0 : ProcedureReturn #False
Else
Protected r = Round(Sqr(v), #PB_Round_Down)
Protected f = 5
While f <= r
If v % f = 0 Or v % (f + 2) = 0 :
ProcedureReturn #False
EndIf
f + 6
Wend
EndIf
ProcedureReturn #True
EndProcedure
Procedure isSG(n.i)
If Not isPrime(n) : ProcedureReturn #False : EndIf
ProcedureReturn isPrime(2*n+1)
EndProcedure
OpenConsole()
c.i = 1
i.i = 1
Print("2 ")
While c < 50
i + 2
If isSG(i):
Print(Str(i) + #TAB$)
c + 1
If c % 10 = 0 : PrintN("") : EndIf
EndIf
Wend
Input()
CloseConsole()
Yabasic
sub isPrime(v)
if v < 2 then return False : fi
if mod(v, 2) = 0 then return v = 2 : fi
if mod(v, 3) = 0 then return v = 3 : fi
d = 5
while d * d <= v
if mod(v, d) = 0 then return False else d = d + 2 : fi
wend
return True
end sub
sub isSG(n)
if not isPrime(n) then return False : fi
return isPrime(2*n+1)
end sub
c = 1
i = 1
print "2 ";
while c < 50
i = i + 2
if isSG(i) then
print i, " ";
c = c + 1
if mod(c, 10) = 0 then print : fi
endif
wend
end
Go
package main
import (
"fmt"
"rcu"
)
func main() {
var sgp []int
p := 2
count := 0
for count < 50 {
if rcu.IsPrime(p) && rcu.IsPrime(2*p+1) {
sgp = append(sgp, p)
count++
}
if p != 2 {
p = p + 2
} else {
p = 3
}
}
fmt.Println("The first 50 Sophie Germain primes are:")
for i := 0; i < len(sgp); i++ {
fmt.Printf("%5s ", rcu.Commatize(sgp[i]))
if (i+1)%10 == 0 {
fmt.Println()
}
}
}
- Output:
The first 50 Sophie Germain primes are:
2 3 5 11 23 29 41 53 83 89
113 131 173 179 191 233 239 251 281 293
359 419 431 443 491 509 593 641 653 659
683 719 743 761 809 911 953 1,013 1,019 1,031
1,049 1,103 1,223 1,229 1,289 1,409 1,439 1,451 1,481 1,499
Haskell
isPrime :: Int -> Bool
isPrime n
|n == 2 = True
|otherwise = all (\i -> mod n i /= 0 ) [2..limit]
where
limit = floor $ sqrt $ fromIntegral n
solution :: [Int]
solution = take 50 $ filter (\i -> isPrime i && isPrime ( 2 * i + 1 ) ) [2 , 3 ..]
- Output:
[2,3,5,11,23,29,41,53,83,89,113,131,173,179,191,233,239,251,281,293,359,419,431,443,491,509,593,641,653,659,683,719,743,761,809,911,953,1013,1019,1031,1049,1103,1223,1229,1289,1409,1439,1451,1481,1499]
J
5 10$(#~ 1 2&p. e. ])p:i.1e5
2 3 5 11 23 29 41 53 83 89
113 131 173 179 191 233 239 251 281 293
359 419 431 443 491 509 593 641 653 659
683 719 743 761 809 911 953 1013 1019 1031
1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
jq
Works with gojq, the Go implementation of jq
See e.g. Find_adjacent_primes_which_differ_by_a_square_integer#jq for suitable implementions of `is_prime/0` and `primes/0` as used here.
limit(50; primes | select(2*. + 1|is_prime))- Output:
2 3 5 ... 1451 1481 1499
Julia
using Primes
for (i, p) in enumerate(filter(x -> isprime(2x + 1), primes(1500)))
print(lpad(p, 5), i % 10 == 0 ? "\n" : "")
end
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
MAD
NORMAL MODE IS INTEGER
INTERNAL FUNCTION(X)
ENTRY TO PRIME.
WHENEVER X.L.2, FUNCTION RETURN 0B
WHENEVER X.E.X/2*2, FUNCTION RETURN X.E.2
WHENEVER X.E.X/3*3, FUNCTION RETURN X.E.3
D = 5
TEST WHENEVER D*D.G.X, FUNCTION RETURN 1B
WHENEVER X.E.X/D*D, FUNCTION RETURN 0B
D = D+2
WHENEVER X.E.X/D*D, FUNCTION RETURN 0B
D = D+4
TRANSFER TO TEST
END OF FUNCTION
INTERNAL FUNCTION(X)
ENTRY TO SOPHIE.
WHENEVER .NOT. PRIME.(X), FUNCTION RETURN 0B
FUNCTION RETURN PRIME.(X*2+1)
END OF FUNCTION
DIMENSION LINE(10)
VECTOR VALUES LINFMT = $10(I6)*$
COL = 0
N = 0
THROUGH LOOP, FOR I=50, -1, I.E.0
SEARCH THROUGH SEARCH, FOR N=N+1, 1, SOPHIE.(N)
LINE(COL) = N
COL = COL + 1
WHENEVER COL.E.10
COL = 0
PRINT FORMAT LINFMT,LINE(0),LINE(1),LINE(2),
0 LINE(3),LINE(4),LINE(5),LINE(6),LINE(7),
1 LINE(8),LINE(9)
LOOP END OF CONDITIONAL
END OF PROGRAM- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Mathematica / Wolfram Language
nextSafe[n_] :=
NestWhile[NextPrime, n + 1, ! (PrimeQ[2 # + 1] && PrimeQ[#]) &]
Labeled[Grid[Partition[NestList[nextSafe, 2, 49], 10],
Alignment -> {Right,
Baseline}], "First 50 Sophie Germain primes:", Top]
- Output:
First 50 Sophie Germain primes:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Maxima
/* Function that generate the pairs below n */
sg_s_pairs(n):=block(
L:makelist([i,2*i+1],i,1,n),
L1:[],
for i from 1 thru length(L) do if map(primep,L[i])=[true,true] then push(L[i],L1),
reverse(L1))$
/* Test case */
/* The first of the pairs is a Sophie Germain pair, first element of the pairs must be extracted */
map(first,sg_s_pairs(1500));
- Output:
[2,3,5,11,23,29,41,53,83,89,113,131,173,179,191,233,239,251,281,293,359,419,431,443,491,509,593,641,653,659,683,719,743,761,809,911,953,1013,1019,1031,1049,1103,1223,1229,1289,1409,1439,1451,1481,1499]
Miranda
main :: [sys_message]
main = [Stdout (table 6 10 (take 50 sophie))]
table :: num->num->[num]->[char]
table cw w ls = lay [concat [rjustify cw (show n) | n<-part] | part<-split w ls]
split :: num->[*]->[[*]]
split n [] = []
split n ls = take n ls:split n (drop n ls)
sophie :: [num]
sophie = [p | p<-primes; (2*p+1) $in primes]
primes :: [num]
primes = sieve [2..] where sieve (p:ns) = p : sieve [n | n<-ns; n mod p>0]
in :: *->[*]->bool
in a [] = False
in a (a:as) = True
in b (a:as) = False, if a>b
in b (a:as) = b $in as- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Modula-2
MODULE SophieGermainPrimes;
FROM InOut IMPORT WriteCard, WriteLn;
VAR n, i: CARDINAL;
PROCEDURE prime(n: CARDINAL): BOOLEAN;
VAR d: CARDINAL;
BEGIN
IF n<2 THEN RETURN FALSE END;
IF n MOD 2=0 THEN RETURN n=2 END;
IF n MOD 3=0 THEN RETURN n=3 END;
d := 5;
WHILE d*d <= n DO
IF n MOD d=0 THEN RETURN FALSE END;
INC(d, 2);
IF n MOD d=0 THEN RETURN FALSE END;
INC(d, 4)
END;
RETURN TRUE
END prime;
PROCEDURE sophieGermain(n: CARDINAL): BOOLEAN;
BEGIN
RETURN prime(n) AND prime(2*n + 1)
END sophieGermain;
BEGIN
n := 0;
FOR i := 1 TO 50 DO
REPEAT INC(n) UNTIL sophieGermain(n);
WriteCard(n, 6);
IF i MOD 10 = 0 THEN WriteLn END
END
END SophieGermainPrimes.
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Nim
import std/strutils
func isPrime(n: Natural): bool =
if n < 2: return false
if (n and 1) == 0: return n == 2
if n mod 3 == 0: return n == 3
var k = 5
var delta = 2
while k * k <= n:
if n mod k == 0: return false
inc k, delta
delta = 6 - delta
result = true
iterator sophieGermainPrimes(): int =
var n = 2
while true:
if isPrime(n) and isPrime(2 * n + 1):
yield n
inc n
echo "First 50 Sophie Germain primes:"
var count = 0
for n in sophieGermainPrimes():
inc count
stdout.write align($n, 4)
stdout.write if count mod 10 == 0: '\n' else: ' '
if count == 50: break
- Output:
First 50 Sophie Germain primes: 2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
PARI/GP
issg(n)=if(isprime(n)&&isprime(1+2*n),1,0)
c = 0
n = 2
while(c<50,if(issg(n),print(n);c=c+1);n=n+1)PascalABC.NET
##
uses school;
var primenumbers := firstprimes(1000);
primes(3000).where(x -> 2 * x + 1 in primenumbers).take(50).println;
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Perl
#!/usr/bin/perl
use strict; # https://rosettacode.org/wiki/Safe_and_Sophie_Germain_primes
use warnings;
use ntheory qw( forprimes is_prime);
my @want;
forprimes { is_prime(2 * $_ + 1) and (50 == push @want, $_)
and print("@want\n" =~ s/.{65}\K /\n/gr) + exit } 2, 1e9;
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Phix
with javascript_semantics function sophie_germain(integer p) return is_prime(2*p+1) end function sequence res = {} integer n = 1 while length(res)<50 do integer p = get_prime(n) if sophie_germain(p) then res &= p end if n += 1 end while printf(1,"First 50: %s\n",{join(shorten(apply(res,sprint),"",5))})
- Output:
First 50: 2 3 5 11 23 ... 1409 1439 1451 1481 1499
PL/I
sophieGermainPrimes: procedure options(main);
prime: procedure(n) returns(bit);
declare (n, d) fixed;
if n<=4 then return(n=2 | n=3);
if mod(n,2)=0 | mod(n,3)=0 then return('0'b);
do d=5 repeat(d+4) while(d*d<=n);
if mod(n,d)=0 then return('0'b);
d = d+2;
if mod(n,d)=0 then return('0'b);
end;
return('1'b);
end prime;
sophie: procedure(n) returns(bit);
declare n fixed;
return(prime(n) & prime(n*2+1));
end sophie;
declare (n, i) fixed;
n = 0;
do i=1 to 50;
do n=n+1 repeat(n+1) while(^sophie(n)); end;
put edit(n) (F(6));
if mod(i, 10)=0 then put skip;
end;
end sophieGermainPrimes;- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
PL/M
100H:
BDOS: PROCEDURE (FN,ARG); DECLARE FN BYTE, ARG ADDRESS; GO TO 5; END BDOS;
EXIT: PROCEDURE; GO TO 0; END EXIT;
PUTCH: PROCEDURE (CHAR); DECLARE CHAR BYTE; CALL BDOS(2, CHAR); END PUTCH;
PUTS: PROCEDURE (STR); DECLARE STR ADDRESS; CALL BDOS(9, STR); END PUTS;
DECLARE FALSE LITERALLY '0', TRUE LITERALLY 'NOT FALSE';
PRINT$NUM: PROCEDURE (N);
DECLARE (N, P) ADDRESS, S (6) BYTE INITIAL ('.....$');
DECLARE DGT BASED P BYTE;
P = .S(5);
DIGIT:
P = P-1;
DGT = N MOD 10 + '0';
IF (N := N/10) > 0 THEN GO TO DIGIT;
CALL PUTS(P);
END PRINT$NUM;
PRIME: PROCEDURE (N) BYTE;
DECLARE (N, D) ADDRESS;
IF N<=4 THEN RETURN N=2 OR N=3;
IF NOT N OR N MOD 3=0 THEN RETURN FALSE;
D = 5;
DO WHILE D*D <= N;
IF N MOD D=0 THEN RETURN FALSE;
D = D+2;
IF N MOD D=0 THEN RETURN FALSE;
D = D+4;
END;
RETURN TRUE;
END PRIME;
SOPHIE$GERMAIN: PROCEDURE (N) BYTE;
DECLARE N ADDRESS;
RETURN PRIME(N) AND PRIME(N+N+1);
END SOPHIE$GERMAIN;
DECLARE I BYTE, N ADDRESS INITIAL (0);
DO I=1 TO 50;
DO WHILE NOT SOPHIE$GERMAIN(N := N+1); END;
CALL PRINT$NUM(N);
IF I MOD 5=0
THEN CALL PUTS(.(13,10,'$'));
ELSE CALL PUTCH(9);
END;
CALL EXIT;
EOF- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Python
print("working...")
row = 0
limit = 1500
Sophie = []
def isPrime(n):
for i in range(2,int(n**0.5)+1):
if n%i==0:
return False
return True
for n in range(2,limit):
p = 2*n + 1
if isPrime(n) and isPrime(p):
Sophie.append(n)
print("Found ",end = "")
print(len(Sophie),end = "")
print(" Safe and Sophie primes.")
print(Sophie)
print("done...")
- Output:
working... Found 50 Safe and Sophie primes. [2, 3, 5, 11, 23, 29, 41, 53, 83, 89, 113, 131, 173, 179, 191, 233, 239, 251, 281, 293, 359, 419, 431, 443, 491, 509, 593, 641, 653, 659, 683, 719, 743, 761, 809, 911, 953, 1013, 1019, 1031, 1049, 1103, 1223, 1229, 1289, 1409, 1439, 1451, 1481, 1499] done...
Quackery
isprime is defined at Primality by trial division#Quackery.
[ temp put [] 0
[ 1+
dup isprime until
dup 2 * 1+ isprime until
dup dip join
over size temp share = until ]
drop
temp release ] is sgprimes ( n --> [ )
50 sgprimes witheach [ echo sp ]- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Raku
put join "\n", (^∞ .grep: { .is-prime && ($_*2+1).is-prime } )[^50].batch(10)».fmt: "%4d";
- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Refal
$ENTRY Go {
= <Table (10 6) <Gen 50 SophieGermain>>;
};
Cell {
s.CW s.N, <Repeat s.CW ' '> <Symb s.N>: e.R,
<Last s.CW e.R>: (e.X) e.C = e.C;
};
Repeat {
0 s.I = ;
s.N s.I = s.I <Repeat <- s.N 1> s.I>;
};
Table {
(s.W s.CW) = ;
(s.W s.CW) e.X, <First s.W e.X>: (e.Row) e.Rest =
<Prout <Each (Cell s.CW) e.Row>>
<Table (s.W s.CW) e.Rest>;
};
Gen {
s.N s.F = <Gen s.N s.F 1>;
0 s.F s.I = ;
s.N s.F s.I, <Mu s.F s.I>: {
True = s.I <Gen <- s.N 1> s.F <+ s.I 1>>;
False = <Gen s.N s.F <+ s.I 1>>;
};
};
SophieGermain {
s.N, <+ 1 <* 2 s.N>>: s.Safe = <And <Prime s.N> <Prime s.Safe>>;
};
Prime {
s.N, <Compare s.N 2>: '-' = False;
s.N = <Prime s.N 2>;
s.N s.D, <Compare s.N <* s.D s.D>>: '-' = True;
s.N s.D, <Mod s.N s.D>: 0 = False;
s.N s.D = <Prime s.N <+ s.D 1>>;
};
And {
True True = True;
s.X s.Y = False;
};
Each {
(e.F) = ;
(e.F) s.I e.X = <Mu e.F s.I> <Each (e.F) e.X>;
};- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Ring
load "stdlib.ring"
see "working..." +nl
row = 0
limit = 1500
Sophie = []
for n = 1 to limit
p = 2*n + 1
if isprime(n) and isprime(p)
add(Sophie,n)
ok
next
see "Found " + len(Sophie) + " Safe and Sophie German primes."+nl
for n = 1 to len(Sophie)
row++
see "" + Sophie[n] + " "
if row % 10 = 0
see nl
ok
next
see "done..." + nl- Output:
working... Found 50 Safe and Sophie primes. 2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499 done...
RPL
≪ DUP + 1 + ISPRIME? ≫ 'SOPHIE?' STO ≪ → function count ≪ { } 2 WHILE OVER SIZE count < REPEAT IF DUP function EVAL THEN SWAP OVER + SWAP END NEXTPRIME END DROP ≫ ≫ 'FIRSTSEQ' STO
≪ SOPHIE? ≫ 50 FIRSTSEQ
{{out}
1: {2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499}
Rust
fn is_prime( num : u32 ) -> bool {
let limit : u32 = (num as f32).sqrt( ).floor( ) as u32 ;
(2..=limit).all( | x | num % x != 0 )
}
fn main() {
let mut sophie_primes : Vec<u32> = Vec::new( ) ;
let mut count : u8 = 0 ;
let mut current : u32 = 2 ;
while count < 50 {
if is_prime( current ) && is_prime( 2 * current + 1 ) {
sophie_primes.push( current ) ;
count += 1 ;
}
current += 1 ;
}
println!("The first 50 Sophie - Germain primes:" ) ;
let mut ct : u8 = 0 ;
for num in sophie_primes {
print!("{:>5}" , num ) ;
ct += 1 ;
if ct % 10 == 0 {
println!("") ;
ct = 0 ;
}
}
}
- Output:
The first 50 Sophie - Germain primes:
2 3 5 11 23 29 41 53 83 89
113 131 173 179 191 233 239 251 281 293
359 419 431 443 491 509 593 641 653 659
683 719 743 761 809 911 953 1013 1019 1031
1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
SETL
program sophie_germain_primes;
amount := 50;
loop until amount = 0 do
if sophie (n +:= 1) then
nprint(lpad(str n, 6));
amount -:= 1;
if amount mod 10=0 then print; end if;
end if;
end loop;
op sophie(n);
return prime n and prime (2*n+1);
end op;
op prime(n);
if n<5 then return n in {2,3}; end if;
if n mod 2=0 or n mod 3=0 then return false; end if;
d := 5;
loop while d*d <= n do
if n mod d=0 then return false; end if;
d +:= 2;
if n mod d=0 then return false; end if;
d +:= 4;
end loop;
return true;
end op;
end program;- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Sidef
^Inf -> lazy.grep{|p| all_prime(p, 2*p + 1) }.first(50).slices(10).each{
.join(', ').say
}
- Output:
2, 3, 5, 11, 23, 29, 41, 53, 83, 89 113, 131, 173, 179, 191, 233, 239, 251, 281, 293 359, 419, 431, 443, 491, 509, 593, 641, 653, 659 683, 719, 743, 761, 809, 911, 953, 1013, 1019, 1031 1049, 1103, 1223, 1229, 1289, 1409, 1439, 1451, 1481, 1499
VTL-2
10 Z=3000
20 I=0
30 I=I+1
40 :I)=0
50 #=I=Z=0*30
60 I=2
70 J=I+I
80 :J)=1
90 J=J+I
100 #=J>Z=0*80
110 I=I+1
120 #=I>Z=0*70
130 N=1
140 I=50
150 N=N+1
160 #=:N)+:N+N+1)>1*150
170 ?=N
180 I=I-1
190 #=I=0*999
200 #=I/5*0=%*230
210 ?=9
220 #=150
230 ?=""
240 #=150- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Wren
import "./math" for Int
import "./fmt" for Fmt
var sgp = []
var p = 2
var count = 0
while (count < 50) {
if (Int.isPrime(p) && Int.isPrime(2*p+1)) {
sgp.add(p)
count = count + 1
}
p = (p != 2) ? p + 2 : 3
}
System.print("The first 50 Sophie Germain primes are:")
Fmt.tprint("$,5d", sgp, 10)
- Output:
The first 50 Sophie Germain primes are:
2 3 5 11 23 29 41 53 83 89
113 131 173 179 191 233 239 251 281 293
359 419 431 443 491 509 593 641 653 659
683 719 743 761 809 911 953 1,013 1,019 1,031
1,049 1,103 1,223 1,229 1,289 1,409 1,439 1,451 1,481 1,499
XPL0
func IsPrime(N); \Return 'true' if N is a prime number
int N, I;
[for I:= 2 to sqrt(N) do
if rem(N/I) = 0 then return false;
return true;
];
int N, Count;
[N:= 2;
Count:= 0;
repeat if IsPrime(N) & IsPrime(2*N+1) then
[IntOut(0, N); ChOut(0, 9\tab\);
Count:= Count+1;
if rem(Count/10) = 0 then CrLf(0);
];
N:= N+1;
until Count >= 50;
]- Output:
2 3 5 11 23 29 41 53 83 89 113 131 173 179 191 233 239 251 281 293 359 419 431 443 491 509 593 641 653 659 683 719 743 761 809 911 953 1013 1019 1031 1049 1103 1223 1229 1289 1409 1439 1451 1481 1499
Categories:
- Draft Programming Tasks
- Prime Numbers
- Mathematics
- 11l
- 8080 Assembly
- ABC
- ALGOL 68
- ALGOL 68-primes
- APL
- Arturo
- AWK
- BQN
- C
- Cowgol
- Delphi
- SysUtils,StdCtrls
- Draco
- EasyLang
- Factor
- Fermat
- BASIC
- FreeBASIC
- GW-BASIC
- BASIC256
- PureBasic
- Yabasic
- Go
- Go-rcu
- Haskell
- J
- Jq
- Julia
- MAD
- Mathematica
- Wolfram Language
- Maxima
- Miranda
- Modula-2
- Nim
- PARI/GP
- PascalABC.NET
- Perl
- Ntheory
- Phix
- PL/I
- PL/M
- Python
- Quackery
- Raku
- Refal
- Ring
- RPL
- Rust
- SETL
- Sidef
- VTL-2
- Wren
- Wren-math
- Wren-fmt
- XPL0