Skip to content

Instantly share code, notes, and snippets.

@ericnormand
Created November 15, 2021 15:30
Show Gist options
  • Save ericnormand/31f4eccd020567c3ecdd8cadd918f725 to your computer and use it in GitHub Desktop.
Save ericnormand/31f4eccd020567c3ecdd8cadd918f725 to your computer and use it in GitHub Desktop.
451 PurelyFunctional.tv Newsletter

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/

@KingCode
Copy link

KingCode commented Nov 17, 2021

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment