# lisp-random / calkin-wilf.lisp

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 ;;;; 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))))