How to Enumerate Pandigital Prime Sets

Question

Project Euler problem 118 reads, "Using all of the digits 1 through 9 and concatenating them freely to form decimal integers, different sets can be formed. Interestingly with the set {2,5,47,89,631} all of the elements belonging to it are prime. How many distinct sets containing each of the digits one through nine exactly once contain only prime elements."

This is the only problem I've encountered so far that doesn't provide an answer to an easier version of the question. As a result, it is difficult for me to check the algorithm I've made. Do you see what could be the problem?

I coded up a pretty straightforward dynamic programming approach. For that, I will answer the question, "How many distinct sets containing each of the digits in S exactly once contain only prime elements", for different values of S, and then combine their results to answer the original question.

First, I generate a list of all the relevant primes. This is every prime with 9 or fewer digits, that has no duplicate digits, and that doesn't have a 0 as one of its digits. I then use this list to create a set count dictionary. The keys to this dictionary are frozensets that represent which digits each prime has, and the values are integers that count the number of relevant primes that reduce to that set:

from itertools import chain, combinations
from functools import cache
from collections import defaultdict
from pyprimesieve import primes


primes = [
    p for p in primes(10**9) if len(set(str(p))) == len(str(p)) and "0" not in str(p)
]
set_counts = defaultdict(int)
for p in primes:
    set_counts[frozenset(map(int, str(p)))] += 1

set_counts will then serve as the base case for our recursion.

Next, we need a way to break down our problem into relevant subproblems. So I wrote a function which will generate all of the ways to split a set into two disjoint nonempty subsets:

@cache
def set_decomps(s):
    """
    Generates all possible ways to partition a set s
        into two non-empty subsets
    """
    l = list(
        map(
            lambda x: (frozenset(x), s - frozenset(x)),
            chain.from_iterable(combinations(s, r) for r in range(1, len(s))),
        )
    )
    return l[: len(l) // 2]

Finally, we should be able to just use a simple memoized approach to use these together to solve the problem, and all of its subvariants.

@cache
def dp(s):
    """
    Returns the number of ways to create a set containing prime numbers
    Such that the set of all the digits in the primes is equal to s
    With no duplicates or extra digits
    """
    base = set_counts[s]
    if len(s) > 1:
        for a, b in set_decomps(s):
            base += dp(a) * dp(b)
    return base
print(dp(frozenset(range(1,10)))) # prints 114566, which is wrong

Obviously there is a flaw in my approach.

Answer 1

The flaw is this. Consider the set {2,5,47,89,631}. Different initial splits can lead to it. One is {1,2,3,6}, {4,5,7,8,9} and another is {1,3,6,8,9} {2,4,5,7}. There are many more. And therefore you are overcounting.

To leave you the fun of the problem, I won't tell you how to fix this overcounting. I'll just tell you that if you are counting multiple ways to get to a particular set, then your number will be too high.

Answer 2

This seemed interesting. To avoid overcounting, I took a recursive approach that relies on a structure that orders the digits of each candidate prime and groups them by the lowest digit, as well as records the count of each equal digit set. For example, 13 and 31 would be listed under the group for lowest-digit 1 as one set, {1, 3} (represented as a bitset by an int), with a count of 2.

Then the search iterates over each digit set, starting with the first group, and tries to build valid, memoised, paths by summing the multiples of the current count with the counts for the growing digit sets, where the next choice is from the group of the next lowest missing digit.

The bottleneck, I found out, was in generating just the primes that are needed rather then checking all primes up to 987,654,321.

Output:

$ python3 prime_sets.py
Creating primes dictionary took 3 seconds.
The answer is 44680.
Recursive search took 0 seconds.

Python code:

# Project Euler 118
 
import collections
import itertools
import random
import time

 
# https://www.codeproject.com/Articles/691200/Primality-test-algorithms-Prime-test-The-fastest-w
def MillerRabinPrimalityTest(number):
    if number == 2:
        return True
    elif number == 1 or number % 2 == 0:
        return False
 
    oddPartOfNumber = number - 1
    timesTwoDividNumber = 0
 
    while oddPartOfNumber % 2 == 0:
        oddPartOfNumber = oddPartOfNumber // 2
        timesTwoDividNumber = timesTwoDividNumber + 1 
 
    for time in range(3):
        while True:
            randomNumber = random.randint(2, number)-1
            if randomNumber != 0 and randomNumber != 1:
                break
 
        randomNumberWithPower = pow(randomNumber, oddPartOfNumber, number)
 
        if (randomNumberWithPower != 1) and (randomNumberWithPower != number - 1):
            iterationNumber = 1
 
            while (iterationNumber <= timesTwoDividNumber - 1) and (randomNumberWithPower != number - 1):
                randomNumberWithPower = pow(randomNumberWithPower, 2, number)
                iterationNumber = iterationNumber + 1
            if (randomNumberWithPower != (number - 1)):
                return False
 
    return True
 

t0 = time.time()
primes = { d:collections.defaultdict(int) for d in range(1, 10) }

 
def add_prime(p, ps):
  ds = 0
  lowest = 9
  while p > 0:
    d = p % 10
    ds = ds | (1 << d)
    lowest = min(lowest, d)
    p = p // 10
 
  primes[lowest][ds] += 1


for i in range(1, 10):
  for comb in itertools.combinations(map(str, range(1, 10)), i):
    for perm in itertools.permutations(comb):
      n = int("".join(perm))
      if MillerRabinPrimalityTest(n):
        add_prime(n, primes)


t1 = time.time()
print(f"Creating primes dictionary took {int(t1 - t0)} seconds.")


def get_lowest_missing(ds):
  for i in range(1, 10):
    if (1 << i) & ds == 0:
      return i
  return None
 

def f(primes, ds, count, lowest, dp):
  if lowest is None:
    return count
 
  key = (ds, count, lowest)
  if key in dp:
    return dp[key]
 
  answer = 0
 
  for other_ds in primes[lowest]:
    if ds & other_ds == 0:
      answer = answer + f(
          primes,
          ds | other_ds,
          count * primes[lowest][other_ds],
          get_lowest_missing(ds | other_ds),
          dp)
 
  dp[key] = answer
  return answer

 
print(f"The answer is {f(primes, 0, 1, 1, {})}")


t2 = time.time()
print(f"Recursive search took {int(t2 - t1)} seconds.")