Calkin-Wilf tree in Python
An infinite binary tree containing every positive rational exactly once in lowest terms, with children of a/b being a/(a+b) and (a+b)/b.
- File
calkin_wilf_tree.py, 93 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/calkin-wilf-tree/. The authors are listed in prov/R06GGM.jsonl.
"""The Calkin-Wilf tree, stored implicitly.
A node is a positive fraction a/b in lowest terms, given as the pair (a, b).
The root is (1, 1); the children of (a, b) are (a, a + b) and (a + b, b).
Nodes are numbered in breadth-first order from 1, so node n has children 2n
and 2n + 1. Every operation works on the pair itself; nothing is stored.
Numerators and denominators are Python integers of any size. Pairs with a
non-positive entry, or not in lowest terms, raise ValueError.
"""
def _check(a: int, b: int) -> None:
if a <= 0 or b <= 0:
raise ValueError("a node is a pair of positive integers")
def children(a: int, b: int) -> tuple[tuple[int, int], tuple[int, int]]:
"""The left and right children of a/b."""
_check(a, b)
return (a, a + b), (a + b, b) # @line L1
def parent(a: int, b: int) -> tuple[int, int] | None:
"""The parent of a/b, or None for the root 1/1."""
_check(a, b)
if a == b: # @line L2
if a != 1:
raise ValueError("not in lowest terms")
return None
if a < b: # @line L3
return a, b - a
return a - b, b # @line L4
def at(n: int) -> tuple[int, int]:
"""The node with breadth-first index n >= 1."""
if n < 1:
raise ValueError("indices start at 1")
a, b = 1, 1 # @line L5
for i in range(n.bit_length() - 2, -1, -1): # @line L6
if (n >> i) & 1 == 0: # @line L7
b = a + b
else:
a = a + b
return a, b # @line L8
def index_of(a: int, b: int) -> int:
"""The breadth-first index of a/b, which must be in lowest terms."""
_check(a, b)
n, length = 0, 0 # @line L9
while a != b: # @line L10
if a > b: # @line L11
q = (a - 1) // b # right steps: a - q b lands in (0, b]
a -= q * b
n |= ((1 << q) - 1) << length
else: # @line L12
q = (b - 1) // a # left steps: zero digits
b -= q * a
length += q # @line L13
if a != 1: # @line L14
raise ValueError("not in lowest terms")
return n | (1 << length) # @line L15
def next_node(a: int, b: int) -> tuple[int, int]:
"""The node after a/b in breadth-first order (Newman's formula)."""
_check(a, b)
q = a // b # @line L16
return b, (2 * q + 1) * b - a # @line L17
def depth(a: int, b: int) -> int:
"""The level of a/b; the root has depth 0."""
return index_of(a, b).bit_length() - 1
def fusc(n: int) -> int:
"""Stern's diatomic sequence: fusc(0) = 0, fusc(1) = 1,
fusc(2m) = fusc(m), fusc(2m + 1) = fusc(m) + fusc(m + 1)."""
if n < 0:
raise ValueError("fusc is defined for n >= 0")
x, y = 1, 0 # invariant: fusc(original n) = x fusc(n) + y fusc(n + 1)
while n > 0:
if n & 1:
y += x
else:
x += y
n >>= 1
return y
test_calkin_wilf_tree.py 166 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/calkin-wilf-tree/. The authors are listed in prov/R06GGM.jsonl.
import math
import random
from fractions import Fraction
import pytest
from calkin_wilf_tree import at, children, depth, fusc, index_of, next_node, parent
def fib(n: int) -> int:
a, b = 0, 1
for _ in range(n):
a, b = b, a + b
return a
def brute_fusc(limit: int) -> list[int]:
f = [0, 1] + [0] * limit
for n in range(2, limit + 1):
f[n] = f[n // 2] if n % 2 == 0 else f[n // 2] + f[n // 2 + 1]
return f
def brute_levels(levels: int) -> list[tuple[int, int]]:
"""Breadth-first list of the first levels, built from the child rule."""
order = [(1, 1)]
i = 0
while len(order) < (1 << levels) - 1:
a, b = order[i]
order.append((a, a + b))
order.append((a + b, b))
i += 1
return order
def cf_terms(a: int, b: int) -> list[int]:
terms = []
while b:
terms.append(a // b)
a, b = b, a % b
return terms
def test_example() -> None:
assert [at(n) for n in range(1, 16)] == [
(1, 1), (1, 2), (2, 1), (1, 3), (3, 2), (2, 3), (3, 1),
(1, 4), (4, 3), (3, 5), (5, 2), (2, 5), (5, 3), (3, 4), (4, 1),
]
assert index_of(3, 5) == 10
assert index_of(5, 2) == 11
assert index_of(11, 7) == 37
assert children(3, 5) == ((3, 8), (8, 5))
assert parent(3, 5) == (3, 2)
assert parent(1, 1) is None
assert next_node(3, 5) == (5, 2)
assert next_node(4, 1) == (1, 5)
def test_against_breadth_first_construction() -> None:
order = brute_levels(13)
f = brute_fusc(len(order) + 2)
for n, node in enumerate(order, start=1):
assert at(n) == node
assert index_of(*node) == n
assert at(n) == (f[n], f[n + 1])
assert fusc(n) == f[n]
if n < len(order):
assert next_node(*node) == order[n]
if n > 1:
assert parent(*node) == order[n // 2 - 1]
assert children(*order[n // 2 - 1])[n % 2] == node
def test_every_reduced_fraction_once() -> None:
seen = set()
for a in range(1, 80):
for b in range(1, 80):
if math.gcd(a, b) != 1: # reference only
with pytest.raises(ValueError):
index_of(a, b)
continue
n = index_of(a, b)
assert n not in seen
seen.add(n)
assert at(n) == (a, b)
assert depth(a, b) == sum(cf_terms(a, b)) - 1
def test_newman_formula_on_fractions() -> None:
x = Fraction(1)
for n in range(1, 3000):
assert (x.numerator, x.denominator) == at(n)
x = 1 / (2 * math.floor(x) - x + 1)
def test_powers_of_two_and_edges() -> None:
for k in range(0, 200):
assert at(1 << k) == (1, k + 1)
assert at((1 << (k + 1)) - 1) == (k + 1, 1)
assert index_of(k + 1, 1) == (1 << (k + 1)) - 1
assert index_of(1, k + 1) == 1 << k
def test_invalid_inputs() -> None:
for a, b in [(0, 1), (1, 0), (0, 0), (-1, 2), (2, -3), (-1, -1)]:
with pytest.raises(ValueError):
index_of(a, b)
with pytest.raises(ValueError):
children(a, b)
with pytest.raises(ValueError):
next_node(a, b)
for n in (0, -1, -(10**20)):
with pytest.raises(ValueError):
at(n)
with pytest.raises(ValueError):
fusc(-1)
with pytest.raises(ValueError):
parent(4, 4)
with pytest.raises(ValueError):
index_of(6, 4)
assert index_of(1, 1) == 1
assert fusc(0) == 0
def test_fibonacci_and_large_values() -> None:
for k in range(2, 300):
a, b = fib(k + 1), fib(k)
n = index_of(a, b)
assert at(n) == (a, b)
assert n.bit_length() - 1 == sum(cf_terms(a, b)) - 1
rng = random.Random(99)
done = 0
while done < 600:
digits = rng.choice([1, 3, 10, 50, 300, 400])
a = rng.randrange(1, 10**digits)
b = rng.randrange(1, 10**digits)
if math.gcd(a, b) != 1: # reference only
continue
terms = cf_terms(a, b)
if sum(terms) > 100000:
continue
n = index_of(a, b)
assert at(n) == (a, b)
p = parent(a, b)
if p is not None:
assert index_of(*p) == n // 2
left, right = children(a, b)
assert index_of(*left) == 2 * n
assert index_of(*right) == 2 * n + 1
assert at(n + 1) == next_node(a, b)
done += 1
def test_odd_length_continued_fraction_rule() -> None:
# The index read from its lowest digit is 1^{a0} 0^{a1} 1^{a2} ... for the
# continued fraction [a0; a1, ..., ak] of odd length.
for n in range(1, 5000):
a, b = at(n)
terms = cf_terms(a, b)
if len(terms) % 2 == 0:
terms = terms[:-1] + [terms[-1] - 1, 1]
digits = "".join(("1" if i % 2 == 0 else "0") * t for i, t in enumerate(terms))
assert int(digits[::-1], 2) == n
What the next level needs
- the tests did not run in this build