Sophie Safe primes

Question

Description

Write a program or function that takes in a positive integer \$n\$ as input and outputs all Sophie Germain primes that are safe primes less than or equal to \$n\$. A prime number \$p\$ is a Sophie Germain prime if \2ドルp+1\$ is also a prime. A prime number \$p\$ is a safe prime if \$p=2q+1\$, where \$q\$ is also a prime. The output should be a list of the Sophie Germain primes that are also safe primes, in ascending order only.

Test Cases

[20] => 5, 11
[10000] => 5, 11, 23, 83, 179, 359, 719, 1019, 1439, 2039, 2063, 2459, 2819, 2903, 2963, 3023, 3623, 3779, 3803, 3863, 4919, 5399, 5639, 6899, 6983, 7079, 7643, 7823
[1000] => 5, 11, 23, 83, 179, 359, 719

Task and Criterion

Output primes \$p\$ such that both \$\frac{p-1}{2}\$ and \2ドルp+1\$ are prime less then, equal to a specific \$n\$ given in input. Shortest bytes answer will win, in case if bytes are equated, a tie will occur.

• Arnauld: This is A059455.
• Arnauld: Please consider using standard sequence I/O rules for future challenges.
• Neil: @Arnauld But they don't allow for the "all terms whose value does not exceed n" outputs...
• Arnauld: @Neil True. Maybe that should be added? (But that question belongs to codegolf.meta.)

Answer 1

Japt -f, 10 bytes

[UUÑÄUz]ej

Try it

Answer 2

Wolfram Language (Mathematica), 46 bytes

Select[Range@#,And@@PrimeQ@{#,(#-1)/2,2#+1}&]&

Try it online!

Answer 3

R, 71 bytes

f=\(n)if(n>4)c(f(n-1),n[all(Map(\(y)sum(!y%%2:y)<2,(n+1)*2^(-1:1)-1))])

Attempt This Online!

Heavily inspired by @Dominic van Essen's answer to Imtiaz Germain Primes.

Answer 4

05AB1E, 10 bytes

Lʒx>y2÷)pP

Try it online.

Explanation:

L # Push a list in the range [1, (implicit) input-integer]
 ʒ # Filter it by:
 x # Double the current value (without popping)
 > # And increase it by 1
 y # Push the current value again
 2÷ # Integer-divide it by 2
 ) # Wrap all three values into a list: [n,2n+1,(n-1)/2]
 p # Check for each whether it's a prime number
 P # Check if all three are truthy
 # (after which the filtered list is output implicitly as result)

• Shaggy: Tested on ATO and it doesn't appear to work without the §.
• Kevin Cruijssen: @Shaggy Ah ok. Thanks for the confirmation. Then I'll report it as a bug later today. :)
• Shaggy: Does 05AB1E have floor division? If so, could you save anything by floor dividing by 2 instead?
• Kevin Cruijssen: @Shaggy I already did when I made my reply to you. ;) I noticed that integer-dividing by 2 was a shorter alternative when I made my Java answer.

Answer 5

Vyxal, 10 bytes

':d›n‹1⁄2WæA

Try it Online!

Ports japt. Quite literally get all numbers where n, 2n incremented, and n decremented halved are all prime. Here's a more fun 11 byter

Vyxal, 11 bytes

~æ~≬‹1⁄2æ'd›æ

Try it Online!

Explained

~æ~≬‹1⁄2æ'd›æ
~æ # keep only from range [1, input] where prime
 ~≬--- # from that, keep where:
 ‹1⁄2 # (n - 1) / 2
 æ # is prime
 ' # from that, keep where:
 d› # 2n + 1
 æ # is prime

Answer 6

Java 8, (削除) 117 (削除ここまで) (削除) 115 (削除ここまで) 110 bytes

i->{for(;i-->5;)if(p(i)+p(i-~i)+p(i/2)<4)System.out.println(i);};int p(int i){for(int k=i;k%--i>0;);return i;}

Outputs each Sophie Safe prime on a separated newline in reversed order.
Based on my answer for the Imtiaz Germain Primes challenge.

Try it online.

Explanation:

i->{ // Method with integer parameter and no return-type
 for(;i-->5;) // Loop `i` in the range (input,5]:
 if(p(i) // If `i` is a prime number
 +p(i-~i) // and `2i+1` is a prime number
 +p(i/2) // and `i` integer-divided by 2 is a prime number
 <4) // (by checking if all three are 0 or 1)
 System.out.println(i);} // Print `i` with trailing newline
// Separated method with integer as both parameter and return-type,
// to check whether the given number (≥2) is a prime number (0 or 1)
int p(int i){
 for(int k=i; // Set `k` to the given integer `i`
 k%--i>0;); // Decrease `i` before every iteration with `--i`
 // Keep looping as long as `k` is NOT divisible by `i`
 return i;} // After the loop, return `i`
 // (if it's 1 or 0, it means it's a prime number)

Answer 7

JavaScript (ES6), 71 bytes

