451 PurelyFunctional.tv Newsletter

00 Ulam sequence.md

Ulam sequence

The Ulam sequence is an interesting mathematical sequence of integers. It starts with [1 2]. At each step, the next element is:

• not already in the sequence
• a sum of two previous elements
• the number must be produced by only one sum
• the smallest in case there are multiple candidates

Let's walk through an example.

As we said before, the first two elements are 1 and 2. Those sum to 3, and there's only one way to get it, so 3 is the next element.

[1 2 3]

Now, there are 3 possible sums: 1+2=3, 1+3=4, and 2+3=5. 3 already exists, so it's out. 4 and 5 both only have one way to produce them, but 4 is smaller, so it's the next element.

[1 2 3 4]

Now: 1+2=3 (done), 1+3=4 (done), 1+4=5, 2+3=5, 2+4=6, 3+4=7. 5 is produced two ways, so it's out. 6<7, so the next element is 6.

[1 2 3 4 6]

Your task is to create a lazy Ulam sequence.

Example

(take 2 (ulam)) ;=> (1 2)
(take 5 (ulam)) ;=> (1 2 3 4 6)
(take 17 (ulam)) ;=> (1 2 3 4 6 8 11 13 16 18 26 28 36 38 47 48 53)

Thanks to this site for the problem idea, where it is rated Expert in Java. The problem has been modified.

Please submit your solutions as comments on this gist.

To subscribe: https://purelyfunctional.tv/newsletter/

