diff --git a/include/vde/sdf/sdf_operations.h b/include/vde/sdf/sdf_operations.h
new file mode 100644
index 0000000..919d594
--- /dev/null
+++ b/include/vde/sdf/sdf_operations.h
@@ -0,0 +1,165 @@
+#pragma once
+#include <cmath>
+#include <algorithm>
+#include "vde/core/point.h"
+
+namespace vde::sdf {
+
+// ═══════════════════════════════════════════════════
+//  Boolean Operations (sharp)
+// ═══════════════════════════════════════════════════
+
+/// Union (min) — combine two shapes
+[[nodiscard]] inline double op_union(double d1, double d2) {
+    return std::min(d1, d2);
+}
+
+/// Intersection (max) — region common to both shapes
+[[nodiscard]] inline double op_intersection(double d1, double d2) {
+    return std::max(d1, d2);
+}
+
+/// Difference — subtract d2 from d1 (d1 \ d2)
+[[nodiscard]] inline double op_difference(double d1, double d2) {
+    return std::max(-d1, d2);
+}
+
+// ═══════════════════════════════════════════════════
+//  Smooth Boolean Operations (Inigo Quilez)
+// ═══════════════════════════════════════════════════
+
+namespace detail {
+    template <typename T>
+    [[nodiscard]] constexpr T mix(T a, T b, T t) noexcept {
+        return a + t * (b - a);
+    }
+}
+
+/// Smooth union with blend radius k
+/// Falls back to sharp union when k <= 0
+[[nodiscard]] inline double op_smooth_union(double d1, double d2, double k) {
+    if (k <= 0.0) return op_union(d1, d2);
+    double h = std::clamp(0.5 + 0.5 * (d2 - d1) / k, 0.0, 1.0);
+    return detail::mix(d2, d1, h) - k * h * (1.0 - h);
+}
+
+/// Smooth intersection with blend radius k
+/// Falls back to sharp intersection when k <= 0
+[[nodiscard]] inline double op_smooth_intersection(double d1, double d2, double k) {
+    if (k <= 0.0) return op_intersection(d1, d2);
+    double h = std::clamp(0.5 - 0.5 * (d2 - d1) / k, 0.0, 1.0);
+    return detail::mix(d1, d2, h) + k * h * (1.0 - h);
+}
+
+/// Smooth difference with blend radius k (d1 \ d2)
+/// Falls back to sharp difference when k <= 0
+[[nodiscard]] inline double op_smooth_difference(double d1, double d2, double k) {
+    if (k <= 0.0) return op_difference(d1, d2);
+    double h = std::clamp(0.5 - 0.5 * (d2 + d1) / k, 0.0, 1.0);
+    return detail::mix(d2, -d1, h) + k * h * (1.0 - h);
+}
+
+// ═══════════════════════════════════════════════════
+//  Shape Modifiers
+// ═══════════════════════════════════════════════════
+
+/// Round a shape — shrink-wraps by radius r, then expands
+/// Equivalent to: shape_distance - r
+[[nodiscard]] inline double op_round(double d, double r) {
+    return d - r;
+}
+
+/// Onion — create a thin shell of thickness t
+/// Returns SDF of the shell region (t/2 inward, t/2 outward from surface)
+[[nodiscard]] inline double op_onion(double d, double thickness) {
+    return std::abs(d) - thickness * 0.5;
+}
+
+// ═══════════════════════════════════════════════════
+//  Domain Deformations (coordinate remapping)
+//  Apply BEFORE evaluating the SDF primitive.
+// ═══════════════════════════════════════════════════
+
+/// Infinite repetition along each axis
+/// c.x/y/z: cell size per axis (0 = no repeat along that axis)
+[[nodiscard]] inline vde::core::Point3D op_repeat(const vde::core::Point3D& p, const vde::core::Point3D& c) {
+    // p - c * round(p / c)
+    // For axes where c_i == 0, skip repetition
+    auto q = p;
+    if (c.x() > 0.0) q.x() = p.x() - c.x() * std::round(p.x() / c.x());
+    if (c.y() > 0.0) q.y() = p.y() - c.y() * std::round(p.y() / c.y());
+    if (c.z() > 0.0) q.z() = p.z() - c.z() * std::round(p.z() / c.z());
+    return q;
+}
+
+/// Mirror across X=0 plane (folds everything to x >= offset)
+[[nodiscard]] inline vde::core::Point3D op_mirror_x(const vde::core::Point3D& p, double offset = 0.0) {
+    return {offset + std::abs(p.x() - offset), p.y(), p.z()};
+}
+
+/// Mirror across Y=0 plane (folds everything to y >= offset)
+[[nodiscard]] inline vde::core::Point3D op_mirror_y(const vde::core::Point3D& p, double offset = 0.0) {
+    return {p.x(), offset + std::abs(p.y() - offset), p.z()};
+}
+
+/// Mirror across Z=0 plane (folds everything to z >= offset)
+[[nodiscard]] inline vde::core::Point3D op_mirror_z(const vde::core::Point3D& p, double offset = 0.0) {
+    return {p.x(), p.y(), offset + std::abs(p.z() - offset)};
+}
+
+/// Rotate point around Y axis by angle (radians)
+[[nodiscard]] inline vde::core::Point3D op_rotate(const vde::core::Point3D& p, double angle_rad) {
+    double s = std::sin(angle_rad);
+    double c = std::cos(angle_rad);
+    return {c * p.x() + s * p.z(), p.y(), -s * p.x() + c * p.z()};
+}
+
+/// Translate point
+[[nodiscard]] inline vde::core::Point3D op_translate(const vde::core::Point3D& p, const vde::core::Point3D& offset) {
+    return {p.x() - offset.x(), p.y() - offset.y(), p.z() - offset.z()};
+}
+
+/// Non-uniform scale — divides coordinates (non-exact SDF scaling)
+/// Use with care: the resulting distance field is approximate.
+[[nodiscard]] inline vde::core::Point3D op_scale(const vde::core::Point3D& p, const vde::core::Point3D& s) {
+    return {p.x() / s.x(), p.y() / s.y(), p.z() / s.z()};
+}
+
+/// Twist around Y axis — rotation amount increases with y coordinate
+[[nodiscard]] inline vde::core::Point3D op_twist(const vde::core::Point3D& p, double amount) {
+    double angle = amount * p.y();
+    double s = std::sin(angle);
+    double c = std::cos(angle);
+    return {c * p.x() + s * p.z(), p.y(), -s * p.x() + c * p.z()};
+}
+
+/// Bend along Y axis — curves the shape in x-y plane
+[[nodiscard]] inline vde::core::Point3D op_bend(const vde::core::Point3D& p, double amount) {
+    double c = std::cos(amount * p.y());
+    double s = std::sin(amount * p.y());
+    return {c * p.x() - s * p.y(), s * p.x() + c * p.y(), p.z()};
+}
+
+/// Elongate along axes — stretches coordinates (divides by factor)
+/// Equivalent to non-uniform scale but conceptually "stretches" the shape.
+[[nodiscard]] inline vde::core::Point3D op_elongate(const vde::core::Point3D& p, const vde::core::Point3D& factors) {
+    return {p.x() / factors.x(), p.y() / factors.y(), p.z() / factors.z()};
+}
+
+/// Displacement — adds sinusoidal perturbation to the SDF value
+/// p: current point, d: base SDF value, amplitude/frequency: displacement params
+[[nodiscard]] inline double op_displace(double d, const vde::core::Point3D& p, double amplitude, double frequency) {
+    return d + amplitude * std::sin(p.x() * frequency)
+                           * std::sin(p.y() * frequency)
+                           * std::sin(p.z() * frequency);
+}
+
+/// Cheap bend — bends around Z axis based on x coordinate
+/// Simpler and faster approximation than full bend.
+[[nodiscard]] inline vde::core::Point3D op_cheap_bend(const vde::core::Point3D& p, double k) {
+    double c = std::cos(k * p.x());
+    double s = std::sin(k * p.x());
+    return {c * p.x() - s * p.y(), s * p.x() + c * p.y(), p.z()};
+}
+
+} // namespace vde::sdf
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index bd352bf..c3a522e 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -106,6 +106,19 @@ target_link_libraries(vde_boolean
     PUBLIC  vde_mesh vde_spatial vde_compile_options
 )
 