f=p=>p>4?f(p-1)+[[,p+' '][(P=x=>p%--x?P(x):x)(p)*P(p-=~p)*P(p>>=2)]]:''

Try it online!

Commented

f = p => // recursive function taking a prime candidate p
p > 4 ? // if p is greater than 4:
 f(p - 1) + // do a recursive call with p - 1
 [[, p + ' '][ // append p followed by a space iff the following
 // result is 1:
 ( P = x => // helper recursive function taking x
 p % --x ? // decrement x; if x is not a divisor of p:
 P(x) // do recursive calls until it is
 : // else:
 x // return x (smallest proper divisor of p)
 )(p) * // first call with p
 P(p -= ~p) * // second call with 2p + 1
 P(p >>= 2) // third call with floor((2p + 1) / 4) which is
 // (p - 1) / 2 if p is odd
 // the product is 1 iff all calls return 1,
 ]] // i.e. all tested numbers are prime
: // else:
 '' // stop

Answer 8

Thunno +, \$ 16 \log_{256}(96) \approx\$ 13.17 bytes

gzt2*1+s2,ZPsTNP

Attempt This Online!

Explanation

gzt2*1+s2,ZPsTNP # Implicit input
g # Filter [1..input] by the following: n
 zt # Triplicate the number n, n, n
 2*1+ # Multiply by two and add one n, n, 2n+1
 s2, # Swap and integer divide by two n, 2n+1, n//2
 ZP # Pair the top two items n, [2n+1, n//2]
 sT # Swap and append to the list [2n+1, n//2, n]
 NP # Are they all prime?

Answer 9

Raku, 41 bytes

{grep {is-prime $_&2*$_+1&$_/2-.5},1..$_}

Try it online!

This filters (greps) the numbers from 1 to the input argument to those where the conjunction (&) of the number ($_), twice that number plus one (2 * $_ + 1) and half the number minus one half ($_ / 2 - .5) is prime. The conjunction is prime if all of its elements are.

Answer 10

J, 28 bytes

[:I.1*/@p:i.>:@+:^:_1 0 1@,]

Attempt This Online!

i.,] Integers from 0 to n.
>:@+: Function f that computes i -> 1+2*i.
f^:_1 0 1 Apply f -1 times (inverse), 0 times, and 1 time and collect results in a 3-row matrix.
1 p: For each integer, is it prime?
*/ Take the product of each column.
I. Indices of 1s in the resulting vector.

Answer 11

Excel, 118 bytes

=LET(
 a,SEQUENCE(2*A1+1),
 b,IF(MMULT(N(MOD(a,TOROW(a))=0),a^0)=2,a),
 FILTER(b,1-ISNA(XMATCH(2*b+1,b)*XMATCH((b-1)/2,b)))
)

Fails where A1>3663.

Answer 12

Retina 0.8.2, 75 bytes

.+
$*
¶11$`$`
4A`
A`^(11+)1円+$
1(1?)
1ドル
A`^(11+)1円+$
A`^1((11+)2円+)1円$
%`1

Try it online! Link is to verbose version of code. Explanation: Based on my Retina answer to Imtiaz Germain Primes.

.+
$*

Convert to unary.

¶11$`$`

Generate all numbers of the form 2n+1 up to the input. Due to golfing, the last entry is not odd, so I have to add 2 to each value so that the final odd entry is 2n+1. (Although, the non-golfy way only goes up to 2n-1 so wouldn't work anyway.)

4A`

Deleting composite numbers won't work when the number to be tested is 0 or 1 so delete the small numbers manually.

A`^(11+)1円+$

Delete the remaining non-prime numbers.

1(1?)
1ドル

Integer divide by 2.

A`^(11+)1円+$

Delete all composite numbers.

A`^1((11+)2円+)1円$

Delete all numbers that are 1 more than 2 times a composite number.

%`1

Convert the results to decimal.

Answer 13

Charcoal, 26 bytes

×ばつ⊕ιE3X2λ⬤...2λ%λν

