Robert Smith avatar Robert Smith committed c171f4d

Calkin-Wilf trees, Stern's diatomic sequence, and all that,

Comments (0)

Files changed (1)

+;;;; calkin-wilf.lisp
+;;;; Copyright (c) 2012 Robert Smith
+
+;;;; Compute stuff related to the Calkin-Wilf tree.
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;; UTILITIES ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+(defmacro with-destructured-rational ((n d) rational &body body)
+  (let ((rational-once (gensym "RATIONAL-ONCE-")))
+    `(let ((,rational-once ,rational))
+       (let ((,n (numerator ,rational-once))
+             (,d (denominator ,rational-once)))
+         ,@body))))
+
+
+;;;;;;;;;;;;;;;;;;;;;;;;; CALKIN-WILF TREES ;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+;;; The Calkin-Wilf tree is a binary tree generated by the following
+;;; rule:
+;;; 
+;;; Left[a/b] = a/(a+b)
+;;; Right[a/b] = (a+b)/b
+;;; 
+;;; with the root node being 1/1.
+
+(defconstant +calkin-wilf-root+ 1/1
+  "The root of the Calkin-Wilf tree.")
+
+(defun child-left (a/b &optional (n 1))
+  "Compute the left Calkin-Wilf child of A/B. Iterate N times, i.e.,
+compute teh Nth composition of CHILD-LEFT."
+  (with-destructured-rational (a b) a/b
+    (/ a (+ (* n a) b))))
+
+(defun child-right (a/b &optional (n 1))
+  "Compute the right Calkin-Wilf child of A/B. Iterate N times, i.e.,
+compute the Nth composition of CHILD-RIGHT."
+  (with-destructured-rational (a b) a/b
+    (/ (+ a (* n b)) b)))
+
+(defun sibling (a/b)
+  "Compute the Calkin-Wilf sibling of A/B. Note that 1/1 has no
+siblings."
+  (with-destructured-rational (a b) a/b
+    (cond
+      ((< a b) (/ b (- b a)))
+      ((> a b) (/ (- a b) a)))))
+
+(defun parent (a/b)
+  "Compute the Calkin-Wilf parent of A/B. Note that 1/1 has no
+parents."
+  (with-destructured-rational (a b) a/b
+    (cond
+      ((< a b) (/ a (- b a)))
+      ((> a b) (/ (- a b) b)))))
+
+(defun children (a/b)
+  "Compute the Calkin-Wilf children of the rational A/B."
+  ;; This will recompute the sum A+B within each function. Oh well...
+  (list (child-left a/b) (child-right a/b)))
+
+;;;;;;;;;;;;;;;;;;;;;;;; CALKIN-WILF SEQUENCE ;;;;;;;;;;;;;;;;;;;;;;;;
+
+;;; The Calkin-Wilf sequence CW is the sequence generated by reading
+;;; off the tree in the following way:
+;;;
+;;;                    1
+;;;                   / \ 
+;;;                  2   3
+;;;                 / \ / \ 
+;;;                4  5 6  7
+;;;                ...   ...
+;;; 
+;;; That is, read off each generation in order. The recurrence below
+;;; can be derived from the recurrences of Stern's diatomic sequence,
+;;; treated below.
+
+(defun calkin-wilf-sequence-next (a/b)
+  "Compute the next element in the Calkin-Wilf sequence after A/B."
+  (/ (1+ (- (* 2 (floor a/b))
+            a/b))))
+
+(defun calkin-wilf-sequence-nth (n)
+  "Compute the Nth element of the Calkin-Wilf sequence."
+  (labels ((rec (n CW)
+             (if (zerop n)
+                 CW
+                 (rec (1- n)
+                      (calkin-wilf-sequence-next CW)))))
+    (rec n 1)))
+
+;;;;;;;;;;;;;;;;;;;;; STERN'S DIATOMIC SEQUENCE ;;;;;;;;;;;;;;;;;;;;;;
+
+;;; Stern's diatomic sequence S can be defined as the sequence of
+;;; numbers generated by taking the numerators of the Calkin-Wilf
+;;; sequence, initally starting with 0. If we denote L[X] the left
+;;; child of X and R[X] the right child of X, then we can see that the
+;;; denominator of L[X] is always the same as the numerator of
+;;; R[X]. Moreover, this carries over to other siblings in the same
+;;; generation. Therefore, we can compute CW[N] by taking pairwise
+;;; quotients of S: CW[N] = S[N]/S[N+1]
+
+;;; We can create a recurrence to generate such numbers by noting that
+;;; if X is the Nth element of the Calkin-Wilf sequence, then the
+;;; children will be 2N and 2N+1. With this numbering, 
+;;; 
+;;;                 S[2N] = S[N] = numerator of X,
+;;; 
+;;; and
+;;; 
+;;;   S[2N+1] = S[N] + S[N+1] = numerator of X + denominator of X.
+;;;
+;;; With some initial conditions, we have the recurrence relations
+;;; 
+;;;    S[0]    = 0
+;;;    S[1]    = 1
+;;;    S[2N]   = S[N]
+;;;    S[2N+1] = S[N] + S[N+1].
+
+(defun fusc (n)
+  "Compute the Nth element of Stern's diatomic sequence. The name FUSC
+  comes from Dijkstra."
+  (cond
+    ((zerop n) 0)
+    ((= n 1)   1)
+    ((evenp n) (fusc (floor n 2)))
+    ((oddp n)  (let ((n/2 (floor n 2)))
+                 (+ (fusc n/2) (fusc (1+ n/2)))))))
+
+;;; As said, we can compute CW[N] by the relation
+;;; 
+;;;           CW[n] = S[N]/S[N+1].
+
+(defun calkin-wilf-from-stern (n)
+  "Compute the Nth element of the Calkin-Wilf sequence using Stern's
+diatomic sequence."
+  (/ (fusc n)
+     (fusc (1+ n))))
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.