# Efficient prime k-tuple search ## Prime Constellations Demystified I made this since I have been struggling to find a complete, from start to finish and more importantly beginner friendly explanation on how to search for prime k-tuples or prime constellations quicker, excluding the obvious brute force solution. To follow along, some mathematical background is required such as modular arithmetic, what is a prime constellation is and basic understanding of algorithm complexity. ### The problem to solve The problem that we want to solve is: Given a constellation of length $$k$$ in the form of $$(o\_1, o\_2, ..., o\_k)$$ and a target number $$T$$ find a prime number $$n$$ such that $$(n+o\_1, n+o\_2, ... ,n+o\_k)$$ are all prime and $$n \geq T$$. In our examples we will use the pattern for a prime triplet. You can find further info about the possible patterns [here](https://riecoin.dev/en/Possible_Prime_Constellation_Patterns). ### Required Mathematical Tools #### Primorial A primorial with primorial number $$m$$ is the multiplication of the first $$m$$ primes and it is denoted by $$p\_m #$$. For example $$p\_1 # = 2$$, $$p\_2 # = 2 \cdot 3 = 6$$, $$p\_3 # = 2 \cdot 3 \cdot 5 = 30$$ and so on. We will use the primorial to quickly eliminate candidates by understanding a recurring pattern around the primorial. #### Fermat Primality Test The most important category of tools that we will require is primality testing algorithms. We have multiple options in our disposal, some exact while others are probabilistic. Since [rieMiner](https://github.com/Pttn/rieMiner) is using the Fermat test we will be using that. It seems to be fast enough as well as simple to implement, with negligible probability of a false positive (i.e. a pseudo-prime) and with ever more decreasing probability for false positives the longer the candidates we test are. The Fermat primality test works as follows: If you want to test if $$n$$ is prime 1. Pick a random number $$a$$ such that $$1 < a < n$$ 2. Check for this equality $$a^{n-1} \equiv 1 ; (mod ; n)$$ 3. If true then $$n$$ is prime In our case the following implementation is used: ```python def fermat_prime_p(n): if n == 2: return True if not n & 1: return False if pow(2, n-1, n) == 1: return True return False ``` #### Sieve of Eratosthenes The Sieve of Eratosthenes is a popular algorithm for "sieving" out composite numbers from a prime table up to a given limit. It is usually used to quickly find all prime numbers up to a given limit. One important property of this algorithm is that it trade's off memory with time. The idea is to initialize a boolean array with 1's where each index represents a number and its value represents if its prime or not. The idea is that after we initialize the array with 1's we use known primes such as 2 and 3 and eliminate their multiples by setting them to zero. Then we move on and implement this for the next non eliminated number which in our case will be 5 (another prime). And so on until we reach the end of the table. When we are done we have a table where each prime number is set to 1 and each composite number is set to false. The following function is used in the python script to quickly generate all prime numbers up to `num` and returns a list of all the primes up to `num`. ```python def sieve_of_eratosthenes(num): prime = [True for i in range(num+1)] # boolean array p = 2 while (p * p <= num): # If prime[p] is not # changed, then it is a prime if (prime[p] == True): # Updating all multiples of p for i in range(p * p, num+1, p): prime[i] = False p += 1 primes = [] # Print all prime numbers for p in range(2, num+1): if prime[p]: primes.append(p) return primes ``` We will use the sieve of Eratosthenes to build a prime table in order to be able to quickly use small primes in our calculations. The use-case will become clearer once we see how to sieve factors in the final version of the algorithm. #### Multiplicative Inverse The multiplicative inverse of a number $$a$$ in modular arithmetic is the number $$x$$ such that $$ax \equiv 1 (mod ; n)$$. It can be computed using the extended Euclidean algorithm but in our case we will use the build in python function `pow()`. We require the multiplicative inverse when we compute the inverses of the primorial `p_m` mod each prime `p` in our prime table. This is done in the following function: ```python def inverses_of_eratosthenes(p_m, primes): inverses = [] for p in primes: try: inverse = pow(p_m, -1, p) except: inverse = 0 # if no multiplicative inverse inverses.append(inverse) return inverses ``` ### Trivial Solution: Brute-force The most obvious solution is to brute-force our way into testing each odd $$n \geq T$$ to see if its the base of a constellation with our given pattern. The following function is used to test if $$n$$ is a base of our constellation pattern. ```python def fermat_is_valid_pow(n, constellation_pattern): for offset in constellation_pattern: if not fermat_prime_p(n+offset): return False return True ``` The script provided doesn't implement a brute force search but its quite trivial to do: ```python n_max = 100000 for i in range(5,n_max,2): if fermat_is_valid_pow(i, constellation_pattern): print(i) ``` This is by far the slowest way to do things but its quite important to understand this first before we move on to a more optimized version ### Optimized Solution: Wheel Factorization The idea of this optimization is based on the fact that we would like to test as few $$n$$ as possible using the `fermat_is_valid_pow()` function which will be in a tight loop taking most of the program's time. To do this we would like to sieve out some candidates $$c$$ that are definitely not going to be the base of a prime constellation. We have already done that even in the brute-force version since we ignore even numbers. If we generalize this idea to sieve not only the even numbers, but also odds that are not primes, and even primes that when adding the offset of the desired constellation, produce composites. To be able to sieve out these numbers systematically we need to establish a pattern. To do that we will build a sieve of Eratosthenes, and sieve out $$2,3$$ and $$5$$ , which are the primes the third primorial or $$p\_3#$$ is comprised of. ``` 012345678901234567890123456789 p = 2 xoxoxoxoxoxoxoxoxoxoxoxoxoxoxo p = 3 xooxooxooxooxooxooxooxooxooxoo p = 5 xooooxooooxooooxooooxooooxoooo p_3# xoxxxxxoxxxoxoxxxoxoxxxoxxxxxo ``` The above table represents our sieve of Eratosthenes, where the first row represents each number from 0-29 (the decimal parts of the numbers are left out for readability purposes). The next two rows represent the sieved out numbers based on primes $$2,3,5$$ which our primorial of choice $$p\_3#$$ is comprised of. The final row represents all the eliminated/sieved numbers when taking the whole primorial into account. What we can see is a pattern that emerges. Specifically the pattern `xoxxxxxoxxxoxoxxxoxoxxxoxxxxxo` of eliminated numbers. Note that `x` is an eliminated number and `o` is a number that has not been sieved out. We can see that this pattern repeats forever. Thus, we can use this observation to our advantage. Since the pattern of primes relative to the primorial repeats, then the pattern of prime k-tuples also repeats (relative to that primorial). The idea is to calculate the offsets based on the above table, up to $$p\_m#$$. In our example we notice that we have a few numbers left until $$p\_3#$$ after we sieve out primes $$2,3$$ and $$5$$. For those remaining numbers we also have to check if $$i+o\_1, i+o\_2, ..., i+o\_k$$ are also divided by $$2,3$$ and $$5$$. We will choose a simpler pattern for this example but this applied to any pattern. Lets use the triplet $$n + (0, 2, 6)$$ as an example. For those left numbers $$i$$ we have to check if $$i+2$$ and $$i+6$$ are divisible by our primes $$2,3,5$$. If so we have to also eliminate them. This makes sense as we don't want to consider prime bases that don't match our desired pattern! After doing so the left numbers are only $$(11, 17)$$. These are the offsets for the primorial $$p\_3#$$. Now we only need to take into account numbers of the form $$p\_3# \cdot f + 11$$ and $$p\_3# \cdot f + 17$$ starting from some $$f$$ such that our numbers are bigger than our target $$T$$. This is done by calculating $$T^{'} = T + p\_m # - (T ; mod ; p\_m#)$$ and checking for numbers of the form $$T\_f = T^{'} + p\_m # \cdot f + o$$ or more concretely for our example $$T\_f = T^{'} + p\_3 # \cdot f + 11$$ and $$T\_f = T^{'} + p\_3 # \cdot f + 17$$ In practice for larger primorials these offsets will be a lot more so reality we do something a bit different. The idea is to quickly take the first offset $$o\_1$$ which in the previous example was $$11$$ but without doing the whole sieve for primorial $$p\_m #$$. This is achieved by using the base prime $$n$$ of a constellation (with the same pattern as our desired pattern) that is larger than the $$m^{th}$$ prime $$p\_m$$ of the primorial of choice $$p\_m #$$. In practice the offsets are pre-calculated and given as input to the algorithm. ### More optimized: Eliminating factors The last step in our algorithm is to eliminate possible $$f$$'s that produce candidates that are definitely composite. If we can also eliminate them based on a pattern, and also efficiently compute this pattern, then we could save a lot of time in the primality testing phase as we skip on these $$f$$'s. In the python implementation a dictionary is used, which in a more real world setting would be replaced by a boolean array. Since we are only considering numbers of the form $$p\_3 # \cdot f + o + c\_i$$ if we require that $$o \leq f \leq p$$ where p is a prime in the prime table we can solve for $$f\_p = ((p - ((T^{'} + o + c\_i) ; mod ; p)) \cdot {p\_m#}^{-1}) ; mod ; p$$ where $$c\_i$$ is the offset in our desired constellation pattern, $${p\_m#}^{-1}$$ is the multiplicative inverse of our selected primorial modulo $$p$$ and $$o$$ is the pre-calculated offset that is given. Using that we can calculate an $$f$$ such that $$T\_{f\_p}$$ is divisible by $$p$$ thus we can eliminate $$f\_p$$ and all of its multiples. The following Python code is the implementation: ```python def get_eliminated_factors(primes, inverses, offset, p_m, m, t, constellation_pattern): T2 = t + p_m - (t % p_m) # T' eliminated_factors = {} k_max = 100 i = 0 for p in primes: p_m_inverse = inverses[i] if p_m_inverse != 0: for c_i in constellation_pattern: f_p = ((p - ((T2 + o + c_i) % p))*p_m_inverse) % p eliminated_factors[f_p] = True for k in range(0, k_max): f_p+=p eliminated_factors[f_p] = True i+=1 return eliminated_factors ``` As we can see, we consider all primes in our prime table that are larger than $$p\_m$$ the largest number in our primorial or where the multiplicative inverse of $$p\_m# ; mod ; p$$ exists. We then calculate $$f\_p$$ for each $$c\_i$$ in our constellation pattern and eliminate $$f\_p$$ and its multiples up to an arbitrarily selected $$k\_{max}$$. Finally, when we do wheel factorization we skip over the $$f$$'s that provably do not produce prime candidates. The following code demonstrates that: ```python for f in range(f_start, f_start+f_max): if f not in eliminated_factors: primality_tests +=1 candidate = p_m*f + o if fermat_is_valid_pow(candidate, constellation_pattern): print(candidate) primes+=1 ``` ### Further Reading Despite the fact that we have implemented some optimizations in our code there is a lot to be desired. Firstly we could be using a more suited programming language like C, C++ or Rust which are all generally faster than Python. Since our goal is to be as efficient as possible, this is a must for the final implementation. Another thing to consider is to use an efficiently implemented boolean array to sieve out the candidates. This is usually done by using an unsigned integer array where each bit of the integer represents a number. Then use bit-maks to be able to set a specific bit with 1 or 0 depending on whether that number is a prime or not. This will be a lot quicker and memory efficient as well. Finally we can compute multiple offsets (as many as the CPU cores in our disposal) and parallelize the search on each core, giving each core its own offset. Finally, you can find the complete code [here](https://pastebin.com/1i9HnyRF). --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://makischristou.gitbook.io/primes/efficient-prime-k-tuple-search.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.