+# ── sdf ────────────────────────────────────────────
+add_library(vde_sdf STATIC
+    sdf/sdf_tree.cpp
+    sdf/sdf_to_mesh.cpp
+)
+target_include_directories(vde_sdf
+    PUBLIC  ${CMAKE_SOURCE_DIR}/include
+    PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}
+)
+target_link_libraries(vde_sdf
+    PUBLIC  vde_core vde_compile_options
+)
+
 # ── collision ───────────────────────────────────────
 add_library(vde_collision STATIC
     collision/gjk.cpp
@@ -131,6 +144,7 @@ target_link_libraries(vde
     INTERFACE vde_spatial
     INTERFACE vde_boolean
     INTERFACE vde_collision
+    INTERFACE vde_sdf
     INTERFACE vde_brep
     INTERFACE vde_sketch
 )
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 816e8e4..a96febe 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -12,4 +12,6 @@ add_subdirectory(spatial)
 add_subdirectory(collision)
 add_subdirectory(boolean)
 add_subdirectory(brep)
+add_subdirectory(sdf)
 add_subdirectory(sketch)
+add_subdirectory(sdf)
diff --git a/tests/sdf/CMakeLists.txt b/tests/sdf/CMakeLists.txt
new file mode 100644
index 0000000..c171a34
--- /dev/null
+++ b/tests/sdf/CMakeLists.txt
@@ -0,0 +1,9 @@
+function(add_vde_sdf_test name)
+    add_executable(${name} ${name}.cpp)
+    target_link_libraries(${name} PRIVATE vde_sdf GTest::gtest GTest::gtest_main)
+    target_include_directories(${name} PRIVATE ${CMAKE_SOURCE_DIR}/include)
+    gtest_discover_tests(${name})
+endfunction()
+
+add_vde_sdf_test(test_sdf_tree)
+add_vde_sdf_test(test_sdf_to_mesh)
diff --git a/tests/sdf/test_sdf_operations.cpp b/tests/sdf/test_sdf_operations.cpp
new file mode 100644
index 0000000..1ca1401
--- /dev/null
+++ b/tests/sdf/test_sdf_operations.cpp
@@ -0,0 +1,414 @@
+#include <gtest/gtest.h>
+#include <cmath>
+#include "vde/sdf/sdf_operations.h"
+
+using namespace vde::sdf;
+using vde::core::Point3D;
+
+// ── Minimal inline primitives for testing (no sdf_primitives dependency) ──
+
+inline double sdf_sphere(const Point3D& p, double r = 1.0) {
+    return std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z()) - r;
+}
+
+inline double sdf_box(const Point3D& p, const Point3D& b) {
+    auto q = Point3D(std::abs(p.x()) - b.x(),
+                     std::abs(p.y()) - b.y(),
+                     std::abs(p.z()) - b.z());
+    return std::max({q.x(), q.y(), q.z(), 0.0}) +
+           std::min(std::max({q.x(), q.y(), q.z()}), 0.0);
+}
+
+// ═══════════════════════════════════════════════════
+//  Sharp Boolean Operations
+// ═══════════════════════════════════════════════════
+
+TEST(SdfOperations, Union) {
+    EXPECT_DOUBLE_EQ(op_union(1.0, 2.0), 1.0);
+    EXPECT_DOUBLE_EQ(op_union(-1.0, 2.0), -1.0);
+    EXPECT_DOUBLE_EQ(op_union(3.0, -1.0), -1.0);
+    EXPECT_DOUBLE_EQ(op_union(-3.0, -1.0), -3.0);
+}
+
+TEST(SdfOperations, Intersection) {
+    EXPECT_DOUBLE_EQ(op_intersection(1.0, 2.0), 2.0);
+    EXPECT_DOUBLE_EQ(op_intersection(-1.0, 2.0), 2.0);
+    EXPECT_DOUBLE_EQ(op_intersection(3.0, -1.0), 3.0);
+    EXPECT_DOUBLE_EQ(op_intersection(-3.0, -1.0), -1.0);
+}
+
+TEST(SdfOperations, Difference) {
+    // d1 is the body, d2 is the cutter
+    // difference = max(-d1, d2)
+    EXPECT_DOUBLE_EQ(op_difference(1.0, 2.0), 2.0);   // -d1=-1, max(-1,2)=2
+    EXPECT_DOUBLE_EQ(op_difference(3.0, 2.0), 2.0);   // max(-3,2)=2
+    EXPECT_DOUBLE_EQ(op_difference(-1.0, 2.0), 2.0);  // max(1,2)=2
+    EXPECT_DOUBLE_EQ(op_difference(2.0, -3.0), -3.0); // max(-2,-3)=-2
+}
+
+TEST(SdfOperations, UnionIdentity) {
+    // Union with positive infinity yields the other value
+    double inf = std::numeric_limits<double>::infinity();
+    EXPECT_DOUBLE_EQ(op_union(5.0, inf), 5.0);
+    EXPECT_DOUBLE_EQ(op_union(inf, 5.0), 5.0);
+}
+
+TEST(SdfOperations, IntersectionIdentity) {
+    double inf = std::numeric_limits<double>::infinity();
+    EXPECT_DOUBLE_EQ(op_intersection(5.0, -inf), 5.0);
+    EXPECT_DOUBLE_EQ(op_intersection(-inf, 5.0), 5.0);
+}
+
+// ═══════════════════════════════════════════════════
+//  Smooth Boolean Operations
+// ═══════════════════════════════════════════════════
+
+TEST(SdfOperations, SmoothUnionFallsBackToSharp) {
+    // k=0 should fall back to sharp union
+    EXPECT_DOUBLE_EQ(op_smooth_union(1.0, 2.0, 0.0), 1.0);
+    EXPECT_DOUBLE_EQ(op_smooth_union(-1.0, 2.0, 0.0), -1.0);
+    // k negative should also fall back
+    EXPECT_DOUBLE_EQ(op_smooth_union(1.0, 2.0, -5.0), 1.0);
+}
+
+TEST(SdfOperations, SmoothUnionBlending) {
+    // With k > 0, result should be less than min(d1,d2) (blend digs in)
+    double k = 0.1;
+    double result = op_smooth_union(1.0, 1.0, k);
+    EXPECT_LT(result, 1.0);           // blend reduces distance
+    EXPECT_GT(result, 1.0 - k * 0.25); // bound: max reduction is k/4 at h=0.5
+
+    // Far apart — should approach sharp union
+    double far_result = op_smooth_union(1.0, 10.0, 0.1);
+    EXPECT_NEAR(far_result, 1.0, 0.01);
+}
+
+TEST(SdfOperations, SmoothIntersectionFallsBackToSharp) {
+    EXPECT_DOUBLE_EQ(op_smooth_intersection(1.0, 2.0, 0.0), 2.0);
+    EXPECT_DOUBLE_EQ(op_smooth_intersection(1.0, 2.0, -1.0), 2.0);
+}
+
+TEST(SdfOperations, SmoothIntersectionBlending) {
+    double k = 0.1;
+    double result = op_smooth_intersection(1.0, 1.0, k);
+    EXPECT_GT(result, 1.0);          // blend increases distance
+    EXPECT_LT(result, 1.0 + k * 0.25);
+}
+
+TEST(SdfOperations, SmoothDifferenceFallsBackToSharp) {
+    EXPECT_DOUBLE_EQ(op_smooth_difference(1.0, 2.0, 0.0), 2.0);
+    EXPECT_DOUBLE_EQ(op_smooth_difference(-1.0, 2.0, 0.0), 2.0);
+}
+
+// ═══════════════════════════════════════════════════
+//  Shape Modifiers
+// ═══════════════════════════════════════════════════
+
+TEST(SdfOperations, Round) {
+    EXPECT_DOUBLE_EQ(op_round(5.0, 2.0), 3.0);
+    EXPECT_DOUBLE_EQ(op_round(-1.0, 0.5), -1.5);
+    EXPECT_DOUBLE_EQ(op_round(0.0, 3.0), -3.0);
+}
+
+TEST(SdfOperations, Onion) {
+    // Surface at d=0 → onion creates shell distance
+    EXPECT_DOUBLE_EQ(op_onion(0.0, 2.0), -1.0);   // abs(0) - 1 = -1
+    EXPECT_DOUBLE_EQ(op_onion(2.0, 2.0), 1.0);    // abs(2) - 1 = 1
+    EXPECT_DOUBLE_EQ(op_onion(-2.0, 2.0), 1.0);   // abs(-2) - 1 = 1
+    EXPECT_DOUBLE_EQ(op_onion(0.5, 2.0), -0.5);   // abs(0.5) - 1 = -0.5
+}
+
+TEST(SdfOperations, OnionSphere) {
+    // Sphere of radius 1 with thickness 0.2 creates a shell
+    // Surface of sphere: d=0 at r=1
+    // Onion shell: abs(d)-0.1 → interior at r=0.9 through r=1.1
+    double d = sdf_sphere(Point3D(0, 0, 0.9));
+    double shell = op_onion(d, 0.2);
+    EXPECT_NEAR(shell, 0.0, 1e-10);  // inner surface
+
+    d = sdf_sphere(Point3D(0, 0, 1.0));
+    shell = op_onion(d, 0.2);
+    EXPECT_NEAR(shell, -0.1, 1e-10); // inside shell
+}
+
+// ═══════════════════════════════════════════════════
+//  Domain Deformations
+// ═══════════════════════════════════════════════════
+
+TEST(SdfOperations, RepeatIdentity) {
+    // repeat with zero cell size = no repetition
+    auto p = Point3D(5.0, 3.0, -2.0);
+    auto result = op_repeat(p, Point3D(0, 0, 0));
+    EXPECT_DOUBLE_EQ(result.x(), p.x());
+    EXPECT_DOUBLE_EQ(result.y(), p.y());
+    EXPECT_DOUBLE_EQ(result.z(), p.z());
+}
+
+TEST(SdfOperations, RepeatMapsToUnitCell) {
+    // With cell size (2,2,2), all points should map to [-1, 1]
+    auto c = Point3D(2.0, 2.0, 2.0);
+
+    // Center of cell → maps to 0
+    auto r1 = op_repeat(Point3D(0, 0, 0), c);
+    EXPECT_NEAR(r1.x(), 0.0, 1e-10);
+    EXPECT_NEAR(r1.y(), 0.0, 1e-10);
+    EXPECT_NEAR(r1.z(), 0.0, 1e-10);
+
+    // Cell boundary → maps near ±1
+    auto r2 = op_repeat(Point3D(1.0, 1.0, 1.0), c);
+    EXPECT_NEAR(r2.x(), 1.0, 1e-10);
+    EXPECT_NEAR(r2.y(), 1.0, 1e-10);
+    EXPECT_NEAR(r2.z(), 1.0, 1e-10);
+
+    // Outside the cell
+    auto r3 = op_repeat(Point3D(3.5, 2.5, 1.5), c);
+    EXPECT_NEAR(r3.x(), -0.5, 1e-10);  // 3.5 - 2*2 = -0.5
+    EXPECT_NEAR(r3.y(), 0.5, 1e-10);   // 2.5 - 1*2 = 0.5
+    EXPECT_NEAR(r3.z(), -0.5, 1e-10);  // 1.5 - 1*2 = -0.5
+}
+
+TEST(SdfOperations, RepeatSingleAxis) {
+    // Only repeat along X
+    auto c = Point3D(3.0, 0.0, 0.0);
+    auto r = op_repeat(Point3D(5.0, 2.0, 1.0), c);
+    EXPECT_NEAR(r.x(), -1.0, 1e-10);   // 5 - 2*3 = -1
+    EXPECT_DOUBLE_EQ(r.y(), 2.0);       // unchanged
+    EXPECT_DOUBLE_EQ(r.z(), 1.0);       // unchanged
+}
+
+TEST(SdfOperations, RepeatSphereGrid) {
+    // Repeated sphere: check that sphere exists at multiple cell positions
+    auto sphere_repeated = [](const Point3D& p) -> double {
+        return sdf_sphere(op_repeat(p, Point3D(3.0, 3.0, 3.0)), 0.5);
+    };
+
+    // Center of any cell should be inside the sphere
+    EXPECT_LT(sphere_repeated(Point3D(0, 0, 0)), 0.0);
+    EXPECT_LT(sphere_repeated(Point3D(3, 0, 0)), 0.0);
+    EXPECT_LT(sphere_repeated(Point3D(-3, 6, -9)), 0.0);
+
+    // Midpoint between cells should be outside
+    EXPECT_GT(sphere_repeated(Point3D(1.5, 0, 0)), 0.0);
+    EXPECT_GT(sphere_repeated(Point3D(0, 1.5, 0)), 0.0);
+}
+
+TEST(SdfOperations, Mirror) {
+    // Mirror across X=0: everything maps to x >= 0
+    auto r1 = op_mirror_x(Point3D(3.0, 1.0, 1.0));
+    EXPECT_DOUBLE_EQ(r1.x(), 3.0);  // positive stays positive
+    EXPECT_DOUBLE_EQ(r1.y(), 1.0);
+    EXPECT_DOUBLE_EQ(r1.z(), 1.0);
+
+    auto r2 = op_mirror_x(Point3D(-3.0, 2.0, 2.0));
+    EXPECT_DOUBLE_EQ(r2.x(), 3.0);  // negative folds to positive
+    EXPECT_DOUBLE_EQ(r2.y(), 2.0);
+    EXPECT_DOUBLE_EQ(r2.z(), 2.0);
+
+    // With offset
+    auto r3 = op_mirror_x(Point3D(-3.0, 0, 0), 1.0);
+    EXPECT_DOUBLE_EQ(r3.x(), 5.0);  // 1 + abs(-3 - 1) = 1 + 4 = 5
+}
+
+TEST(SdfOperations, RotateY) {
+    // 90° rotation around Y: (1,0,0) → (0,0,-1)
+    auto r1 = op_rotate(Point3D(1, 0, 0), M_PI / 2.0);
+    EXPECT_NEAR(r1.x(), 0.0, 1e-10);
+    EXPECT_NEAR(r1.y(), 0.0, 1e-10);
+    EXPECT_NEAR(r1.z(), -1.0, 1e-10);
+
+    // 180° rotation: (1,0,0) → (-1,0,0)
+    auto r2 = op_rotate(Point3D(1, 0, 0), M_PI);
+    EXPECT_NEAR(r2.x(), -1.0, 1e-10);
+
+    // 360° → identity
+    auto r3 = op_rotate(Point3D(3.0, 2.0, 1.0), 2.0 * M_PI);
+    EXPECT_NEAR(r3.x(), 3.0, 1e-10);
+    EXPECT_NEAR(r3.y(), 2.0, 1e-10);
+    EXPECT_NEAR(r3.z(), 1.0, 1e-10);
+}
+
+TEST(SdfOperations, Translate) {
+    auto r = op_translate(Point3D(5, 3, 1), Point3D(2, 1, 0));
+    EXPECT_DOUBLE_EQ(r.x(), 3.0);  // 5 - 2
+    EXPECT_DOUBLE_EQ(r.y(), 2.0);  // 3 - 1
+    EXPECT_DOUBLE_EQ(r.z(), 1.0);  // 1 - 0
+}
+
+TEST(SdfOperations, Scale) {
+    auto s = Point3D(2.0, 4.0, 0.5);
+    auto r = op_scale(Point3D(4, 8, 1), s);
+    EXPECT_DOUBLE_EQ(r.x(), 2.0);
+    EXPECT_DOUBLE_EQ(r.y(), 2.0);
+    EXPECT_DOUBLE_EQ(r.z(), 2.0);
+}
+
+TEST(SdfOperations, Twist) {
+    // At y=0, twist has no effect (angle = 0)
+    auto r1 = op_twist(Point3D(1, 0, 0), 1.0);
+    EXPECT_NEAR(r1.x(), 1.0, 1e-10);
+    EXPECT_NEAR(r1.z(), 0.0, 1e-10);
+
+    // At y>0, twist rotates in xz plane
+    auto r2 = op_twist(Point3D(1, 0, 0), 0.5);  // still y=0 so no twist
+    EXPECT_NEAR(r2.x(), 1.0, 1e-10);
+    EXPECT_NEAR(r2.z(), 0.0, 1e-10);
+
+    // Twist should be invertible: twist of -amount should cancel
+    auto p = Point3D(1.5, 2.0, 0.5);
+    auto twisted = op_twist(p, 0.3);
+    // Length should be preserved (rotation preserves norm in xz plane)
+    double orig_xz_len = std::sqrt(p.x()*p.x() + p.z()*p.z());
+    double twisted_xz_len = std::sqrt(twisted.x()*twisted.x() + twisted.z()*twisted.z());
+    EXPECT_NEAR(orig_xz_len, twisted_xz_len, 1e-10);
+    EXPECT_DOUBLE_EQ(twisted.y(), p.y());  // y unchanged
+}
+
+TEST(SdfOperations, Bend) {
+    // No bend: amount=0 → identity
+    auto p = Point3D(1.0, 2.0, 3.0);
+    auto r = op_bend(p, 0.0);
+    EXPECT_NEAR(r.x(), p.x(), 1e-10);
+    EXPECT_NEAR(r.y(), p.y(), 1e-10);
+    EXPECT_NEAR(r.z(), p.z(), 1e-10);
+
+    // Bend with small amount produces valid coordinates
+    auto r2 = op_bend(p, 0.1);
+    EXPECT_TRUE(std::isfinite(r2.x()));
+    EXPECT_TRUE(std::isfinite(r2.y()));
+    EXPECT_TRUE(std::isfinite(r2.z()));
+}
+
+TEST(SdfOperations, Elongate) {
+    auto f = Point3D(2.0, 1.0, 0.5);
+    auto r = op_elongate(Point3D(4, 3, 2), f);
+    EXPECT_DOUBLE_EQ(r.x(), 2.0);
+    EXPECT_DOUBLE_EQ(r.y(), 3.0);
+    EXPECT_DOUBLE_EQ(r.z(), 4.0);
+}
+
+TEST(SdfOperations, Displace) {
+    // At origin, all sin terms = 0 → no displacement
+    double d = op_displace(1.0, Point3D(0, 0, 0), 0.5, 2.0);
+    EXPECT_DOUBLE_EQ(d, 1.0);
+
+    // Displacement is bounded by amplitude
+    double d2 = op_displace(1.0, Point3D(1.0, 2.0, 3.0), 0.5, 1.0);
+    EXPECT_GE(d2, 0.5);
+    EXPECT_LE(d2, 1.5);
+}
+
+TEST(SdfOperations, CheapBend) {
+    // No bend: k=0 → identity
+    auto p = Point3D(2.0, 3.0, 4.0);
+    auto r1 = op_cheap_bend(p, 0.0);
+    EXPECT_NEAR(r1.x(), p.x(), 1e-10);
+    EXPECT_NEAR(r1.y(), p.y(), 1e-10);
+    EXPECT_NEAR(r1.z(), p.z(), 1e-10);
+
+    // At x=0, cheap bend has no effect (sin(0)=0, cos(0)=1)
+    auto r2 = op_cheap_bend(Point3D(0, 3, 1), 1.0);
+    EXPECT_NEAR(r2.x(), 0.0, 1e-10);
+    EXPECT_NEAR(r2.y(), 3.0, 1e-10);
+    EXPECT_NEAR(r2.z(), 1.0, 1e-10);
+}
+
+// ═══════════════════════════════════════════════════
+//  Composition Tests
+// ═══════════════════════════════════════════════════
+
+TEST(SdfOperations, CompositionRepeatTwistSphere) {
+    // Repeat + twist + sphere: sphere repeated in grid, then twisted
+    auto composed = [](const Point3D& p) -> double {
+        auto q = op_repeat(p, Point3D(2.0, 2.0, 2.0));
+        q = op_twist(q, 0.5);
+        return sdf_sphere(q, 0.3);
+    };
+
+    // Near origin (center of a cell) should be inside sphere
+    EXPECT_LT(composed(Point3D(0, 0, 0)), 0.0);
+
+    // Further away in the grid should also have spheres
+    // (repeat maps cell centers to origin)
+    EXPECT_LT(composed(Point3D(2, 0, 0)), 0.0);
+    EXPECT_LT(composed(Point3D(0, 2, 0)), 0.0);
+
+    // Away from cell centers should be outside
+    EXPECT_GT(composed(Point3D(1.0, 1.0, 1.0)), 0.0);
+}
+
+TEST(SdfOperations, CompositionBooleanWithDeformations) {
+    // Union of two translated spheres
+    auto scene = [](const Point3D& p) -> double {
+        double d1 = sdf_sphere(op_translate(p, Point3D(1.5, 0, 0)), 1.0);
+        double d2 = sdf_sphere(op_translate(p, Point3D(-1.5, 0, 0)), 1.0);
+        return op_union(d1, d2);
+    };
+
+    // Between the two spheres (origin) should be outside
+    EXPECT_GT(scene(Point3D(0, 0, 0)), 0.0);
+
+    // Center of each sphere should be inside
+    EXPECT_LT(scene(Point3D(1.5, 0, 0)), 0.0);
+    EXPECT_LT(scene(Point3D(-1.5, 0, 0)), 0.0);
+}
+
+TEST(SdfOperations, CompositionSmoothUnionSpheres) {
+    auto scene = [](const Point3D& p) -> double {
+        double d1 = sdf_sphere(op_translate(p, Point3D(1.0, 0, 0)), 1.0);
+        double d2 = sdf_sphere(op_translate(p, Point3D(-1.0, 0, 0)), 1.0);
+        return op_smooth_union(d1, d2, 0.3);
+    };
+
+    // Midpoint between spheres should be within the blend region
+    double d = scene(Point3D(0, 0, 0));
+    // Both spheres at origin: distance = sqrt(1^2)-1 = 0
+    // h = clamp(0.5 + 0.5*(0-0)/0.3, 0, 1) = 0.5
+    // mix(0,0,0.5) - 0.3*0.5*0.5 = 0 - 0.075 = -0.075
+    EXPECT_LT(d, 0.0);  // should be inside blend region
+    EXPECT_GT(d, -0.3); // but not too deep
+}
+
+TEST(SdfOperations, CompositionDifferenceWithTwist) {
+    // Sphere minus a smaller twisted sphere
+    auto scene = [](const Point3D& p) -> double {
+        double body = sdf_sphere(p, 2.0);
+        auto q = op_twist(p, 1.0);
+        double cutter = sdf_sphere(q, 1.0);
+        return op_difference(body, cutter);
+    };
+
+    // At origin, both spheres contain the point
+    // body: sqrt(0) - 2 = -2; cutter: sqrt(0) - 1 = -1
+    // difference = max(-(-2), -1) = max(2, -1) = 2 → outside!
+    // Wait: difference = max(-d1, d2). So max(-(-2), -1) = max(2, -1) = 2
+    EXPECT_GT(scene(Point3D(0, 0, 0)), 0.0);
+
+    // At r=1.5 from origin: body=-0.5, cutter=0.5 (at y=0 twist=0 → sphere at r=1)
+    // difference = max(0.5, 0.5) = 0.5 → outside
+    // Actually let's test a clear case: far outside both
+    EXPECT_GT(scene(Point3D(3.0, 0, 0)), 0.0);  // far outside
+}
+
+TEST(SdfOperations, MultipleDeformationsChain) {
+    // Translate → rotate → twist → repeat → sphere
+    auto composed = [](const Point3D& p) -> double {
+        auto q = op_translate(p, Point3D(1.0, 0.0, 0.0));
+        q = op_rotate(q, M_PI / 4.0);
+        q = op_repeat(q, Point3D(3.0, 3.0, 3.0));
+        return sdf_sphere(q, 0.5);
+    };
+
+    // Should not crash or produce NaN
+    double d = composed(Point3D(0, 0, 0));
+    EXPECT_TRUE(std::isfinite(d));
+
+    // Check at multiple points
+    for (double x = -4.0; x <= 4.0; x += 2.0) {
+        for (double y = -4.0; y <= 4.0; y += 2.0) {
+            for (double z = -4.0; z <= 4.0; z += 2.0) {
+                double val = composed(Point3D(x, y, z));
+                EXPECT_TRUE(std::isfinite(val));
+            }
+        }
+    }
+}