## Project Euler Problem 51 Solution

## Question

By replacing the 1st digit of *3, it turns out that six of the nine possible values: 13, 23, 43, 53, 73, and 83, are all prime.

By replacing the 3rd and 4th digits of 56**3 with the same digit, this 5-digit number is the first example having seven primes among the ten generated numbers, yielding the family: 56003, 56113, 56333, 56443, 56663, 56773, and 56993. Consequently 56003, being the first member of this family, is the smallest prime with this property.

Find the smallest prime which, by replacing part of the number (not necessarily adjacent digits) with the same digit, is part of an eight prime value family.

## Haskell

```
import Data.List (sort, group)
primes :: [Int]
primes = 2 : sieve primes [3,5..] where
sieve (p:ps) xs = h ++ sieve ps [x | x <- t, rem x p /= 0]
where (h, t) = span (< p*p) xs
isPrime :: Int -> Bool
isPrime n | n < 1 = False
| otherwise = not $ or [n `rem` x == 0 | x <- [2..floor $ sqrt $ fromIntegral n]]
groupFrequency :: Ord a => [a] -> [(Int, a)]
groupFrequency xs = reverse $ sort [(length cs, head cs) | cs <- group $ sort xs]
candidates :: [(Char, String)]
candidates = [(snd $ head g, s) | p <- primes,
let s = show p,
let g = groupFrequency s,
any ((== 3) . fst) g]
replace :: Eq a => a -> a -> [a] -> [a]
replace old new = foldr (\x acc -> (if x == old then new else x) : acc) []
expand :: (Char, String) -> [String]
expand (d, ds) = [new | n <- ['1'..'9'], let new = replace d n ds, isPrime (read new)]
main :: IO ()
main = putStrLn $ snd . head $ filter ((== 8) . length . expand) candidates
```

```
$ ghc -O2 -o prime-replacement prime-replacement.hs
$ time ./prime-replacement
real 0m0.131s
user 0m0.124s
sys 0m0.000s
```

## Python

```
#!/usr/bin/env python
from itertools import *
def primes(n):
if n==2: return [2]
elif n<2: return []
s=list(range(3,n+1,2))
mroot = n ** 0.5
half=(n+1)//2-1
i=0
m=3
while m <= mroot:
if s[i]:
j=(m*m-3)//2
s[j]=0
while j<half:
s[j]=0
j+=m
i=i+1
m=2*i+3
return [2]+[x for x in s if x]
from bisect import bisect_left
# sqrt(1000000000) = 31622
__primes = primes(31622)
def is_prime(n):
# if prime is already in the list, just pick it
if n <= 31622:
i = bisect_left(__primes, n)
return i != len(__primes) and __primes[i] == n
# Divide by each known prime
limit = int(n ** .5)
for p in __primes:
if p > limit: return True
if n % p == 0: return False
# fall back on trial division if n > 1 billion
for f in range(31627, limit, 6): # 31627 is the next prime
if n % f == 0 or n % (f + 4) == 0:
return False
return True
def number_families(num):
digits = [d for d in str(num)]
products = list(product((True, False), repeat=len(digits)))[1:-1]
for p in products:
pattern = ''
for i, x in enumerate(p):
if x:
pattern += digits[i]
else:
pattern += 'X'
yield [int(pattern.replace('X', str(n))) for n in range(10)]
def main():
for prime in primes(1000000):
for number_family in number_families(prime):
prime_family = [n for n in number_family if is_prime(n) and len(str(n)) == len(str(prime))]
if len(prime_family) == 8 and prime in prime_family:
print(prime)
return
if __name__ == "__main__":
main()
```

```
$ time python3 changing-primes.py
real 0m13.523s
user 0m13.416s
sys 0m0.020s
```

Questions? Comments? Send me an email: [email protected]