diff --git a/.gitignore b/.gitignore
new file mode 100644
index 00000000..d4ed3165
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,28 @@
+# Python bytecode
+__pycache__/
+*.py[cod]
+*.pyo
+*.pyd
+
+# Distribution / packaging
+dist/
+build/
+*.egg-info/
+*.egg
+MANIFEST
+
+# Virtual environments
+.venv/
+venv/
+env/
+
+# Testing / coverage
+.pytest_cache/
+.coverage
+htmlcov/
+
+# Type checkers / editors
+.mypy_cache/
+.ruff_cache/
+.idea/
+.vscode/
diff --git a/README.md b/README.md
index 48ce88c8..678e6288 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,194 @@
-# ucns
-unit circle number system - the unit circle is a Mobius disk with recursive epicycles
+# ucns — Unit Circle Number System
+
+> **The unit circle is a Möbius disk with recursive epicycles.**
+
+A zero-dependency Python library for creating **compact, efficient embeddings**
+using a novel **Unit Circle Number System (UCNS)**.
+
+---
+
+## Why UCNS?
+
+Traditional dense embeddings (word2vec, BERT, OpenAI Ada, …) store each
+dimension as a 32-bit float.  UCNS embeddings encode every dimension as an
+**angle** θ ∈ \[0, 2π\) on the unit circle:
+
+| Property | float32 embedding | UCNS embedding |
+|---|---|---|
+| Bytes per dimension | 4 | 2 (uint16) |
+| Similarity computation | dot product + two L2 norms | mean cos(Δθ) – no normalisation |
+| Geometric space | Euclidean Rⁿ | Unit torus (S¹)ⁿ |
+| Hierarchical structure | no | yes (Möbius / Poincaré disk) |
+| External dependencies | numpy / torch / … | **none** |
+
+The compression and speed gains come from two structural properties:
+
+1. **All embeddings already live on the unit sphere** – the inner product
+   `cos(θᵢ − φᵢ)` never needs a length normalisation step.
+2. **Angles fit in 16 bits** – 0.0001 rad resolution with half the storage of
+   float32.
+
+---
+
+## Architecture
+
+```
+input data
+    │
+    ▼
+real-valued signal  (ordinals / floats / bytes …)
+    │
+    ▼  FFT  (Cooley–Tukey O(n log n), pure Python)
+Epicycle decomposition  ──►  amplitudes  +  phases
+                                              │
+                                              ▼
+                                    UCNS embedding vector
+                                    (list of n angles in [0, τ))
+```
+
+The **Möbius disk** (Poincaré disk model of the hyperbolic plane) is available
+as a companion geometry for encoding *hierarchical* relationships.  Points deep
+in the tree live near the boundary of the disk (high hyperbolic radius); root
+nodes sit near the centre.
+
+---
+
+## Installation
+
+```bash
+pip install ucns          # from PyPI (no dependencies)
+# or from source:
+pip install .
+```
+
+Python ≥ 3.8 required.  No third-party packages needed.
+
+---
+
+## Quick start
+
+```python
+from ucns import UCNEmbedding
+
+emb = UCNEmbedding(dim=64)
+
+# Encode any data to a list of 64 angles
+v1 = emb.encode("hello world")
+v2 = emb.encode("hello world")
+v3 = emb.encode("completely different")
+
+print(emb.similarity(v1, v2))   # 1.0  (identical)
+print(emb.similarity(v1, v3))   # < 1.0
+
+# Compact storage: 64 × 2 bytes = 128 bytes (vs 256 bytes for float32)
+packed = emb.encode_packed("hello world")
+print(len(packed))              # 128
+restored = UCNEmbedding.unpack(packed)
+
+# Nearest-neighbour search
+corpus = [emb.encode(w) for w in ["cat", "dog", "fish", "bird"]]
+idx, score = emb.nearest(emb.encode("cat"), corpus)
+print(idx, score)               # 0  1.0
+```
+
+---
+
+## API reference
+
+### `UCN` — core unit-circle number
+
+```python
+from ucns import UCN, TAU
+
+u = UCN(1.23)           # angle in radians, normalised to [0, τ)
+v = UCN.from_real(0.5)  # map float → UCN
+w = u * v               # rotation (angle addition)
+d = u.arc_distance(v)   # geodesic distance on S¹ ∈ [0, π]
+s = u.dot(v)            # cos(θ_u − θ_v) ∈ [−1, 1]
+b = u.to_bytes()        # 2-byte compact serialisation
+```
+
+### `EpicycleDecomposition` — FFT on the unit circle
+
+```python
+from ucns import EpicycleDecomposition
+
+d = EpicycleDecomposition([1, 2, 3, 4, 5, 6, 7, 8])
+print(d.amplitudes)          # per-frequency radii
+print(d.phases)              # per-frequency UCN angles
+print(d.reconstruct())       # lossless signal reconstruction
+sim = d.phase_similarity(d2) # amplitude-weighted phase cosine
+packed = d.pack()            # uint16 serialisation (2 bytes/freq)
+```
+
+### `MobiusTransform` — Möbius disk automorphisms
+
+```python
+from ucns import MobiusTransform, poincare_distance
+
+T = MobiusTransform(a=0.3 + 0.1j, phi=0.5)   # a ∈ open unit disk
+w = T(0.2 + 0j)                                # apply transform
+T_inv = T.inverse()                             # T_inv(T(z)) == z
+d = poincare_distance(0.1 + 0j, 0.5 + 0j)     # hyperbolic metric
+```
+
+### Similarity metrics
+
+```python
+from ucns.similarity import phase_cosine, arc_distance, hyperbolic_cosine, top_k_overlap
+
+phase_cosine(a, b)          # mean cos(θᵢ − φᵢ)  ∈ [−1, 1]
+arc_distance(a, b)          # mean arc distance,  ∈ [0, 1]
+hyperbolic_cosine(a, b)     # Poincaré-disk based ∈ [−1, 1]
+top_k_overlap(amps_a, amps_b, k=8)  # dominant-frequency Jaccard ∈ [0, 1]
+```
+
+---
+
+## Running the tests
+
+```bash
+python -m pytest tests/ -v
+# or without pytest:
+python -m unittest discover tests/
+```
+
+---
+
+## Mathematical background
+
+### Unit Circle Number System
+
+Every UCN is a point on the unit circle S¹ ⊂ ℂ:
+
+    z = e^(iθ),   θ ∈ [0, 2π)
+
+S¹ is a compact abelian group under multiplication.  Representing data as
+angles exploits this group structure: similarity becomes circular correlation,
+arithmetic becomes rotation, and conjugation becomes reflection.
+
+### Möbius disk
+
+The open unit disk D = {z ∈ ℂ : |z| < 1} with the Poincaré metric is a
+model of the hyperbolic plane.  Every conformal automorphism has the form
+
+    T_{a,φ}(z) = e^(iφ) · (z − a) / (1 − ā·z)
+
+These transformations preserve the circular boundary ∂D = S¹ and the
+hyperbolic metric  d(z,w) = 2 arctanh(|(z−w)/(1−w̄z)|).
+
+### Epicycles
+
+Any periodic signal can be written as a sum of circular motions:
+
+    x(t) = Σ_k  Aₖ · e^(i(2πkt/N + φₖ))
+
+This is the Fourier series interpreted geometrically.  The FFT computes the
+amplitudes *Aₖ* and phases *φₖ* in O(N log N) time.  UCNS stores only the
+phases (the angular part), giving a compact multi-scale fingerprint.
+
+---
+
+## License
+
+Apache 2.0 — see [LICENSE](LICENSE).
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 00000000..87dd054e
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,23 @@
+[build-system]
+requires = ["setuptools>=61"]
+build-backend = "setuptools.backends.legacy:build"
+
+[project]
+name = "ucns"
+version = "0.1.0"
+description = "Unit Circle Number System – zero-dependency compact embeddings via epicycle decomposition on the Möbius disk"
+readme = "README.md"
+license = { file = "LICENSE" }
+requires-python = ">=3.8"
+# No runtime dependencies – pure Python standard library only
+dependencies = []
+
+[project.optional-dependencies]
+dev = []
+
+[tool.setuptools.packages.find]
+where = ["."]
+include = ["ucns*"]
+
+[tool.setuptools.package-data]
+ucns = ["py.typed"]
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/test_core.py b/tests/test_core.py
new file mode 100644
index 00000000..6bc5d31c
--- /dev/null
+++ b/tests/test_core.py
@@ -0,0 +1,184 @@
+"""Tests for ucns.core – Unit Circle Number arithmetic."""
+
+import math
+import struct
+import unittest
+
+from ucns.core import UCN, TAU
+
+
+class TestUCNConstruction(unittest.TestCase):
+    def test_theta_normalised_positive(self):
+        u = UCN(3 * TAU + 1.0)
+        self.assertAlmostEqual(u.theta, 1.0, places=10)
+
+    def test_theta_normalised_negative(self):
+        u = UCN(-0.5)
+        self.assertAlmostEqual(u.theta, TAU - 0.5, places=10)
+
+    def test_theta_zero(self):
+        self.assertAlmostEqual(UCN(0.0).theta, 0.0)
+
+    def test_theta_tau_wraps_to_zero(self):
+        self.assertAlmostEqual(UCN(TAU).theta, 0.0, places=10)
+
+    def test_from_complex_unit_circle(self):
+        import cmath
+        z = cmath.exp(1j * 1.23)
+        u = UCN.from_complex(z)
+        self.assertAlmostEqual(u.theta, 1.23, places=10)
+
+    def test_from_complex_non_unit(self):
+        # phase is preserved regardless of magnitude
+        u = UCN.from_complex(2.0 + 2.0j)
+        self.assertAlmostEqual(u.theta, math.pi / 4, places=10)
+
+    def test_from_real_midpoint(self):
+        u = UCN.from_real(0.0, -1.0, 1.0)
+        self.assertAlmostEqual(u.theta, TAU / 2, places=10)
+
+    def test_from_real_lo(self):
+        u = UCN.from_real(-1.0, -1.0, 1.0)
+        self.assertAlmostEqual(u.theta, 0.0, places=10)
+
+    def test_from_real_hi(self):
+        u = UCN.from_real(1.0, -1.0, 1.0)
+        # t=1.0 → theta=τ, which normalises to 0
+        self.assertAlmostEqual(u.theta, 0.0, places=10)
+
+    def test_from_real_clamping(self):
+        u_lo = UCN.from_real(-999.0, -1.0, 1.0)
+        u_hi = UCN.from_real(999.0, -1.0, 1.0)
+        self.assertAlmostEqual(u_lo.theta, 0.0, places=10)
+        self.assertAlmostEqual(u_hi.theta, 0.0, places=10)  # 1.0 wraps to 0
+
+    def test_from_real_equal_lo_hi_raises(self):
+        with self.assertRaises(ValueError):
+            UCN.from_real(0.0, lo=1.0, hi=1.0)
+
+
+class TestUCNProperties(unittest.TestCase):
+    def test_real_imag(self):
+        theta = math.pi / 3
+        u = UCN(theta)
+        self.assertAlmostEqual(u.real, math.cos(theta))
+        self.assertAlmostEqual(u.imag, math.sin(theta))
+
+    def test_complex_on_unit_circle(self):
+        u = UCN(1.1)
+        self.assertAlmostEqual(abs(u.complex), 1.0, places=12)
+
+
+class TestUCNArithmetic(unittest.TestCase):
+    def test_multiplication_adds_angles(self):
+        u = UCN(1.0)
+        v = UCN(2.0)
+        self.assertAlmostEqual((u * v).theta, 3.0, places=10)
+
+    def test_multiplication_wraps(self):
+        u = UCN(TAU - 0.1)
+        v = UCN(0.2)
+        expected = (TAU - 0.1 + 0.2) % TAU
+        self.assertAlmostEqual((u * v).theta, expected, places=10)
+
+    def test_division_subtracts_angles(self):
+        u = UCN(3.0)
+        v = UCN(1.0)
+        self.assertAlmostEqual((u / v).theta, 2.0, places=10)
+
+    def test_conjugate(self):
+        u = UCN(1.0)
+        c = u.conjugate()
+        self.assertAlmostEqual(c.theta, TAU - 1.0, places=10)
+
+    def test_conjugate_of_zero(self):
+        u = UCN(0.0)
+        self.assertAlmostEqual(u.conjugate().theta, 0.0, places=10)
+
+    def test_mul_div_roundtrip(self):
+        u = UCN(1.5)
+        v = UCN(0.7)
+        self.assertAlmostEqual(((u * v) / v).theta, u.theta, places=10)
+
+    def test_identity_element(self):
+        u = UCN(1.23)
+        identity = UCN(0.0)
+        self.assertAlmostEqual((u * identity).theta, u.theta, places=10)
+
+
+class TestUCNMetrics(unittest.TestCase):
+    def test_dot_identical(self):
+        u = UCN(1.0)
+        self.assertAlmostEqual(u.dot(u), 1.0, places=10)
+
+    def test_dot_opposite(self):
+        u = UCN(0.0)
+        v = UCN(math.pi)
+        self.assertAlmostEqual(u.dot(v), -1.0, places=10)
+
+    def test_dot_quarter_turn(self):
+        u = UCN(0.0)
+        v = UCN(math.pi / 2)
+        self.assertAlmostEqual(u.dot(v), 0.0, places=10)
+
+    def test_arc_distance_same(self):
+        u = UCN(1.0)
+        self.assertAlmostEqual(u.arc_distance(u), 0.0, places=10)
+
+    def test_arc_distance_half_circle(self):
+        u = UCN(0.0)
+        v = UCN(math.pi)
+        self.assertAlmostEqual(u.arc_distance(v), math.pi, places=10)
+
+    def test_arc_distance_symmetry(self):
+        u = UCN(0.3)
+        v = UCN(2.7)
+        self.assertAlmostEqual(u.arc_distance(v), v.arc_distance(u), places=10)
+
+    def test_arc_distance_short_arc(self):
+        u = UCN(0.1)
+        v = UCN(TAU - 0.1)
+        self.assertAlmostEqual(u.arc_distance(v), 0.2, places=10)
+
+
+class TestUCNSerialisation(unittest.TestCase):
+    def test_int16_roundtrip(self):
+        for theta in [0.0, 1.0, math.pi, TAU * 0.75]:
+            u = UCN(theta)
+            v = UCN.from_int16(u.to_int16())
+            self.assertAlmostEqual(u.theta, v.theta, delta=TAU / 65535 + 1e-9)
+
+    def test_bytes_roundtrip(self):
+        u = UCN(2.718)
+        v = UCN.from_bytes(u.to_bytes())
+        self.assertAlmostEqual(u.theta, v.theta, delta=TAU / 65535 + 1e-9)
+
+    def test_bytes_length(self):
+        self.assertEqual(len(UCN(1.0).to_bytes()), 2)
+
+
+class TestUCNDunder(unittest.TestCase):
+    def test_repr(self):
+        r = repr(UCN(1.23456))
+        self.assertIn("UCN", r)
+        self.assertIn("1.23456", r)
+
+    def test_equality(self):
+        self.assertEqual(UCN(1.0), UCN(1.0 + TAU))
+
+    def test_inequality(self):
+        self.assertNotEqual(UCN(1.0), UCN(2.0))
+
+    def test_hash_consistent_with_eq(self):
+        self.assertEqual(hash(UCN(1.0)), hash(UCN(1.0 + TAU)))
+
+    def test_float_conversion(self):
+        u = UCN(1.23)
+        self.assertAlmostEqual(float(u), 1.23, places=10)
+
+    def test_lt(self):
+        self.assertLess(UCN(0.5), UCN(1.5))
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_embedding.py b/tests/test_embedding.py
new file mode 100644
index 00000000..2032533d
--- /dev/null
+++ b/tests/test_embedding.py
@@ -0,0 +1,202 @@
+"""Tests for ucns.embedding – UCNEmbedding API."""
+
+import math
+import struct
+import unittest
+
+from ucns.embedding import UCNEmbedding
+
+_TAU = 2.0 * math.pi
+
+
+class TestUCNEmbeddingConstruction(unittest.TestCase):
+    def test_dim_power_of_two(self):
+        emb = UCNEmbedding(dim=64)
+        self.assertEqual(emb.dim, 64)
+
+    def test_dim_rounded_up(self):
+        emb = UCNEmbedding(dim=50)
+        self.assertEqual(emb.dim, 64)
+
+    def test_dim_one(self):
+        emb = UCNEmbedding(dim=1)
+        self.assertEqual(emb.dim, 1)
+
+    def test_dim_zero_raises(self):
+        with self.assertRaises(ValueError):
+            UCNEmbedding(dim=0)
+
+    def test_repr_contains_dim(self):
+        emb = UCNEmbedding(dim=32)
+        self.assertIn("32", repr(emb))
+
+
+class TestUCNEmbeddingEncode(unittest.TestCase):
+    def setUp(self):
+        self.emb = UCNEmbedding(dim=16)
+
+    def test_encode_returns_list(self):
+        v = self.emb.encode("hello")
+        self.assertIsInstance(v, list)
+
+    def test_encode_correct_length(self):
+        v = self.emb.encode("hello")
+        self.assertEqual(len(v), self.emb.dim)
+
+    def test_encode_angles_in_range(self):
+        for data in ["hello", [1.0, 2.0, 3.0], 42.0, b"abc"]:
+            v = self.emb.encode(data)
+            for angle in v:
+                self.assertGreaterEqual(angle, 0.0)
+                self.assertLess(angle, _TAU + 1e-9)
+
+    def test_encode_deterministic(self):
+        v1 = self.emb.encode("hello world")
+        v2 = self.emb.encode("hello world")
+        self.assertEqual(v1, v2)
+
+    def test_encode_string(self):
+        v = self.emb.encode("abc")
+        self.assertEqual(len(v), self.emb.dim)
+
+    def test_encode_float(self):
+        v = self.emb.encode(3.14)
+        self.assertEqual(len(v), self.emb.dim)
+
+    def test_encode_int(self):
+        v = self.emb.encode(42)
+        self.assertEqual(len(v), self.emb.dim)
+
+    def test_encode_list(self):
+        v = self.emb.encode([1.0, 2.0, 3.0])
+        self.assertEqual(len(v), self.emb.dim)
+
+    def test_encode_bytes(self):
+        v = self.emb.encode(b"\x00\xff\x80")
+        self.assertEqual(len(v), self.emb.dim)
+
+    def test_encode_unsupported_type_raises(self):
+        with self.assertRaises(TypeError):
+            self.emb.encode({"key": "value"})
+
+    def test_encode_long_signal_truncated(self):
+        long_signal = list(range(self.emb.dim * 3))
+        v = self.emb.encode(long_signal)
+        self.assertEqual(len(v), self.emb.dim)
+
+    def test_encode_empty_string(self):
+        v = self.emb.encode("")
+        self.assertEqual(len(v), self.emb.dim)
+
+
+class TestUCNEmbeddingPackedSerialization(unittest.TestCase):
+    def setUp(self):
+        self.emb = UCNEmbedding(dim=16)
+
+    def test_packed_length(self):
+        packed = self.emb.encode_packed("hello")
+        self.assertEqual(len(packed), self.emb.dim * 2)
+
+    def test_pack_unpack_roundtrip(self):
+        v = self.emb.encode("hello world")
+        packed = self.emb.encode_packed("hello world")
+        unpacked = UCNEmbedding.unpack(packed)
+        for orig, rec in zip(v, unpacked):
+            self.assertAlmostEqual(orig, rec, delta=_TAU / 65535 + 1e-6)
+
+    def test_packed_is_bytes(self):
+        self.assertIsInstance(self.emb.encode_packed("test"), bytes)
+
+
+class TestUCNEmbeddingSimilarity(unittest.TestCase):
+    def setUp(self):
+        self.emb = UCNEmbedding(dim=32)
+
+    def test_similarity_self(self):
+        v = self.emb.encode("hello")
+        self.assertAlmostEqual(self.emb.similarity(v, v), 1.0, places=10)
+
+    def test_similarity_in_range(self):
+        v1 = self.emb.encode("hello")
+        v2 = self.emb.encode("world")
+        s = self.emb.similarity(v1, v2)
+        self.assertGreaterEqual(s, -1.0)
+        self.assertLessEqual(s, 1.0)
+
+    def test_similarity_symmetric(self):
+        v1 = self.emb.encode("foo")
+        v2 = self.emb.encode("bar")
+        self.assertAlmostEqual(
+            self.emb.similarity(v1, v2),
+            self.emb.similarity(v2, v1),
+            places=12,
+        )
+
+    def test_similarity_different_lengths_raises(self):
+        with self.assertRaises(ValueError):
+            self.emb.similarity([1.0, 2.0], [1.0])
+
+    def test_similar_inputs_higher_score(self):
+        v1 = self.emb.encode("the cat sat")
+        v2 = self.emb.encode("the cat sat")  # identical
+        v3 = self.emb.encode("completely different xyz 999")
+        s_same = self.emb.similarity(v1, v2)
+        s_diff = self.emb.similarity(v1, v3)
+        self.assertGreater(s_same, s_diff)
+
+
+class TestUCNEmbeddingNearest(unittest.TestCase):
+    def setUp(self):
+        self.emb = UCNEmbedding(dim=32)
+
+    def test_nearest_finds_identical(self):
+        query = self.emb.encode("apple")
+        corpus = [
+            self.emb.encode("banana"),
+            self.emb.encode("apple"),
+            self.emb.encode("cherry"),
+        ]
+        idx, score = self.emb.nearest(query, corpus)
+        self.assertEqual(idx, 1)
+        self.assertAlmostEqual(score, 1.0, places=10)
+
+    def test_nearest_empty_corpus_raises(self):
+        query = self.emb.encode("hello")
+        with self.assertRaises(ValueError):
+            self.emb.nearest(query, [])
+
+    def test_nearest_returns_valid_index(self):
+        query = self.emb.encode("test")
+        corpus = [self.emb.encode(str(i)) for i in range(10)]
+        idx, score = self.emb.nearest(query, corpus)
+        self.assertGreaterEqual(idx, 0)
+        self.assertLess(idx, 10)
+
+
+class TestUCNEmbeddingToSignal(unittest.TestCase):
+    def test_int_gives_one_element(self):
+        sig = UCNEmbedding._to_signal(5)
+        self.assertEqual(sig, [5.0])
+
+    def test_float_gives_one_element(self):
+        sig = UCNEmbedding._to_signal(3.14)
+        self.assertAlmostEqual(sig[0], 3.14)
+
+    def test_str_encodes_ordinals(self):
+        sig = UCNEmbedding._to_signal("AB")
+        self.assertAlmostEqual(sig[0], 65.0)
+        self.assertAlmostEqual(sig[1], 66.0)
+
+    def test_bytes_encodes_byte_values(self):
+        sig = UCNEmbedding._to_signal(b"\x00\x01\xff")
+        self.assertAlmostEqual(sig[0], 0.0)
+        self.assertAlmostEqual(sig[1], 1.0)
+        self.assertAlmostEqual(sig[2], 255.0)
+
+    def test_list_passthrough(self):
+        sig = UCNEmbedding._to_signal([1.5, 2.5])
+        self.assertEqual(sig, [1.5, 2.5])
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_epicycle.py b/tests/test_epicycle.py
new file mode 100644
index 00000000..b149a9b7
--- /dev/null
+++ b/tests/test_epicycle.py
@@ -0,0 +1,165 @@
+"""Tests for ucns.epicycle – FFT and EpicycleDecomposition."""
+
+import cmath
+import math
+import unittest
+
+from ucns.epicycle import EpicycleDecomposition, fft, ifft, _next_pow2
+
+_TAU = 2.0 * math.pi
+
+
+class TestNextPow2(unittest.TestCase):
+    def test_exact_power_of_two(self):
+        for p in [1, 2, 4, 8, 16, 32, 64]:
+            self.assertEqual(_next_pow2(p), p)
+
+    def test_non_power_of_two(self):
+        self.assertEqual(_next_pow2(3), 4)
+        self.assertEqual(_next_pow2(5), 8)
+        self.assertEqual(_next_pow2(100), 128)
+
+
+class TestFFT(unittest.TestCase):
+    def _slow_dft(self, x):
+        n = len(x)
+        return [
+            sum(x[j] * cmath.exp(-1j * _TAU * k * j / n) for j in range(n))
+            for k in range(n)
+        ]
+
+    def test_fft_equals_dft_n4(self):
+        signal = [1.0, 2.0, 3.0, 4.0]
+        fast = fft(signal)
+        slow = self._slow_dft(signal)
+        for a, b in zip(fast, slow):
+            self.assertAlmostEqual(a, b, places=10)
+
+    def test_fft_equals_dft_n8(self):
+        import math
+        signal = [math.sin(_TAU * k / 8) for k in range(8)]
+        fast = fft(signal)
+        slow = self._slow_dft(signal)
+        for a, b in zip(fast, slow):
+            self.assertAlmostEqual(a, b, places=10)
+
+    def test_fft_zero_pads(self):
+        # Length 3 → padded to 4
+        result = fft([1.0, 2.0, 3.0])
+        self.assertEqual(len(result), 4)
+
+    def test_fft_single_element(self):
+        result = fft([5.0])
+        self.assertAlmostEqual(result[0], 5.0, places=12)
+
+    def test_fft_dc_component(self):
+        """DC component (k=0) should equal sum of the signal."""
+        signal = [1.0, 2.0, 3.0, 4.0]
+        result = fft(signal)
+        self.assertAlmostEqual(result[0].real, sum(signal), places=10)
+
+    def test_ifft_fft_roundtrip(self):
+        signal = [1.0, -2.0, 3.5, -0.5, 0.0, 1.0, 2.0, -3.0]
+        restored = ifft(fft(signal))
+        for orig, rec in zip(signal, restored):
+            self.assertAlmostEqual(orig, rec.real, places=10)
+            self.assertAlmostEqual(rec.imag, 0.0, places=10)
+
+    def test_ifft_empty(self):
+        self.assertEqual(ifft([]), [])
+
+    def test_parseval_theorem(self):
+        """Energy in time domain == energy in frequency domain / N."""
+        signal = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+        spectrum = fft(signal)
+        n = len(spectrum)
+        energy_time = sum(x * x for x in signal)
+        energy_freq = sum(abs(X) ** 2 for X in spectrum) / n
+        self.assertAlmostEqual(energy_time, energy_freq, places=8)
+
+
+class TestEpicycleDecomposition(unittest.TestCase):
+    def test_empty_raises(self):
+        with self.assertRaises(ValueError):
+            EpicycleDecomposition([])
+
+    def test_lengths(self):
+        d = EpicycleDecomposition([1.0, 2.0, 3.0, 4.0])
+        self.assertEqual(len(d.amplitudes), d.n)
+        self.assertEqual(len(d.phases), d.n)
+        self.assertEqual(len(d.frequencies), d.n)
+
+    def test_n_is_power_of_two(self):
+        d = EpicycleDecomposition([1.0, 2.0, 3.0])
+        self.assertEqual(d.n, 4)
+
+    def test_amplitudes_nonnegative(self):
+        d = EpicycleDecomposition([1.0, -2.0, 3.0, -4.0])
+        for a in d.amplitudes:
+            self.assertGreaterEqual(a, 0.0)
+
+    def test_phases_in_range(self):
+        d = EpicycleDecomposition([1.5, 2.5, 0.5, -1.0])
+        for p in d.phases:
+            self.assertGreaterEqual(p, 0.0)
+            self.assertLess(p, _TAU + 1e-9)
+
+    def test_reconstruction_lossless(self):
+        signal = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]
+        d = EpicycleDecomposition(signal)
+        rec = d.reconstruct()
+        for orig, r in zip(signal, rec):
+            self.assertAlmostEqual(orig, r, places=8)
+
+    def test_reconstruction_padded(self):
+        """Non-power-of-2 input: first n_orig values should reconstruct."""
+        signal = [3.0, 1.0, 4.0, 1.0, 5.0]  # length 5 → padded to 8
+        d = EpicycleDecomposition(signal)
+        rec = d.reconstruct()
+        self.assertEqual(len(rec), len(signal))
+        for orig, r in zip(signal, rec):
+            self.assertAlmostEqual(orig, r, places=8)
+
+    def test_phase_vector_length(self):
+        d = EpicycleDecomposition([1.0, 2.0, 3.0, 4.0])
+        self.assertEqual(len(d.phase_vector), d.n)
+
+    def test_phase_similarity_self(self):
+        d = EpicycleDecomposition([1.0, 2.0, 3.0, 4.0])
+        self.assertAlmostEqual(d.phase_similarity(d), 1.0, places=10)
+
+    def test_phase_similarity_range(self):
+        d1 = EpicycleDecomposition([1.0, 0.0, -1.0, 0.0])
+        d2 = EpicycleDecomposition([0.0, 1.0, 0.0, -1.0])
+        sim = d1.phase_similarity(d2)
+        self.assertGreaterEqual(sim, -1.0)
+        self.assertLessEqual(sim, 1.0)
+
+    def test_pack_unpack_roundtrip(self):
+        signal = [1.0, 2.0, -1.0, 0.5, 0.0, -0.5, 1.5, -2.0]
+        d = EpicycleDecomposition(signal)
+        packed = d.pack()
+        self.assertEqual(len(packed), d.n * 2)
+        unpacked = EpicycleDecomposition.unpack_phases(packed)
+        for orig, rec in zip(d.phases, unpacked):
+            self.assertAlmostEqual(orig, rec, delta=_TAU / 65535 + 1e-6)
+
+    def test_dominant_frequency(self):
+        """Pure sine at frequency k should give dominant_frequency ≈ k."""
+        n = 16
+        k = 3
+        signal = [math.sin(_TAU * k * j / n) for j in range(n)]
+        d = EpicycleDecomposition(signal)
+        self.assertIn(d.dominant_frequency, {k, n - k})
+
+    def test_repr(self):
+        d = EpicycleDecomposition([1.0, 2.0, 3.0, 4.0])
+        self.assertIn("EpicycleDecomposition", repr(d))
+
+    def test_len(self):
+        d = EpicycleDecomposition([1.0, 2.0])
+        self.assertEqual(len(d), d.n)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_mobius.py b/tests/test_mobius.py
new file mode 100644
index 00000000..ec3c2c79
--- /dev/null
+++ b/tests/test_mobius.py
@@ -0,0 +1,145 @@
+"""Tests for ucns.mobius – Möbius disk transformations."""
+
+import cmath
+import math
+import unittest
+
+from ucns.mobius import (
+    MobiusTransform,
+    circle_to_disk,
+    disk_to_circle,
+    poincare_distance,
+)
+
+_TAU = 2.0 * math.pi
+
+
+class TestMobiusConstruction(unittest.TestCase):
+    def test_valid_construction(self):
+        t = MobiusTransform(0.5 + 0j)
+        self.assertEqual(t.a, 0.5 + 0j)
+
+    def test_a_on_boundary_raises(self):
+        with self.assertRaises(ValueError):
+            MobiusTransform(1.0 + 0j)
+
+    def test_a_outside_disk_raises(self):
+        with self.assertRaises(ValueError):
+            MobiusTransform(1.5 + 0j)
+
+    def test_phi_normalised(self):
+        t = MobiusTransform(0j, phi=_TAU + 1.0)
+        self.assertAlmostEqual(t.phi, 1.0, places=10)
+
+
+class TestMobiusEvaluation(unittest.TestCase):
+    def test_maps_a_to_zero(self):
+        a = 0.3 + 0.2j
+        t = MobiusTransform(a)
+        result = t(a)
+        self.assertAlmostEqual(abs(result), 0.0, places=12)
+
+    def test_maps_zero_to_minus_a_rotated(self):
+        a = 0.4 + 0j
+        phi = 0.5
+        t = MobiusTransform(a, phi=phi)
+        expected = cmath.exp(1j * phi) * (-a)
+        self.assertAlmostEqual(t(0j), expected, places=12)
+
+    def test_preserves_unit_circle(self):
+        """Points on ∂D should map to ∂D."""
+        a = 0.3 + 0.1j
+        t = MobiusTransform(a, phi=0.7)
+        for angle in [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]:
+            z = cmath.exp(1j * angle)
+            w = t(z)
+            self.assertAlmostEqual(abs(w), 1.0, places=10,
+                                   msg=f"Unit circle not preserved at angle={angle}")
+
+    def test_maps_disk_to_disk(self):
+        """Interior points should stay inside D."""
+        a = 0.2 - 0.3j
+        t = MobiusTransform(a, phi=1.2)
+        for z in [0.1 + 0j, 0.5j, -0.4 - 0.4j]:
+            self.assertLess(abs(t(z)), 1.0)
+
+
+class TestMobiusInverse(unittest.TestCase):
+    def test_inverse_roundtrip(self):
+        a = 0.35 + 0.15j
+        t = MobiusTransform(a, phi=0.9)
+        t_inv = t.inverse()
+        for z in [0j, 0.1 + 0.2j, -0.3 + 0.0j]:
+            self.assertAlmostEqual(t_inv(t(z)), z, places=10)
+
+    def test_identity_transform(self):
+        t = MobiusTransform(0j, phi=0.0)
+        for z in [0.3 + 0j, -0.5j, 0.1 + 0.1j]:
+            self.assertAlmostEqual(t(z), z, places=12)
+
+
+class TestPoincareDistance(unittest.TestCase):
+    def test_distance_zero_with_itself(self):
+        z = 0.3 + 0.2j
+        self.assertAlmostEqual(poincare_distance(z, z), 0.0, places=12)
+
+    def test_distance_positive(self):
+        d = poincare_distance(0.1 + 0j, -0.1 + 0j)
+        self.assertGreater(d, 0.0)
+
+    def test_distance_symmetric(self):
+        z = 0.2 + 0.1j
+        w = -0.3 + 0.15j
+        self.assertAlmostEqual(poincare_distance(z, w), poincare_distance(w, z), places=12)
+
+    def test_distance_grows_near_boundary(self):
+        """Points near the boundary should be "far" from the origin."""
+        d_near = poincare_distance(0j, 0.1 + 0j)
+        d_far = poincare_distance(0j, 0.99 + 0j)
+        self.assertGreater(d_far, d_near)
+
+    def test_outside_disk_raises(self):
+        with self.assertRaises(ValueError):
+            poincare_distance(1.5 + 0j, 0j)
+
+    def test_on_boundary_raises(self):
+        with self.assertRaises(ValueError):
+            poincare_distance(1.0 + 0j, 0j)
+
+    def test_mobius_isometry(self):
+        """Möbius transform must preserve hyperbolic distance."""
+        z = 0.3 + 0.1j
+        w = -0.2 + 0.4j
+        t = MobiusTransform(0.1 + 0.05j, phi=0.3)
+        d_before = poincare_distance(z, w)
+        d_after = poincare_distance(t(z), t(w))
+        self.assertAlmostEqual(d_before, d_after, places=8)
+
+
+class TestDiskCircleProjection(unittest.TestCase):
+    def test_disk_to_circle_angle(self):
+        theta = 1.2
+        z = 0.7 * cmath.exp(1j * theta)
+        self.assertAlmostEqual(disk_to_circle(z), theta, places=10)
+
+    def test_disk_to_circle_zero(self):
+        self.assertAlmostEqual(disk_to_circle(0j), 0.0, places=10)
+
+    def test_circle_to_disk_modulus(self):
+        z = circle_to_disk(1.0, r=0.5)
+        self.assertAlmostEqual(abs(z), 0.5, places=12)
+
+    def test_circle_to_disk_angle(self):
+        theta = 2.5
+        z = circle_to_disk(theta, r=0.3)
+        self.assertAlmostEqual(cmath.phase(z), theta, places=10)
+
+    def test_circle_to_disk_invalid_radius(self):
+        with self.assertRaises(ValueError):
+            circle_to_disk(1.0, r=0.0)
+        with self.assertRaises(ValueError):
+            circle_to_disk(1.0, r=1.0)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_similarity.py b/tests/test_similarity.py
new file mode 100644
index 00000000..e8b089d5
--- /dev/null
+++ b/tests/test_similarity.py
@@ -0,0 +1,143 @@
+"""Tests for ucns.similarity – metric functions."""
+
+import math
+import unittest
+
+from ucns.similarity import arc_distance, hyperbolic_cosine, phase_cosine, top_k_overlap
+
+_TAU = 2.0 * math.pi
+
+
+class TestPhaseCosine(unittest.TestCase):
+    def test_identical_gives_one(self):
+        a = [0.1, 0.5, 1.2, 2.0]
+        self.assertAlmostEqual(phase_cosine(a, a), 1.0, places=12)
+
+    def test_opposite_gives_minus_one(self):
+        a = [0.0]
+        b = [math.pi]
+        self.assertAlmostEqual(phase_cosine(a, b), -1.0, places=12)
+
+    def test_quarter_turn_gives_zero(self):
+        a = [0.0]
+        b = [math.pi / 2]
+        self.assertAlmostEqual(phase_cosine(a, b), 0.0, places=12)
+
+    def test_symmetric(self):
+        a = [0.1, 0.5, 1.0]
+        b = [0.3, 0.7, 0.2]
+        self.assertAlmostEqual(phase_cosine(a, b), phase_cosine(b, a), places=12)
+
+    def test_empty_returns_zero(self):
+        self.assertEqual(phase_cosine([], []), 0.0)
+
+    def test_different_lengths_raises(self):
+        with self.assertRaises(ValueError):
+            phase_cosine([1.0], [1.0, 2.0])
+
+
+class TestArcDistance(unittest.TestCase):
+    def test_identical_gives_zero(self):
+        a = [0.5, 1.0, 2.0]
+        self.assertAlmostEqual(arc_distance(a, a), 0.0, places=12)
+
+    def test_half_circle_gives_one(self):
+        a = [0.0]
+        b = [math.pi]
+        self.assertAlmostEqual(arc_distance(a, b), 1.0, places=12)
+
+    def test_short_arc_used(self):
+        """Going 0 → 2π-0.1 should use short arc 0.1, not 2π-0.1."""
+        a = [0.0]
+        b = [_TAU - 0.1]
+        result = arc_distance(a, b)
+        self.assertAlmostEqual(result, 0.1 / math.pi, places=10)
+
+    def test_symmetric(self):
+        a = [0.3, 1.1, 2.5]
+        b = [1.0, 0.5, 0.1]
+        self.assertAlmostEqual(arc_distance(a, b), arc_distance(b, a), places=12)
+
+    def test_result_in_range(self):
+        import random
+        rng = random.Random(42)
+        a = [rng.uniform(0, _TAU) for _ in range(20)]
+        b = [rng.uniform(0, _TAU) for _ in range(20)]
+        result = arc_distance(a, b)
+        self.assertGreaterEqual(result, 0.0)
+        self.assertLessEqual(result, 1.0)
+
+    def test_empty_returns_zero(self):
+        self.assertEqual(arc_distance([], []), 0.0)
+
+    def test_different_lengths_raises(self):
+        with self.assertRaises(ValueError):
+            arc_distance([1.0], [1.0, 2.0])
+
+
+class TestHyperbolicCosine(unittest.TestCase):
+    def test_identical_gives_one(self):
+        a = [0.1, 0.5, 1.2]
+        self.assertAlmostEqual(hyperbolic_cosine(a, a), 1.0, places=10)
+
+    def test_result_in_range(self):
+        a = [0.0, 1.0, 2.0, 3.0]
+        b = [math.pi, 0.5, 1.5, 2.5]
+        result = hyperbolic_cosine(a, b)
+        self.assertGreaterEqual(result, -1.0)
+        self.assertLessEqual(result, 1.0)
+
+    def test_symmetric(self):
+        a = [0.3, 1.1]
+        b = [1.0, 0.5]
+        self.assertAlmostEqual(
+            hyperbolic_cosine(a, b),
+            hyperbolic_cosine(b, a),
+            places=12,
+        )
+
+    def test_invalid_radius_raises(self):
+        with self.assertRaises(ValueError):
+            hyperbolic_cosine([0.0], [0.0], radius=0.0)
+        with self.assertRaises(ValueError):
+            hyperbolic_cosine([0.0], [0.0], radius=1.0)
+
+    def test_empty_returns_zero(self):
+        self.assertEqual(hyperbolic_cosine([], []), 0.0)
+
+    def test_different_lengths_raises(self):
+        with self.assertRaises(ValueError):
+            hyperbolic_cosine([1.0], [1.0, 2.0])
+
+
+class TestTopKOverlap(unittest.TestCase):
+    def test_identical_amplitudes_gives_one(self):
+        a = [3.0, 1.0, 2.0, 0.5]
+        self.assertAlmostEqual(top_k_overlap(a, a, k=2), 1.0, places=12)
+
+    def test_completely_different_gives_zero(self):
+        a = [1.0, 0.0, 0.0, 0.0]
+        b = [0.0, 0.0, 0.0, 1.0]
+        self.assertAlmostEqual(top_k_overlap(a, b, k=1), 0.0, places=12)
+
+    def test_result_in_range(self):
+        a = [3.0, 1.0, 2.0, 4.0]
+        b = [1.0, 4.0, 0.5, 2.0]
+        result = top_k_overlap(a, b, k=2)
+        self.assertGreaterEqual(result, 0.0)
+        self.assertLessEqual(result, 1.0)
+
+    def test_k_larger_than_n(self):
+        """k is clamped to len(amplitudes)."""
+        a = [1.0, 2.0]
+        b = [2.0, 1.0]
+        # with k=100 (clamped to 2), both sets are the same → overlap = 1
+        self.assertAlmostEqual(top_k_overlap(a, b, k=100), 1.0, places=12)
+
+    def test_different_lengths_raises(self):
+        with self.assertRaises(ValueError):
+            top_k_overlap([1.0], [1.0, 2.0], k=1)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/ucns/__init__.py b/ucns/__init__.py
new file mode 100644
index 00000000..49026c3a
--- /dev/null
+++ b/ucns/__init__.py
@@ -0,0 +1,76 @@
+"""
+ucns – Unit Circle Number System
+=================================
+A zero-dependency Python library for creating compact, efficient embeddings
+using a novel **Unit Circle Number System (UCNS)**.
+
+Overview
+--------
+Every number in UCNS is an *angle* θ ∈ [0, 2π) that identifies a point on the
+unit circle e^(iθ).  A sequence of such angles forms a UCNS embedding:
+
+* **Compact** – angles stored as ``uint16`` use only 2 bytes per dimension
+  (2× smaller than float32).
+* **Fast similarity** – inner product = mean of cos(θᵢ − φᵢ); no length
+  normalisation required.
+* **Hierarchical** – the recursive epicycle (FFT) structure captures
+  multi-scale patterns; the Möbius disk geometry supports hyperbolic
+  (tree-like) relationships.
+* **Zero dependencies** – pure Python standard library only.
+
+Quick start
+-----------
+>>> from ucns import UCNEmbedding
+>>> emb = UCNEmbedding(dim=64)
+>>> v1 = emb.encode("hello world")
+>>> v2 = emb.encode("hello world")
+>>> emb.similarity(v1, v2)
+1.0
+>>> packed = emb.encode_packed("hello world")
+>>> len(packed)   # 64 angles × 2 bytes
+128
+
+Building blocks
+---------------
+``UCN``
+    Single unit-circle number (angle + arithmetic).
+``EpicycleDecomposition``
+    Decompose any signal into weighted unit-circle rotations via FFT.
+``MobiusTransform``
+    Conformal automorphism of the Poincaré disk.
+``UCNEmbedding``
+    High-level embedding API.
+Similarity functions
+    ``phase_cosine``, ``arc_distance``, ``hyperbolic_cosine``,
+    ``top_k_overlap``.
+"""
+
+from .core import UCN, TAU
+from .epicycle import EpicycleDecomposition, fft, ifft
+from .embedding import UCNEmbedding
+from .mobius import MobiusTransform, poincare_distance, disk_to_circle, circle_to_disk
+from .similarity import phase_cosine, arc_distance, hyperbolic_cosine, top_k_overlap
+
+__all__ = [
+    # Core number type
+    "UCN",
+    "TAU",
+    # Epicycle / FFT
+    "EpicycleDecomposition",
+    "fft",
+    "ifft",
+    # Möbius disk
+    "MobiusTransform",
+    "poincare_distance",
+    "disk_to_circle",
+    "circle_to_disk",
+    # Embedding
+    "UCNEmbedding",
+    # Similarity metrics
+    "phase_cosine",
+    "arc_distance",
+    "hyperbolic_cosine",
+    "top_k_overlap",
+]
+
+__version__ = "0.1.0"
diff --git a/ucns/core.py b/ucns/core.py
new file mode 100644
index 00000000..046eb828
--- /dev/null
+++ b/ucns/core.py
@@ -0,0 +1,165 @@
+"""
+ucns.core
+=========
+Unit Circle Number (UCN) – the fundamental numeric primitive.
+
+Every UCN is an angle θ ∈ [0, 2π) that identifies a point on the unit circle
+e^(iθ) ∈ ℂ.  Because |e^(iθ)| = 1 for all θ, the set of all UCNs forms a
+compact abelian group under multiplication (rotation), making them a natural
+substrate for periodic / cyclic data and for efficient angular embeddings.
+
+Key properties
+--------------
+* **Closure**: multiplying two UCNs (adding angles) always stays on the unit
+  circle.
+* **Compact storage**: an angle fits in a 16-bit integer (65 536 steps vs 32
+  bits for a single-precision float).
+* **Fast inner product**: dot(u, v) = cos(θ_u − θ_v) – no square root needed.
+* **No external dependencies**: pure Python / math stdlib only.
+"""
+
+from __future__ import annotations
+
+import math
+import cmath
+import struct
+
+__all__ = ["UCN", "TAU"]
+
+TAU: float = 2.0 * math.pi  # full turn = τ
+
+
+class UCN:
+    """Unit Circle Number – a real number encoded as an angle on the unit circle.
+
+    Parameters
+    ----------
+    theta:
+        Angle in radians.  Automatically reduced modulo τ = 2π so that
+        ``self.theta`` is always in ``[0, τ)``.
+    """
+
+    __slots__ = ("_theta",)
+
+    def __init__(self, theta: float) -> None:
+        self._theta: float = float(theta) % TAU
+
+    # ------------------------------------------------------------------
+    # Properties
+    # ------------------------------------------------------------------
+
+    @property
+    def theta(self) -> float:
+        """Angle in radians, normalised to ``[0, τ)``."""
+        return self._theta
+
+    @property
+    def real(self) -> float:
+        """Real part of the corresponding unit-circle point: cos θ."""
+        return math.cos(self._theta)
+
+    @property
+    def imag(self) -> float:
+        """Imaginary part of the corresponding unit-circle point: sin θ."""
+        return math.sin(self._theta)
+
+    @property
+    def complex(self) -> complex:
+        """The unit-circle point as a Python ``complex``: e^(iθ)."""
+        return cmath.exp(1j * self._theta)
+
+    # ------------------------------------------------------------------
+    # Constructors
+    # ------------------------------------------------------------------
+
+    @classmethod
+    def from_complex(cls, z: complex) -> "UCN":
+        """Project any complex number *z* onto the unit circle (keep phase)."""
+        return cls(cmath.phase(z))
+
+    @classmethod
+    def from_real(cls, x: float, lo: float = -1.0, hi: float = 1.0) -> "UCN":
+        """Map a real number *x* ∈ [lo, hi] uniformly onto [0, τ).
+
+        Values outside ``[lo, hi]`` are clamped before mapping.
+        """
+        if hi == lo:
+            raise ValueError("lo and hi must differ")
+        t = (max(lo, min(hi, x)) - lo) / (hi - lo)  # ∈ [0, 1]
+        return cls(t * TAU)
+
+    # ------------------------------------------------------------------
+    # Group arithmetic  (unit circle = ℝ/τℤ)
+    # ------------------------------------------------------------------
+
+    def __mul__(self, other: "UCN") -> "UCN":
+        """Rotation: θ₁ ⊗ θ₂ ≡ θ₁ + θ₂  (mod τ)."""
+        return UCN(self._theta + other._theta)
+
+    def __truediv__(self, other: "UCN") -> "UCN":
+        """Inverse rotation: θ₁ ⊘ θ₂ ≡ θ₁ − θ₂  (mod τ)."""
+        return UCN(self._theta - other._theta)
+
+    def conjugate(self) -> "UCN":
+        """Conjugate (reflection): θ* ≡ −θ  (mod τ)."""
+        return UCN(-self._theta)
+
+    # ------------------------------------------------------------------
+    # Metric / similarity
+    # ------------------------------------------------------------------
+
+    def dot(self, other: "UCN") -> float:
+        """Angular inner product: cos(θ_self − θ_other) ∈ [−1, 1]."""
+        return math.cos(self._theta - other._theta)
+
+    def arc_distance(self, other: "UCN") -> float:
+        """Geodesic (arc-length) distance on the unit circle ∈ [0, π]."""
+        diff = abs(self._theta - other._theta) % TAU
+        return min(diff, TAU - diff)
+
+    # ------------------------------------------------------------------
+    # Compact serialisation
+    # ------------------------------------------------------------------
+
+    def to_int16(self) -> int:
+        """Quantise angle to an unsigned 16-bit integer (0 … 65 535).
+
+        Provides ~0.0001 rad (≈ 0.006°) angular resolution with only 2 bytes.
+        """
+        return min(65535, int(self._theta * 65535.0 / TAU))
+
+    @classmethod
+    def from_int16(cls, v: int) -> "UCN":
+        """Restore a UCN from a 16-bit quantised integer."""
+        return cls(int(v) * TAU / 65535.0)
+
+    def to_bytes(self) -> bytes:
+        """Serialise to 2 bytes (little-endian unsigned short)."""
+        return struct.pack("<H", self.to_int16())
+
+    @classmethod
+    def from_bytes(cls, data: bytes) -> "UCN":
+        """Deserialise from 2 bytes."""
+        (v,) = struct.unpack("<H", data[:2])
+        return cls.from_int16(v)
+
+    # ------------------------------------------------------------------
+    # Dunder helpers
+    # ------------------------------------------------------------------
+
+    def __repr__(self) -> str:
+        return f"UCN({self._theta:.6f})"
+
+    def __eq__(self, other: object) -> bool:
+        if not isinstance(other, UCN):
+            return NotImplemented
+        return abs(self._theta - other._theta) < 1e-9
+
+    def __hash__(self) -> int:
+        return hash(round(self._theta, 9))
+
+    def __float__(self) -> float:
+        return self._theta
+
+    def __lt__(self, other: "UCN") -> bool:
+        return self._theta < other._theta
diff --git a/ucns/embedding.py b/ucns/embedding.py
new file mode 100644
index 00000000..0585956d
--- /dev/null
+++ b/ucns/embedding.py
@@ -0,0 +1,212 @@
+"""
+ucns.embedding
+==============
+High-level embedding API built on the Unit Circle Number System.
+
+**Why UCNS embeddings?**
+
+Traditional neural embeddings (e.g. word2vec, BERT, OpenAI Ada) produce
+*dense float32 vectors* of 384–3072 dimensions.  UCNS embeddings offer three
+concrete advantages:
+
+1. **Compact storage** – angles are quantised to ``uint16`` (2 bytes), giving
+   a 2× space saving over ``float32`` with negligible information loss
+   (≈0.0001 rad angular resolution).
+
+2. **Fast similarity** – the inner product ``cos(θᵢ − φᵢ)`` is computed with
+   a single subtraction and cosine lookup per dimension.  No square root for
+   normalisation is needed because all embeddings already live on the unit
+   torus (‖embedding‖ = 1 by construction).
+
+3. **Zero dependencies** – implemented entirely in the Python standard library
+   (``math``, ``cmath``, ``struct``).
+
+Architecture
+------------
+``UCNEmbedding`` uses ``EpicycleDecomposition`` under the hood:
+
+    input data  →  real-valued signal  →  FFT  →  phases  =  embedding
+
+The embedding dimension is always a power of two (the next power-of-two ≥
+*dim*) because the FFT requires it.  Extra dimensions are zeroed out.
+
+Supported input types
+---------------------
+* ``float`` / ``int`` – single-element signal.
+* ``str`` – ordinal encoding of Unicode code points.
+* ``list[float]`` / ``tuple[float]`` – arbitrary real-valued signal.
+* ``bytes`` – unsigned byte values as signal.
+"""
+
+from __future__ import annotations
+
+import math
+import struct
+from typing import Union
+
+from .epicycle import EpicycleDecomposition, _next_pow2
+
+__all__ = ["UCNEmbedding"]
+
+_TAU = 2.0 * math.pi
+
+# Type accepted by UCNEmbedding.encode
+Encodable = Union[int, float, str, bytes, list, tuple]
+
+
+class UCNEmbedding:
+    """Generate and compare Unit Circle Number System embeddings.
+
+    Parameters
+    ----------
+    dim:
+        Desired embedding dimension.  The actual dimension used is the next
+        power of two ≥ *dim* (because the FFT requires it).
+
+    Examples
+    --------
+    >>> emb = UCNEmbedding(dim=16)
+    >>> v1 = emb.encode("hello")
+    >>> v2 = emb.encode("hello")
+    >>> emb.similarity(v1, v2)
+    1.0
+    >>> v3 = emb.encode("world")
+    >>> -1.0 <= emb.similarity(v1, v3) <= 1.0
+    True
+    """
+
+    def __init__(self, dim: int = 64) -> None:
+        if dim < 1:
+            raise ValueError("dim must be at least 1")
+        self._dim_requested: int = dim
+        self._dim: int = _next_pow2(dim)
+
+    @property
+    def dim(self) -> int:
+        """Actual embedding dimension (next power of two ≥ the requested dim)."""
+        return self._dim
+
+    # ------------------------------------------------------------------
+    # Encoding
+    # ------------------------------------------------------------------
+
+    def encode(self, data: Encodable) -> list[float]:
+        """Encode *data* as a UCNS embedding vector.
+
+        Returns a list of ``dim`` angles in ``[0, τ)``.  Identical data always
+        produces identical embeddings.
+
+        Parameters
+        ----------
+        data:
+            Input to encode.  See module docstring for supported types.
+        """
+        signal = self._to_signal(data)
+        # Pad / truncate to self._dim
+        if len(signal) < self._dim:
+            signal = signal + [0.0] * (self._dim - len(signal))
+        else:
+            signal = signal[: self._dim]
+        decomp = EpicycleDecomposition(signal)
+        return decomp.phase_vector
+
+    def encode_packed(self, data: Encodable) -> bytes:
+        """Encode and immediately serialise to compact ``uint16`` bytes.
+
+        Each of the ``dim`` angles is stored as a 16-bit unsigned integer,
+        giving ``2 * dim`` bytes total (vs. ``4 * dim`` for float32).
+        """
+        phases = self.encode(data)
+        scale = 65535.0 / _TAU
+        ints = [min(65535, int(p * scale)) for p in phases]
+        return struct.pack(f"<{len(ints)}H", *ints)
+
+    @staticmethod
+    def unpack(data: bytes) -> list[float]:
+        """Unpack ``uint16`` bytes back to a list of float angles."""
+        n = len(data) // 2
+        ints = struct.unpack(f"<{n}H", data)
+        scale = _TAU / 65535.0
+        return [v * scale for v in ints]
+
+    # ------------------------------------------------------------------
+    # Similarity
+    # ------------------------------------------------------------------
+
+    def similarity(self, a: list[float], b: list[float]) -> float:
+        """Mean phase-cosine similarity between two embeddings ∈ [−1, 1].
+
+        This is the canonical UCNS inner product:
+
+            sim(a, b) = (1/dim) · Σᵢ cos(aᵢ − bᵢ)
+
+        All embedding vectors have unit "norm" under this metric, so the
+        result is a pure cosine without any length normalisation step.
+        """
+        if len(a) != len(b):
+            raise ValueError(
+                f"Embeddings must have equal length; got {len(a)} and {len(b)}"
+            )
+        if not a:
+            return 0.0
+        return sum(math.cos(ai - bi) for ai, bi in zip(a, b)) / len(a)
+
+    def nearest(
+        self,
+        query: list[float],
+        corpus: list[list[float]],
+    ) -> tuple[int, float]:
+        """Find the index and score of the most similar embedding in *corpus*.
+
+        Parameters
+        ----------
+        query:
+            Query embedding (list of angles).
+        corpus:
+            List of candidate embeddings to compare against.
+
+        Returns
+        -------
+        (index, score):
+            Index of the best match and its similarity score in ``[−1, 1]``.
+        """
+        if not corpus:
+            raise ValueError("corpus is empty")
+        best_idx = 0
+        best_score = self.similarity(query, corpus[0])
+        for i, candidate in enumerate(corpus[1:], start=1):
+            score = self.similarity(query, candidate)
+            if score > best_score:
+                best_score = score
+                best_idx = i
+        return best_idx, best_score
+
+    # ------------------------------------------------------------------
+    # Internal helpers
+    # ------------------------------------------------------------------
+
+    @staticmethod
+    def _to_signal(data: Encodable) -> list[float]:
+        """Convert supported input types to a list of floats."""
+        if isinstance(data, (int, float)):
+            return [float(data)]
+        if isinstance(data, str):
+            return [float(ord(c)) for c in data]
+        if isinstance(data, bytes):
+            return [float(b) for b in data]
+        if isinstance(data, (list, tuple)):
+            return [float(x) for x in data]
+        raise TypeError(
+            f"Unsupported type {type(data).__name__!r}. "
+            "Expected int, float, str, bytes, list, or tuple."
+        )
+
+    # ------------------------------------------------------------------
+    # Dunder
+    # ------------------------------------------------------------------
+
+    def __repr__(self) -> str:
+        return (
+            f"UCNEmbedding(dim={self._dim}, "
+            f"bytes_per_embedding={self._dim * 2})"
+        )
diff --git a/ucns/epicycle.py b/ucns/epicycle.py
new file mode 100644
index 00000000..5ba24ca5
--- /dev/null
+++ b/ucns/epicycle.py
@@ -0,0 +1,278 @@
+"""
+ucns.epicycle
+=============
+Epicycle decomposition via the Fast Fourier Transform (FFT).
+
+**What are epicycles?**
+In the Ptolemaic model of the solar system, planets move on small circles
+(epicycles) whose centres trace larger circles.  Mathematically this is
+equivalent to a finite Fourier series: any smooth periodic signal can be
+written as a superposition of circular motions at integer multiples of a
+fundamental frequency.
+
+**Why epicycles for UCNS?**
+The unit circle supports exactly this structure.  A signal of length *n* is
+fully characterised by *n* complex Fourier coefficients.  Each coefficient
+describes *one epicycle*:
+
+    amplitude  = radius of the epicycle
+    frequency  = how many full turns per signal period
+    phase      = UCN angle at which that epicycle starts
+
+Storing only the *phases* (and discarding amplitude information) gives a
+compact angular "fingerprint" of the signal.  Including amplitudes allows
+exact reconstruction (the full DFT is invertible).
+
+**Efficiency**
+This module implements a pure-Python Cooley–Tukey radix-2 FFT so that large
+signals are handled in O(n log n) time rather than the O(n²) naive DFT.
+Inputs whose length is not a power of two are zero-padded automatically.
+"""
+
+from __future__ import annotations
+
+import cmath
+import math
+import struct
+
+__all__ = [
+    "fft",
+    "ifft",
+    "EpicycleDecomposition",
+]
+
+_TAU = 2.0 * math.pi
+
+
+# ------------------------------------------------------------------
+# Low-level FFT (radix-2 Cooley–Tukey, in-place)
+# ------------------------------------------------------------------
+
+
+def _next_pow2(n: int) -> int:
+    """Return the smallest power of two ≥ *n*."""
+    p = 1
+    while p < n:
+        p <<= 1
+    return p
+
+
+def fft(signal: list[float | complex]) -> list[complex]:
+    """Compute the 1-D Discrete Fourier Transform using the Cooley–Tukey
+    radix-2 FFT algorithm.
+
+    Parameters
+    ----------
+    signal:
+        Sequence of real or complex samples.  If its length is not a power of
+        two it is **zero-padded** to the next power of two.
+
+    Returns
+    -------
+    list[complex]
+        Length-*N* list of complex frequency-domain coefficients where
+        ``N = next_pow2(len(signal))``.  The *k*-th entry equals
+
+            X[k] = Σ_{j=0}^{N-1}  x[j] · e^{−2πi·j·k/N}.
+    """
+    n = _next_pow2(len(signal))
+    x: list[complex] = [complex(v) for v in signal]
+    x += [0j] * (n - len(x))  # zero-pad
+    _fft_inplace(x, inverse=False)
+    return x
+
+
+def ifft(spectrum: list[complex]) -> list[complex]:
+    """Compute the inverse DFT (normalised so that ``ifft(fft(x)) ≈ x``).
+
+    Parameters
+    ----------
+    spectrum:
+        Length-*N* sequence of complex frequency coefficients
+        (``N`` must be a power of two).
+
+    Returns
+    -------
+    list[complex]
+        Length-*N* list of complex time-domain samples.
+    """
+    n = len(spectrum)
+    if n == 0:
+        return []
+    x = [complex(v) for v in spectrum]
+    _fft_inplace(x, inverse=True)
+    inv_n = 1.0 / n
+    return [v * inv_n for v in x]
+
+
+def _fft_inplace(x: list[complex], *, inverse: bool) -> None:
+    """Cooley–Tukey iterative FFT (bit-reversal permutation + butterfly).
+
+    Operates in-place on list *x* whose length **must** be a power of two.
+    """
+    n = len(x)
+    if n <= 1:
+        return
+
+    # Bit-reversal permutation
+    j = 0
+    for i in range(1, n):
+        bit = n >> 1
+        while j & bit:
+            j ^= bit
+            bit >>= 1
+        j ^= bit
+        if i < j:
+            x[i], x[j] = x[j], x[i]
+
+    # Butterfly passes
+    sign = 1.0 if inverse else -1.0
+    length = 2
+    while length <= n:
+        half = length >> 1
+        w_n = cmath.exp(1j * sign * _TAU / length)
+        for i in range(0, n, length):
+            w = 1.0 + 0j
+            for k in range(half):
+                u = x[i + k]
+                v = x[i + k + half] * w
+                x[i + k] = u + v
+                x[i + k + half] = u - v
+                w *= w_n
+        length <<= 1
+
+
+# ------------------------------------------------------------------
+# High-level epicycle decomposition
+# ------------------------------------------------------------------
+
+
+class EpicycleDecomposition:
+    """Represent a real signal as a set of weighted unit-circle rotations.
+
+    After construction, the signal is fully described by three parallel arrays:
+
+    * ``amplitudes[k]`` – radius of the *k*-th epicycle (≥ 0).
+    * ``phases[k]``     – UCN angle of the *k*-th epicycle ∈ [0, τ).
+    * ``frequencies[k]`` – integer frequency index (0, 1, …, N−1).
+
+    Parameters
+    ----------
+    signal:
+        1-D sequence of real (or complex) numbers.
+    """
+
+    __slots__ = ("_n_orig", "_n", "amplitudes", "phases", "frequencies")
+
+    def __init__(self, signal: list[float | complex]) -> None:
+        if not signal:
+            raise ValueError("signal must be non-empty")
+        self._n_orig: int = len(signal)
+        spectrum = fft(signal)
+        n = len(spectrum)
+        self._n: int = n
+        self.frequencies: list[int] = list(range(n))
+        self.amplitudes: list[float] = [abs(c) / n for c in spectrum]
+        self.phases: list[float] = [
+            cmath.phase(c) % _TAU if abs(c) > 1e-12 else 0.0
+            for c in spectrum
+        ]
+
+    # ------------------------------------------------------------------
+    # Properties
+    # ------------------------------------------------------------------
+
+    @property
+    def n(self) -> int:
+        """Length of the (possibly padded) transform."""
+        return self._n
+
+    @property
+    def phase_vector(self) -> list[float]:
+        """All epicycle phases as a flat list of angles in ``[0, τ)``."""
+        return list(self.phases)
+
+    @property
+    def dominant_frequency(self) -> int:
+        """Index of the epicycle with the largest amplitude."""
+        return max(range(self._n), key=lambda k: self.amplitudes[k])
+
+    # ------------------------------------------------------------------
+    # Reconstruction
+    # ------------------------------------------------------------------
+
+    def reconstruct(self) -> list[float]:
+        """Reconstruct the original signal (first ``n_orig`` samples).
+
+        Reconstruction is lossless when the input length was already a power
+        of two; otherwise only the first ``n_orig`` values are meaningful.
+        """
+        n = self._n
+        spectrum = [
+            self.amplitudes[k] * n * cmath.exp(1j * self.phases[k])
+            for k in range(n)
+        ]
+        samples = ifft(spectrum)
+        return [s.real for s in samples[: self._n_orig]]
+
+    # ------------------------------------------------------------------
+    # Similarity
+    # ------------------------------------------------------------------
+
+    def phase_similarity(self, other: "EpicycleDecomposition") -> float:
+        """Amplitude-weighted phase-cosine similarity ∈ [−1, 1].
+
+        Computes
+
+            sim = Σ_k  (A_k · B_k · cos(φ_k − ψ_k))  /  (‖A‖₂ · ‖B‖₂)
+
+        where *A_k*, *B_k* are the amplitudes and *φ_k*, *ψ_k* are the phases
+        of the two decompositions.  Uses only the shared frequency bands when
+        the two transforms have different lengths.
+        """
+        n = min(self._n, other._n)
+        total = sum(
+            self.amplitudes[k] * other.amplitudes[k]
+            * math.cos(self.phases[k] - other.phases[k])
+            for k in range(n)
+        )
+        norm_a = math.sqrt(sum(a * a for a in self.amplitudes[:n])) or 1.0
+        norm_b = math.sqrt(sum(b * b for b in other.amplitudes[:n])) or 1.0
+        return total / (norm_a * norm_b)
+
+    # ------------------------------------------------------------------
+    # Compact serialisation
+    # ------------------------------------------------------------------
+
+    def pack(self) -> bytes:
+        """Serialise phases as 16-bit unsigned integers (2 bytes each).
+
+        Provides ~0.0001 rad resolution with 2 bytes per dimension, giving a
+        **2× compression** vs. 32-bit floats for the same angular data.
+        """
+        scale = 65535.0 / _TAU
+        ints = [min(65535, int(p * scale)) for p in self.phases]
+        return struct.pack(f"<{len(ints)}H", *ints)
+
+    @classmethod
+    def unpack_phases(cls, data: bytes) -> list[float]:
+        """Deserialise 16-bit integer phases back to floats in ``[0, τ)``."""
+        n = len(data) // 2
+        ints = struct.unpack(f"<{n}H", data)
+        scale = _TAU / 65535.0
+        return [v * scale for v in ints]
+
+    # ------------------------------------------------------------------
+    # Dunder
+    # ------------------------------------------------------------------
+
+    def __repr__(self) -> str:
+        dom = self.dominant_frequency
+        return (
+            f"EpicycleDecomposition(n={self._n}, "
+            f"dominant_freq={dom}, "
+            f"dominant_amplitude={self.amplitudes[dom]:.4f})"
+        )
+
+    def __len__(self) -> int:
+        return self._n
diff --git a/ucns/mobius.py b/ucns/mobius.py
new file mode 100644
index 00000000..7bdcf5b0
--- /dev/null
+++ b/ucns/mobius.py
@@ -0,0 +1,191 @@
+"""
+ucns.mobius
+===========
+Möbius (bilinear) transformations of the unit disk.
+
+The **Poincaré disk model** represents the hyperbolic plane as the open unit
+disk  D = {z ∈ ℂ : |z| < 1}.  Its boundary ∂D is the unit circle – the home
+of every UCN.  Conformal automorphisms of D are Möbius transformations of the
+form
+
+    T_{a,φ}(z) = e^(iφ) · (z − a) / (1 − ā·z),   a ∈ D, φ ∈ ℝ.
+
+These transformations:
+
+* **preserve** the unit circle (boundary maps to boundary),
+* **preserve** the hyperbolic (Poincaré) metric,
+* compose to form the group Aut(D) ≅ PU(1,1).
+
+Geometric intuition
+-------------------
+Think of the disk as a "rubber sheet" that can be stretched or compressed while
+keeping the circular boundary fixed.  Embedding data on the interior of the
+disk naturally encodes *hierarchical* relationships: nearby points are close in
+hyperbolic distance; points near the boundary are conceptually "far out" (low
+frequency, coarse-grained).  Combining this with the recursive epicycle
+structure (see ``ucns.epicycle``) yields a multi-scale embedding space that
+is simultaneously compact (unit circle) and hierarchical (Möbius disk).
+
+References
+----------
+* Poincaré disk model – Wikipedia
+* "Poincaré Embeddings for Learning Hierarchical Representations" – Nickel & Kiela 2017
+"""
+
+from __future__ import annotations
+
+import cmath
+import math
+
+__all__ = ["MobiusTransform", "poincare_distance", "disk_to_circle", "circle_to_disk"]
+
+_TAU = 2.0 * math.pi
+
+
+class MobiusTransform:
+    """A conformal automorphism of the open unit disk.
+
+    Parameters
+    ----------
+    a:
+        Translation parameter.  Must satisfy ``|a| < 1`` (interior of disk).
+        The transformation maps ``a ↦ 0``.
+    phi:
+        Rotation angle in radians.  Applied after the translation.
+    """
+
+    __slots__ = ("a", "phi")
+
+    def __init__(self, a: complex, phi: float = 0.0) -> None:
+        if abs(a) >= 1.0:
+            raise ValueError(
+                f"|a| must be strictly less than 1; got |a| = {abs(a):.6f}"
+            )
+        self.a: complex = complex(a)
+        self.phi: float = float(phi) % _TAU
+
+    # ------------------------------------------------------------------
+    # Evaluation
+    # ------------------------------------------------------------------
+
+    def __call__(self, z: complex) -> complex:
+        """Apply the transform: T(z) = e^(iφ) · (z − a) / (1 − ā·z)."""
+        denom = 1.0 - self.a.conjugate() * z
+        if abs(denom) < 1e-15:
+            raise ValueError("z is the image of infinity under this transform")
+        return cmath.exp(1j * self.phi) * (z - self.a) / denom
+
+    # ------------------------------------------------------------------
+    # Group structure
+    # ------------------------------------------------------------------
+
+    def inverse(self) -> "MobiusTransform":
+        """Return T⁻¹ such that T⁻¹(T(z)) = z for all z ∈ D."""
+        return MobiusTransform(
+            -self.a * cmath.exp(1j * self.phi),
+            -self.phi,
+        )
+
+    def compose(self, other: "MobiusTransform") -> "MobiusTransform":
+        """Return the composed transform self ∘ other (apply *other* first)."""
+        # Numerically stable composition via a sample point
+        # We derive the new (a, phi) from where each transform sends 0.
+        # T_composed sends 0 → self(other(0))
+        a_new = self(other(0j))
+        # Determine rotation by evaluating at a second point
+        p = other(0.5 + 0j)
+        q = self(p)
+        # q = e^(i*phi_new) * (q_unnorm) → phi_new = arg(q / T_a_new(a_sample))
+        t_check = MobiusTransform(a_new)
+        sample = 0.5 + 0j
+        out = q
+        raw = (out - a_new) / (1.0 - a_new.conjugate() * out)
+        phi_new = cmath.phase(raw) if abs(raw) > 1e-15 else 0.0
+        return MobiusTransform(a_new, phi_new)
+
+    # ------------------------------------------------------------------
+    # Hyperbolic geometry
+    # ------------------------------------------------------------------
+
+    def hyperbolic_distance(self, z: complex, w: complex) -> float:
+        """Poincaré disk metric  d(z, w) = 2·arctanh(|T_z(w)|).
+
+        This is the intrinsic distance in the hyperbolic plane modelled by D.
+        """
+        return poincare_distance(z, w)
+
+    # ------------------------------------------------------------------
+    # Dunder
+    # ------------------------------------------------------------------
+
+    def __repr__(self) -> str:
+        return f"MobiusTransform(a={self.a:.4f}, phi={self.phi:.4f})"
+
+    def __eq__(self, other: object) -> bool:
+        if not isinstance(other, MobiusTransform):
+            return NotImplemented
+        return (
+            abs(self.a - other.a) < 1e-9
+            and abs(self.phi - other.phi) < 1e-9
+        )
+
+
+# ------------------------------------------------------------------
+# Standalone geometric helpers
+# ------------------------------------------------------------------
+
+
+def poincare_distance(z: complex, w: complex) -> float:
+    """Hyperbolic distance between two points in the Poincaré disk.
+
+    Parameters
+    ----------
+    z, w:
+        Points in the open unit disk (``|z|, |w| < 1``).
+
+    Returns
+    -------
+    float
+        ``d(z, w) = 2·arctanh(|(z − w) / (1 − w̄·z)|)`` ≥ 0.
+    """
+    for name, p in (("z", z), ("w", w)):
+        if abs(p) >= 1.0:
+            raise ValueError(
+                f"Point {name} = {p} lies outside or on the unit disk boundary"
+            )
+    denom = 1.0 - w.conjugate() * z
+    if abs(denom) < 1e-15:
+        return float("inf")
+    rho = abs((z - w) / denom)
+    rho = min(rho, 1.0 - 1e-15)  # guard against numerical overshoot
+    return 2.0 * math.atanh(rho)
+
+
+def disk_to_circle(z: complex) -> float:
+    """Project a point *z* in the unit disk to its angle on ∂D (the unit circle).
+
+    Uses the Cayley-like radial projection z ↦ z/|z| for z ≠ 0; for z = 0
+    returns 0.  The returned angle θ ∈ [0, τ) gives the UCN associated with
+    the boundary limit of the radial ray through *z*.
+    """
+    if abs(z) < 1e-15:
+        return 0.0
+    return cmath.phase(z) % _TAU
+
+
+def circle_to_disk(theta: float, r: float = 0.5) -> complex:
+    """Embed a unit-circle point (angle) into the interior of the disk at radius *r*.
+
+    This is useful when you want to treat UCN angles as interior hyperbolic
+    points (r < 1 keeps them strictly inside D).
+
+    Parameters
+    ----------
+    theta:
+        Angle on the unit circle (radians).
+    r:
+        Radial depth in (0, 1).  Defaults to 0.5.
+    """
+    if not (0.0 < r < 1.0):
+        raise ValueError("r must be in (0, 1)")
+    return r * cmath.exp(1j * theta)
diff --git a/ucns/similarity.py b/ucns/similarity.py
new file mode 100644
index 00000000..b410af48
--- /dev/null
+++ b/ucns/similarity.py
@@ -0,0 +1,181 @@
+"""
+ucns.similarity
+===============
+Similarity and distance metrics for Unit Circle Number embeddings.
+
+All functions operate on plain Python lists of angles (floats in [0, τ)) as
+returned by ``UCNEmbedding.encode``.  No external libraries required.
+
+Metric catalogue
+----------------
+``phase_cosine``
+    Mean of cos(θᵢ − φᵢ) – the natural "dot product" on the torus.
+    Range: [−1, 1].  Value 1 means identical embeddings.
+
+``arc_distance``
+    Mean minimum arc length |θᵢ − φᵢ|_circle, normalised to [0, 1].
+    Value 0 means identical; value 1 means diametrically opposite.
+
+``hyperbolic_cosine``
+    Uses the Poincaré disk: maps each angle to an interior point of the
+    unit disk and computes hyperbolic cosine similarity.  Sensitive to
+    hierarchical structure (high-frequency vs low-frequency components).
+
+``top_k_overlap``
+    Jaccard-like overlap of the *k* dominant frequency indices.  Good for
+    sparse or categorical data.
+"""
+
+from __future__ import annotations
+
+import math
+import cmath
+from typing import Sequence
+
+__all__ = [
+    "phase_cosine",
+    "arc_distance",
+    "hyperbolic_cosine",
+    "top_k_overlap",
+]
+
+_TAU = 2.0 * math.pi
+
+
+def _check_same_length(a: Sequence[float], b: Sequence[float]) -> int:
+    if len(a) != len(b):
+        raise ValueError(
+            f"Embeddings must have the same length; got {len(a)} and {len(b)}"
+        )
+    return len(a)
+
+
+def phase_cosine(a: Sequence[float], b: Sequence[float]) -> float:
+    """Mean angular cosine similarity between two UCNS embeddings.
+
+    Parameters
+    ----------
+    a, b:
+        Lists of angles (radians) of equal length.
+
+    Returns
+    -------
+    float
+        Value in ``[−1, 1]``.  A value of 1 indicates identical phase
+        patterns; −1 indicates perfectly anti-phase patterns.
+    """
+    n = _check_same_length(a, b)
+    if n == 0:
+        return 0.0
+    return sum(math.cos(ai - bi) for ai, bi in zip(a, b)) / n
+
+
+def arc_distance(a: Sequence[float], b: Sequence[float]) -> float:
+    """Mean normalised arc distance between two UCNS embeddings.
+
+    Each per-dimension distance is the shorter arc between the two angles,
+    normalised by π so the result lies in ``[0, 1]``.
+
+    Parameters
+    ----------
+    a, b:
+        Lists of angles (radians) of equal length.
+
+    Returns
+    -------
+    float
+        Value in ``[0, 1]``.  0 means identical; 1 means maximally different.
+    """
+    n = _check_same_length(a, b)
+    if n == 0:
+        return 0.0
+    total = 0.0
+    for ai, bi in zip(a, b):
+        diff = abs(ai - bi) % _TAU
+        total += min(diff, _TAU - diff)
+    return total / (n * math.pi)
+
+
+def hyperbolic_cosine(
+    a: Sequence[float],
+    b: Sequence[float],
+    *,
+    radius: float = 0.5,
+) -> float:
+    """Similarity via per-dimension hyperbolic cosine in the Poincaré disk.
+
+    Each angle θ is embedded at ``r·e^(iθ)`` inside the unit disk.  The
+    hyperbolic distance between the two disk points is converted to a cosine:
+
+        sim = mean_k( cos( d_hyp(r·e^{iθ_k}, r·e^{iφ_k}) ) )
+
+    This metric is more sensitive to *low-frequency* (large-amplitude)
+    components than ``phase_cosine`` and captures hierarchical relationships.
+
+    Parameters
+    ----------
+    a, b:
+        Lists of angles (radians) of equal length.
+    radius:
+        Radial depth in ``(0, 1)`` for the disk embedding.  Smaller values
+        compress the hyperbolic scale; larger values expand it.
+
+    Returns
+    -------
+    float
+        Value in ``[−1, 1]``.
+    """
+    if not (0.0 < radius < 1.0):
+        raise ValueError("radius must be in (0, 1)")
+    n = _check_same_length(a, b)
+    if n == 0:
+        return 0.0
+    total = 0.0
+    for ai, bi in zip(a, b):
+        za = radius * cmath.exp(1j * ai)
+        zb = radius * cmath.exp(1j * bi)
+        denom = 1.0 - zb.conjugate() * za
+        if abs(denom) < 1e-15:
+            d = 0.0
+        else:
+            rho = abs((za - zb) / denom)
+            rho = min(rho, 1.0 - 1e-15)
+            d = 2.0 * math.atanh(rho)
+        total += math.cos(d)
+    return total / n
+
+
+def top_k_overlap(
+    amplitudes_a: Sequence[float],
+    amplitudes_b: Sequence[float],
+    *,
+    k: int = 8,
+) -> float:
+    """Jaccard-like overlap of the *k* most energetic frequency components.
+
+    Parameters
+    ----------
+    amplitudes_a, amplitudes_b:
+        Amplitude arrays (as from ``EpicycleDecomposition.amplitudes``).
+    k:
+        Number of top components to compare.
+
+    Returns
+    -------
+    float
+        Value in ``[0, 1]``.  1 means the top-*k* frequency sets are
+        identical; 0 means completely disjoint.
+    """
+    _check_same_length(amplitudes_a, amplitudes_b)
+    k = max(1, min(k, len(amplitudes_a)))
+
+    def top_k_indices(amps: Sequence[float]) -> set:
+        return set(
+            sorted(range(len(amps)), key=lambda i: amps[i], reverse=True)[:k]
+        )
+
+    sa = top_k_indices(amplitudes_a)
+    sb = top_k_indices(amplitudes_b)
+    intersection = len(sa & sb)
+    union = len(sa | sb)
+    return intersection / union if union else 0.0