I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Safe and Sophie Germain primes

From Rosetta Code
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[edit]

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 

ALGOL 68[edit]

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

Arturo[edit]

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[edit]

 
# 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

Factor[edit]

Works with: Factor version 0.99 2022-04-03
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[edit]

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[edit]

FreeBASIC[edit]

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

GW-BASIC[edit]

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[edit]

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[edit]

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[edit]

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[edit]

Translation of: Wren
Library: Go-rcu
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 

jq[edit]

Works with: 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[edit]

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

Mathematica / Wolfram Language[edit]

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

PARI/GP[edit]

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)

Perl[edit]

Library: ntheory
#!/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[edit]

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

Python[edit]

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

Raku[edit]

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

Ring[edit]

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

Wren[edit]

Library: Wren-math
Library: Wren-seq
Library: Wren-fmt
import "./math" for Int
import "./seq" for Lst
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:")
for (chunk in Lst.chunks(sgp, 10)) Fmt.print("$,5d", chunk)
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[edit]

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