Source

lisp-random / calkin-wilf.lisp

;;;; 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.