# Random pair of primes

Given a positive input $$\n\$$, output a random pair of primes whose difference is $$\n\$$. It's fine if there's another prime between them.

Every pair should possibly appear and the program should have zero possibility to fall into infinite loop.

It's promised that such pair exist, but not always infinite, so beware. You also shouldn't assume it's infinite for even $$\n\$$ until a proof exist.

Test cases:

2 -> [3,5], [11,13], [29,31], etc.
3 -> [2,5]
4 -> [3,7], [19,23], etc.
7 -> Undefined behavior


Shortest code win.

• I'm confused, are we supposed to output a single pair $[x, x+n]$, such that both are prime, or should we output all such pairs? Apr 9, 2021 at 12:41
• @l4m2 That limitation seems... arbitrary? If you have an odd gap then one of the primes must be even. Two is the only even prime so odd gaps can exist only in pairs where one member is 2. Thus the only odd gaps that exist are pairs [2, gap + 2] where gap + 2 is prime. Apr 9, 2021 at 13:56
• You also shouldn't assume it's infinite for even $n$ until a proof exist. This new rule is basically invalidating all existing answers, right? Apr 9, 2021 at 14:05
• It's promised that such pair exist, but not always infinite Can you clarify why that means, please? Apr 9, 2021 at 15:02
• I am retracting my VTC after reading Delfad0r's answer, which will work. Would be good to add this as a potential approach. Apr 9, 2021 at 17:36

# R, 1111051037370 69 bytes

Edit: -2 bytes thanks to Robin Ryder (as well as -5 more bytes for the now-discarded previous version)

function(s){while(sum(!(z=T+s*0:1)%%rep(2:T,e=2))>1)T=1+rpois(1,9)
z}


Try it online!

Picks random numbers n until (n,n+s) are both prime.

• c(0,s) can be s*0:1 for -1 byte. Additionally, I think this works 105 bytes. Apr 9, 2021 at 14:57
• @RobinRyder - Lovely! That is much, much nicer. Thankyou. Apr 9, 2021 at 15:03
• !sum()<2 can be sum>1, and you run no risk allowing T to be 0, leading to 67 bytes Apr 9, 2021 at 16:39
• @RobinRyder - Thanks for the >1 tip! Apr 9, 2021 at 21:23
• @RobinRyder - The problem when T is zero is that %%0 is undefined and gives NA. This doesn't happen often when the lambda is set to 9 in rpois, but if you reduce it to 1 like this the problem becomes obvious... Apr 9, 2021 at 21:25

# 05AB1E, 13 bytes

Generates a random integer $$\k\$$ until both $$\k\$$ and $$\k+n\$$ are prime.

