Project Euler Problem 51

Ivar Thorson bio photo By Ivar Thorson

Problem 51 asks us to find a tricky class of prime numbers generated by replacing certain digits of a pattern with the same digit.

My first approach to solving this problem was to try a simple algorithm:

  1. Incrementally search through primes in increasing order.
  2. For each prime, generate a list of possible digits to mask off and replace with the same number
  3. Map through the masked-off numbers, testing each for primality
  4. If we find enough (8) prime numbers for a particular mask, we are done.

Using clojure.contrib.combinatorics sounds tempting here to generate permutations and combinations, but it’s just as easy to create our own binary-masking template, which I called n-digit-masks.

(use '[clojure.contrib.lazy-seqs :only (primes)]
     '[clojure.contrib.math :only (expt)])

(defn as-int [coll] (Integer/parseInt (apply str coll)))

(defn prime? [n]
  (and (< 1 n)
       (not-any? #(zero? (rem n %)) (take-while #(<= (* % %) n) primes)))) 

(defn n-digit-masks
  "Returns all n-digit binary numbers as strings."
  (map #(rest (Integer/toBinaryString (bit-or (bit-shift-left 1 n) %)))
       (range (expt 2 n))))

(defn apply-mask
  "Replaces digits of p with value n according to nonzero locations in mask m."
  [m p n]
  (map #(if (= \0 %1) %2 n) m p))

(defn euler-51 [family]
   (for [p (map str primes)
         d [(count p)]
         m (rest (butlast (n-digit-masks d)))
         f [(filter #(prime? (as-int %))
                    (map #(apply-mask m p %) "0123456789"))]
         g [(filter #(not (= \0 (first %))) f)]
         :when (<= family (count g))]
     (map #(apply str %) g))))

(time (euler-51 8)) ;; "Elapsed time: 29890.717326 msecs"

Whoa, that’s kind of slow! Let’s see if instead of the following a generative approach like above a reduction-oriented approach might be faster. That is, if we start with the list of prime numbers and just apply filters we can avoid primality tests completely…

(defn n-digit-primes [d]
  (->> primes
       (drop-while #(< % (expt 10 (dec d))))
       (take-while #(< % (expt 10 d)))
       (map str)))

(defn all-unmasked-same? [m n]
  (let [z (map #(if (= \0 %1) %2 \x) m n)]
    (= 1 (count (distinct (filter #(not (= \x %)) z))))))

(defn mask-off [m a]
  (apply str (map #(if (= \0 %1) \x %2) m a)))

(defn euler-51-alternate [family]
   (for [d (iterate inc 2)
         ps [(n-digit-primes d)]
         m (reverse (rest (butlast (n-digit-masks d))))
         f [(filter #(all-unmasked-same? m %) ps)]
         g [(frequencies (map #(mask-off m %) f))]
         k (sort (keys g))
         :when (<= family (g k))]

(time (euler-51-alternate 8)) ;; "Elapsed time: 21633.206711 msecs"

That’s mildly better, but certainly not a showstopper, and I think it may be subtly incorrect in the order that masks are applied, even though it gives the correct answer in this case.

Hmm…I guess I’ll just leave this one here for now. Suggestions on a better algorithm are welcome!