Euler #51 Prime Digit Replacements (Spoiler)

euler51.py

from math import log10, ceil, sqrt
import numpy as np
from itertools import combinations
import operator
import array

def basicNumpySieve(bound):
    isPrime = np.ones((bound), dtype=np.bool)
    for i in range(2, ceil(sqrt(bound)) + 1):
        if isPrime[i]:
            isPrime[i*i:bound:i] = False
    numPrimes = np.count_nonzero(isPrime) - 2
    result = np.empty((numPrimes), dtype=np.int64)
    size = 0
    for i in isPrime[2:bound].nonzero()[0]:
        result[size] = i + 2
        size += 1
    assert size == numPrimes
    return result

def digitsOf(num, targetSize=0):
    if num == 0:
        return [0]
    result = array.array('B')
    resultSize = 0
    while num > 0:
        num, digit = divmod(num, 10)
        result.append(digit)
    remainingDigits = targetSize - len(result)
    if remainingDigits > 0:
        result.extend((0,) * remainingDigits)
    result.reverse()
    return result

def fromDigits(digits):
    result = 0
    for digit in digits:
        result *= 10
        result += digit
    return result

def primeDigitMatrix(amount):
    primes =  basicNumpySieve(10**amount)
    primeDigits = []
    result = np.zeros((10,) * amount, dtype=bool)
    for prime in primes:
        digits = tuple(digitsOf(prime, amount))
        primeDigits.append(digits)
        assert len(digits) == amount
        result[digits] = True
    return result, primes, primeDigits

def isSorted(l):
    return all(e >= l[index] for index, e in enumerate(l[1:]))

def digitReplacementPrimeFamilies(boundDigits, minimumSize=1):
    assert boundDigits > 1, f"Invalid bound digits: {boundDigits}"
    digitMatrix, primes, primeDigits = primeDigitMatrix(boundDigits)
    digitReplacementCombinations = []
    for numReplaced in range(1, boundDigits):
        digitReplacementCombinations.extend(combinations(range(boundDigits), r=numReplaced))
    assert boundDigits != 5 or (2, 3) in digitReplacementCombinations, ""
    primeDigitMap = {digits: prime for prime, digits in zip(primes, primeDigits)}
    for prime, digits in zip(primes, primeDigits):
        if digits[0] == 0:
            continue
        originalDigits = np.array(digits, dtype=np.int8)
        # Try replacing parts of the digits
        for replacementIndexes in digitReplacementCombinations:
            digits = originalDigits.copy()
            primeFamily = []
            for value in range(10):
                digits.put(replacementIndexes, value)
                digitTuple = tuple(digits)
                if digitMatrix[digitTuple] and digitTuple[0] != 0:
                    primeFamily.append(primeDigitMap[digitTuple])
            if len(primeFamily) >= minimumSize:
                assert all(family in primes for family in primeFamily)
                return min(primeFamily), primeFamily

for b, minimumSize in ((2, 5),(5, 7)):
    smallestReplacementFamily = digitReplacementPrimeFamilies(b, minimumSize)
    print(f"Smallest digit replacement prime family less than {10**b}: {smallestReplacementFamily}")
    
smallestEightFamily = digitReplacementPrimeFamilies(6, minimumSize=8)
print(f"Smallest family of eight prime: {smallestEightFamily}")