finding the maximum product of 2 primes below a given number

Question

Given a number N, how do we find maximum P*Q < N, such that P and Q are prime numbers?

My (brute force) attempt:

• find a list {P, N/P} for all primes P < √N
• find a list of primes Q, such that Q is the largest prime just below N/P in the list above
• Determine the maximum product P*Q from above

While this brute force approach will work, is there a formal (more sensible) solution to this question?

Example: N=27

√N = 5.196

Candidate primes: 2,3,5 --> [{2,13.5},{3,9},{5,5.4}] ->[{2,13},{3,7},{5,5}]

Solution: Max([2*13, 3*7, 5*5]) = 2*13 = 26

Hence, the brute force solution works.

Taking this one step further, we see that Q_max <= N/2 and if indeed we agree that P < Q, then we have Q >= √N.

We can refine our solution set to only those values {P, N\2}, where N\2 >= √N.

I have opted for integer division "\", since we are only interested in integer values, and integer division is indeed much faster than regular division "/"

The problem reduces to:

Example: N=27

√N = 5.196

Candidate P: 2,3 --> [{2,13},{3,9}] -->[{2,13},{3,7}]

(we drop {5,5} since N\P < √N i.e. 5 < 5.196)

Solution set: max([2*13, 3*7]) = 2*13 = 26

It might look trivial, but it just eliminated 1/3 of the possible solution set.

Are there other clever procedures we can add to reduce the set further?

Answer 1

Another brute force attempt:

1. Find all prime numbers p <= N/2
2. Iterate over the array from the smallest p as long as p < √N and multiply it with the largest q < N/p, retaining the largest product.

If enough memory is available (N/2 bit), one could make a bitarray of that size. Initialize it with all TRUE but the first position. Iterating over the bit array, calculate the multiples of the position you are at and set all multiples to false. If the next position is false already, you do not need to recalculate all its multiples, they are already set to false.

Finding all primes therefore is < O(N^2).

a[1] := false;
m := n \ 2; // sizeof(a)
for i := 2 to m do
   a[i] := true;
for i := 2 to m do
   if a[i] then
     for j := 2*a[i] to m step a[i] do
       a[i*j] := false;

Step 2) is < O(n^2) as well:

result := 0;
for i := 2 to √N do
  if not a[i] then continue; // next i;
  for j := (n \ i) downto i do
     if not a[j] then continue; // next j
     if a[j] * a[i] < N
        result :=  max(result, a[j] * a[i]);
        break; // next i;
  if result = N then break; // you are finished

This can be optimized further, I guess. You can keep (i,j) to know the two prime numbers.

Answer 2

This is similar to what @RalphMRickenback describes, but with tighter complexity bounds.

The prime finding algorithm he describes is the sieve of Erathostenes, which needs space O(n), but has time complexity O(n log log n), you may want to see the discussion on Wikipedia if you want to be more careful about this.

After finding a list of primes smaller than n // 2, you can scan it a single time, i.e. with O(n) complexity, by having a pointer start at the beginning and another at the end. If the product of those two primes is larger than your value, reduce the high pointer. If the product is smaller, compare it to a stored maximum product, and increase the low pointer.

EDIT As mentioned in the comments, the time complexity of the scan of the primes is better than linear on n, since it is only over the primes less than n, so O(n / log n).

Rather than pseudo-code, here's full implementation in Python:

def prime_sieve(n):
    sieve = [False, False] + [True] * (n - 1)

    for num, is_prime in enumerate(sieve):
        if num * num > n:
            break
        if not is_prime:
            continue
        for not_a_prime in range(num * num, n + 1, num):
            sieve[not_a_prime] = False

    return [num for num, is_prime in enumerate(sieve) if is_prime]

def max_prime_product(n):
    primes = prime_sieve(n // 2)

    lo, hi = 0, len(primes) - 1
    max_prod = 0
    max_pair = None

    while lo <= hi:
        prod = primes[lo] * primes[hi]
        if prod <= n:
            if prod > max_prod:
                max_prod = prod
                max_pair = (primes[lo], primes[hi])
            lo += 1
        else:
            hi -= 1
    return max_prod, max_pair

With your example this produces:

>>> max_prime_product(27)
(26, (2, 13))