[Î‚[>₄Ω#}DpP#


Try it online! Runs the program 1000 times and shows the unique results.

Commented:

[              # loop:
Î‚            # push [0, input]
[    }      # loop:
>          #   increment both integers
₄Ω        #   draw a random digit of 1000
#       #   if this is 1, stop the loop
D     # duplicate the current pair
p    # for both integers: are they prime?
P   # take the product
#  # if both are prime, end the main loop


Old answer that always just generates one random integer, 22 bytes:

∞+∞ø.ΔpP}D[₄Ω#>}DpP_i\


Try it online! Runs the program 100 times to demonstrate different outputs.

∞+∞ø.ΔpP} finds the first valid pair.
[₄Ω#>} is a loops which increments the pair and has a $$\50\%\$$ probability to stop at every iteration.
DpP_i\ removes the new pair if any of the numbers is not prime.

• How do your description deal finite case?
– l4m2
Apr 9, 2021 at 13:15
• This seems to time out on TIO for the case 13. Apr 9, 2021 at 13:38
• It always finds the first valid pair, then adds a random integer from the geometric distribution to each integer and checks whether it found a new pair. [TΩ#>} will terminate with probability $1$, I hope that is fine? I'm not quite sure I understood your question.
– ovs
Apr 9, 2021 at 13:39
• @Etheryte there is no pair of two primes with distance 13, so I don't have to handle it. (It's promised that such pair exist)
– ovs
Apr 9, 2021 at 13:41
• @l4m2 I see that this deals with the finite case like so: it will find the lowest pair, then it will test a bigger pair of integers and terminate with 50% probability regardless of If it's a valid result, and yield the largest valid pair found so far. Thus the lowest pair is output with a high probability, but there is a non zero probability of any valid pair being output. Apr 9, 2021 at 17:50

f$[2..]>>= \x->[2..x] f(x:y)n=do b<-randomIO;last$f y n:[pure(x,x+n)|b,p x,p$x+n] p n=all((>0).mod n)[2..n-1] import System.Random  Try it online! The relevant function is f$[2..]>>= \x->[2..x], which takes n as input and returns a random pair of primes (x,x+n).

Definitely not the shortest possible answer, but what I like about it is that, when given enough time, it will actually generate every possible pair of primes at distance $$\n\$$, no matter how large. Also, unlike my two other answers below, it has the nice advantage of actually taking less than a googol years to produce some output.

## How?

The idea is to iterate over the list

[2..]>>= \x->[2..x]≡[2,2,3,2,3,4,2,3,4,5,2,3,4,5,6,2,3,4,5,6,7,...]


which contains every integer $$\\ge 2\$$ infinitely many times. Whenever we encounter a number $$\x\$$ such that $$\x\$$ and $$\x+n\$$ are prime, we output the pair $$\(x,x+n)\$$ with probability $$\\frac{1}{2}\$$ (b<-randomIO picks a Bool at random), otherwise we keep going. The program will halt with probability $$\1\$$, however the probability of outputting anything other than the first pair is very low, especially for larger values of $$\n\$$. The TIO link above runs the function $$\1,\!000,\!000\$$ times, and collect the unique results.

# Haskell, 115 102 104 bytes

f n=do x<-randomIO;last$f n:[pure(x,x+n)|x>1,p x,p$x+n]
p n=all((>0).mod n)[2..n-1]
import System.Random


Try it online!

The relevant function here is f.

This code will almost never run in a reasonable amount of time, but theoretically, if given enough time, it should output all the prime pairs $$\(x,x+n)\$$ with $$\x<2^{63}\$$ (the System.Random random generator can't generate numbers larger than this).

## How?

This is the naive approach: generate a random number $$\x\$$, check if $$\x\$$ and $$\x+n\$$ are prime, and in this case return $$\(x,x+n)\$$. Otherwise try again.

(2?)
x?n=do d<-randomIO;last$(x+d)?n:[pure(x,x+n)|d>0,x>1,p x,p$x+n]
p n=all((>0).mod n)[2..n-1]
import System.Random


Try it online!

The relevant function is (2?).

A mix of the two solutions above. It has the theoretical guarantees of the first (given enough time, it should generate all the prime pairs with difference $$\n\$$, even if the primes are larger than $$\2^{63}\$$). However, the probability of generating some output within a reasonable amount of time is close to $$\0\$$, except if $$\(2+n)\$$ is prime.

## How?

Start with x=2. If x and x+n are both primes, then return (x,x+n) with probability $$\\frac{1}{2}\$$. Otherwise, add a random number d (which can be positive or negative) to x and repeat.

• For the second approach, if you generated x randomly, why do you only need to output (x,x+n) with probability 1/2? Can't you always output it? Apr 9, 2021 at 15:49
• @DominicvanEssen Ah yes, very true. I didn't think too much about the second approach because I liked the first one better anyway. Thanks! Apr 9, 2021 at 15:51
• For the new third approach, surely there's a problem if fn is odd: from your explanation, it seems that after you find the first (and, in this case, only) solution, there's a 1/2 chance that you'll add a random number and never be able to find another solution...? Apr 9, 2021 at 21:31
• @DominicvanEssen The trick here is that d<-randomIO can also be negative. I should probably edit my explanation to include that. Apr 9, 2021 at 22:18
• Ah - That makes sense now. Thanks for expaining. Apr 9, 2021 at 22:22

# MATL, 16 15 bytes

,1r/kYq]-G-}2M


Outputs the two primes in decreasing order. Try it online!

### Explanation

Let $$\q(k)\$$ denote the $$\k\$$-th prime: $$\q(1)=2\$$, $$\q(2)=3\$$, $$\q(3)=5\$$, $$\\ldots\$$

The code generates two random primes $$\r_1, r_2\$$ independently with distribution

$$\\quad \displaystyle \Pr[r_i = q(k)] = \frac 1 {k\cdot(k+1)}; \quad i = 1,2; \quad k = 1,2,3,\ldots\$$

The distribution of $$\k\$$ is easily obtained as the inverse of a normalized uniform random variable rounded down. Each pair of primes has nonzero probability of occuring.

If the two generated primes happen to differ by the input $$\n\$$, they are displayed and the program ends. Else the process is repeated.

        % Do while
,      %   Do twice
1    %     Push 1
r    %     Push random number between 0 and 1
/    %     Divide
k    %     Round down. Gives a random positive integer, k
Yq   %     k-th prime
]      %   End
-      %   Subtract. Gives the difference between the two primes (*)
G      %   Push program input, n
-      %   Subtract. (This will used as the loop condition) (**)
}        % Finally (execute on loop exit)
2M     %   Push inputs to second-to-last normal function (*): latest two primes
% End (implicit). A new iteration will be run if (**) is nonzero
% Display stack (implicit)


# Husk, 15 bytes

ḟΛṗmoSe+⁰L↓x1İπ


Try it online!

Note: requires input of a random 'seed' to generate different output on different runs

Husk is completely deterministic, and has no built-in random number generator. However, some googling indicates that the digits of Pi have been shown to be 'random' (in a certain statistical sense), if a set of constraints about the distribution of chaotic sequences holds true (1). Although this conjecture is still unproven, it suggests that we can use the digits of Pi as the basis for a random number generator, until proven otherwise (2).

