Project Euler Problem 51 Solution
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.
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
#!/usr/bin/env python from itertools import * def primes(n): if n==2: return  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 +[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: firstname.lastname@example.org