(defn next-ulam-set [current-set]
  (->> (+ a b)
       (for [a current-set
             b current-set :when (< a b)])
       (remove current-set)
       frequencies
       (filter #(= 1 (second %)))
       (map first)
       (apply min)
       (conj current-set)))

(defn ulam []
  (->> (sorted-set 1 2)
       (iterate next-ulam-set)
       (map (comp first rseq))
       (cons 1)))

(defn- tails [xs]
  (->> xs (iterate rest) (take-while seq)))

(defn- keep-unique [xs]
  (->> (frequencies xs)
       (keep (fn [[x n]] (when (= 1 n) x)))))

(defn ulam-next [seen]
  (->> (for [[a & bs] (tails seen)
             b bs
             :let [c (+ a b)]
             :when (not (seen c))]
         c)
       keep-unique
       (apply min)
       (conj seen)))

(defn ulam []
  (let [init (sorted-set 1 2)]
    (concat init (->> init (iterate ulam-next) (drop 1) (map last)))))

(require '[clojure.test :refer [deftest are is]])
(deftest ulam-test
  (are [range expected]
       (is (= expected (->> (ulam) (take (second range)) (drop (first range)))))
    [0 2]     [1 2]
    [0 5]     [1 2 3 4 6]
    [0 17]    [1 2 3 4 6 8 11 13 16 18 26 28 36 38 47 48 53]
    [205 207] [1856 1858]
    [451 453] [5018 5020]))
;; https://oeis.org/A002858/b002858.txt

@jonasseglare that some interesting trick with threading expression into for loop. Don't know if I love it or hate it 🤣

@nbardiuk Thanks! That code is meant to troll you 😆

There are some optimized algorithms for this but they seem like a lot of work for this gist.

Not as pretty as the above solutions, I tried to exploit the ordering for efficient updating of candidates to add to the Ulam sequence:

(defn merge-sorted [xs ys]
  (lazy-seq
    (if-let [fx (first xs)]
      (if-let [fy (first ys)]
        (if (< fx fy)
          (cons fx (merge-sorted (rest xs) ys))
          (cons fy (merge-sorted xs (rest ys))))
        xs)
      ys)))

(defn trim-leading-dups [xs]
  (if-let [x (first xs)]
    (if (= x (second xs))
      (recur (drop-while #{x} xs))
      xs)
    ()))

(defn ulam* [ulam-vec sums]
  (lazy-seq
    (let [ulam-next (first sums)
          new-sums (map + ulam-vec (repeat ulam-next))
          sums' (trim-leading-dups (merge-sorted (rest sums) new-sums))]
      (cons ulam-next (ulam* (conj ulam-vec ulam-next) sums')))))

(defn ulam []
  (cons 1 (ulam* [1] [2])))

(defn sums-without-duplicates [coll]
  (loop [coll            coll
         first-iteration coll
         sums            []]
    (cond
      (= (count coll) 2)
      (as-> sums v
        (conj v (+ (first coll) (second coll)))
        (frequencies v)
        (group-by val v)
        (get v 1)
        (keys v)
        (sort v))
      (= (count first-iteration) 2)
      (recur (rest coll) (rest coll) (conj sums (+ (first first-iteration) (second first-iteration))))
      :else (recur coll (butlast first-iteration) (conj sums (+ (first first-iteration) (last first-iteration)))))))

(defn ulam
  ([] (ulam [1]))
  ([coll]
   (lazy-seq
     (let [number (count coll)]
       (if (< number 2)
         (cons (last coll) (ulam (conj coll (inc number))))
         (cons (last coll) (ulam (conj coll (some #(when-not (.contains coll %) %) (sums-without-duplicates coll))))))))))

Edit: refactored from clojure.math.combinatorics/combinations to for loop ala @jonasseglare. The code is shorter than the name, lol. Btw, I know mine is strikingly similar - but the for loop is the only thing I borrowed. The rest was just the end point of several code golf refactorings.

(defn nexul [us]
  (->> (for [a us b us :when (< a b)] (+ a b))
       frequencies
       (filter #(= 1 (val %)))
       (map first)
       (remove (set us))
       (apply min)
       (conj us)))

(defn ulam []
  (->> (iterate nexul [1 2])
       (map last)
       (cons 1)))

This was a fun exercise - I made many edits, ending with two side-by-side variations of the basic algorithm,
which looks at those existing Ulam numbers below half of the current candidate value, mentioned in this paper.

One variation is a hand-rolled lazy-sequence using a recursive function taking the state, consing the computed value from it,
and feeding it to the next call). The other one uses iterate and a transducing (map first).

In both cases the approach is to pick a candidate c (c0 = U(n) + 1) , then go through each number U(k) in the
Ulam sequence search space in descending order, and test that diff(c,k) = c - U(k) is also in the sequence:

• if it is and no other test with a U(j), j < k yields a member diff(c,j) then c is a comfirmed Ulam number;
• else, continue with c + 1 and repeat.

To quickly extract the bottom portion of the existing sequence as the search space argument to next-ulam
I tried two implementatons of a fast binary-search function (big-O(logN)): one hand-rolled and the other using
java.util.Collections.binarySearch (bin-search, bin-search2 resp.).

Surprisingly, the bench tests show that the hand-rolled versions (ulam and bin-search) perform slightly better,
especially bin-search which I expected the java version to be much faster.

(defn next-ulam 
  "Yields the next-ulam number given a search space (a sorted 
   sub/vector of exising ulam numbers), a set of all previously
   known ulam numbers, a binary predicate validating distinct sum 
   operands, and the minimal starting candidate."
  [ulies uly? start]

   ;; Starting with minimal candidate kand = U-n + 1,
   ;; pluck each sum partner p in the search space of (some) 
   ;; previous Ulam numbers starting with the highest.
   ;; If kand - p = p', keep [p,p'] if p' is a member.
   ;; If no other pair is found, kand is the next Ulam;
   ;; otherwise incrment kand and repeat. 
   
   ;; Traversal order enables O(N) decision for the call,
   ;; and also helps detect a previously seen pair
   ;; with its operands reversed (depending on new-sum? pred.

  (loop [kand start 
         sums-count 0 
         search ulies]
    (let [[sum-pal diff] (when-let [top (peek search)]
                           [top, (- kand top)])]
      (cond
        (< 1 sums-count)
        (recur (inc kand) 0 ulies)

        (empty? search) 
        (if (pos? sums-count)
          kand
          (recur (inc kand) 0 ulies))

        (and (uly? diff) (< sum-pal diff))
        (recur kand (inc sums-count) (pop search))

        :else
        (recur kand sums-count (pop search))))))

(defn bin-search
  "Store is a sorted sequential. Yields the index of the value v 
   closest to n s.t. v < n.

   If no element is found satisfying the constraint (when
   n <= first element) nil is returned."
  [store n]
  (let [k (count store)]
    (cond 
      (<= n (nth store 0)) nil, ;; boundary cases
      (< (peek store) n) (dec k),
      (< k 3) 1
      :else ;; main
      (loop [lo 0 hi k]
        (let [mid (-> hi (+ lo) (quot 2))
              v (nth store mid)
              v' (nth store (inc mid))]
          (cond 
            (and (< v n) (<= n v'))
            mid ;;found

            (< v n) ;; too low
            (recur mid hi)

            :else ;; too high
            (recur lo mid)))))))

(defn bin-search2 [store n]
  (when (< 1 n)
    (let [k (java.util.Collections/binarySearch store n)]
      (if (pos? k)
        (dec k)
        (- k)))))

(def ^:dynamic *bin-search* bin-search)

(defn bottom-ulies [ulies limit-twice+]
  (let [cutoff (-> ulies 
                   (*bin-search* (inc (quot limit-twice+ 2)))
                   (or (-> ulies count dec)))]
    (subvec ulies 0 (inc cutoff))))

(defn ulam1 
  ([] (ulam1 1 2))
  ([u0 u1]
   (let [ufn
         (fn u1 [xs xset]
           (let [min-kand (inc (peek xs)) 
                 ret (next-ulam (bottom-ulies xs min-kand),
                                xset, min-kand)]
             (cons ret (lazy-seq (u1 (conj xs ret)
                                     (conj xset ret))))))]
     (cons u0 (cons u1 (ufn [u0 u1] #{u0 u1}))))))

(defn ulam2 
  ([] (ulam2 1 2))
  ([u0 u1]
   (let [kont
         (->> [u1 [u0 u1] #{u0 u1}] 
              (iterate 
               (fn [[prev useq uset]]
                 (let [min-kand (inc prev) 
                       ret (-> useq (bottom-ulies min-kand)
                               (next-ulam uset, min-kand))]
                   [ret (conj useq ret) (conj uset ret)])))
              rest
              (sequence (map first)))]
     (cons u0 (cons u1 kont)))))

;; TESTING CODE
(defn ->useq [ufn-or-useq seeds]
         (if (fn? ufn-or-useq)
           (apply ufn-or-useq (take 2 seeds))
           ufn-or-useq))

(defn useq-samples 
  ([ending-with-nth ufn-or-useq & seeds]
   (->> (range 1 (inc ending-with-nth))
        (map #(take % (->useq ufn-or-useq seeds))))))

(defn urange-samples
  ([start-nth n ufn-or-useq & seeds]
   (->> (->useq ufn-or-useq seeds)
        (drop (dec start-nth))
        (take n))))

;;Accuracy? 
(useq-samples 20 ulam1)
;;=> ((1)
 ;;   (1 2)
 ;;   (1 2 3)
 ;;   ...
 ;;   (1 2 3 4 6 8 11 13 16 18 26 28 36 38 47 48 53 57 62 69))

(= (useq-samples 100 ulam1)
   (useq-samples 100 ulam2))
;;=> true

;; Test further against
;; 5 entries starting with n1670 = 21083 
;; in Jud McCranie's table referenced at https://oeis.org/A002858, see:
;; https://oeis.org/A002858/b002858.txt
(urange-samples 1670 5 ulam1)
;;=> (21083 21105 21108 21110 21130)

(->> [ulam1 ulam2] (map #(urange-samples 1670 5 %)) =)
;;=> true

;; Performance?
(require '[criterium.core :refer [bench quick-bench]])


(quick-bench (->> (ulam1) (drop 1000) (take 10) doall))
;;=> Evaluation count : 6 in 6 samples of 1 calls.
             ;; Execution time mean : 397.219837 ms
    ;; Execution time std-deviation : 3.857484 ms

;; (bench (->> (ulam1) (drop 1000) (take 10) doall))
;;=> Evaluation count : 180 in 60 samples of 3 calls.
             ;; Execution time mean : 393.787551 ms
    ;; Execution time std-deviation : 2.792384 ms

(binding [*bin-search* bin-search2] 
  (quick-bench (->> (ulam1) (drop 1000) (take 10) doall)))
;;=> Evaluation count : 6 in 6 samples of 1 calls.
             ;; Execution time mean : 403.468283 ms
    ;; Execution time std-deviation : 2.252653 ms

(quick-bench (->> (ulam2) (drop 1000) (take 10) doall))
;;=> Execution time mean : 406.314845 ms
    ;; Execution time std-deviation : 6.954487 ms
   ;; Execution time lower quantile : 401.708926 ms ( 2.5%)

;; (bench (->> (ulam2) (drop 1000) (take 10) doall))

;;=> Evaluation count : 180 in 60 samples of 3 calls.
             ;; Execution time mean : 410.378755 ms
    ;; Execution time std-deviation : 11.957382 ms

(binding [*bin-search* bin-search2] 
  (quick-bench (->> (ulam2) (drop 1000) (take 10) doall)))
;;=> Evaluation count : 6 in 6 samples of 1 calls.
             ;; Execution time mean : 411.399288 ms
    ;; Execution time std-deviation : 2.756681 ms

(def ulams-1 (ulam1))
(def ulams-2 (ulam2))
;;Force realized access
(->> ulams-1 (drop 1000) (take 10) doall)
(->> ulams-2 (drop 1000) (take 10) doall)

(quick-bench (->> ulams-1 (drop 1000) (take 10) doall))
;;=> Evaluation count : 16086 in 6 samples of 2681 calls.
             ;; Execution time mean : 37.383059 µs
    ;; Execution time std-deviation : 273.345596 ns

(quick-bench (->> ulams-2 (drop 1000) (take 10) doall))
;;=> Evaluation count : 15576 in 6 samples of 2596 calls.
             ;; Execution time mean : 39.063882 µs
    ;; Execution time std-deviation : 342.094132 ns