Collatz Primes

It’s has been a busy week in between getting myself out of jet lag, teaching, preparing seminar lectures, and learning about making music on computers and having fun.

As I was winding down and preparing to sleep, I started reading Programming Praxis. Ah! the Collatz Sequence. Can’t resist.

Description of the problem

Given \(n\), write a function which returns the smallest initial value of a Collatz sequence which contains at least \(n\) primes. Recall that the Collatz sequence is defined recursively as

\[ c_{n+1} = \begin{cases} 1 & \text{ if } c_n = 1 \\ c_n/2 & \text{ if $c_n$ is even } \\ 3 c_n + 1 & \text{ otherwise } \end{cases} \]

Solution

First I need a general purpose filtering function:

(defun filter (pred xs &optional carry)
   (cond
      ((null xs) (nreverse carry))
      ((funcall pred (car xs)) (filter pred (cdr xs) (cons (car xs) carry)))
      (t (filter pred (cdr xs) carry))))
FILTER

Next, a predicate to check primeness. I will memoize this:

(let ((table-a (make-hash-table)))
  (defun primep (n)
    (multiple-value-bind
          (val foundp) (gethash n table-a)
      (if foundp
          val
          (setf (gethash n table-a)
                (or (= 2 n)
                    (and (oddp n)
                         (let ((m (floor (sqrt n))))
                           (do ((i 3 (+ 2 i)))
                               ((or (zerop (mod n i)) (> i m))
                                (> i m)))))))))))

And finally the Collatz Sequence generator. I will modify the sequence as it is clear that we don’t care about the even part of the sequence.

(defun collatz (n &optional carry)
   (cond
      ((= n 1) carry)
      ((evenp n) (collatz (/ n 2) carry))
      (t (collatz (1+ (* 3 n)) (cons n carry)))))
COLLATZ

Now, let us search the case for which we get 65 primes:

(let ((carry nil))
   (do ((n 2 (1+ n)))
       ((>= (length carry) 65) (list n (length carry) carry))
      (setf carry (filter #'primep (collatz n)))))
(1805312 66
(5 53 23 61 433 577 1367 911 1619 719 479 283 251 167 593 263 233 103 137 107
 71 47 1777 1579 2213 983 2621 1747 6211 117773 93053 248141 348553 275399
 122399 580259 386839 2063141 2445203 2173513 2898017 1932011 381631 508841
 150767 100511 59561 70589 47059 1784753 352543 2228417 660271 880361 2473217
 1648811 1737017 1158011 14639539 52051693 46268171 30845447 9139391 6092927
 4061951 2707967))

Addendum

@simendiferlerin on twitter pointed out that my primep returns t for even numbers, and it did. Not by intention, of course. But it worked in this case as I never send an even number to its path :) So, I changed the code above accordingly.

While at it, one can change the last part of the code to make it more like a lisp code should be, recursive :)

(defun collatz-primes (n &optional (m 2))
  (let ((res (filter #'primep (collatz m))))
    (if (>= (length res) n)
        (list m (length res) res)
        (collatz-primes n (1+ m)))))

And now you can call collatz-primes with the number of prime numbers you want and it returns the initial value for the sequence, the number of primes in the Collatz sequence, and the list of primes themselves.

And finally: here is the same code (without the memoization) translated to clojure in case common lisp is non-fashionable for you :-p

(defn prime? [n]
  (cond
    (= n 2) true
    (even? n) false
    :else (let [m (Math/floor (Math/sqrt n))]
            (loop [i 3]
              (if (or (= 0 (mod n i))
                      (> i m))
                (> i m)
                (recur (inc i)))))))

(defn collatz [n]
  (loop [m n
         carry '()]
    (cond
      (= m 1) carry
      (even? m) (recur (/ m 2) carry)
      :else (recur (+ 1 (* 3 m)) (conj carry m)))))

(defn collatz-prime [n]
  (loop [m 2]
    (let [res (filter prime? (collatz m))]
      (if (>= (count res) n)
        {:arg m :len (count res) :res res}
        (recur (inc m))))))