Safe and Sophie Germain primes

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.

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

Works with: Dyalog 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(11+↕)|)
Sophie  Prime  (Prime 1+2×⊢)
•Show 510⥊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

Works with: Delphi version 6.0


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

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

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

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

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 

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

Perl

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

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

Works with: HP version 49g
≪ 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}

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

Library: Wren-math
Library: Wren-fmt
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