Try it online! Link is to verbose version of code. Explanation:

 ... Range from
 2 Literal integer `2 to
 N Input integer
 ⊕ Incremented
 ⊘ Halved
 Φ Filtered where
 3 Literal integer `3`
 E Map over implicit range
 2 Literal integer `2`
 X Raised to power
 λ Inner value
 ×ばつ Vectorised multiply by
 ι Outer value
 ⊕ Incremented
 ⊖ Vectorised decrement
 ⬤ All values satisfy
 ... Range from
 2 Literal integer `2`
 λ To inner value
 ⬤ All values satisfy
 λ Inner value
 % Modulo i.e. is not a multiple of
 ν Innermost value
 ⊗ Vectorised double
 ⊕ Vectorised increment
I Cast to string
 Implicitly print

Answer 14

Excel (ms365), 146 bytes

enter image description here

=LET(x,SEQUENCE(A1*2+1),y,MAP(x,LAMBDA(z,SUM(--((z/x)=INT(z/x)))))=2,q,FILTER(x,y),FILTER(q,MAP(q,LAMBDA(p,SUM(--(q=(p-1)/2))*SUM(--(q=2*p+1))))))

---

EDIT:

Since the above is rather slow (thought should work for all given samples), I was looking for ways to gain more speed. Just for fun, I created my own LAMBDA() based function called 'a' within Excel's name-manager based on one of the more simplistic ways to generate prime numbers; the Sieve of Eratosthenes:

=LAMBDA(x,y,z,IF(x=MAX(y),SORT(y),a(SMALL(y,z),FILTER(y,MOD(y,x)+(y=x)),z+1)))

Now, since we can output all primenumbers, we can call this on the worksheet and do some filtering on the output:

=LET(q,a(2,SEQUENCE(A1*2,,2),1),FILTER(q,MAP(q,LAMBDA(p,SUM(--(q=(p-1)/2))*SUM(--(q=2*p+1))))))

The combined byte-count == 173, but it is working faster then the original answer of 146 bytes.

This makes me wonder if I, or someone, else can come up with ways to implement other known ways of generating primes and gain more speed.

• Jos Woolley: =LET(a,SEQUENCE(2*A1+1),b,IF(MMULT(N(MOD(a,TOROW(a))=0),a^0)=2,a),FILTER(b,1-ISNA(XMATCH(2*b+1,b)*XMATCH((b-1)/2,b)))) is 118 bytes, though your version might be able to cope with an input of 10000 (although looking doubtful on my machine).
• JvdV: @JosWoolley, I did get a result using my (unfortunately more verbose) version. Your version gave me an error pop-up.
• Jos Woolley: Nice implementation of the Sieve of Eratosthenes! I just wonder whether the recursion limit for LAMBDAs will mean that this set-up is too restrictive? I understand that limit is 1024/x, where x is the number of parameters passed to LAMBDA. Since you're passing three parameters here, I presume that means a maximum of 341 iterations. Not sure how that translates into an upper bound on the largest integer parsable, though.
• JvdV: @JosWoolley would be interesting to test, I've only tried the upperbound as 10000*2, which worked relatively fast.
• JvdV: @JosWoolley, FYI the ubound happened to be n == 19432 =). I did expect more of this..

Answer 15

Factor + math.primes, 63 bytes

[ primes-upto [ dup 2/ swap 2 * 1 + [ prime? ] both? ] filter ]

Try it online!

• primes-upto get a list of primes up to the input
• [ ... ] filter select \$p\$ where...
• dup 2/ \$\lfloor{p/2}\rfloor\$
• swap 2 * 1 + and \2ドルp+1\$
• [ prime? ] both? are both prime

Answer 16

Arturo, 45 bytes

$=>[select..1&'n[every?@[1+n*2n/2n]=>prime?]]

Try it

$=>[ ; a function
 select..1&'n[ ; select numbers from 1 to <input>, assign current elt to n
 every? ; is every element in the block...
 @[1+n*2n/2n] ; ... [n/2, n, n*2+1] (where / is integer division)
 =>prime? ; prime?
 ] ; end select
] ; end function

Answer 17

Jelly, 16 bytes

ÆRḤ‘$ÆP$Ƈ’H$ÆP$Ƈ

Try it online!

ÆR - list of prime numbers up to the input
Ḥ‘ - 2x+1
$ÆP$Ƈ - evaluates as primes 
’H - (x-1)/2
$ÆP$Ƈ - evaluates as primes (again)

Answer 18

SageMath, 76 bytes

Run it on SageMathCell!

g=lambda n:[x for x in[1..n-1]if all(is_prime(y)for y in[x,(x-1)//2,2*x+1])]

Answer 19

Uiua, 63 Bytes

×ばつ3∩∩∵(¬∊0↘2◿⇡.)×ばつ2⍥.4↘3+1⇡

Try It

×ばつ3∩∩∵(¬∊0↘2◿⇡.)×ばつ2⍥.4↘3+1⇡
 +1⇡ # generate all numbers from [1, n]
 ↘3 # drop first three since the prime number function will flag 3 (it returns true for 1)
 ⍥.4 # make four copies, 3 for checking 2p+1, p, and floor(p/2), 
 # one extra so that we can use a double both (∩) later
 ×ばつ2 # calculate 2p+1
 ⌊÷2: # calculate floor(p/2)
 ∩∩∵(¬∊0↘2◿⇡.) # apply prime number function to the top four stack elements
 . # copy so we can divide by the number later
 ⇡ # generate numbers [0, x] where x is the number we are checking
 ◿ # calculate mod, to see if it evenly divides
 ↘2 # drop first two so we discard mod by 0 and 1
 ∊0 # check if there is remainder of 0
 ¬ # if there exists no divisors which yield a remainder of zero, it is prime
 ×ばつ3 # use multiplication to and conditions together: prime(floor(p/2)) & prime(2p+1) & prime(p) & prime(p)
▽ # only keep numbers that matched condition