Obviously, Pi is constant, so we will need to 'seed' the random number generator if we wish to obtain different results for successive runs of the program. However, one could argue that this is not substantially different from the 'seeding' of other pseudo-random number generators in other programming languages.

So: our Pi-based random number generator (mL↓x1İπ: try it!) outputs the separation between occurrences of the digit 1 in the decimal representation of Pi, after skipping the first s outputted numbers with s given as the 'seed'. I believe that this series should include all non-negative integers with non-zero frequency. Please correct me if I'm wrong!

With that out of the way, here's the 'random pair of primes' program, which simply outputs the first pair it finds of randomly-picked numbers with spacing n that are both primes:

             İπ     # Get the (decimal) digits of pi as an infinite list;
x1       # split into sublists on every occurrence of the digit '1';
↓         # discard as many initial sublists as the random 'seed';
L          # and get the length of each sublist:
# this is our random sequence.
moSe             # Now, for each random number, make a 2-element list
+⁰           # by combining with itself plus the input value of 'n';
ḟ                   # and output the first 2-element list that satisfies:
Λṗ                 # both elements are primes.


This link tests n=4 for random seeds of 1..100 to confirm that various different outputs can be obtained when there is more-than-one solution;
This link tests n=3 for random seeds of 1..100 to confirm that when there is only one solution, this is always outputted.

(1) https://www.nersc.gov/news-publications/nersc-news/science-news/2001/are-the-digits-of-pi-random/
(2) And if that's not Ok, then this is probably not Ok either: https://codegolf.stackexchange.com/a/223136/95126

# APL (Dyalog Unicode), 40 bytes

{⍵{⍵[⊃(?∘⍴⌷⊢)⍸⍺=-∘.-⍨⍵]}(⊢~∘.×⍨)1↓⍳3300}


Not entirely certain that I understood the problem statement. I run out of memory with primes larger than around 3300. Not great, but it's my first stab at it.

⍵{...} ⍝ passing user input to nested function as left argument
⍵[...] ⍝ retrieving primes for given index values
⊃(?∘⍴⌷⊢) ⍝ randomly selecting pair of primes from list
⍸⍺= ⍝ index value of every difference equal to the user input n
- ⍝ negation
∘.-⍨⍵ ⍝ matrix of all possible differences between primes
(⊢~∘.×⍨) ⍝ list of primes in range (inclusive) found by not being members of the composite matrix (created by outer product)
1↓⍳3300 ⍝ range 2–3300 (determined by tio memory limit, though TryAPL allowed for double that)


Try it online!

# JavaScript (Node.js), 68 bytes

f=n=>(g=d=>k%--d?g(d):d>1)(k=2/Math.random()|0)|g(k+=n)?f(n):[k-n,k]


Try it online!

Still shorter than Arnauld's wrong answer

• @DominicvanEssen No. If I knew such solution exist likely I won't ask
– l4m2
Apr 9, 2021 at 14:35

# Jelly,  16  15 bytes

2‘XḂ$¿,³Ä¹ßẒẠ$?


A full program that prints a random, valid pair - if given enough time and recursion stack/memory (recursion depth is set to $$\2^{30}\$$ in TIOs interpreter/Dennis' GitHub; but we seg-fault before then generally).

Try it online! Note that there is a seg-fault for high n (or always for those invalid n which are odd while n+2 is not prime).

### How?

Chooses a random integer greater than 1, pairs it with the integer that's n greater. Then, if both are prime stops and outputs, otherwise repeats the whole procedure again.

The procedure for choosing a random integer, while golfy, will find low numbers exponentially more often than high numbers, so the lowest possible pair is produced extremely often with even modest n.

2‘XḂ$¿,³Ä¹ßẒẠ$? - Main Link - called with 1 program argument, n
2               - set the left argument, say v, to two
¿          - while...
$- ...this is true: last two links as a monad - f(v): X - random integer in [1,v] Ḃ - is odd? ‘ - ...do: increment ,³ - pair with n -> [v, n] Ä - cumulative sums -> [v, n + v] ? - if...$  - ...condition:
Ẓ    -   is prime (vectorises)
Ạ   -   all?
¹      - ...then: no-op
ß     - ...else: call Main Link again
- implicit print


# Japt, 25 bytes

Warning: This will hang your browser if the random value is unfavorable. It does however not have an infinite loop, just a very, very long runtime in bad cases.

This relies on the given assumption that a satisfactory pair is guaranteed to exist. Starting from a random whole number ensures that any pair can be returned, after that we just move -1, +1, -2, +2 etc until we find a number that satisfies the input.

@T=[X,X+U];Tej}c(1/Mr)c;T
(1/Mr)c   // Starting from a random whole number,
@             }c          // try that number -1, +1, -2, +2 etc until we find one where
T=[X,X+U];               // that number as well as that number plus input
Tej            // are both prime.
;T // Return the resulting array.



Try it here.

# Japt, 14 13 bytes

Takes input as a single integer array. Weighted very heavily towards the earlier pairs in the sequence.

iT Èej}a@m+Xö


Try it